Fresh mature leaves were plucked from two to three plant individuals sampled for Thunia alba (L01 and L02) and Thuniopsis cleistogama (L03, L04 and L05). The living plants were introduced from Pu’er, Yunnan Province in Southewest China and cultivated in the greenhouse of South China Botanical Garden, Chinese Academy of Science (SCBG, CAS). Total genomic DNA was extracted from young leaves using a Trelief TM Plant Genomic DNA Kit (TsingKe Biotechnology Co., Ltd., Beijing, China). After quality testing, DNA was fragmented and used to set up 200 bp short-insert libraries and the qualified libraries were sequenced with PE150 bp on the BGISEQ-500 sequencer according to the manufacturer’s instructions. The sequencing depth was ∼3.0 Gb of 150-bp paired-end reads for each species. Voucher specimens were deposited at the herbarium of South China Botanical Garden, CAS (IBSC).

First, all raw reads were trimmed using fastp (Chen et al., 2018). Subsequently, high-quality reads were mapped to the reference chloroplast genomes of Orchidaceae obtained from GenBank through Bowtie2 v.2.3.4.3 (Langmead and Salzber, 2012). The sequence of the coding gene having the maximum coverage was utilized as a seed sequence for de-novo assembly by NOVOPlasty v4.2.1 (Dierckxsens et al., 2017). To ensure the accuracy of the genome assembly, raw sequencing reads were remapped to the candidate chloroplast genome. The final assembly included contiguous sequence without gaps. The assembled cp genomes were annotated with GeSeq (Tillich et al., 2017) and tRNAscan (Chan and Lowe, 2019), then manual adjusted and confirmed using Geneious 9.1.8 (Kearse et al., 2012). The circular chloroplast genome map was drawn by OrganellarGenomeDRAW tool (OGDRAW) v.1.3.1 (Greiner et al., 2019) for further comparison of gene order and content.

The codon usage patterns of protein-coding genes in the Thuniopsis cleistogama genomes (L03, L04 and L05) and Thunia alba genomes (L01 and L02) were estimated by the software CodonW v.1.4.2. The Relative Synonymous Codon Usage (RSCU) values and the effective number of codon (ENC) values were determined to quantify the extent of the codon usage bias. RSCU was calculated for every codon in each genome according to the published equation (Sharp and Li, 1987). The ENC values indicate the codon bias for individual genes, over a range of values from 20 to 61. Lower ENC values indicate higher codon bias (Wright et al., 1990). The overall GC content and GC content at the first, second and third codon positions (GC1, GC2, and GC3, respectively) as well as the average GC content of the first and second codon positions (GC12) of the genomes were calculated using EMBOSS software suite (Rice et al., 2000). The neutrality plot (GC12 vs. GC3) and ENC-plot analysis were performed to determine the relative contributions of mutational pressure and natural selection in shaping the codon usage patterns. In neutral plot analysis, a plot regression with a slope of zero indicates no effect of directional mutation pressure (Sueoka, 1988). The ENC-plot (ENC vs. GC3s) was mapped with the scatterplot. If the point lies exactly on the standard curve (corresponding to the ENC values), the determinant of codon preference is mutation pressure (Wright et al., 1990). Simple sequence repeats (SSRs) were searched via MISA v1.01 (Beier et al., 2017) with the following criteria: 10, 5, 4, 3, 3, and 3 repeat units for mono-, di-, tri-, tetra-, penta-, and hexa-nucleotides, respectively. Chloroplast genome similarity was assessed using BLAST Atlas on the GView server (https://server.gview.ca/) with 10 kbp connection windows (Petkau et al., 2010) with Thuniopsis cleistogama _L04 genome as a reference. Circos plot of genomic landscape for small genomic variants based on SNP and InDel data was created using the TBtools (Chen et al., 2020). The junction regions between the IR, SSC, and LSC of these plastomes were compared using IRscope online program (Amiryousefi et al., 2018). The divergent regions were visualized using Shuffle-LAGAN mode (Brudno et al., 2003) included in mVISTA v.2.0 (Frazer et al., 2004). To identify polymorphic regions with substantial variability, the aligned sequences were imported in DnaSP v6.12.03 (DNA Sequences Polymorphism) using the sliding window method with a step size of 200 bp and a window length of 600 bp (Rozas et al., 2017).

We reconstructed the phylogenetic relationships among 18 species representing eight genera in subtribe Coelogyninae based on 24 complete chloroplast (cp) genomes, including five newly sequenced plastomes for each species of Thuniopsis, Thunia and 19 previously reported plastomes of 16 species downloaded from GenBank database. Arundina graminifolia (GenBank accession No. MN171408) was included as outgroup based on previous study (Li et al., 2015). All the genome sequences were aligned using MAFFT v7.313 (Katoh and Standley, 2013) and adjusted manually by BioEdit (Hall, 1999). We used maximum likelihood (ML) and Bayesian inference (BI) methods for phylogenetic analyses. The ML tree was conducted using IQ-TREE v.1.6.12 (Nguyen et al., 2015) and web server (http://iqtree.cibiv.univie.ac.at). The best-fitting nucleotide substitution model TVM + F + R2 was determined using the Akaike Information Criterion (AIC) by ModelFinder (Kalyaanamoorthy et al., 2017) in the IQ-TREE package and 1,000 bootstrap replicates. The Bayesian inference was performed with MrBayes v.3.2.7 (Ronquist et al., 2012), employing the TVM + F + R2 model of nucleotide substitution, as determined by ModelTest-NG 0.1.6 (Darriba et al., 2019). Two independent Markov Chain Monte Carlo (MCMC) runs were performed for 1,000,000 steps with a random starting tree and sampled every 1,000 generations. Each chain started with a random tree, and the first 25% sampled trees were discarded as burn-in to construct a majority-rule consensus tree and to estimate posterior probabilities (PP).

We obtained the whole chloroplast genomes of three samples for Thuniopsis cleistogama (GenBank accession Nos. OL809657, OL809660 and OL809661) and two samples for Thunia alba (GenBank accession Nos. OL809658, and OL809659). As most angiosperms plastid genomes, the newly assembled plastomes exhibited the classical quadripartite circular structure, with two inverted repeat (IRA and IRB) regions separated by a Large Single Copy (LSC) and a Small Single Copy (SSC) by OGDRAW (Figure 1). The average assemblies of the five cp genomes varied in size from 159,550 to 159,950 bp (Table 1). The LSC regions ranged in size from 87,223 to 87,533 bp; the SSC regions varied between 18,858 and 18,870 bp; and the IR regions varied from 26,697 bp to 26,778 bp. These cp genomes of the two taxa show highly syntenic nature in gene content and order. Each of the genomes encoded a total of 135 unique genes, including 89 protein-coding genes (65.93%), 38 tRNA genes (28.15%) and eight rRNA genes (5.29%). A total of 19 genes contained two exons, including 8 tRNA genes (two trnA, trnG, two trnI, trnK, trnL, trnV) and 11 protein coding genes (atpF, ndhA, two ndhB, petB, petD, two rpl2, rpl16, rpoC1, rps16), while the other 4 protein coding genes (clpP, ycf3, two rps12) each contained three exons. Eight protein-coding genes (ycf2, ycf15, rps19, rpl2, rpl23, rps7, rps12, and ndhB), eight tRNAs (trnH-GUG, trnI-CAU, trnI-GAU, trnL-CAA, trnN-GUU, trnR-ACG, trnV-GAC, and trnA-UGC), and all four rRNAs (rrn4.5, rrn5, rrn16 and rrn23) included two copies because of their location at the IR regions. In addition, the genes ycf1, ndhF were located at the junctions of SSC/IR and rpl22 was located at the junction of LSC/IR. The rps12 gene was found to be trans-spliced, consisting of three exons, with a single copy rps12. e1 located in the LSC region, whereas two copies of rps12. e2 and rps12. e3 located in the IRs. The Guanine-Cytosine (GC) content of the IR regions (43.18–43.24%) is higher than that of the LSC (35.03–35.09%) and SSC regions (30.19–30.28%).

Codon usage bias (CUB) is an essential feature of genome which provides important information for understanding species evolution. As described by Sharp and Li (1987), relative synonymous codon usage (RSCU) refers to the ratio of its actual usage frequency of a particular codon to expected frequency in the absence of codon usage bias. The RSCU analysis identified totally 80,805–80,913 bp protein-coding genes based on the five cp genomes of Thuniopsis cleistogama and Thunia alba, accounting for 50.58–50.64% of the entire genome sequence. These genes are encoded in 26,935–26,971 codons. The RSCU value of each codon for their amino acid was calculated (Supplementary Table S1). Among the 21 amino acids (Supplementary Figure S1), 19 amino acids are encoded by two, three, four or six codons with the exception of methionine (Met) and tryptophan (Trp). There are three amino acids: Arginine (Arg), leucine (Leu), and serine (Ser), each are encoded by six different codons. On average, leucine (Leu, encoded by UUA, UUG, CUU, CUC, CUA and CUG) was the most frequent amino acid, comprising 2,808–2,809 (10.4%) of the total number of codons whereas cysteine (Cys, encoded by UGU and UGC) was the least frequently encoded amino acid, with 322–324 codons (1.2%). The results indicated that the majority of preferred codons (RSCU >1) ended in A or U, with the exception of UUG (RSCU = 1.21).

To investigate the extent of codon usage bias in the cp genomes of Thunia alba (L01 and L02) and Thuniopsis cleistogama (L03, L04 and L05), the effective number of codons (ENC) values were calculated. The detailed information of ENC values was shown in Supplementary Table S2. The ENC values varied from 30.698 to 60.224, with the highest in Thunia alba; in contrast, lowest in Thuniopsis cleistogama, displaying different trends in codon preferences between the species. Furthermore, we identified only three genes with an ENC value smaller than 35, indicating a low codon usage bias. More specifically, the ENC values of individual genes differ significantly between the two species. For gene petL, the ENC value in Thunia alba was 42.429, whereas the ENC value in Thuniopsis cleistogama was 61. To further explore their differences, the distribution of the ENC values of the coding genes in the genomes was shown in Supplementary Figure S2.

The overall GC content of the genomes was 37.22–37.26%, indicating nearly identical levels among the five chloroplast genomes. As expected, we found that the GC1, GC2 and GC3 contents varied significantly across species and also among genes in the genomes. The greatest difference of GC content was found in GC3 (Supplementary Table S3), which was widely used to better elucidate the codon usage variation. The neutrality plot revealed a weak correlation between GC3 and GC12. As shown in Figure 2, the correlation coefficient was very low, indicating that the GC composition for the three positions of the codon differed. In the neutrality plot of all the genes generated, the slope of the regression line was close to zero, and most plotted points did not lie on or along the diagonal line. These data gave evidence that the codon preference was dominated by natural selection (Sueoka, 1988). Compared to Thuniopsis cleistogama, the slope of the regression line for Thunia alba was smaller, with the data points almost form a horizontal line. This observation suggested that Thunia alba codon usage was more affected by natural selection.

The ENC-plot analysis (Supplementary Figure S3) revealed that only a few points were located on or close to the expected curve, most ENC values lower than expected values were observed to lie below the curve. These results thus confirmed that codon usage preference of the five cp genomes was mainly influenced by natural selection and other factors, while mutation pressure played only a minor role (Wright et al., 1990).

Simple sequence repeats (SSRs), also known as microsatellites, consist of short tandemly repeated DNA sequences of 1–6 base pair units. Their rapid evolution makes them useful for tracing the evolutionary history of populations and investigating patterns of selection (Pauwels et al., 2012; Asaf et al., 2016), In this investigation, we detected a total of 186 simple sequence repeats (SSRs) with four types (i.e., mono-, di-, tri-, and tetra-nucleotide repeats) for the chloroplast genomes of each species sampled. The details of all cpSSRs identified in these plastomes are represented in Supplementary Table S4. The MISA analysis (Beier et al., 2017) identified 36 to 37 SSRs in each Thuniopsis cp genome, whereas 38 SSRs in each Thunia cp genome. The tetra-nucleotide SSR was only identified in Thunia whereas tri-nucleotide SSR only exists in Thuniopsis (Figure 3). In each case, The SSRs varied with the number of repeats for each type depending on the species and individuals examined. The results indicated that the selected SSRs can detect a relatively wide genetic diversity in genomes of different individuals or populations.

Among all SSRs identified across the genomes, mononucleotide repeats were the most common SSRs (29–34, the average percentage of 84.41%), followed by di-nucleotide SSRs, accounted for 12.9% of the total repeats, while tri-nucleotide and tetra-nucleotide SSRs occurred less frequently (1.61, 1.08% of all SSRs, respectively). The majority of SSRs were located in the LSC regions (25–29, the average percentage of 72.58%), followed by SSC regions (7–9, the average percentage of 22.04%). Only two SSRs were identified in IR regions, made up 5.38% of all SSRs. Only one tri-nucleotide SSR (for species Thuniopsis cleistogama) and one tetra-nucleotide SSR (for species Thunia alba) were identified in IR, made up 5.38% of all SSRs. In addition, SSRs were composed primarily of A or T with an obvious A/T bias in these plastomes. The highest mononucleotide SSRs (89.47%) composed of A/T units. Meanwhile, AT/TA repeats were the most common among di-nucleotide SSRs (A+ T: 7.89%).

To investigate intra-generic and interspecific divergence, a graphical genome map was generated using the Gview tool (Petkau et al., 2010). The overall genome features and structural variations were assessed for the plastomes sequenced and assembled here and previously reported plastomes with Thuniopsis cleistogama _L04 genome as a reference. All the nine cp plastomes from six species showed high similarity in genomic structure. However, much higher genetic inconsistency occurred in the LSC and SSC regions compared to the IR regions. No significant difference was detected between individuals or populations of the same species. The variability among species within the same genera was generally not significant than the diversification across species in different genera (Figure 4).

For estimating genetic variations across the genomes of five newly sequenced samples and four plastomes available in GenBank, we identified single nucleotide polymorphisms (SNPs) and insertions/deletions (InDels) in relation to the reference genome. The Circos plot represents the distribution pattern of SNPs and InDels in the genic regions of the nine plastomes (Figure 5). The SNPs and InDels that detected for each plastome alignment were statistically analyzed (Supplementary Table S5). The polymorphisms identified from our alignment were not evenly distributed across the segments and varied clearly across different taxonomic levels (intergeneric and intrageneric). Compared with the reference genome of Thuniopsis cleisostogama_L04, genetic diversity among individuals of the same species are relatively low. We found only 3 SNPs and 8 InDels in the sample of T. cleisostogama L05. In comparison, the number of SNPs detected in two individuals of Thunia alba ranged from 2,006 to 2,044, InDels ranged from 1,358 to 1,446. The number of SNPs detected in Bletilla taxa ranged from 1,958 to 2,026; InDels ranged from 1,195 to 1,291. The number of SNPs detected in Pleione taxa ranged from 1,917 to 1,972; InDels ranged from 1,363 to 1,518. Totally 11,923 SNPs and 8,171 InDels were identified among these plastomes (rings “2–10”). The average numbers of SNPs and InDels are 1,987 and 1,362, respectively, with the ratio 1.459: 1, indicating that SNPs represent the most common type of DNA polymorphism in these genomes. In general, SNPs and InDels were spread over the genomes with a similar distribution pattern. The majority of these polymorphisms were found in the noncoding and intergenic regions of the genomes. Furthermore, genes at the IR and SSC junctions showed comparably higher variations. Among more than 7,000 variant sites, only more than 1,400 were located in the coding regions. Therefore, the overall variants of non-coding regions are higher than those of coding regions.

The LSC/IR and SSC/IR borders of the three cp genomes of Thuniopsis cleistogama were compared with the two cp genomes of Thunia alba, and of the published cp genomes of four other closely related species (Bletilla formosana, B. striata, Pleione bulbocodioides and P. maculata). As shown in Figure 6, the rps3 gene was situated in the LSC region. The ycf1 and ndhF genes were located at the junction of the SSC/IR borders, while rpl22 was mainly located in the boundary of the LSC/IRb junction. Two copies of the rps19 genes were present in the IRa and IRb regions, near the LSC/IR borders. In all the cp genomes at the SSC/IR borders, the IR extended into the ycf1 gene to create a long ycf1 gene of 1,028–1,034 bp. Two copies of the ycf1 genes crossed the SSC/IRa and SSC/IRb borders, respectively. As a whole, the junctions between the IR and LSC/SSC regions slightly varied among these cp genomes. The rpl22 gene extended from the LSC regions into the IRb regions by 35 bp in Thuniopsis cleistogama, Thunia alba, Bletilla formosana and B. striata. This distance was 37 bp away from the LSC/IRb border in Pleione bulbocodioides and P. maculata. The functional ycf1 gene spanned the IRa and SSC regions, with 4,557 to 4,560 bp sequences situated at SSC region in Thuniopsis cleistogama and Thunia alba, and 4,551 to 4,563 bp in Bletilla formosana and B. striata, 4,482 to 4,546 bp in Pleione bulbocodioides and P. maculata. The sizes of the fragments located in the IRa regions were 1,011 bp in Thuniopsis cleistogama and Thunia alba, 1,011 to 1,014 bp in Bletilla formosana and B. striata, 1,014 to 1,022 bp in Pleione bulbocodioides and P. maculata, respectively. The ndhF gene was situated at the IRb/SSC boundary, with 55 bp sequences situated at the IRb region in Thuniopsis cleistogama, Thunia alba and Bletilla formosana, the comparable region in Bletilla striata and Pleione bulbocodioides is 58 bp long, whereas in P. maculata, with 66 bp apart from the SSC/IRb border.

To elucidate the level of the plastome divergence, the complete chloroplast genomes of Thuniopsis were compared and plotted using mVISTA (Frazer et al., 2004) by aligning the nine cp genomes with the reference T. cleistogama L04 genome (Figure 7). The noncoding regions showed a significantly higher degree of sequence divergence than that observed in the coding regions. The LSC and SSC regions were more divergent between species than the two IR regions. Most sequence variations were found concentrated in the intergenic regions, which exhibited noticeably higher divergence than the other regions.

In order to detect highly variable regions, polymorphic sites and nucleotide variability (polymorphism information, Pi) were calculated using a sliding window analysis (Figure 8). Among the nine plastomes, Pi values ranged from 0 to 0.13278 (petN-petM). The sequences were more conserved in the IR regions than in the LSC and SSC regions in all plastomes. The nucleotide variability was higher in the SSC (Pi = 0.0190) and LSC (Pi = 0.0135) regions, when compared to IR regions (Pi = 0.0031), which had a much lower nucleotide diversity. The Pi values in non-coding regions (with an average value of 0.019448) showed comparably higher divergence levels than the coding regions (with an average value of 0.01008).

In the LSC region, one intergenic region (petN-petM) showed the highest Pi value of 0.13278, followed by the atpB-rbcL with Pi value of 0.12639. In the SSC region, ndhA showed a higher evolution rate compared with other genes, with Pi value of 0.12972. Eight hypervariable regions among the cp genomes were recognized as hotspot regions with nucleotide diversity >0.08. These regions were petN-petM (0.13278), ndhA (0.12972), atpB-rbcL (0.12639), ndhF-rpl32 (0.09556), psaC-ndhE (0.08667), ndhJ (0.08278), trnW-trnP (0.08194) and rps16-trnQ (0.08111). Among them, ndhF-rpl32 spanned the SSC/IRb boundary. ndhA and psaC-ndhE were situated in the SSC region, whereas five of eight were located in the LSC region (petN-petM, atpB-rbcL, ndhJ, trnW-trnP and rps16-trnQ).

Plastome sequences varied most in the intergenic spacer (IGS) regions. Six regions showed Pi values higher than 0.08 and all of these regions were located in the IGS region, e.g., petN-petM, atpB-rbcL, ndhF-rpl32, psaC-ndhE, trnW-trnP and rps16-trnQ. Two genes: ndhJ of the LSC region and ndhA of the SSC region present higher nucleotide variability than the mean gene. These divergent hotspots show potential value for the development of molecular markers for phylogenetic and phylogeographic analysis.

The ML and BI analyses of the complete chloroplast genomes generated consistent topologies. The ML topology was selected for discussion, with ML bootstrap (MLBS) and posterior probabilities (PP) values are given near nodes (Figure 9). The sister relationship between genera Bletilla and Thunia was recovered (PP = 1.00; BS = 99). Thunia represented by T. alba, and Bletilla represented by B. formosana, B. striata and B. ochracea. The same clade was seen in previous analyses (Li et al., 2015). All phylogenetic analyses consistently indicated that three sampled specimens of Thuniopsis cleistogama formed a well-supported monophyletic lineage (PP = 1.00; BS = 100), sister to the clade containing species of the rest genera examined within the subtribe Coelogyninae, including Bulleyia, Coelogyne, Panisea, Pholidota and Pleione. These results indicated that Thuniopsis formed an independent lineage that genetically separated from the rest genera in Coelogyninae.

Generally, our results showed d that the chloroplast genomes from six species, representing the four genera (Bletilla, Thunia, Thuniopsis and Pleione) were rather similar in terms of structural organization, gene content and gene order. Genome annotation identified a total of 135 unique genes, consisting of 89 protein-coding genes, 38 tRNA genes, and eight rRNA genes in each genome (Table 1). Nevertheless, the cp genome size varied within each genus, species, and even within individual from the same species, ranging from 158,394 in P. maculata to 159,269 bp to 159,950 bp in Thunia alba_L02. The LSC regions ranged in size from 86,603 bp in P. maculata to 87,533 bp in Thunia alba_L02; the SSC regions ranged from 18,499 bp in P. maculata to 18,870 bp in Thuniopsis cleistogama_L03 and L04, and the IR regions ranged from 26,646 bp in P. maculata to 26,801 bp in Bletilla formosana. In comparison with other species of related genera, Thuniopsis cleistogama possesses a relatively larger SSC region (18,868–18,870 bp).

Genomic DNA base composition (GC content) is predicted to significantly affect genome functioning (Šmarda et al., 2014). Our analysis showed that the GC contents varied among different genomes and different regions within a genome. When compared with the LSC (34.96–35.12%) and SSC (30.19–30.41%) regions, we observed significantly higher GC contents in the IR regions (43.18–43.24%). The high GC content could be caused by the presence of eight rRNA genes with high GC content in these regions. The result was consistent with previous studies (Darshetkar et al., 2019; Wu et al., 2020). GC content generally showed a quadratic relationship with genome size. Interestingly, with respect to its smallest genome (158,394 bp), Pleione maculata genome (MW699846) has a relatively higher GC content (37.32%) than expected. According to Šmarda et al. (2014), in monocots, increased GC content was associated with increased tolerance and ability to grow in seasonally cold and/or dry climates.

Although the cp genomes of Thuniopsis and its related genera showed extremely conserved structure, slight variations were found at the LSC/IR and IR/SSC boundaries of the plastomes. Compared with the LSC/IR border, the IR/SSC border displayed higher variations among species within the four genera (Figure 6). The relatively higher divergence was seen in genes ycf1 and ndhF, which might play an important role in the evolutionary of the cp genome of these species.

Consistent with the diversity patterns found in most angiosperms (Lu et al., 2017; Lee et al., 2019; Zheng et al., 2020), mVISTA analysis (Figure 7) has revealed that sequence variability in non-coding regions was greater than that in coding regions. Overall, sequence divergence was concentrated in the non-coding region, and the sequence divergence in the protein-coding regions (CDSs) was limited. Further, the cp genome sequence comparison has detected large variability in the conserved noncoding sequences regions (CNS). The diversity in CNS is nonrandom, with variants shared across different genera. Many studies have found that CNSs are enriched significantly in regulatory sequence elements. CNSs may have roles in the evolution and some critical biological function (Burgess and Freeling, 2014; Hettiarachchi et al., 2014; Xie et al., 2018).

The nucleotide diversity (Pi) values also indicated the conserved nature of the IR regions (Figure 8). The sequences in the IR regions were more conserved as compared to the LSC and SSC regions. In particular, intergenic regions (IGS) showed significantly higher variability than that in genic regions. Eight highly variable regions were identified, which show high potential value for future DNA barcoding and phylogenetic reconstruction. They are six intergenic markers (petN-petM, atpB-rbcL, ndhF-rpl32, psaC-ndhE, trnW-trnP and rps16-trnQ) and two genic makers (ndhA, ndhJ).

As for the cp genomes newly sequenced in this study, three examined accessions of Thuniopsis cleistogama were resolved into a highly supported monophyletic lineage. All four accessions of Thunia alba including two accessions obtained from GenBank clustered together as expected. The phylogenetic relationship of Thuniopsis and its closely related genera was basically congruent with those of previous studies but found strong support.

In the previous phylogenies by Li et al. (2015) based on partial DNA markers (ITS, matK and trnL), Thuniopsis formed a poorly supported sister clade to the Bletilla, Dilochia and Thunia group, with rather low ML bootstrap value of 62%. In this study inferred from the complete plastomes, Thuniopsis and Thunia-Bletilla were strongly supported as successively diverging lineages. Bletilla and Thunia were grouped into one clade in both trees with high bootstrap values. The main lineages investigated here are clearly classified with a better resolution. Currently, we were unable to obtain the sequences for genus Dilochia. Dilochia represented the sister group of Thunia in former study (weakly supported). More extensive cp genome sampling is necessary to further illustrate the relationships of species in subtribe Coelogyninae.

In addition to its significant genetic divergence, several studies have suggested that Thuniopsis exhibits unique morphological features that differentiate it from other members of Coelogyninae (Li et al., 2015; Kang et al., 2019; Kurzweil et al., 2020). Thuniopsis is superficially similar to Thunia in the elongate leafy stem, carrying distichous, soft and plicate leaves, but can be easily distinguished from the latter by having prominent subterranean corms; much smaller and spurless flowers; prominent stigma and bilobed rostellum. By contrast, Thunia is characterized by its distinctly fleshy, cane-like stems; large and showy flowers. Thuniopsis grows in comparatively dry and hot environment. Normally, its stem and leaves of the plant die annually and become dormant during the winter months. The taxonomic status of the genus Thuniopsis was confirmed by both molecular analyses and morphological characters. The complete chloroplast genomes were proved to be more informative than cp DNA fragments in revealing the phylogeny of Coelogyninae.