The program DOGMA was used to annotate the O. minuta cp genome (Wyman et al., 2004). The annotation results were checked manually, and codon positions were adjusted by comparison to homologs from the cp genomes of O. australiensis and O. sativa ssp. indica in the database. All transfer RNA sequences were verified using tRNAscan-SE version 1.21 (Schattner et al., 2005) with the default settings. OGDRAW (Lohse et al., 2007) was applied to illustrate the structural features of the O. minuta cp genome. To examine deviations in synonymous codon usage by avoiding the influence of amino acid composition, the relative synonymous codon usage (RSCU) was determined using MEGA 6 software (Kumar et al., 2008). mVISTA software was used in the Shuffle-LAGAN mode to compare the complete variation in the O. minuta cp genome with eleven other cp genomes using the O. minuta annotation as a reference (Frazer et al., 2004).

We employed REPuter to identify repeat sequences, including palindromic, reverse, and direct repeats, within the cp genome (Kurtz et al., 2001). The following settings for repeat identification were used: (1) Hamming distance of 3; (2) 90% or greater sequence identity; (3) a minimum repeat size of 30 bp. Phobos version 3.3.12 (Leese et al., 2008) was used to detect (SSRs) within the cp genome, with the search parameters set at ten repeat units ≥10 for mononucleotides, eight repeat units ≥8 for dinucleotides, four repeat units ≥4 for trinucleotides and tetranucleotides, and three repeat units ≥3 for pentanucleotide and hexanucleotide SSRs. Tandem repeats in the O. minuta cp genome were identified using Tandem Repeats Finder version 4.07 b (Benson, 1999) with the default settings.

Complete cp genomes as well as a separate partition using only 65 shared genes were employed to analyze the average pairwise sequence divergence for 11 Oryza species: O. australiensis, O. nivara, O. rufipogon, O. sativa L. ssp. indica, O. sativa L. ssp. japonica, O. barthii, O. glumipatula, O. longistaminata, O. meridionalis, O. officinalis, and O. punctata. Missing and ambiguous gene annotations were confirmed by comparative sequence analysis after a multiple sequence alignment and gene order comparison. These regions were aligned using MAFFT (version 7.222) (Katoh and Standley, 2013) with the default parameters. Kimura's two-parameter (K2P) model was selected to calculate pairwise sequence divergences (Kimura, 1980). To resolve the O. minuta phylogenetic position within the rice tribe (Oryzeae), 13 published cp genomes were downloaded from the NCBI database for analyses. First, multiple alignments were performed using the complete cp genomes based on the conserved structure and gene order of the chloroplast genomes (Wicke et al., 2011). Four methods were employed to construct phylogenetic trees, including Bayesian inference (BI) implemented with MrBayes 3.12 (Ronquist and Huelsenbeck, 2003), maximum parsimony (MP) with PAUP 4.0 (Swofford, 1993), and maximum likelihood (ML) and neighbor-joining (NJ) with MEGA 6 (Kumar et al., 2008) using described settings (Wu et al., 2015; Asaf et al., 2016a). In the second phylogenetic analysis, 65 shared genes from the cp genomes of 12 Oryza species and two Zizania outgroup species were aligned in ClustalX using the default settings, followed by manual adjustment to preserve reading frames. The above four phylogenetic-inference methods were used to infer trees from the 65 concatenated genes using the same settings (Wu et al., 2015; Asaf et al., 2016a).

The O. minuta cp genome was assembled by mapping all Illumina reads to the draft cp genome sequence using CLC Genomics Workbench v7.0. A total of 1,577,251 reads were obtained, with an average length of 120 bp, for 504.211X coverage of the cp genome. The consensus sequence for a specific position was generated by assembling reads mapped with at least 875 reads per position and was used to construct the complete sequence of the O. minuta cp genome. The complete O. minuta cp genome is 135,094 bp in size (GenBank: KU179220), which is similar to the already reported cp genome sizes of related Oryza species and is within the range of other angiosperms (Yang et al., 2010). The cp genome possesses a typical quadripartite structure, which includes a pair of inverted repeats (IRa and IRb 20,836 bp) and separate SSC (12,446 bp) and LSC (80,974 bp) regions (Table 1, Figure 1). The GC content (39%) of the O. minuta cp genome is very similar to that of other Oryza species cp genomes (Table 1) (Wu et al., 2015). However, the GC content is unequally distributed in the O. minuta cp genome: it is highest in the IR regions (44.3%), moderate in the LSC regions (37.1%) and lowest in the SSC regions (33.3%). This high IR GC percentage is due to the presence of eight ribosomal RNA (rRNA) sequences in these regions. These results are similar to a previously reported high GC percentage in IR regions (Qian et al., 2013).

A total of 139 genes were found in the O. minuta cp genome, of which 110 are unique, including 91 protein-coding genes, 40 tRNA genes, and 8 rRNA genes (Figure 1, Table 2). Of these, 11 protein-coding, four rRNA, and eight tRNA genes are duplicated in the IR regions. The LSC region comprises 62 protein-coding and 24 tRNA genes, whereas the SSC region comprises 11 protein-coding genes and one tRNA gene. The protein-coding genes present in the O. minuta cp genome include nine genes encoding large ribosomal proteins (rpl2, 14, 16, 20, 22, 23, 32, 33, 36), 12 genes encoding small ribosomal proteins (rps2, 3, 4, 7, 8, 11, 12, 14, 15, 16, 18, 19), five genes encoding photosystem I components (psaA, B, C, I, J), 10 genes related to photosystem II (Table 2), and six genes (atpA, B, E, F, H, I) encoding ATP synthase and electron transport chain components (Table 2). A similar pattern of protein-coding genes is also present in O. sativa (Zhang et al., 2012) and O. glaberrima (Wambugu et al., 2015). There are 11 intron-containing genes, 10 of which contain one intron, with only ycf3 genes having two introns (Table S2). The ndhA gene has the longest intron (965 bp).

Protein, rRNAs, and tRNAs are encoded by 45.1, 6.83, and 2.2% of the entire cp genome, respectively, and the remaining 45.8% is composed of non-coding regions (Table 3). The total protein-coding sequences (CDSs) are 60,948 bp in length and consist of 91 genes encoding 20,354 codons (Tables 1, 4). The O. minuta cp genome codon usage frequency was determined based on tRNA and protein-coding gene sequences (Table 5). Leucine (10.7%) and cysteine (1.2%) are the maximum and minimum commonly encoded amino acids, and isoleucine, serine, glycine, arginine and alanine are encoded by 7.9, 7.5, 7.4, 6.5, and 6.1% of CDSs, respectively (Figure S1). Similar ratios for amino acids are present in previously reported cp genomes (Qian et al., 2013; Chen et al., 2015).

Among these, the maximum and minimum codons used are ATT (820), encoding isoleucine, and TTG and ATT (1, 1), encoding methionine. The AT content is 52.5, 60.0, and 68.7% at the 1st, 2nd, and 3rd codon positions, respectively, within CDS regions (Table 4). The preference for a high AT content at the 3rd codon position is similar to the A and T concentrations reported in various terrestrial plant cp genomes (Morton, 1998; Nie et al., 2012; Qian et al., 2013). In total, 42.65 and 57% of all types of preferred synonymous codons (RSCU>1) ending with A and U and C and G, respectively, were found. Non-preferred synonymous codons (RSCU <1) are 42.40 and 57.50% for C and G and A and U. Usage of the start codon AUG and UGG, the latter encoding tryptophan, has no bias (RSCU = 1) (Table 5).

Repeat sequences, which play a role in genome rearrangements, are very helpful in phylogenetic studies (Cavalier-Smith, 2002; Nie et al., 2012). Furthermore, analyses of various cp genomes revealed that repeat sequences are essential to induce indels and substitutions (Yi et al., 2013). Repeat analysis of the O. minuta cp genome showed 20 palindromic repeats, 30 forward repeats, and 28 tandem repeats (Figure 2A). Among these, 17 forward repeats are 30–44 bp in length, with only three tandem repeats of the same length and 18 15–29 bp in length (Figures 2A–D). Similarly, 11 palindromic repeats are 30–44 bp, and 6 repeats are 45–59 bp in length (Figure 2B). Overall, 78 repeats were found in the O. minuta cp genome. Similarly, 73, 73, 76, 71 72, 78, 72, 71, 73, 77, and 74 repeat pairs were found in previously reported O. australiensis, O. nivara, O. rufipogon, O. sativa L. ssp. indica, O. sativa L. ssp. japonica, O. barthii, O. glumipatula, O. longistaminata, O. meridionalis, O. officinalis and O. punctata genomes, respectively (Figure 2A). This suggests that O. minuta is more similar to O. barthii and O. officinalis in terms of repeats. Approximately 29.4% of these repeats are distributed in protein-coding regions. Previous reports suggest that sequence variation and genome rearrangement occur due to the slipped-strand mispairing and improper recombination of these repeat sequences (Cavalier-Smith, 2002; Asano et al., 2004; Timme et al., 2007). Furthermore, the presence of these repeats indicates that the locus is a crucial hotspot for genome reconfiguration (Gao et al., 2009; Nie et al., 2012). Additionally, these repeats are an informative source for developing genetic markers for phylogenetic and population studies (Nie et al., 2012).

Simple sequence repeats (SSRs), or microsatellites, are repeating sequences of typically 1–6 bp that are distributed throughout the genome. In this study, we detected perfect SSRs in O. minuta together with 11 other Oryza species cp genomes (Figure 3A). Certain parameters were set because SSRs of 10 bp or longer are prone to slipped-strand mispairing, which is believed to be the main mechanism for SSR polymorphisms (Rose and Falush, 1998; Raubeson et al., 2007; Huotari and Korpelainen, 2012). A total of 419 perfect microsatellites were found in the O. minuta cp genome (Figure 3A). Similarly, 418, 413, 416, 416, 419, 420, 419, 419, 421, 429, and 422 SSRs were detected in O. australiensis, O. nivara, O. rufipogon, O. sativa L. ssp. indica, O. sativa L. ssp. japonica, O. barthii, O. glumipatula, O. longistaminata, O. meridionalis, O. officinalis and O. punctata, respectively (Figure 3A). The majority of SSRs in these cp genomes possess a dinucleotide repeat motif, varying in quantity from 269 in O. sativa ssp. indica to 276 in O. officinalis. Mononucleotide SSRs are the second most common, ranging from 92 in O. nivara to 100 in O. officinalis. Using our search criterion, only one pentanucleotide SSR was found in O. nivara, O. rufipogon, O. indica and O. officinalis (Figure 3A). In O. minuta, most mononucleotide SSRs are A (97%) and T (2.12.30%) motifs, with the majority of dinucleotide SSRs being A/G (47.05%) and A/T (38.60%) motifs (Figure 3B). Approximately 62% of SSRs are located in non-coding regions; approximately 4.3% are present in rRNA sequences and 2.3% in tRNA genes (Figure 3C). Further analysis revealed that approximately 66.82% of SSRs occur in the LSC region, whereas 24.34 and 8.83% were found in IR and SSC regions, respectively (Figure 3D). These results are similar to previous reports that SSRs are unevenly distributed in cp genomes, and the findings might provide more information for selecting effective molecular markers for detecting intra- and interspecific polymorphisms (Powell et al., 1995a,b; Provan et al., 1997; Pauwels et al., 2012). Furthermore, most mononucleotides and dinucleotides are composed of A and T, which may contribute to bias in base composition, consistent with other cp genomes (Li et al., 2013). Our findings are comparable to previous reports that SSRs found in cp genome are generally composed of polythymine (polyT) or polyadenine (polyA) repeats and infrequently contain tandem cytosine (C) and guanine (G) repeats (Kuang et al., 2011). Therefore, these SSRs identified contribute to the AT richness of the O. minuta cp genome, as previously reported for various species (Kuang et al., 2011; Chen et al., 2015).

Eleven complete cp genomes within the Oryza genus (O. australiensis, O. nivara, O. rufipogon, O. sativa L. ssp. indica, O. sativa L. ssp. japonica, O. barthii, O. glumipatula, O. longistaminata, O. meridionalis, O. officinalis, and O. punctata) were selected for comparison with that of O. minuta (135,094 bp). O. australiensis has the largest genome, and this difference is mostly attributed to variation in the length of the LSC region (Table 1). Analysis of genes with known functions showed that O. minuta shares 65 protein-coding genes with eleven other Oryza species. The number of unique genes found in O. australiensis, O. nivara, O. rufipogon, O. sativa L. ssp. indica, O. sativa L. ssp. japonica, O. barthii, O. glumipatula, O. longistaminata, O. meridionalis, O. officinalis, and O. punctata was 110, 100, 101, 108, 80, 104, 104, 104, 100, 104 and 114, respectively (Table S3). Furthermore, the O. minuta cp genome has a gene content and organization that are similar to other Oryza species and members of Poaceae (Wicke et al., 2011); however, as for other grasses, it lacks a ycf1 gene, and the accD gene is a truncated pseudogene. Because these genes are essential for the survival of photosynthetic plants (Drescher et al., 2000; Kode et al., 2005), they were most likely functionally transferred to the nucleus or functionally replaced by a eukaryotic gene, as observed for the accD plastid gene in other plant families (Babiychuk et al., 2011; Rousseau-Gueutin et al., 2011).

Pairwise cp genomic alignment between O. minuta and the 11 other genomes showed a high degree of synteny. The O. minuta cp genome annotation was used as a reference for plotting the overall sequence identity of the cp genomes of the 11 Oryza species in mVISTA (Figure 4), and the results revealed high sequence identity with all 11 Oryza species. However, except for O. australiensis, relatively lower identity was also observed with these species in various comparable genomic regions, particularly the rps3, rpl22, rpl23, rpl2, and rps19 regions (Figure 4). In addition, the LSC and SSC regions show less similarity than the two IR regions in all Oryza species. In addition, non-coding regions exhibit greater divergence than coding regions. These highly divergent regions include rbcL, rps16-trnQ, trnfM-trnM, psbM-petN, rpoC2, atpI-atpH, ndhA rpl33, petA-psbJ, ccsA, ndhF-rpl32, and ycf3. Similar results related to these genes were also reported by Qian et al. (2013). Our results also confirm similar differences among various coding regions in the analyzed species, as suggested by Kumar et al. (2009).

We compared the cp genomes and calculated the average pairwise sequence divergence among the 12 species (Table S4). Of these, the O. minuta genome has 0.005 average sequence divergence, and high divergence was found for O. australiensis (0.00725); O. officinalis has the lowest average sequence divergence (0.0044). Furthermore, the twelve most divergent genes among these genomes are petG, matK, infA, ccsA, rpoC2, clcP, psbE, rbcL, psbN, rps18, rpl36, and ndhF. The highest average sequence distance was found for rpoC2 (0.01983), followed by petG (0.0154) (Figure 5). Both these genes are located in LSC regions and display a trend toward more rapid evolution.

Expansion and contraction at the borders of IR regions are the main reason for size variations in the cp genome and play a vital role in its evolution (Raubeson et al., 2007; Wang et al., 2008; Yang et al., 2010, 2014). A detailed comparison on four junctions (J_LA, J_LB, J_SA, and J_SB) between the two IRs (IRa and IRb) and the two single-copy regions (LSC and SSC) was performed among O. australiensis, O. nivara, O. rufipogon, O. sativa L. ssp. indica, O. sativa L. ssp. japonica, O. barthii, O. glumipatula, O. longistaminata, O. meridionalis, O. officinalis and O. punctata with regard to O. minuta by carefully analyzing the exact IR border positions and adjacent genes (Figure 6). Despite the similar length of the O. minuta IR region with the other eleven Oryza species, from 20,836 bp to 20,840 bp, some IR expansion and contraction was observed. J_LA is located between rps19 and psbA, and variation in distances between rps19 and J_LA range from 40 to 49 bp across all species; the distance in O. minuta is 46 bp. The distance between psbA and J_LA is 81 bp in O. minuta, which is similar to the other genomes (81 bp). The distance between rpl22 and J_LB varies from 23 bp to 29 bp. In O. minuta, 1-bp variations exist in the J_SA border region compared to the other cp genomes. The ndhH gene traverses the SSC and IRa regions, with approximately 164 bp located in the IR region for O. minuta. Furthermore, there are 16-bp variations observed compared with O. officinalis for ndhF, ndhH and rps15 in the SSC and IRb regions, located 41 bp, 164 bp and 302 bp from the J_SB and J_SA border regions, respectively.

The Oryza genus is composed of 23 species distributed in different regions of America, Africa, Asia, and Australia (Ge et al., 1999). Continued efforts have expanded our ability to differentiate among and to understand the genomic structure and phylogenetic relationships of rice species (Khush, 1997). Taxonomy and phylogeny of the rice genus have been extensively investigated at genus level (Ge et al., 1999; Zhu and Ge, 2005; Jacquemin et al., 2013). Previous evolutionary relationships among different rice genomes and species were estimated by nuclear and chloroplast DNA restriction fragment-length polymorphisms (Ge et al., 1999; Zou et al., 2015), but complete genome sequencing provides more detailed insight (Wambugu et al., 2015; Wu et al., 2015; Asaf et al., 2016b). In this regard, O. minuta has been poorly investigated. In this study, the phylogenetic position of O. minuta within Oryza was established by utilizing complete cp genomes and 65 shared genes among 12 Oryza members (Figures 7A,B). Two species, Zizania aquatic and Zizania latifolia were set as outgroups. Phylogenetic analysis using Bayesian inference (BI), maximum parsimony (MP), maximum likelihood (ML) and neighbor-joining (NJ) methods were performed. The results showed same phylogenetic signals for the complete cp genomes and 65 shared genes of O. minuta. The complete genome sequences (Table S5) and 65 shared genes (Tables S3, S6) from all species generated phylogenetic trees with same topologies (Figures 7A,B). In these phylogenetic trees based on the entire genome data set and 65 shared genes, O. minuta formed a single clade with O. punctata, with high BI and bootstrap support using four different methods (Figures 7A,B). Furthermore, the tree topology confirmed the relationship inferred from the phylogenetic work conducted by Ge et al. (1999) and Zou et al. (2015). This position of O. minuta confirms the previously published phylogeny described by Ge et al. (1999). Ge et al. (1999) reported that O. minuta BBCC shares a clade with O. punctata BB with regard to Adh1, whereas it forms a clade with O. officinalis CC in the Adh2 phylogenetic analysis. Similar resuls was suggested by Zou et al. (2015), whereby phylogenetic analysis of the four nuclear loci and three meternally interited chloroplast fragments from different Oryza species grouped O. minuta in a clade with maternal parent O. punctata BB (Zou et al., 2015). As the phylogenetic tree based on the matK gene represents the maternal genealogy of rice species, which can offer an opportunity to identify maternal parents of allotetraploid species, we performed an additional phylogenetic analysis of O. minuta using the matK gene from related species (Figure S2). The results revealed a single clade for O. minuta with parental O. punctata. Similar results was also suggested by Ge et al. (1999), whereby phylogenetic analysis of the matK gene from different Oryza species grouped O. minuta in a clade with the maternal parent O. punctata BB instead of O. officinalis CC. Furthermore, the result suggests that there is no conflict between the entire genome data set and 65 shared genes of these cp genomes.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.