Fresh leaves of the four species of Oryza were collected from the Kunming Institute of Botany, Chinese Academy of Sciences and the Rice Research Institute, Guangdong Academy of Agricultural Sciences. Fresh leaves from each accession were immediately dried with silica gel before DNA extraction. Total genomic DNA was extracted and purified following the method of Li et al. (2013). DNA samples were randomly fragmented into 400–600 bp fragments using an ultrasonicator. An Illumina paired-end DNA library with 500-bp insert size was constructed using a NEBNext^® Ultra^TMDNA Library Prep Kit following the manufacturer’s instructions. Paired-end sequencing (2 × 150 bp) was conducted on an Illumina HiSeq X Ten platform.

The paired-end reads were qualitatively assessed and assembled with SPAdes 3.6.1 (Bankevich et al., 2012). Chloroplast genome sequence contigs were selected from SPAdes software by performing a BLAST search using the Oryza sativa chloroplast genome sequence as a reference (GenBank accession number: JN861110). The selected contigs were assembled with Sequencher 5.4.5 (Gene Codes, Ann Arbor, MI, United States). Small gaps in the assemblies were bridged with specific primers designed for PCR based on their flanking sequences and then by Sanger sequencing. Based on the reference chloroplast genome, the four junctions between the inverted repeats (IRs) and single-copy regions were checked by amplification with specific primers followed by Sanger sequencing (Dong et al., 2013). Chloroplast genome annotation was performed with Plann (Huang and Cronk, 2015) using the Oryza sativa reference sequence from Genbank. A chloroplast genome map was drawn using Genome Vx software (Conant and Wolfe, 2008).

MIcroSAtellite (MISA¹) was used to detect single sequence repeats (SSR) within the chloroplast genome, with the parameters set at >10 for mononucleotide, >5 for dinucleotide, >4 for trinucleotide, and >3 for tetranucleotide, pentanucleotide, and hexanucleotide SSRs. The web-based software REPuter was used to visualize the dispersed repeat sequences in Oryza by forward vs. reverse complement (palindromic) alignment (Kurtz et al., 2001). The following settings for repeat identification were used: (1) a minimum repeat size of 30 bp; (2) Hamming distance of 3; and (3) 90% or greater sequences identity. Tandem repeats were identified using web-based Tandem Repeats Finder², with 2, 7, and 7 set for the alignment parameters match, mismatch, and indel, respectively.

Complete chloroplast genomes were used to analyze the average pairwise sequence divergence for the 18 Oryza species (Table 1). All Oryza sequenced chloroplast genomes were aligned using MAFFT v7 (Katoh and Standley, 2013) assuming collinear genomes for the full alignment and then adjusted manually using Se-Al 2.0 (Rambaut, 1996). Variable and parsimony-informative base sites were calculated using MEGA 6.0 software (Tamura et al., 2013). The p-distances among Oryza chloroplast genomes were calculated to evaluate the divergence of Oryza species using MEGA software.

Phylogenetic trees were constructed by maximum parsimony (MP), maximum likelihood (ML), and Bayesian inference (BI) using the entire chloroplast genome. The chloroplast genome sequences of Leersia tisserantii, Zizania latifolia, Zizania aquatic, Rhynchoryza subulata, Chikusichloa aquatic, Potamophila parviflora, and Microlaena stipoides were used as out-groups.

Maximum parsimony analyses were conducted using PAUP v4b10, which was performed using heuristic searches with the ‘MulTrees’ option followed by tree bisection–reconnection (TBR) branch swapping. To assess node support, bootstrap analyses were performed using 1000 replicates with 10 random taxon additions and heuristic search options. All characters were equally weighted, gaps were treated as missing, and character states were treated as unordered. ML analyses were performed using RAxML v.8.1.24 with general time reversible (GTR) + G model. Supports for nodes were assessed with 1,000 rapid bootstrapping replicates. BI was conducted with Mrbayes v3.2 (Ronquist et al., 2012). The analysis was run for 1,000,000 generations and sampled every 1,000 generations. The first 25% of the trees were discarded as burn-in, the remaining trees were used to build a 50% majority-rule consensus tree and estimate the Bayesian posterior probabilities.

Hypervariable chloroplast regions were identified using slideAnalyses of SPIDER (Brown et al., 2012) version 1.2-0. This function extracts all the passable windows of a chosen size in a DNA alignment and performs pairwise distance (K2P) analyses of each window. Two factors were considered for the definition of hypervariable regions: first, the mean distance of each window, and second, the proportion of zero pairwise distances for each species in the matrix. The step size was set to 200 bp, with an 800-bp window length. Two data sets were created for this analysis: (1) all 18 species data set and (2) only AA genome species data set (included nine species).

We analyzed the hypervariable barcodes and compared the chloroplast genes rbcL, matK, and trnH-psbA using two different methods. To assess the barcoding resolution for all barcodes, the distance and the tree-building methods were the two analytical approaches. Distance may be the most commonly used method for DNA sequences classification, which was used to evaluate the barcode performances of the newly identified highly variable regions. The function nearNeighbour of SPIDER was used for barcoding analyses (Brown et al., 2012).

Tree building analyses provide a convenient and visualized method for evaluating discriminatory performance by calculating the proportion of monophyletic species. MP trees were constructed for each hypervariable marker and the different marker combinations using MEGA 6.0. Relative support for the branches of the MP tree was assessed via 1,000 bootstrap replicates. Species were considered successfully discriminated when all individuals of a species formed a single and exclusive clade in a MP tree with bootstrap value above 50%.

Using an Illumina HiSeq X ten system, the four species of Oryza were sequenced, producing from 9,642,763 to 10,287,100 paired-end raw reads (150 bp for average read length). After mapping the paired-end reads with the reference chloroplast genome of Oryza sativa, 105,824 to 499,727 reads were extracted, with 117× to 556× coverage (Table 2). PCR-based sequencing validated four junction regions in each of the Oryza chloroplast genomes. The four Oryza chloroplast genome sequences were deposited in GenBank (accession numbers, MF401450–MF401453).

The chloroplast genome of the four species of Oryza ranged from 134,748 bp (O. rhizomasis) to 135,239 bp (O. meyeriana). Chloroplast genomes showed a typical quadripartite structure, consisting of a pair of IRs (20,817–20,822 bp) separated by the LSC (80,788–81,212 bp) and SSC (12,326–12,497 bp) regions (Figure 1 and Table 2). For all four chloroplast genomes, the average GC content was 39.0%. The annotation chloroplast genomes of these four species were represented in one circular map, because their gene number, order and names were the same. With the duplicated genes in IR regions only counted once, the chloroplast genomes of the four species of Oryza harbored 110 different genes, including 77 protein-coding genes, 29 tRNA genes and 4 rRNA genes. The gene number and the genome organization are very similar to the chloroplast genomes of other Oryza species.

In this study, we detected SSRs and dispersed repeats in all published Oryza species chloroplast genomes. Each Oryza chloroplast genome contained 22–31 SSRs (Figure 2). Among those SSRs, most were located in the LSC/SSC regions (92.0%). For each species, mono-, di-, tri-, tetra-, penta-, and hexanucleotide SSRs were all detected, with the average of mono-, di-, tri-, and tetranucleotide SSRs 38.69%, 15.64%, 11.63%, and 32.77%, respectively. Additionally, we did not find many penta- and hexanucleotide SSRs in Oryza, which is similar to the chloroplast genome of angiosperm species (Choi et al., 2016; Xu et al., 2017). SSRs were particularly rich in AT in the Oryza chloroplast genomes. The majority of SSRs in all species were A/T mononucleotides.

Dispersed repeat sequences, which play a role in genome rearrangement, have been used as a source to understand the phylogenetic relationships of species. Repeat sequences with a repeat unit longer than 30 bp were analyzed (Figure 3). Each Oryza chloroplast genome contained 50–63 repeat sequences, including 30–42 forward repeats, 19–22 palindromic repeats, and 22–34 tandem repeats. The repeats were primarily located in non-coding regions of chloroplast genomes.

Eighteen Oryza chloroplast genomes were fully aligned, and the alignment matrix of 139,988 bp was obtained. The alignment revealed high sequence similarity across the Oryza chloroplast genomes, suggesting that they are highly conserved. We retrieved 3,593 variable sites (2.57%), including 2,017 singletons and 1,576 parsimony-informative sites (1.13%), in the total chloroplast genome.

The number of nucleotide substitutions and p-distance were used to estimate the divergence among the eighteen Oryza species. The number of nucleotide substitutions among seven species varied from 13 to 1,953, and the p-distance varied from 0.0001 to 0.0149 (Supplementary Table S1). The overall sequence divergence estimated by p-distance among the four genomes was only 0.0052. The largest sequence divergence was observed between O. minuta and O. brachyantha and the lowest divergence was between O. barthii and O. glaberrima. At the level of the genome type, the largest sequence divergence occurred between the FF genome group and AA or BB genome group (Supplementary Table S2), with the sequence distance (0.0147). The smallest sequence divergence was observed between the CC genome group and CCDD genome group, with the sequence distance (0.0012).

Phylogenetic relationships within the Oryza were reconstructed using MP, ML, and BI analyses. The topologies based on the three methods were highly supported and largely congruent. Figure 4 illustrates the phylogeny generated by the ML analysis, including three types of support values: MP bootstrap values (MP-BS), ML bootstrap values (ML-BS), and BI posterior probabilities (BI-PP). All analyses fully supported the monophyly of the genus Oryza (ML-BS/BI-PP/MP-BS = 100/1.0/100). Within Oryza, the eighteen species divided into seven major groups with different genome types (AA, BB, CC, CCDD, EE, FF, and GG). The CCDD genome species formed a monophyletic clade with the CC genome, which suggested that the CC genome served as the maternal parent of the CCDD genome species. Two BBCC genome species, O. minuta and O. eichingeri, had different maternal origins, with a maternal parent of the BB genome for O. minuta and a maternal parent of the CC genome for O. eichingeri. O. longistaminata, and O. glumaepatula formed the basal clade in the AA group, which are found in Africa and South America, respectively.

To identify the hypervariable chloroplast regions, pairwise distance values within 800 bp in the eighteen Oryza chloroplast genomes were calculated with slideAnalyses of SPIDER (Figure 5). Three regions showed remarkably high distances and low proportions of zero pairwise distances for each species in the 18 species data set. One region was the coding region of rpoC2 and two were the intergenic regions (rps16-trnQ and rbcL-accD). For the data set of AA genome species, four variable chloroplast regions (rps16-trnQ, trnTEYD, psbE-petL and rbcL-accD) were identified.

The variability of the five developed variable regions was tested together with three conventional candidate DNA barcodes (matK, rbcL and trnH-psbA). Features of the eight barcode data set are shown in Table 3. The conventional candidate DNA barcodes had lower variability than that of the newly developed markers. The highest variability was in the rpoC2 region (6.80%), followed by the regions psbE-petL region (6.75%), rps16-trnQ (6.73%), trnTEYD (6.00%), and rbcL-accD (5.22%).

Of the single-region barcodes, rpoC2 and rbcL-accD had the highest rate of correct identifications (72.22%), followed by rps16-trnQ and psbE-petL (66.67%) and trnTEYD (50.00%). The conventional candidate DNA barcodes had less discriminatory power, e.g., trnH-psbA had only a 16.67% success rate. Of the two-region barcodes, the best performing was rpoC2 + psbE-petL (94.44%). With the two core DNA barcodes rbcL and matK combined, success was only 33.33%. When analyzing multi-region barcodes, the highest correct identifications (100%) were with rps16-trnQ + trnTEYD + rbcL-accD, rps16-trnQ + trnTEYD + rpoC2, and rpoC2 + trnTEYD + psbE-petL. The MP tree method generated a graphical representation of the results, and the results of the tree-based method were the same (Figure 6).