Comparative Plastomes and Phylogenetic Analysis of Cleistogenes and Closely Related Genera (Poaceae)

Cleistogenes (Orininae, Cynodonteae, Chloridoideae, Poaceae) is an ecologically important genus. The phylogenetic placement of Cleistogenes and phylogenetic relationships among Cleistogenes taxa remain controversial for a long time. To resolve the intra- and inter-generic relationships of Cleistogenes, the plastomes of 12 Cleistogenes taxa (including 8 species and 4 varieties), one Orinus species, 15 Triodia species, two Tripogon species, and two Aeluropus species were included in the present study. All the taxa showed a similar pattern in plastome structure, gene order, gene content, and IR boundaries. The number of simple sequence repeats ranged from 145 (O. kokonorica) to 161 (T. plurinervata and T. schinzii). Moreover, 1,687 repeats were identified in these taxa, including 1,012 forward, 650 palindromic, 24 reverse, and one complement. Codon usage analysis revealed that these plastomes contained 16,633 (T. stipoides) to 16,678 (T. tomentosa) codons. Sequence divergence analysis among Cleistogenes and closely related genera identified five non-coding regions (trnS-UGA-psbZ, rpl32-trnL-UAG, trnQ-UUG-psbK, trnD-GUC-psbM, trnT-GGU-trnE-UUC). Phylogenetic analysis of complete plastomes indicated that Cleistogenes is sister to a clade composed of Orinus and Triodia, whereas it did not support the sister relationship between the recently proposed subtribe Orininae (Cleistogenes and Orinus) and Triodia. The subtribe Orininae was not supported by our complete plastome data. The split between Cleistogenes and Orinus-Triodia clade go back to 14.01 Ma. Besides, our findings suggested that C. squarrosa and C. songorica are the successive early diverging groups in the phylogenetic analysis. The other 10 taxa are divided into two groups: a monophyletic group composed of Cleistogenes sp. nov. and C. caespitosa var. ramosa is sister to other eight Cleistogenes taxa. Cleistogenes was estimated to have experienced rapid divergence within a short period, which could be a major obstacle in resolving phylogenetic relationships within Cleistogenes. Collectively, our results provided valuable insights into the phylogenetic study of grass species.


INTRODUCTION
As an ecologically important genus in the grass family (Poaceae) (Liang et al., 2002;Wang et al., 2003), Cleistogenes is composed of about 13 species, which are mainly distributed in South Europe to Turkey and eastward through central Asia, China, Pakistan, Northwest India, and Japan (Chen, 1990;Chen et al., 2006;Clayton et al., 2006). A large proportion of such species are found in semi-arid regions, where they are excellent forage grass and sand-binding grass (Shan et al., 2010;Shi, 2011;Tao et al., 2017). It is well documented that Cleistogenes has the C 4 photosynthetic pathway and can adapt to a dry environment (Redmann et al., 1995;Su et al., 2011). Two Cleistogenes species including C. squarrosa and C. songorica are commonly used as plant materials to study grass competition and drought stress (Gao et al., 2005;Niu and Nan, 2017). This genus is remarkable for the cleistogamous spikelets, which often exist in the leaf sheaths to ensure the seed production even under severe conditions. These cleistogamous spikelets usually have fewer florets, which are smaller and narrower lemmas with longer awns compared with the chasmogamous spikelets (Chen et al., 2006).
Several phylogenetic studies have studied Cleistogenes and its closely related genera, while the inter-generic relationships remain uncertain. Cleistogenes is either classified within Eragrostideae or Cynodonteae based on morphological characters in traditional classification (Tutin and Al, 1980;Tzvelev, 1989;Chen, 1990;Chen et al., 2006). These two tribes are belonging to the subfamily Chloridoideae. Cleistogenes has been suggested to be within Cynodonteae in molecular phylogenetic analyses Soreng et al., 2015Soreng et al., , 2017. Cleistogenes is considered as an incertae sedis genus due to its different placement in ITS, plastid, and combined trees Soreng et al., 2015). Recent studies show that Aeluropodinae (Aeluropus and Odyssea), Triodiinae (Triodia), and Orininae (Cleistogenes and Orinus) form a monophyletic clade Soreng et al., 2017), while the relationships among them are incongruence. Cleistogenes is sister to Orinus in an ITS tree, and in a combined tree (Peterson et al., , 2016. The clade composed of Cleistogenes and Orinus is sister to the clade consisting of Aeluropodinae and Triodiinae in the ITS tree  and a successive sister to Aeluropodinae and Triodiinae in the combined tree (Peterson et al., 2016). Cleistogenes and Orinus are difficult to seperate morphologically, and Orinus kokonorica was once described as a Cleistogenes species (Hao, 1938). Cleistogamous spikelets are morphological characters used to distinguish Cleistogenes from Orinus (Chen, 1990;Chen et al., 2006). Based on these studies, Peterson et al. (2016) have proposed a new subtribe Orininae, including Cleistogenes and Orinus. However, inter-generic relationships of Cleistogenes and closely related genera remain poorly resolved in the plastid tree . The phylogenetic relationships among Cleistogenes and closely related genera are controversial. Moreover, the sister relationship between Cleistogenes and Orinus was merely moderately supported. Poor phylogenetic resolution is likely due to the lack of informative sites in the selected molecular markers, which appears to be a common problem in multi-locus phylogenetic studies. A robust phylogenetic tree between Cleistogenes and closely related genera was urgently need to understand their phylogenetic relationships among Cleistogenes and closely related genera.
Only very few studies have studied the intra-generic relationships of Cleistogenes. Species delimitation is explored using morphological characters. Cleistogenes was previously considered to be invalid because the name is coincident with a technical term. Kengia was originally published to replace Cleistogenes Keng (Packer, 1960). However, Art. 20.2 of the Vienna Code retroactively permits the formation of names based on non-Latin terms, and, thereby, validates Cleistogenes and retroactively makes Kengia illegitimate and superfluous. Yu and Zhao (2005) have suggested that there are 18 species and eight varieties of Kengia in China. In Flora of China, there are 10 species and a variety, including C. squarrosa, C. songorica, C. ramiflora, C. mucronate, C. festucacea, C. caespitosa var. caespitosa, C. kitagawae, C. polyphylla, C. hackelii var. hackelii, C. hancei, and C. hackelii var. nakaii (Chen et al., 2006). Moreover, 13 records have been accepted in the plant list 1 , including C. caespitosa var. caespitosa, C. festucacea, C. gatacrei, C. hackelii var. hackelii, C. hancei, C. kitagawae, C. mucronate, C. nedoluzhkoi, C. polyphylla, C. ramiflora, C. serotine, C. songorica, and C. squarrosa. Two species including C. songorica and C. polyphylla have been shown as successive early diverging groups based on seven plastid regions and ITS regions, while the relationships among the remaining taxa remain largely unexplored (Peterson et al., 2016). Therefore, more evidence is required to clarify intrageneric relationships of Cleistogenes.
As an important organelle in plant cells, chloroplast participates in many biological processes (Leister, 2003;Lancien et al., 2006;Leister and Kleine, 2008). Chloroplast genome (plastome) of angiosperm usually has a typical quadripartite structure with two inverted repeat (IR) regions separated by a large single-copy (LSC) region and a small single-copy (SSC) region (Jansen and Ruhlman, 2012). Plastomes usually contain 110-130 genes, including ∼80 protein-coding genes (PCGs), ∼30 transfer RNA (tRNAs) genes, and four ribosomal RNA (rRNA) genes (Guisinger et al., 2010). Although plastomes are conserved in grass species, they exhibit different levels of sequence divergence in different regions of plastomes, such as IR regions, PCGs, and intergenic regions. Plastome sequences can offer valuable insights into phylogenetic relationships in plants (Dong et al., 2012;Curci et al., 2016). They have been increasingly adopted in phylogenetic studies of Poaceae species due to the rapid development of next-generation sequencing technology (Saarela and Graham, 2010;Burke et al., 2016;Duvall et al., 2016Duvall et al., , 2020Guo et al., 2019;Orton et al., 2019;Liu et al., 2020;Welker et al., 2020;Hardion et al., 2021). Plastomes have been used to study phylogenetic relationships within Eragrostis and Triodia of Chloridoideae (Anderson et al., 2019;Somaratne et al., 2019). Relationships among tribes and some genera of chloridoid grasses have been clarified based on the complete plastome (Duvall et al., 2016).

Taxon Sampling, DNA Extraction, and Sequencing
In the present study, 32 plastomes representing 32 taxa were included in the phylogenetic analysis. Among them, 15 taxa including two outgroup taxa (Tripogon bromoides and T. chinensis) were newly sequenced, and the other 17 taxa were downloaded from GenBank [Aeluropus lagopoides (NC_042858), Orinus kokonorica (NC_042859), Triodia basedowii (NC_042860), T. chichesterensis (NC_042861), . Fresh leaves of these 15 taxa were collected in the field, followed by the drying process with silica gel. Voucher specimens and silica-dried leaves were stored at the College of Life Sciences, Shandong Normal University (SDNU), Jinan, China. Table 1 lists the sampling locality and GenBank accession numbers of these taxa. Total genome DNA was extracted using a modified CTAB method (Doyle and Doyle, 1987). The isolated DNA was subjected to agarose gel electrophoresis to assess the DNA quality. The concentration of total DNA was determined using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, United States). A final DNA concentration of >30 ng/µL was selected for library preparation. Sequencing library was generated using NEB Next R Ultra TM DNA Library Prep Kit (Illumina, NEB, United States). Genomic DNA were fragmented by sonication to a size of 350 bp. Then DNA fragments were end-polished, A-tailed, and ligated with the full-length adapter for Illumina sequencing, followed by further PCR amplification. DNA libraries were sequenced through the Illumina NovaSeq 6000 platform and 150 bp paired-end reads were generated at Novogene (Beijing, China).

Plastome Assembly and Annotation
After sequencing data were obtained, plastomes were assembled as described in Qu et al. (2019a). Both the Organelle Genome Assembler (OGA) pipeline and Spades v3.13.0 (Bankevich et al., 2012) were used in plastome assembly. The latter one was with the "careful"-option and k-mers of 61, 81, 101, and 121. Annotation was carried out using PGA (Qu et al., 2019b). A manual correction was conducted using Geneious v9.1.4. Plastome sequences downloaded from the NCBI database were re-annotated using PGA, followed by manual correction. The circular maps for newly sequenced plastomes were generated using the OGDRAW tool (Stephan et al., 2019).

Codon Usage Analysis
Referring to the unequal use of synonymous codons in an organism, codon usage can be used to determine the gene expression level, etc. Codon usage analysis of 32 taxa was performed on PCGs. The PCGs shorter than 300 bp were removed to avoid sampling bias (Wright, 1990;Rosenberg et al., 2003). The frequency of the nucleotides G + C at the third position (GC3s), the frequency of each base A, T, G, and C at the third position of codons (A3s, T3s, G3s, and C3s), the relative values of synonymous codon usage (RSCU) and GC content were determined using CodonW v1.4.2 3 . RSCU is a simple measure of non-uniform usage of synonymous codons in a coding sequence.
If RSCU value of a codon >1, that codon is frequently used than expected whereas RSCU value < 1, means that the codon is less frequently used than expected. If RSCU equals 1, it means that the codon is used randomly and equally with other synonymous codons (Sharp et al., 1986).

Phylogenetic Analysis and Divergence Time Estimations
One copy of IR was removed before the phylogenetic analysis. All 32 taxa were used for phylogenetic analyses of complete plastomes, PCGs, and NCRs. Sequences were aligned using MAFFT v7.313 (Katoh and Standley, 2013). The maximum likelihood (ML) tree was constructed using RAxML v8.2.10 (Alexandros, 2014) with 1,000 bootstrap replicates and the GTR + I + G model. The best model of evolution was determined using jModelTest2 (Darriba et al., 2012). The Bayesian inference (BI) tree was constructed using MrBayes v3.2.6 with The Markov Chain Monte Carlo (MCMC) was run for 1,000,000 steps with a random starting tree, birth-death default priors, and sampled one tree every 1,000 steps. We discarded the first 25% steps as burn-in. A relaxed clock method and penalized likelihood were involved in dating analyses using treePL (Sanderson, 2002;Smith and O'Meara, 2012). We generated 1,000 ML bootstrap trees with branch lengths by using RAxML. The minimum age of the Tripogon crown node was constrained at 20 Ma (Anderson et al., 2019). The maximum age of Aeluropus, Cleistogenes, Orinus, Triodia, and Tripogon crown node was assigned as 25 million years ago (Ma) (Anderson et al., 2019). The minimum age of Aeluropus, Cleistogenes, Orinus, Triodia, and Tripogon crown was assigned as 7.9 Ma (Anderson et al., 2019). The minimum and maximum age for the internal nodes were calculated from dating 1,000 bootstrap trees by using treePL and TreeAnnotator v1.8.495 (Drummond et al., 2012).

Repeat Sequence Analysis
A total of 1,687 repeats were identified in these 32 plastomes by Reputer, including 1,012 forward, 650 palindromic, 24 reverse, and one complement. The sole complement repeat was only identified in C. chinensis. The number of repeats ranged from 45 (T. stipoides) to 76 (T. lanigera). Most taxa had 45-60 pairs of repeats. The repeat lengths in 30-39 bp were most common, accounting for 64.1%. Repeats with the lowest proportion (7.1%) ranged in size from 50 to 59 bp. Repeat lengths in 50-59 bp were not identified in Cleistogenes taxa. Figure 2 shows that the number of SSRs for all taxa ranged from 145 (Orinus kokonorica) to 161 (T. plurinervata and T. schinzii). The number of SSRs found in Cleistogenes taxa ranged from 155 to 159. There were six categories of SSRs detected in these taxa, including mono-nucleotide, di-nucleotide, tri-nucleotide, tetra-nucleotide, penta-nucleotide, and hexanucleotide repeats. In all taxa, mono-nucleotides were the most common SSRs, accounting for 74.4%, followed by di-nucleotide (19.1%). There were some types of SSRs found in all plastomes, including mono-(A/C/G/T), di-(AG/AT/CT/GA/TA/TC), tri-(AAT), and tetra-(AACG/AATA/GTAG/YCGT/TTCT) (Supplementary Table 2). Penta-nucleotide repeats were only found in T. longiceps and T. plurinervata, and hexa-nucleotide repeats (AAATAT) were only found in Cleistogenes sp. nov. and C. caespitosa var. ramosa.   Table 3). The GC contents of the entire plastomes, the frequency of the nucleotides G + C at the third position (GC3s), and the frequency of each base A, T, G, and C at the third position of codons (A3s, T3s, G3s, and C3s) were similar in all taxa. The GC content was approximately 39%, suggesting an abundance of AT. A3s and T3s were larger than G3s and C3s. Leucine was the most frequent amino acid in the plastomes (∼10.8%), while cysteine was the least frequent amino acid (∼1.0%) (Supplementary Table 4). The relative values of synonymous codon usage (RSCU) were determined to estimate the preference for the use of synonymous codons. In all taxa, nearly all the amino acid codons had a bias (RSCU > 1 or RSCU < 1) except for methionine (AUG, RSCU = 1) and tryptophan (UGG, RSCU = 1).

IR Expansion and Contraction of Ael-Cle-Ori-Trio Clade Plastomes
The IR/SSC and IR/LSC boundary regions of these taxa were very conserved (Figure 3). The junctions of LSC/IRb and IRa/LSC were located in intergenic regions. LSC/IRb was located between rpl22 and rps19. The fragment size of rpl22-rps19 located in the IRb region was 34 bp in O. kokonorica, 36 bp in T. plurinervata,
The mVISTA results showed that all aligned plastome sequences were highly similar, and the NCRs showed a higher divergence compared with the PCGs (Figure 5).

Phylogenetic Analysis and Divergence Time Estimations
Maximum likelihood and Bayesian inference phylogenetic analyses using three datasets (complete plastomes, PCGs, and NCRs) generated identical topologies for Cleistogenes and its closely related genera in the present study (Figure 6,  Supplementary Figures 1-6). Aeluropus lagopoides and A. sinensis form a monophyletic clade with a strong support (BS = 100, PP = 1.00). Cleistogenes is strongly supported as a monophyletic group in ML and BI analysis (BS = 100, PP = 1.00).

Morphological Comparison of O. kokonorica and C. songorica
Spikelet, lemma and underground part of O. kokonorica and C. songorica were compared in the present study (Figure 8). Spikelets of O. kokonorica and C. songorica were laterally compressed, fewer florets in O. kokonorica. O. kokonorica is charactered by lemma thin, dorsally black-brown but yellowbrown at base and apex, with loosely pilose at margins and lower keel. Lemma of C. songorica was wider than O. kokonorica and glabrous. Underground part of these two taxa were significant different. Long scaly rhizomes can be observed in this part of O. kokonorica, while they were not present in C. songorica.

Plastome Evolution
Plastomes of land plants have greatly conserved genome size, gene content, gene order, and organization (Palmer, 1991;Bock, 2007). Previous studies have assessed the organization and evolution of Poaceae plastomes (Katayama and Ogihara, 1993;Guisinger et al., 2010). In the present study, the plastomes had high structural similarity among all taxa, including genome size, gene content, and gene order. The plastomes of these taxa had an average genome size of approximately 135 kb. Plastomes of these taxa showed the typical quadripartite structure as previously reported Poaceae species, consisting of an LSC region and an SSC region, and separated by a pair of IR regions (Rousseau-Gueutin et al., 2015;Huang et al., 2017). Each plastome contained the same number of genes (86 PCGs, 38 tRNAs, and eight rRNAs), with a similar GC content and conserved intron positions. There were some unusual features of Poaceae plastomes compared with other angiosperm species, including three gene loss (accD, ycf1, and ycf2) and intron loss from two genes (clpP and rpoC1) (Katayama and Ogihara, 1996;Guisinger et al., 2010).
In the present study, accD and ycf1 were not annotated in all taxa. These two genes might be useful molecular markers for phylogenetic analysis of land plants, playing essential roles in leaf development (Kode et al., 2005) and plant viability (Kikuchi et al., 2013). In the case of ycf2 gene loss, different lengths of ycf2 were found, suggesting a progressive degradation. The ycf2 gene retained segments in different lengths, ranging from 105 to 792 bp. The ycf2 gene plays a significant undetermined function in higher plants (Drescher et al., 2000). This gene can provide effective variable sites for reconstructing a generally wellsupported phylogeny (Huang et al., 2010). Intron losses of clpP and rpoC1 were detected in all 32 taxa. Gene and intron losses can lead to a decrease in grass plastome size . The expansion and contraction of the IR region have been proposed as an important source of length variation in plastomes (Palmer et al., 1987;Aii et al., 1997;Liu et al., 2018). However, in our present study, all taxa exhibited a highly conserved pattern of IR boundaries with only slight structural variations. In Poaceae, ndhH and ndhF are located near opposite ends of the SSC region and can extend into IR regions (Maier et al., 1990;Davis and Soreng, 2010). In the early diverged grasses and closest relatives, ndhH extends 175 to 200 nucleotides into the IR region, while ndhF is confined to the SSC region (Davis and Soreng, 2010). Davis and Soreng (2010) have reported that a part of ndhH and ndhF less than 30 nucleotides migrates into IR regions in PACMAD clade of the Poaceae family. A similar pattern has been observed in Eragrostideae species except for Eragrostis tenellula (Somaratne et al., 2019). In our present study, this phenomenon was observed in most taxa, except for Triodia stipoides, in which ndhF extended 34 nucleotides into the IRb region. The different migrations of gene termini relative to the two SSC-IR junctions are also of significance in terms of the gene overlap (Davis and Soreng, 2010).
As tandemly repetitive DNA sequences, SSRs show high levels of polymorphism. They are abundant across plastomes. As they are reliable and highly polymorphic, they have been increasingly used in phylogenetics, species identification and population genetic studies (Kalia et al., 2011). Similar to the results of other Poaceae species, the predominant type of SSRs is mononucleotides, of which A or T repeats account for the majority (Tanaka et al., 2016;Somaratne et al., 2019). The abundance and distribution pattern of penta-nucleotide and hexa-nucleotide repeats were different in the present study. Penta-nucleotide repeats and hexa-nucleotide repeats are the most common types among all the Oryza species (Tripathy et al., 2019), while they were species-specific in the present study. Penta-nucleotide repeats were only found in T. longiceps and T. plurinervata, and hexa-nucleotide repeats (AAATAT) were only detected in Cleistogenes sp. nov. and C. caespitosa var. ramosa. In the present study, penta-nucleotide repeats were distributed in both coding regions and NCRs, while hexa-nucleotide repeats were distributed in NCRs. Repeats identified by Reputer were assigned into four categories, including forward, reverse, complement, and palindromic. Most of the repeats were forward and reverse repeats. All the identified repeats in the present study might be useful in future studies on population genetics of these 32 taxa.
Synonymous codon usage bias represents the differences in the relative frequency of synonymous codons for individual amino acids in protein-coding sequences, and many biological factors, such as gene expression, amino acid composition of protein, and protein structure, are correlated to such bias (Angellotti et al., 2007). Base composition is a balance between G + C nucleotide pairs and mutational pressure, which is a pervasive effect on codon usage (Sueoka, 1962). In the present study, the GC content was highly conserved, suggesting that these plastomes achieved evolutionary equilibrium (Zhang Y. R. et al., 2012). Our findings were consistent with previously reported Poaceae plastomes that all the 32 plastomes had similar codon usage patterns, and the ATrich bias was the strongest in the third codon position (Somaratne et al., 2019). Preferred codons could be useful in designing degenerate primers, the introduction of point mutation, and evolutionary studies. In the standard genetic codes, there are 64 codons, 61 of which encode 20 different amino acids, and 3 of which are stop codons. All amino acids are encoded by two to six synonymous codons except for Met and Trp. Among the stop codons (UAA, UAG, and UGA), the predominantly used one was UAA (RSCU > 1.63) in the present study. Similar favorite stop codons have been found in some Poaceae species, e.g., Zizania latifolia (Zhang et al., 2016) and Pennisetum glaucum (Raveendar et al., 2019).

Phylogenetically Informative Markers
Plastid regions are widely used in the phylogenetic study of Poaceae (Johnson, 1997;Bouchenak-Khelladi et al., 2008;Rua et al., 2010;Barbera et al., 2020). Previous studies of Cleistogenes and its closely related genera mainly include 11 plastid regions (ndhA, ndhA intron, ndhF, rps16-trnK, rps16 intron, rps3, rpl32-trnL, rps16, rpoC2, ccsA, and matK), while these regions are not effective in resolving the intra-and intergeneric relationships of Cleistogenes. Some of them (e.g., rps16, rpoC2, ndhA intron, and rps16 intron) had a low percentage of parsim-info sites (Figure 4). Therefore, phylogenetic trees based on these regions generated low resolution and different topologies. In the present study, rps3 was the most divergent PCG, while Pi% was only 4.44%. Insufficient sequence polymorphism of the markers hindered their use for phylogenetic study. Therefore, it is important to identify more informative markers to study the intra-and inter-genus phylogeny for Cleistogenes.
In Poaceae plastomes, PCGs and IR regions are more conserved at the sequence divergence level compared with NCRs and single-copy (SC) regions, respectively (Bortiri et al., 2008;Qiu et al., 2019;Somaratne et al., 2019). In the present study, the whole plastome comparative analysis showed that IR regions were more conserved compared with the SC regions, and comparison of the percentage of parsim-info sites confirmed that NCRs had higher Pi% than PCGs. NCRs evolve via accumulating nucleotide substitutions and microstructural mutations (insertions, deletions, and inversions) (Yamane et al., 2006). Currently, NCRs of plastomes have become an important tool in studying phylogenetic relationships in Poaceae (Liu et al., 2008;Oliveira et al., 2014). In the present study, five regions displayed remarkably higher Pi% (>5%), including four regions (trnS-UGA-psbZ, trnQ-UUG-psbK, trnD-GUC-psbM, trnT-GGU-trnE-UUC) located in the LSC region and one region (rpl32-trnL-UAG) located in the SSC region. These regions might be potential molecular markers in the phylogenetic analysis of Cleistogenes and its related genera. Among the five phylogenetically informative markers, two regions (rpl32-trnL-UAG and trnT-trnE) have been reported as highly variable markers to study phylogenetic relationships of Poaceae species (Rousseau-Gueutin et al., 2015;Qiu et al., 2019).

Whether Orininae Should Be Established?
The monophyly of the Ael-Cle-Ori-Trio clade has been strongly supported in previous studies (Soreng et al., 2015(Soreng et al., , 2017. However, there are some incongruences about the phylogenetic relationships among these four genera (Chen et al., 2006;Shan et al., 2010;Soreng et al., 2015Soreng et al., , 2017Peterson et al., 2016). In the present study, the analysis of complete plastomes strongly supported the highly congruent phylogenetic topology for Ael-Cle-Ori-Trio clade (Figure 6). Aeluropus species including A. lagopoides and A. sinensis form a monophyletic clade, which is located at the basal position in the tree. The Cle-Ori-Trio clade is detected with a high bootstrap value. Cleistogenes is strongly supported as a monophyletic group and sister to the clade composed of Orinus and Triodia. The position of Cleistogenes recovered here is incongruent with other studies. The placement of Cleistogenes is not resolved in the plastid tree based on seven plastid regions. However, Cleistogenes is sister to Orinus in the ITS tree  and combined tree (Peterson et al., , 2016. Subtribe Orininae has been newly proposed in a phylogenetic analysis of 213 species based on seven plastid regions and ITS (Peterson et al., 2016). The present study showed that Cleistogenes is sister to the clade composed of Orinus and Triodia, which indeed did not support the establishment of Orininae (Peterson et al., 2016). These three genera are different in morphology. Inflorescences of Cleistogenes and Orinus are sparse panicle (Chen et al., 2006). Inflorescence of Triodia is usually a panicle of solitary spikelets, and sometimes it is a spike or raceme (Lazarides, 1997). Cleistogenes is a genus of Eurasian flowering plants in the grass family. It is remarkable for cleistogamous spikelets hidden in the upper leaf sheaths (Chen et al., 2006), while cleistogamous spikelets are not found in Orinus and Triodia. Cleistogamy flowering assures plant reproduction under variable environmental conditions, and its development is known to be affected by drought, chilling, salinity, and light (Morinaga et al., 2008). Orinus is a genus of Asian plants in the grass family (Chen et al., 2006). Triodia is a large genus of hummock-forming grass endemic to Australia (Lazarides, 1997;Gamage et al., 2012). Eurasia and Australia have been separated by oceans at ca. 90 Ma (Davis et al., 2002), earlier than the origin of Cle-Ori-Trio clade. Therefore, long-distance dispersal must have played a dominant role in the transoceanic distribution of this clade. The split between these three genera go back to 14.01 Ma (Miocene). The Malay Archipelago probably facilitated biotic dispersal between Asia and Australia during the Miocene (Welzen et al., 2005). The habitat conditions of Orinus and Triodia were similar. They all distribute in arid regions on sandy or stony soils. Morphologically, leaf blades of Orinus and Triodia species are linear to involute and rigidly straight. Leaf involution is often hypothesized to confer survival value in xeric habitats by reducing stomatal and cuticular conductance and postponing desiccation. We compared the morphology of O. kokonorica and C. songorica (Figure 8). Lemmas of Cleistogenes are almost glabrous. Lemmas of Orinus are hairy all over or only on margins. Lemmas of Triodia species are entire or minutely to deeply 2-or 3-lobed and hairy. Trichomes (also called leaf hairs) are specialized cell types in the epidermal layer, that enhance the protection of plant tissues from the external factors both mechanically and chemically (Hauser, 2014). Long scaly rhizomes are recognized as an important character of Orinus (Chen, 1990;Chen et al., 2006). The scaly rhizomes allow the plants to survive unfavorable dry spells or seasons in dormancy and to re-sprout from the axils of the scale-like leaves when conditions improved. Orinus-Triodia and Cleistogenes may represent two different ways to adapt to extreme environments. Meanwhile, there are some similar micromorphological characters recognized in Orinus and Triodia. Papillae occurring on long cells and dumbbell-shaped silica cells are observed as two similar lemma micromorphological characters in Triodia and Orinus, while papillae and silica cells are absent in Cleistogenes taxa (Liu et al., 2010). From the above-mentioned results, some conflicts existed between plastome-and ITS-based trees. These conflicts could mostly result from rapid divergence, hybridization and incomplete lineage sorting. The Cle-Ori-Trio clade was estimated to have experienced rapid divergence in a short period (12.54 to 14.01 Ma), which could be an obstacle in resolving phylogenetic relationships within the clade. Natural hybridization and incomplete lineage sorting are very common in numerous plants (Ellstrand et al., 1996;Maddison, 1997). Hybridization would lead to alleles with mixed histories, therefore resulting in high character conflict (Folk et al., 2016). Incomplete lineage sorting may be defined as the failure for two or more allelic lineages to coalesce within a population (Degnan and Rosenberg, 2009). Hybridization and incomplete lineage sorting can lead to discord in phylogenetic analyses (Kubatko and Degnan, 2007). In conclusion, the establishment of Orininae was not supported by plastomes data.

Phylogenetic Relationship Within Cleistogenes
Cleistogenes is a monophyletic genus in the tree based on complete plastomes, which reinforce the previous study using ITS and plastid sequences (Peterson et al., 2016). Cleistogenes began to diversify in Pliocene or Pleistocene (1.34-5.2 Ma). It is suggested that there is a climate shift in the Pliocene across large area of East Asia, as well as an even colder and drier Quaternary compared with the late Pliocene . As a kind of cool-and dry-adapted grass, Cleistogenes was probably diversified caused by climate change. The phylogenetic relationship among Cleistogenes taxa have been rarely studied. The lack of a well-supported intrageneric phylogenetic framework has posed a problem for the classification of Cleistogenes. Phylogenetic relationships within Cleistogenes recovered here are not congruent with the previous study. It is shown that C. songorica and C. polyphylla are early diverging groups in the phylogenetic analysis of Cynodonteae based on seven plastid regions and ITS regions (Peterson et al., 2016). In the present study, C. songorica and C. squarrosa are resolved as a highly supported successive sister group to the remaining 10 taxa. It is reported that C. squarrosa is a dominant grass species in grasslands (Wang et al., 2000;Liang et al., 2002). Dense foliage is on the low mounds of the culm, and upper internodes of the dried culm are curling and elongate. Awned spikelets are always with few flowers. These above-mentioned characteristics distinguish it from the remaining taxa of Cleistogenes (Chen et al., 2006). C. songorica is a dominant species of desert grassland, and it has excellent drought tolerance Zhang et al., 2014;Muvunyi et al., 2018). This species has dense tufts of graygreen leaves. Spikelets are purple and awnless, usually with broader lemmas. These two species including C. songorica and C. squarrosa can adapt to dry environment. In the complete plastome phylogenetic tree, C. caespitosa var. ramosa was sister to Cleistogenes sp. nov., while C. caespitosa var. caespitosa clusters with the other seven Cleistogenes taxa. The difference between this variety and its original variety was that the upper stem was much branched, and the panicle was extremely narrow. Similar morphological characters of C. caespitosa var. ramose and C. caespitosa var. caespitosa could be attributed to parallel evolution. Palea of Cleistogenes sp. nov. significantly shorter compared with the other Cleistogenes taxa. The clade composed of C. hackelii var. hackelii, C. polyphylla, C. festucacea, C. hancei, C. hackelii var. nakaii, C. chinensis, C. caespitosa var. caespitosa, and C. gracilis is found with short internal branch length. The short branches most likely represent a case of rapid radiation with lineage-specific heterogeneity in rates of substitution. Cleistogenes was estimated to have experienced rapid divergence within a short period, which could be a main reason for short internal branches between these 10 Cleistogenes taxa in phylogenetic tree. Lemma micromorphological characters of C. caespitosa var. caespitosa, C. chinensis, and C. hackelii var.
hackelii are highly conserved (Liu et al., 2010). Straight outline of intercostal long cells, semi-circle-shaped cork cell, triangularshaped stomata subsidiary cell, and papillate-based macro hair have been observed in C. caespitosa var. caespitosa, C. chinensis, and C. hackelii var. hackelii. In the present study. the phylogeny of complete plastomes revealed the intra-generic framework of Cleistogenes. However, future studies with a larger sample size at the population level or the next-generation sequencing approach may help resolve their phylogenetic relationships.

CONCLUSION
Collectively, we sequenced, assembled, and annotated 15 plastomes of Cleistogenes and its closely related genera. Plastomes of 12 Cleistogenes taxa (including 8 species and 4 varieties), 15 Triodia species, one Orinus species, two Tripogon species, and two Aeluropus species were included in the present study. All the plastomes were highly conserved in plastome structure, gene content, gene order, and IR boundaries. The type and number of repeat sequences and SSRs were calculated and could be used in studies of phylogenetics, species identification, and population genetics. By examining variable sites and parsim-info sites, five highly variable regions were identified, which could be used for phylogenetics, population genetics, and biogeography of Cleistogenes and its closely related genera. Phylogenetic analysis revealed that Cleistogenes is sister to a clade composed of Orinus and Triodia, whereas recently proposed Orininae (Cleistogenes and Orinus) is not supported in the present study. It is suggested that C. squarrosa and C. songorica are the successive early diverging groups. The other 10 taxa are divided into two groups. Cleistogenes sp. nov. and C. caespitosa var. ramosa form a monophyletic group sister to the other eight Cleistogenes taxa. Our current findings provided valuable insights into the phylogenetic study of grass species. Moreover, more nuclear data will be necessary for further study on the intra-and inter-generic relationships of Cleistogenes.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.

AUTHOR CONTRIBUTIONS
X-JZ, X-JQ, and S-JF designed the experiments. RW and KL carried out the experiment and analyzed the data. RW wrote the first draft of the manuscript. X-JZ, W-LC, X-JQ, and S-JF supervised and completed the writing. All authors collected the samples, revised, and approved the final manuscript.