The materials of Coleanthus subtilis were collected from Harbin, China, in June 2021, and subsequently deposited in the Herbarium of the Wuhan Botanical Garden (HIB), Chinese Academy of Sciences (China), with herbarium number ZXX21129. For drying and long-term preservation of molecular samples, fresh leaves were preserved in silica gel (Chase and Hills, 2019). The complete genomic DNA of C. subtilis chloroplast was extracted using a modified CTAB procedure (Allen et al., 2006) and then sequenced at Novogene Co., Ltd. (Beijing, China) with Illumina paired-end technology platform. Purified high-quality genomic DNA was broken into short fragments of approximately 350 bp, and paired-end (PE) libraries were constructed by adding A-tails, PCR amplification and other steps, followed by sequencing in 150 bp paired-end mode on an Illumina HiSeq 2500 platform. The final number of raw reads obtained was 36,062,743 and that of clean reads after filtering was 35,335,540. The raw data has been uploaded to the NCBI database (BioProject ID: PRJNA802068).

Get Organelle v1.7.4 (Jin et al., 2020) was used to assemble the chloroplast genome with default parameters. The low-quality reads and adapters were first filtered, then a de novo assembly performed, and the results were further purified to generate the complete chloroplast genomes. The results were visualized with Bandage (Wick et al., 2015). The Plastid Genome Annotator (PGA) software (Qu et al., 2019) was used to perform the annotation of the entire chloroplast genome, and in addition to using Amborella trichopoda as the reference genome, some Poaceae species were also selected to enhance the credibility of the annotation results. Furthermore, to ensure the accuracy of the annotation results, the genome was also annotated simultaneously with the help of GeSeq online tool² (Tillich et al., 2017).

The check of annotated genes was implemented in the software Geneious-v10.2.3 (Kearse et al., 2012), which was used to further verify and refine the annotation results and to manually correct errors detected in gene annotation. Special attention was paid to some genes located at the boundaries and the highly variable genes, such as ndhF, ndhK, ycf2, accD, etc. The circular chloroplast genome map of Coleanthus subtilis was drawn and visualized using OGDraw online tool³ (Greiner et al., 2019). Lastly, the annotated sequence was submitted to GenBank on the NCBI website, with an accession number OL692806.

The chloroplast genome characteristics of Coleanthus subtilis were analyzed in Geneious-v10.2.3 software by comparing chloroplast genomes with those of Poaceae species downloaded from the NCBI database (Supplementary Table 1). A total of 24 species representing 10 subtribes (5 tribes) were used for the comparative analysis of chloroplast genomes. Additionally, to determine genomic divergence among these species, genomic similarity analysis was performed using the Glocal alignment program (shuffle-LAGAN mode) in mVISTA (Brudno et al., 2003; Frazer et al., 2004) with C. subtilis as the reference. The SC/IR boundary analysis was done using the IRscope (Amiryousefi et al., 2018) to observe the contraction and/or the expansion of the genes at the borders. For the codon usage bias analysis, MEGA 7.0 software (Kumar et al., 2016) was chosen to calculate relative synonymous codon usage (RSCU) values based on the coding sequences (CDS regions).

The REPuter tool⁴ (Kurtz et al., 2001) was used to identify repeats including forward, reverse, palindrome, and complement sequences. When the Hamming distance is equal to 3, the length and identity of repeats are limited to ≥30 bp and >90%, respectively. The simple sequence repeats (SSRs) were analyzed using the MISA (Beier et al., 2017) with the basic repeat setting: a threshold of 10, 5, 4, 3, 3, and 3 for mono-, di-, tri-, tetra-, penta-, and hexa-nucleotides, respectively. The DnaSP-v5.10 software (Librado and Rozas, 2009) was used to calculate nucleotide variability (Pi) values and variable sites using the aligned chloroplast genome sequences with a window length of 600 bp and a step size of 200 bp.

The EasyCodeML program in PAML package (Gao et al., 2019) was utilized to identify positive sites in protein-coding genes to quantify selection pressure. This software provided four site models (M0 vs. M3, M1a vs. M2a, M7 vs. M8, and M8a vs. M8), Bayes Empirical Bayes (BEB) analysis (Yang et al., 2005) and Naive Empirical Bayes (NEB) analysis were performed in each model to measure the loci with positive selection pressure.

To understand the phylogenetic position of Coleanthus subtilis in the family Poaceae and its affinities with other species, a phylogenetic tree was reconstructed using the Maximum Likelihood (Felsenstein, 1981) and Bayesian Inference analysis (Huelsenbeck et al., 2001). This was based on 76 shared protein-coding genes of the Chloroplast genome from a total of 53 species from 26 genera in Poaceae, with Acidosasa purpurea as the outgroup (Supplementary Table 1). Each protein-coding sequence was first aligned in the software MAFFT-v7.409 (Katoh and Standley, 2013), followed by removing the stop codon and discarding the bad fragment with the Gblock program (Talavera and Castresana, 2007) and later concatenated using the concatenated in-built PhyloSuite program (Zhang et al., 2020). ML analysis in IQ-tree and BI analysis in MrBayes were used to infer phylogenetic relationships. The best-fit models for each of the two analyses were found in Model Finder (Kalyaanamoorthy et al., 2017) according to the Bayesian Information Criterion (BIC), and the most suitable model for Bayesian analysis was detected as GTR + F + I + G4, while GTR + F + R3 was used for the Maximum Likelihood analysis. Subsequently, the BI tree was constructed by the software MrBayes-3.2.6 (Ronquist et al., 2012) for 1,000,000 generations, sampling every 1000 generations, and the software IQ–TREE was implemented to construct the ML tree with bootstrap replications of 1000 (Lam-Tung et al., 2015). The phylogenetic trees were visualized in the software Figtree-v1.4.4⁵. Both phylogenetic trees were combined manually using AI software based on consistent topological structures. The results were imported into the software Figtree-v1.4.4 to view the generated phylogenetic trees and to enhance their visualization. Considering the consistent topology, the phylogenetic trees constructed by both methods were manually combined in the AI software.

The chloroplast genome of Coleanthus subtilis is 135,915 bp in size and consists of four regions that together form a loop structure. These four regions are the large single copy region (LSC) of 80,100 bp, a small single copy region (SSC) of 12,757 bp, and two inverted repeat regions (IR) of 21,529 bp in length, respectively. In addition, a pair of inverted repeat regions separate the two single-copy regions (Figure 1 and Table 1). GC content varies in different regions of the chloroplast genome. The highest GC content of 43.9% was found in the IR regions of C. subtilis, while the two single copy regions had 36.3% (LSC) and 32.4% (SSC) (Table 1).

A total of 129 genes were annotated in the chloroplast genome of C. subtilis, with 83 protein-coding genes (PCGs), 38 tRNA genes, and 8 rRNA genes. In addition, the accD gene was found missing in the chloroplast genome, while ycf1, ycf2, ycf15, and ycf68 were pseudogenes (Table 2). These genes were divided into three groups based on their different functions. Nineteen genes were observed to replicate in the inverted repeat regions, seven of which were PCGs (ndhB, rpl2, rpl23, rps7, rps12, rps15, rps19), eight were tRNA genes (trnA-UGC, trnI-CAU, trnI-GAU, trnH-GUG, trnL-CAA, trnN-GUU, trnR-ACG, and trnV-GAC) and the remaining four genes were rRNA (rrn4.5, rrn5, rrn16, and rrn23). In addition, the largest number of genes in the LSC region was 82, while only 11 genes were located in the SSC region. More interestingly, all rRNA genes were distributed in the IR regions (Supplementary Table 2). We identified 15 genes containing one intron in C. subtilis, with six being tRNAs and nine being PCGs. It’s important to highlight that trnK-UUU had the longest intron with 2480 bp, which completely wrapped the MatK gene. Meanwhile, the ycf3 gene contained two codons with lengths of 774 bp and 726 bp (Supplementary Table 3).

The chloroplast genome of C. subtilis showed high similarities with other Poaceae species in terms of genome length and structure, GC content and gene number. The complete genomes length varied from 133608 bp (Colpodium humile) to 137370 bp (Stipa purpurea), LSC from 78636 bp (Colpodium humile) to 81252 bp (Brachypodium stacei), SSC from 12390 bp (Zingeria biebersteiniana) to 12842 bp (Stipa purpurea), and IR from 20831 bp (Melica mutica) to 22917 bp (Briza maxima) (Table 1). The overall GC content was around 38.5%, and each of the four regions also differed only insignificantly. In particular, gene number and composition were almost identical in 24 species, with only the trnL-UAA gene missing in Bromus vulgaris. In addition, no structural rearrangements were found in any of them (Figure 2).

To observe the variation of IR boundaries, we did a comparative analysis of the junction structure based on the chloroplast genomes of Coleanthus subtilis and 23 other Poaceae species (Figure 3). The results showed that their boundary features were similar, the genes found at the nodes were mainly rpl22, rps19, rps15, ndhF, ndhH, and psbA. The rps19 and rps15 genes were replicated and fully embedded in the IR region, with lengths of 13–46 bp and 293–479 bp from the two IR/LSC boundaries, respectively. The ndhF genes were located entirely on the left of the IRb/SSC and were 27 to 122 bp from this boundary. Also, the ndhH gene occupied the IRa/SSC junction and was overwhelmingly located within the SSC region, with only a small portion of 156 to 316 bp extending into the IRa region. It should be noted that the ndhH gene of Colpodium humile was slightly shorter in length and was therefore completely encapsulated in the SSC. In addition, Brachypodium stacei and Briza maxima showed significant differences in boundary characteristics from the other species. It was clearly observed that the IR regions of Brachypodium stacei were contracted, resulting in the distribution of the rps19 gene originally located in this region to the LSC. However, the IR region of Briza maxima expanded, wrapping the rpl22 that should have been located in the LSC.

Whole sequence alignment of the chloroplast genomes of 24 Pooideae species was performed to detect the differences that exist in their structures (Figure 4). The annotation of Coleanthus subtilis were used as a reference. The chloroplast genomes of these species were largely identical in terms of the number and arrangement of genes. However, some highly variable regions were still detected, such as rbcL-psaI, psbE-petL, trnD-GUC-psbM, trnG-UCC-trnT-GGU, rpl32-trnL-UAG and other intergenic regions. Overall, the non-coding regions showed a higher potential for variation compared to the coding regions. Although the protein-coding regions were relatively conserved, larger variants were observed in the rpoC2, infA, cemA and matK genes. Besides, variations were also presented in some genes located at the IR/SC boundary, such as rps19 and ndhF. However, the rRNA and tRNA sequences were highly conserved, where genes such as rrn16, rrn23, trnV-GAC, and trnR-ACG were almost unchanged. At the same time, IR regions of these species were minimally altered and significantly more conserved than the two single-copy regions.

There were 19838 codons eventually found in chloroplast genome of Coleanthus subtilis. Methionine and Tryptophan amino acids were encoded by a single codon, AUG and CGG, respectively. However, the remaining amino acids were encoded by two to six codons and showed a clear preference for codon usage (Figure 5). The most abundant amino acid in the C. subtilis was leucine 2135 (10.76%). Conversely, the least abundant amino acid was cysteine 218, which accounted for only 1.10% of the total. Meanwhile, among the six codons encoding leucine, UUA had the highest RSCU value of 2.10, which indicated that it had a high preference and was the most commonly used codon. Interestingly, most of the codons with RSCU values greater than 1 had A/U as the terminal codon, while those with C/G as the terminal codon usually had RSCU values less than 1.

The RSCU values of the five species were compared in order to understand the differences in their codon usage (Figure 6). For one amino acid, the sum of the RSCU values of all codons involved in its encoding was almost equal. Also, the RSCU values of the same codons were almost identical in these species, indicating that their codon usage habits were more stable and hardly change (Figure 6 and Supplementary Tables 4, 5).

We detected only palindromic and forward repeats in chloroplast genomes of Coleanthus subtilis and its related species, where the proportion of forward repeats was higher than that of palindromic repeats (Figure 7A and Supplementary Table 6). Most of repeats were 30–34 bp in length and were mainly distributed in the LSC region (Figures 7B,D and Supplementary Tables 7, 8). Also, the CDS regions contained most of the repeats, followed by the IGS regions (Figure 7C and Supplementary Table 9). Some repeats were also shared between IGS, CDS, tRNA, and intron regions.

A total of 26 SSRs were detected in C. subtilis, while 28, 28, 30, and 33 microsatellites were found in Phippsia algida, Puccinellia nuttalliana, Sclerochloa dura, and Zingeria biebersteiniana, respectively (Figure 8A and Supplementary Table 10). These SSRs were classified into five types, namely Mono-, di-, tri-, tetra-, and penta-nucleotides repeats. The mono-nucleotide repeats accounted for 50.34% of the 145 microsatellites and were the most abundant SSR types in the five species, followed by tetra-nucleotide repeats (26.21%). Most of the microsatellites were distributed in the LSC region and consisted of A/T motifs (Figures 8B,C and Supplementary Tables 11, 12).

To comprehensively understand the sequence divergence of the chloroplast genomes of Coleanthus subtilis and its related species, we calculated Pi values for nucleotide diversity. Pi values fluctuated between 0 and 0.0697, with a mean value of 0.02172 (Figure 9). We identified 13 polymorphic regions (matK, trnK-UUU/rps16, rps16/trnQ-UUG, trnG-UCC/trnT-GGU, trnT-GGU/trnE-UUC, petN/trnC-GCA, trnC-GCA/rpoB, rps4/trnL-UAA, trnL-UAA/ndhJ, ndhC/trnV-UAC, ndhF, ndhF/rpl32, and ndhA) with nucleotide diversity >0.05, 10 of which were intergenic spacer regions and the remaining three were protein-coding regions. Meanwhile, no highly variable loci were detected in the IR regions and the nucleotide diversity values were significantly lower than those in the single copy regions (Figure 9 and Supplementary Table 13).

In this study, dN/dS values were calculated based on 76 CDS regions with site models in EasyCodeML. According to the M8 model, only the atpF gene possessed a significant positive site in the BEB approach (Table 3). Meanwhile, a total of 45 loci corresponding to 21 genes were identified in the NEB method, of which 13 genes (atpA, atpF, atpI, ccsA, clpP, infA, ndhA, ndhD, ndhK, rbcL, rpoA, rps16, and rps3) had a significant positive site. In addition, the ndhF, psaA, psaB, psbC, and rpoC1 genes contained two significant positive selection loci, while the cemA and matK genes were detected with three and eight loci under positive selection, respectively. Moreover, the rpoC2 gene was found to have the highest number of positive selection sites, including 11 significant positive sites.

In the current study, we utilized the protein-coding regions of chloroplast genomes for the first time to explore the phylogenetic position of Coleanthus subtilis. The topologies of the phylogenetic trees generated with maximum likelihood (ML) and Bayesian analysis (BI) were identical, with generally high branch bootstrap values and posterior probabilities. Based on consistent topologies, we showed the phylogenetic tree represented by the ML method (Figure 10). The 53 species representing 26 genera were divided into ten subtribes and six tribes. Among them, Coleanthus was placed in the big clade containing Phippsia, Puccinellia, Sclerochloa, and Zingeria, which were components of the subtribe Coleanthinae. In addition, C. subtilis formed a sister branch with the genus Phippsia, while this branch was also sister to other taxa of this subtribe (BS = 100, PP = 1). The genus Colpodium was nested in the subclade Loliinae and had a sister relationship with the genus Castellia (BS = 100, PP = 1).

The chloroplast genome of Coleanthus subtilis exhibited a tetrad structure of 135915 bp in length, which is similar to the length and structural characteristics of cp genomes of other higher plants (Jansen et al., 2005; Daniell et al., 2016). We found that the GC content in the cp genome of Pooideae species was unevenly distributed, with the IR regions having a higher GC content than the two single copy regions. This may be attributed to the fact that four rRNA genes with high GC content were located in the IR regions, which supported the speculation of previous studies (Mardanov et al., 2008; Gao et al., 2009; Wanga et al., 2021). The accD gene has been lost within the cp genomes of Pooideae species, while ycf1, ycf2, ycf15, and ycf68 were pseudogenes, which is a relatively common phenomenon in Poaceae (Huang et al., 2017). There is a correlation between gene loss and evolution, and some studies suggest that it may be an adaptive strategy with positive effects on survival and reproduction (Xu and Guo, 2020). In addition, we also found trnL-UAA gene loss in Bromus vulgaris. Pseudogenization of tRNA (trnT-GGU) has also been observed in the Asteraceae family (Abdullah et al., 2021a). Sixteen intron-containing genes were detected in 24 species in which introns of rpoC1 and clpP genes were lost. Besides, the trnK-UUU has the longest intron that completely wraps the matK gene, a result that has been reported in other studies (Li X. et al., 2019; Souza et al., 2020). The rpoC1 gene has been reported to contain introns in most land plants (Ohyama et al., 1986; Kugita et al., 2003). However, deletion of the rpoC1 intron was observed in some angiosperm lineages, such as most Poaceae and some species of the families Fabaceae, Cactaceae, and Aizoaceae (Downie et al., 1996; Wallace and Cota, 1996; Huang et al., 2017). Our study on the subfamily Pooideae further confirm that the absence of the rpoC1 intron is universal in the Poaceae. Similarly, the clpP gene usually contained two introns. Nevertheless, both introns have been lost in Pinus and some species from the genera Oenothera, Silene, and Menodora (Lee et al., 2007; Huang et al., 2017). Also, it was demonstrated that the loss of clpP introns were present in all Poaceae species (Guisinger et al., 2010), which was supported by our findings. This study revealed that genomic structure, gene content and total GC content were significantly similar or identical within 24 genera from Pooideae, which were consistent with the genus Blumea and the families Solanaceae, Malvaceae, and Araceae (Abdullah et al., 2020b,c, 2021b).

Length variation in the IR region of the chloroplast genome was a common phenomenon during the evolution of land plants, which has led to the formation of diverse boundary features (Yang et al., 2010; Wang et al., 2017; Ding et al., 2021). The study demonstrated that boundary genes in the species of the subfamily Pooideae were mainly rpl22, rps19, rps15, ndhF, ndhH, and psbA, which differ from the boundaries of Clethra and Blumea species (Abdullah et al., 2021b; Ding et al., 2021). In general, the subfamily Pooideae shared many similarities at the nodes, which further endorsed the idea that the boundary features were relatively stable among closely related species (Liu et al., 2018). This phenomenon has also been observed in the subfamily Asteroideae (Abdullah et al., 2021b). However, distinct junction characteristics also existed in related species, such as Brachypodium stacei and Briza maxima. The present study found that although both were species of the subfamily Pooideae, they formed different boundary features due to noticeable contraction or expansion of the IR Regions, respectively. The same findings were also noted in the genera Pelargonium and Psilotum (Chumley et al., 2006; Grewe et al., 2013; Sun et al., 2013).

The results of the mVISTA analysis showed that the coding regions were more conserved than the non-coding regions in the cp genomes of the subfamily Pooideae, and the two single copy regions showed higher variation potential than the IR regions. These two findings agreed with previous studies in other plant taxa (Gu et al., 2016; Xu et al., 2017; Alzahrani et al., 2020). We detected some highly variable non-coding regions, such as rbcL-psaI, psbE-petL, trnD-GUC-psbM, and rpl32-trnL-UAG. Despite the relative conservation of the protein-coding regions, variations were also observed in rpoC2, infA, cemA, and matK genes. The highly variable regions detected in this study were promising to be developed as specific DNA barcodes for the subfamily Pooideae, which has positive implications for the identification of species. In addition, the high GC content might be one of the reasons for less variation in tRNA sequences and IR regions, which further demonstrates the significance of GC content in maintaining sequence stability (Necsulea and Lobry, 2007; Kim et al., 2019).

The codon usage preference is closely related to gene expression and affects protein and mRNA levels in the genome (Zhou et al., 2013; Lyu and Liu, 2020). The most abundant amino acid in the C. subtilis was leucine 2135 (10.76%), which has also been frequently reported in the chloroplast genomes of other angiosperms (Jian et al., 2018; Somaratne et al., 2019). More interestingly, most codons ending in A/U have RSCU values greater than 1, while those ending in C/G are less than 1. This pattern also applies to the preference of codon usage in other plants (Wang et al., 2018; Liu X.Y. et al., 2020).

Oligonucleotide repeats are very common in plastid genome and are thought to be a proxy for identifying mutational hotspots (Ahmed et al., 2012; Lee et al., 2014; Abdullah et al., 2020a,d; Liu Q. et al., 2020). In the present study, we detected both forward and palindromic repeats, mostly distributed in the LSC region. Additionally, most of the repeats were 30–40 bp in length, which was similar to those found in other species (Chen et al., 2018; Li D.M. et al., 2019; Wu et al., 2020). Simple sequence repeats (SSRs) were often used as a molecular marker to explore population relationships and evolutionary history due to its polymorphism, co-dominance and reliability (Oliveira et al., 2006; Sonah et al., 2011; Gao et al., 2018). A total of five types of SSRs were detected in the cp genomes of C. subtilis and its related species, of which mono-nucleotide repeats were the most common. Similarly, the most abundant SSR type in the genus Quercus was also mono-nucleotide repeats (Yang et al., 2016). However, there are other possibilities, such as tri-nucleotide repeats occurring most frequently in Urophysa (Xie et al., 2018). Furthermore, this study not only found that most SSR types were mono-nucleotide repeats, but they had A/T preference. This phenomenon can also be observed in numerous other taxa (Wheeler et al., 2014; Munyao et al., 2020).

We identified 13 polymorphic regions (matK, trnK-UUU/rps16, rps16/trnQ-UUG, trnG-UCC/trnT-GGU, trnT-GGU/trnE-UUC, petN/trnC-GCA, trnC-GCA/rpoB, rps4/trnL-UAA, trnL-UAA/ndhJ, ndhC/trnV-UAC, ndhF, ndhF/rpl32, and ndhA) with nucleotide diversity >0.05, mainly located in the LSC region. In addition, the nucleotide diversity values within the IR regions were significantly lower than those in the single copy regions, which is consistent with the pattern found in previous studies (Li D.M. et al., 2019; Ding et al., 2021). The dN/dS analysis was regarded as one of the most popular and reliable measures to quantify selective pressure (Kryazhimskiy and Plotkin, 2008; Mugal et al., 2014). We performed a selection pressure analysis on different genera of the subfamily Pooideae, and the result indicated that there are some genes under positive selective pressure, which was crucial for understanding the evolutionary history of these genera. The positively selected genes identified were nearly identical to those previously reported for other species in the family Poaceae, and our findings further support the plausibility of these loci (Piot et al., 2018). Furthermore, these genes are associated with photosynthesis, self-expression and regulatory activity (Piot et al., 2018), which has a positive effect on understanding the mechanisms of selection pressure generation.

In the current study, the 76 protein-coding regions of the chloroplast genome were used for the first time to explore the phylogenetic position of Coleanthus subtilis. The reconstructed phylogenetic tree divided the 53 species into ten subtribes and six tribes, which coincided with the broad framework of the Poaceae phylogeny (Soreng et al., 2017; Saarela et al., 2018; Tkach et al., 2020). Phylogenetic analysis strongly demonstrated that C. subtilis formed a sister branch with the genus Phippsia (BS = 100, PP = 1), which further justified the results of previous morphological treatments and phylogenetic studies based on chloroplast fragments (Tzvelev, 1976; Soreng et al., 2015; Gnutikov et al., 2020). Moreover, our data revealed that Colpodium was nested in the subtribe Loliinae and was particularly closely related to the genus Castellia, while Zingeria was located in the subtribe Coleanthinae (BS = 100, PP = 1). This finding differed from that of earlier studies and provided a new perspective on the relationships between Colpodium, Zingeria and Coleanthinae. Some previous studies suggested that the genera Zingeria and Colpodium are sister groups and rather distantly related to the subtribe Coleanthinae, forming a branch known as the two-chromosome grasses (Rodionov et al., 2008; Kim et al., 2009). At the same time, these two genera were considered as constituent members of Coleanthinae (Soreng et al., 2015). However, apart from the fact that Zingeria belongs to the subtribe Coleanthinae, our results do not support the previously reported relationship between Colpodium, Zingeria and Coleanthinae. This work will not only contribute to further insight into the phylogenetic position of C. subtilis and the composition of the subtribe Coleanthinae, but also provide valuable chloroplast genomic information for future exploration of the origin and differentiation between C. subtilis and its related species at the cp genome level.