1. Introduction

Echinochloa phyllopogon (= Echinochloa oryzicola) proliferation seriously threatens rice production worldwide. As a C4-photosynthetic weed, E. phyllopogon is highly adapted to rice (C3-photosynthesis type) planting environments, where it causes significant rice yield loss (Holm et al., 1979; Rao et al., 2007; Yamasue, 2001). Furthermore, E. phyllopogon has evolved resistance to various herbicides in different areas (Heap, 2015). Understanding the genetic diversity of agricultural pests, such as E. phyllopogon, is important for both evolutionary and population biology, and critical for agricultural management (Sun et al., 2015).

Microsatellite markers (simple sequence repeats, SSR) and single-nucleotide polymorphisms (SNP) are useful tools for studying genetic diversity and evolution (Zhang et al., 2011), and for developing high density genetic maps (Zhang et al., 2012). SSRs are short tandem repetitive sequences, which are co-dominant, abundant, multi-allelic, uniformly distributed, and can be detected by simple reproducible assays (Wang et al., 2015). SNPs are usually bi-allelic and characterized by low mutation rates; and thus, SNPS are stable from generation to generation across the genome (Kruglyak, 1997). This stability coupled with the abundance of SNPs makes them very useful both for linkage and genetic diversity studies (Talukder et al., 2014). To date, there are only eight SSR markers available for E. phyllopogon (Osuna et al., 2011; Lee et al., 2015), and an even more limited number of SNPs.

One promising approach to reduced-representation genomics is restriction site-associated DNA (RAD) sequencing, which sequences short DNA fragments flanking restriction enzyme cut sites, allowing orthologous sequences to be targeted across multiple samples to identify and score thousands of genetic markers (Miller et al., 2007). Therefore, a RAD sequencing approach can be successfully used to identify genome-wide SSRs (Gupta et al., 2015; Orjuela et al., 2010) and SNPs (Baird et al., 2008; Talukder et al., 2014; Vandepitte et al., 2013) in different species. In this study, we describe the generation of genomic RAD tags from E. phyllopogon plants. The RAD tags were sequenced using the Illumina platform and then annotated/categorized. These data allowed the discovery of a large number of SSR and SNP markers.

2. Material and methods 2.1. DNA isolation

Seeds from E. phyllopogon individuals were collected and cultivated to fruiting stage in a greenhouse at Nanjing Agricultural University. Two E. phyllopogon plants with typical characteristics were used for SSR identification. Total genomic DNA was extracted from young leaves using DNeasy Plant Mini Kits (Qiagen, USA) according to the manufacturer's protocol.

2.2. RAD library preparation, sequencing and assembly

The RAD library was constructed at Hengchuang Inc. (China), according tothe protocol described by Baird et al. (2008). Briefly, genomic DNA (300 ng) was digested for 60 min at 37 ℃ in a 50 mL reaction containing 20 U each of SgrAI and PstI (New England Biolabs, Beverly MA, USA). Reactions were stopped by incubating at 65 ℃ for 20 min. The P1 adapter (a modified Illumina adapter, see Baird et al., 2008) was ligated to the products of the restriction reaction, and the "barcoding" of the various samples was achieved with a set of index nucleotides in the P1 adapter sequence. A 2.5 mL aliquot of 100 nM P1 adapter was added to each sample, along with 1 mL 10 mM ATP (Promega), 1 mL 10 × NEBBuffer4, 1mL (equivalent to 1000 U) T4 DNA ligase (Enzymatics, Inc) and 5 mL water, then incubated at room temperature for 20 min, before heat-inactivated (20 min at 65 ℃). The reactions were then pooled and the products randomly sheared to a mean size of 500 bp using a Bioruptor (Diagenode). The material was electrophoresed through a 1.5% agarose gel, and the DNA in the range 300–800 bp isolated using a MinElute Gel Extraction Kit (Qiagen). dsDNA ends were treated with end blunting enzymes (Enzymatics, Inc) to remove overhangs, and the samples purified using a MinElute column (Qiagen). 3'-adenine overhangs were then added by the addition of 15 U Klenow exo-(Enzymatics), followed by incubation at 37 ℃ for 10 min. Following repurification, 1 mL 10 mM P2 adapter (a modified Illumina adapter, see Baird et al., 2008) was ligated, as described above for P1. The samples were then purified as above, and eluted in a volume of 50 mL. Following quantification (Qubitfluorimeter), 20 ng were taken as thetemplatefor a 100 mL PCR containing 20mL Phusion Master Mix (NEB), 5mL 10 mM P1 adapter primer (Illumina), 5 mL 10 mM P2 adapter primer (Illumina) and water. The Phusion PCR settings followed product guidelines (NEB) over 18 cycles. The amplicons were gel purified, the size range 300–700 bp was excised from the gel, with the DNA content adjusted to 3 ng/mL. The constructed RAD libraries were sequenced on the NGS Illumina platform PE150 at Hengchuang Inc. (China), following the manufacturer's protocol.

To obtain clean, high quality reads, we discarded low quality raw sequences with adapter contamination or N content >10%. We used Stacks software for RAD tag clustering for each sample (ustacks). The Reads group (Read1 and Read2) at a same enzyme loci RAD were assembled by using the ABYSS software (Catchen et al., 2011).

2.3. SSR identification

SSR motifs were identified by SSRIT software (http://www.gramene.org/db/markers/ssrtool) using default parameters (Temnykh et al., 2001). Both perfect and imperfect di-, tri-, tetra-, penta-and hexa-nucleotide motifs were targeted. Di-nucleotide motifs with at least 4 repeats and other motifs with at least 3 repeats were selected. We used Primer3 software (http://sourceforge.net/projects/primer3/) to design primers in the flank regions of SSR sequences (SSR sequences were not contained in the primers), the replicated primers were removed and unique primers and relative loci were retained.

To analyze the frequency of SSR motifs, SSRs were first standardized (Wang et al., 2015). For example, SSRs with motifs of AT and TA were analyzed as AT, and motifs of ATG, TGA, GAT, TAC, ACT and CAT are analyzed as ATG.

2.4. Sequence annotation

For the contigs with SSR loci, sequence annotation and Gene Ontology analyses were further conducted. BlastN searches were performed against the Gene Ontology database (http://www.geneontology.org/), using 90% identity and a minimum alignment of 100 bp as cut-off parameters. A threshold E-value of e^-15 was adopted for each annotation. The annotated sequences were assigned a function based on the Gene Ontology database (http://www.geneontology.org/); GO terms were determined with respect to cellular component, biological process and molecular function (Barchi et al., 2011).

2.5. SNP discovery

SNPs were detected by Stacks pipeline, ustacks software was used to build loci, cstacks software was used to create a catalog of loci, and sstacks software was used to match samples back against the catalog (Catchen et al., 2011). Default settings were used in Stacks.

2.6. Microsatellites amplification

To test the validity of the SSRs identified by RAD sequencing here, we used eight SSRs (Table 1) to study the genetic diversity of four E. phyllopogon populations collected from rice fields in China. We extracted total genomic DNA from four-leaf stage plants using a DNeasy PlantMini Kit(TiangenBiotech, Beijing, China) according tothe manufacturer's instructions. Isolated DNA concentration and relative purity were checked using Nanodrop ND-1000 (Thermo Scientific), and adjusted to 30–40 ng/mL. Forward primers of SSRs were labeled with fluorescent tags (Table 1). PCR amplification was conducted in a total volume of 10 mL. The PCR mixture contained 0.2 mL of DNA, 0.4 mL of each primer (10 mM), 5 mL of 2 × PCR Taq Mix (Dongsheng Biotech, China), and ddH₂O to a final volume of 10 mL. The amplifications were performed using the following cycling program: initial denaturation at 94 ℃ for 4 min, followed by 35 cycles of 94 ℃ for 30 s, relative annealing temperatures for 30 s, and 72 ℃ for 1 min, with a final extension step at 72 ℃ for 10 min. The amplification products were combined with formamide and a size standard GeneScan-500 LIZ (Applied Biosystems, Foster City, California, USA), and separated on a 3730 ABI automated sequencer (Applied Biosystems). Sample profiles were scored manually using GeneMarker v. 2.4 (Applied Biosystems).

2.7. Data analysis

The multilocus data were transformed to a binary matrix of presence/absence of each allele for each individual, which was used for further analysis with GenAlex 6.5 (Peakall and Smouse, 2012; Teixeira et al., 2014). Total number of alleles and the number of private alleles for each population were determined using GenAlex 6.5, and genetic diversity was determined using GenoDive2.0b23 (Teixeira et al., 2014), according to the tutorials (www.patrickmeirmans.com/software/GenoDive.html). GenoDive allows analyzing polyploids with unknown dosage of alleles (Meirmans and VanTienderen, 2004).

3. Results 3.1. Sequencing and contig assembly

The sequencing procedure generated 71.45 million reads for the two E. phyllopogon samples (Table 2). After editing/trimming, 10, 440.6 Mb of high quality sequences were available, which were assembled into 37, 662 contigs. Average contig lengths for the two samples were 334 and 346 bp. The GC content of E. phyllopogon was 45.8%.

3.2. Identification of SSRs

A screen of the dataset resulted in the identification of 4710 putative SSRs that permitted PCR primer design for E. phyllopogon. Tables S1 and S2 show motifs, number of repeats, sequence of 5'-and 3'-flanking, sequences and annealing temperatures of primers, sequence of PCR products and the potential relative genes for each SSR loci.

The majority of motifs among the RAD SSRs were dinucleotide (>82%) for both samples, and 14%–15% of the SSR motifs were trinucleotide (Table 3). The majority of SSRs were four motifrepeats. The abundance of SSRs decreased significantly (P < 0.01) with increasing motif-repeats for E. phyllopogon (Fig. 1).

Nearly all (97.3%) E. phyllopogon SSR motifs consisted of dinucleotide plus trinucleotide repeats. Thus, we further analyzed dinucleotide and trinucleotide motifs. Before the analysis, SSRs were standardized. For example, SSRs with motifs of AT and TA were analyzed as AT, and motifs of ATG, TGA, GAT, TAC, ACT and CAT were analyzed as ATG. AT was the most frequent, accounting for 36.3%–37.5%, followed by AG and AC (Table 4). Among the four kinds of dinucleotide motifs, CG dinucleotide repeats represented the lowest percentage of all SSRs ( < 6%). CCG was the most frequent kind of trinucleotide motif for both samples (Table 4), accounting for about 4% of the total SSRs for E. phyllopogon. The predicted length of PCR products amplified by SSR primers designed in this study are shown in Table 4.

In total, 78 putative polymorphic SSR loci were found by RAD sequencing (Table 5). These 78 SSRs include 65 SSRs with dinucleotide motifs, 10 SSRs with trinucleotide motifs, two with tetranucleotide motifs and one with a pentanucleotide motif. The AT dinucleotide repeat, which accounts for 49.4% of all motifs, was the most frequent kind.

To test the validity of the SSRs identified by RAD sequencing here, we used eight SSRs to study the genetic diversity of four E. phyllopogon populations collected from rice fields in China. We amplified 66 alleles from the eight microsatellite loci. The primer sequence EG_305 amplified 14 alleles, EG_302 amplified 12 alleles, and EG_320 and EG_1 amplified five alleles (Table 1). EG_305 amplified three to six alleles per sample, while EG_307 and EG_320 amplified two to five and two to four alleles per sample, respectively. Moreover, EG_305 amplified the most alleles on average (4.1). On average, 3.1–4.8 alleles were amplified from one locus per population (Table 6). All four populations showed private alleles, among which the populations EP13 and EP50 showed 13 and eight private alleles, respectively. The heterozygosity values of these populations ranged from 0.064 to 0.091, and their Shannon's information indices ranged from 0.087 to 0.381. Analysis of molecular variance (AMOVA) indicated that 39% of diversity occurs among populations, while 61% of diversity occurs within populations (Table 7).

3.3. Annotation of contigs with SSR loci

Using two E. phyllopogon individuals, we identified 4710 SSR loci in 4132 contigs, and annotated 643 contigs (Table S2). Among these 643 contigs, 8631 annotations, potentially referring to 2155 unigenes, were searched (a given gene product can be associated with more than one annotation). Annotated E. phyllopogon sequences with SSR loci were functionally assigned and arranged into Gen Ontology (GO) slim categories (Fig. 2). GO analyses suggested that contigs with SSR loci were mostly related to metabolic processes (12.1% of the total 2155 unigenes) and cellular processes (10.5%) among biological processes; cell (9.8%), cell part (9.7%) and organelle (8.6%) among cellular components; and binding (12.6%) and catalytic activity (8.7%) among molecular functions.

3.4. SNP discovery

In total, 49, 179 SNPs were discovered between the two samples of E. phyllopogon. Table S3 shows the kind, sequence and location of 49, 179 SNPs discovered between two samples of E. phyllopogon. Among these SNPs, transversions (67.1% of total SNPs) were much more frequent than transitions (Fig. 3).

4. Discussion 4.1. High GC content of E. phyllopogon genome

Higher GC content in plant genomes possibly contributes to an increased ability to adapt to various arable lands that are mainly maintained and regulated by human disturbance. Šmarda et al. (2014) studied GC content in 239 different plant genomes, finding that the GC content of monocots varied between 33.6% and 48.9%, and increased GC content was documented in species able to grow in seasonally cold and/or dry climates, which possibly indicates GCrich DNA may confer more stability during cell freezing and desiccation. The GC content of E. phyllopogon was higher than those of many monocots such as Juncus inflexus (33.7%), Luzula badia (33.6%), Carex acutiformis (35.6%), Schoenoplectus lacustris (35.8%), Canna indica (39.7%), Oryza sativa (43.6%) and Triticum aestivum (44.7%); and only lower than those of a few Poaceae species such as Stipa calamagrostis (47.5%) and Zea mays (47.4%) (Raats et al., 2013).

4.2. Characteristics on SSR motifs of E. phyllopogon

The majority of RAD SSR motifs were dinucleotide and with four motif-repeats. Gupta et al. (2015) identified SSR motifs in peanut (Arachis hypogaea) through RAD sequencing, and found that 67.6% of the motifs were dinucleotide, 14.6% were trinucleotide, 12.5% were tetranucleotide, 3.2% were pentanucleotide and 2.2% were hexanucleotide. Nevertheless, in eggplant (Solanum melongena), the percentages among total motifs with two to six nucleotides of dinucleotide, trinucleotide, tetranucleotide, pentanucleotide and hexanucleotide were 20.4%, 37.9%, 12.8%, 18.1% and 10.9% (Barchi et al., 2011). Using RAD sequencing in eggplant, Barchi et al. (2011) found that AAC was the most frequent kind of motif, accounting for 19.0% of the total SSRs, followed by AT (9.6%). Wang et al. (2015) analyzed the genomes of nine plant species from the Poaceae family, and found that among the genome SSRs of O. sativa ssp. indica, O. sativa ssp. japonica, Phyllostachys heterocycla, Sorghum bicolor and Z. mays, AT was the most frequent motif, and also very frequent in other Poaceae plants.

To test the validity of the SSRs identified by RAD sequencing here, we used eight SSRs to study the genetic diversity of four E. phyllopogon populations collected from rice fields in China. All eight loci were polymorphic, particularly when compared with the five SSRs that have been used for Echinochloa since 2002 (Danquah et al., 2002; Nozawa et al., 2006; Lee et al., 2015).

4.3. Potential usage of the SSRs and SNPs identified

A great number of Echinochloa species are aggressive invaders and managing crop lands requires unique strategies for each (Holm et al., 1979; Tabacchi et al., 2006). Thus, correctly identifying Echinochloa spp. is of agronomical and economic importance. The genus Echinochloa contains about 35 species that are widespread in both tropical and temperate regions and in dry or water-flooded soils (Flora of China, 2015). The taxonomy of this genus is complex, and Echinochloa species show wide variability in morphological, biological and physiological features (Danquah et al., 2002; Tabacchi et al., 2006; Vidotto et al., 2007). Conventionally, the identification of Echinochloa species has been attempted taxonomically using morphological assessment of plants, which has frequently been found to be difficult and uncertain (Tabacchi et al., 2006). Moreover, there are different taxonomic key systems for Echinochloa species, which may lead to misidentification (Flora of China, 2015; Tabacchi et al., 2006). Molecular identification of the Echinochloa species is not yet reliable and requires further study (Danquah et al., 2002; Kaya et al., 2014; Tabacchi et al., 2006). In addition, molecular markers may be very useful in studying the origin and distribution of herbicide-resistant populations (Okada et al., 2013; Osuna et al., 2011). SNPs and SSRs are ideal molecular tools for gene location and molecular breeding (Danquah et al., 2002; Gupta et al., 2015; Vandepitte et al., 2013; Zhang et al., 2011).

Acknowledgements

This research was supported by China Postdoctoral Science Foundation (2015M571763) and the Special Fund for Agroscientific Research in the Public Interest of China (201303022). We thank Kui Wu and Zhi-hui Yan (Nanjing Agricultural University, China) for providing helps on plant cultivation. Thanks are also due to the reviewers and editors for their helpful comments and English polish on earlier drafts of the manuscript.

Appendix A. Supplementary data

Supplementary data related to this article can be found at https://doi.org/10.1016/j.pld.2017.08.004.