1. Introduction

Bulbophyllum Thouars is probably the largest pantropical genus of Orchidaceae, with ca. 2200 species (Pridgeon et al., 2014). The Paleotropical region is the richest in species of this genus, with ca. 1700 species occurring in Asian (sub) tropical area, followed by Africa and then the neotropics (Pridgeon et al., 2014). It is well known that Bulbophyllum represents one of the most significant diversifications genera in the orchid family (Gamisch and Comes, 2019). Bulbophyllum species are notoriously difficult to classify taxonomically, due to highly morphological variation (Seidenfaden, 1979). Many independent and closely related genera established then combined. Numerous sectional classifications remained controversial (Pridgeon et al., 2014).

Several studies have been reported in recent years examining the phylogeny of Bulbophyllum using molecular data (Chase et al., 2015; Hosseini et al., 2012, 2016; Gravendeel et al., 2014; Fischer et al., 2007; Smidt et al., 2011). The finding of these works well supported the monophyly of continental groups, e.g., Asian, African and Neotropical clade, but contributed little to elucidating the relationships at infrageneric, interspecific level. Very recently, Hu et al. (2020) combined DNA and floral characters to reconstruct the phylogenetic relationship in Asian Cirrhopetalum alliance, and further indicated that the majority of traditional sectional classifications were shown to be highly artificial. The elucidation of the phylogenic relationships within Bulbophyllum still requires more detailed study.

The chloroplast (cp) genomes are considered to be useful and informative genetic resources for studies on evolutionary relationships in the plant kingdom at various taxonomic levels (Asaf et al., 2017; Wu et al., 2020), because of their uniparentally inherited and relatively conserved genome structure and higher evolutionary rate (Shaw et al., 2014; Wicke et al., 2011). Especially in Orchidaceae, the size and gene content of cp genome vary greatly. The size of the cp genome within orchids is ranging from 10 to 170 kb (Dong et al., 2018; Lin et al., 2015; Kim et al., 2015). The most reduced chloroplast genomes have been found in several mycoheterotrophic genera, e.g., Epipogium, Gastrodia, and Aphyllorchis (Feng et al., 2016; Schelkunov et al., 2015; Yuan et al., 2018), which have lost many photosynthetic genes (Graham et al., 2017). There are more than 100 orchid cp genomes, including Dendrobium 40 spp., Holcoglossum 11 spp., Corallorhiza 10 spp., Cymbidium 9 spp., Neottia 7 spp., have been reported and available in the recent years (Barrett and Davis, 2012; Kim et al., 2020a; Li et al., 2019; Logacheva et al., 2011; Niu et al., 2017; Shao et al., 2020; Yang et al., 2013; Zhu et al., 2018). However, the cp genomes are still little known in Bulbophyllum. Only nine cp genomes of Bulbophyllum species (eight neotropical and one Asian) were reported until now (Kim et al., 2020a; Zavala-Páez et al., 2020). Comparative cp genome studies at different taxonomical level are needed to clarify the diversification pattern in this species-rich group, and then to develop suitable molecular markers to study speciation, evolutionary adaptation and phylogeny. For taxonomists, uncertainty in taxonomic status of the Bulbophyllum species is quite often, which also requires suitable molecular markers for taxonomic revision (Hosseini et al., 2016; Zavala-Páez et al., 2020).

The Asian section Macrocaulia (Bl.) Aver. is one of acceptable natural sections in Bulbophyllun supported both by morphological characters and molecular evidence, and it is a medium section with more or less 68 species, mainly distributed in Southeastern Asia (Averyanov and Komarov, 1994; Vermeulen, 2014). There are six species recorded in (sub) tropical regions of China, which are the most northern parts for section Macrocaulia (Chen and Vermeulen, 2009; Jin et al., 2019; Huang et al., 2020). The section Macrocaulia is characterized by mini-miniature sized, epiphytes with creeping, close set pseudobulbs carrying a single, apical, usually persistent leaf that blooms on a single flower that have morphological diverse labellum and sepals. Based on the nrITS, the section Macrocaulia is monophyletic and regarded as the basal group sister to the rest of Asian Bulbophyllum clade (Gravendeel et al., 2014).

In this study, the complete chloroplast genomes of three species from section Macrocaulia were sequenced, assembled and characterized, then to investigate chloroplast genome diversification pattern and examine species divergence and evolutionary adaptation. The main objectives of the present study were to (1) understand the interspecific variation pattern in the cp genomes of closely related species in the genus of Bulbophyllum; (2) test the monophyly of section Macrocaulia using chloroplast genome; (3) identify genes underlying positive selection as potentially genes for adaptive evolution, radiation and speciation; (4) examine hotspot regions as candidate sequences for species delimitation and phylogenetic analysis in the genus Bulbophyllum.

2. Materials and methods 2.1. Plant material

Three species, Bulbophyllum lingii M.Z. Huang, G.S. Yang & J.M. Yin, B. pantenurum Seidenf., and B. menghaiense Z.H. Tsi, are closely related species from section Macrocaulia. B. lingii is newly found (Huang et al., 2020) and B. pantenurum is a new record (unpublished) in China recently. Bulbophyllum gedangense Y. Luo, J. P. Deng & Jian W. Li, which is also newly found in China recently, might belonging to section Desmosanthes (Bl.) J. J. Smith (Luo et al., 2020). Fresh leaves of four Bulbophyllum species were collected from the living collections at the greenhouse of Xishuangbanna Tropical Botanical Garden, Chinese Academy of Sciences and Yunnan Fengchunfang Biotechnology Co. Ltd. The voucher specimens were deposited in the herbarium of Xishuangbanna Tropical Botanical Garden, Chinese Academy of Sciences (HITBC) (Table 1). In addition, we downloaded the available complete cp genomes of nine other Bulbophyllum species from GenBank (Table 1) to analyze phylogenetic relationships.

2.2. DNA extraction and sequencing

Total genomic DNA from fresh leaves was extracted by using TiangenDNA kit (TIANGEN, China). An Illumina paired-end DNA library was constructed using the IlluminaTruSeq Library Preparation Kit (San Diego, CA, USA) following the manufacturer's instructions. The library was sequenced by the Illumina Hiseq 2500 sequencing platform (Illumina, CA, USA) at Personal Biotechnology Co., Ltd (Shanghai, China).

2.3. Chloroplast genome assembly and annotation

The total 3 G data of sequencing reads were trimmed to high quality reads (Table 2), and then GetOrganelle pipeline were used for de novo chloroplast genome assembly (Jin et al., 2020). CPGAVAS2 (Shi et al., 2019) and GeSeq (Tillich et al., 2017) were used to annotate the assembled chloroplast genome using default parameters to predict protein coding, rRNA, and tRNA genes. The annotation results were checked plus manual corrections and codon positions were adjusted by comparing them with the more similar six Bulbophyllum reference species present in the database using Geneious prime. The genome map of the species was illustrated with the help of OGDRAW (Greiner et al., 2019), and the annotated chloroplast genome sequences were submitted to NCBI (Table 1). GC content was calculated using Geneious prime.

2.4. Simple repeat sequence analysis

We identified simple sequence repeats (SSRs) of Bulbophyllum with MISA (Beier et al., 2017) by setting the minimum number of repeats to 10, 5, 4, 3, 3, and 3 for mono-, di-, tri-, tetra-, penta-, and hexa-nucleotides, respectively.

2.5. Comparative genome analysis

The complete chloroplast genomes of four Bulbophyllum species were compared using the program mVISTA (Frazer et al., 2004). We also compared the borders of large single copy (LSC), small single copy (SSC), and inverted repeat (IR) regions in the genomes of the four species. The nucleotide diversity value (Pi) of the cp genomes of all four Bulbophyllum species was calculated using DnaSP v5.10 (Librado and Rozas, 2009). The number of mutations and indels events was calculated by DnaSP v5.0 (Rozas et al., 2004), and their rate was calculated per 100 bp. Due to indels having higher levels of homoplasy than mutations, indels were considered as events instead of sites in the alignment (Ingvarsson et al., 2003).

2.6. Selective pressure analysis

In order to detect the protein coding genes under selection in four Bulbophyllum species, the synonymous (dS), non-synonymous (dN) nucleotide substitution rates and the dN/dS ratio were calculated using the model in DnaSP v.5.0 (Librado and Rozas, 2009). To explore the selection patterns and identify positive selection on the protein coding genes, six models, model M0, M3, M1a, M2a, M7 and M8, implemented in the PAMLX Codeml program (Xu and Yang, 2013). The codon frequencies were determined by the F3 × 4 model. 79 common genes were examined. The Bayes empirical Bayes (BEB) method was used to statistically identify sites under positive selection with posterior probabilities≥0.95. Pairwise distances and branch lengths were computed using the maximum likelihood (ML) criterion under a codon model. The likelihood ratio test (LRT) was used to estimate the quality of each model (Yang and Nielsen, 2002; Yang et al., 2005).

2.7. Phylogenetic analysis

The complete chloroplast genomes of the four newly sequenced Bulbophyllum species, together with one Asian (from Japan) and eight neotropical Bulbophyllum species, and two species from Dendrobium (Table 1), which download from GenBank, were selected to reconstruct phylogenetic tree. Phylogenetic analyses were performed for two datasets: complete plastid genome sequences and the 68 protein coding sequences. The chloroplast DNA sequences were aligned with MAFFT v.7 (auto strategy) (Katoh and Standley, 2013). Maximum parsimony (MP), Maximum Likelihood (ML) and Bayesian (BI) analyses were run for each data set. Parsimony analyses were performed with MegaX (Kumar et al., 2008). Heuristic searches were conducted with tree bisection-reconnection (TBR) branch swapping, and random stepwise addition with 1000 replications. ML trees with 1000 BS were produced with IQ-tree (Trifinopoulos et al., 2016). Bootstrap analyses were used to evaluate the support for individual clades with 1000 replicates. BI tree used MrBayes 3.2.4 (Ronquist et al., 2012) in PhyloSuite (Zhang et al., 2020). The Markov chain Monte Carlo (MCMC) algorithm was run with two independent chains with a random starting tree and default priors for 2 × 10⁶ generations, with trees sampled every 100 generations. The first 25% of all generations were excluded, and a consensus phylogenetic tree was obtained based on majority rule from the remaining trees.

3. Results 3.1. Characteristics of Bulbophyllum chloroplast genome

The cp genomes of three species of section Macrocaulia ranged in size from 156, 182 bp (Bulbophyllum pentaneurum) to 156, 689 bp (B. lingii), and that of B. gedangense is 158, 524 bp (Fig. 1, Table 2). The genome length differences between the largest and smallest genomes were the 2342 bp. All of the cp genomes had a typical quadripartite structure, similar to those of most land plants, consisting two inverted repeats (IRs, 26, 838 to 26, 919 bp) separated by large (LSC, 84, 240 to 86, 200 bp) and small (SSC, 18, 105 to 18, 632 bp) single copy regions, respectively (Fig. 1, Table 2). GC content was similar (36.6–36.8%; Table 2). In addition, the gene composition and order in the cp genomes of four Bulbophyllum species are nearly identical (Fig. 1, Table 2). Four each cp genome, 113 unique genes were annotated, including 79 protein-coding genes, 30 tRNAs, and 4 rRNAs (Table 2). Most of genes occur as a single copy in LSC or SSC, while 20 genes are duplicated in the IRs, including eight protein coding genes, eight tRNAs and four rRNAs (Table 3). Eight protein coding genes and six tRNA genes included one intron, while three genes (ycf3, clpP, rps12) included two introns (Table 3).

3.2. Characters of simple sequence repeats

A total of 45 (Bulbophyllum pentanurum) to 58 (B. menghaiense) SSRs were detected in the cp genome of four Bulbophyllum species, and six categories of SSRs (mono-, di-, tri-, tetra-, penta-, hexa-nucleotide repeats) were found with lengths ranged from 10 to 26 bp (Fig. 2, Table 4). In the four Bulbophyllum species, the most abundant SSRs were mono-nucleotide repeats, accounting for 63.64%–78.57%, and number varies from 28 to 44, followed by 4 (6.9%) to 8 (14.29%) di-nucleotide repeats, 0 to 4 (6.9%) tri-nucleotide repeats, 2 (3.57%) to 10 (17.24%) tetra-nucleotides repeats, 0 to 2 (3.45%) penta-nucleotide repeats, and 0 to 2 (3.45%) hexa-nucleotide repeats (Fig. 2A, Table 4). Interestingly, all mono-nucleotide SSRs belonged to A or T type and the majority of di-, tri-, tetra-, penta-, and hexa-nucleotide SSRs were especially rich in A or T (Table S1). Of the SSR regions, 31 to 40 SSRs were identified in intergenic spaces, while 5 to 13 SSRs were located in the coding DNA sequence and 3 to 11 SSRs in the coding sequence introns (Fig. 2B).

According to the criteria of SSRs, i.e., presence of the same repeat units, being located in the homologous regions, and having different number of the repeat units, 35 SSRs were identified as polymorphic SSRs between Bulbophyllum species (Table 5). Most of SSR loci were specific for each species. B. menghaiense had 10 exclusive SSRs, e.g., (TTA)₄ in trnS-trnG and atpB-rbcL regions; B. lingii had six exclusive SSRs, e.g., (AT)₅ in trnL-UAA and (TTATA)₃ in clpP genes; B. gedangense had six exclusive SSRs e.g., (AT)₆ in trnT-trnL and (TA)₁₃ in trnE-trnT regions; B. pentanurum had five exclusive SSRs, e.g., (TA)₅ in psbM-trnD and (CTT)₄ in rpl14-rpl16 regions (Table 5).

3.3. Comparative genomic analysis 3.3.1. Sequence divergence

The divergence of sequences in the cp genomes of four Bulbophyllum species was plotted by mVISTA program using the annotated B. gedangense sequence as a reference (Fig. 3). A whole genome alignment revealed an overall high sequence similarity among the four Bulbophyllum species, and also showed lower sequence divergence between three species from section Macrocaulia. Furthermore, LSC and SSC regions were more divergent compared with IR regions. The regions with highest level of divergence were rps16-trnQ, rrn16, ndhF-rpl32, psbI-trnS, psaA-ycf3, trnS-rps4, ndhF-rpl32, and rps15-ycf1 (Fig. 3).

3.3.2. Mutations and indels

The number of mutations and indels events were compared among these four Bulbophyllum cp genomes. There were 4600 mutations and 943 indels calculated. The intergenic regions had the most 2320 mutations and 710 indels respectively. The intergenic regions had the highest mutation rate, 4.65 mutations per 100 bp, following by introns and the coding regions with 3.32 and 2.11, respectively (Table 6). We further analyzed the details of indels (indel size > 10 bp) of the four Bulbophyllum species based on multiple whole cp genome alignment. 189 positions showed DNA indels in at least one of the cp genomes. The three species from section Macrocaulia shared 67 deletions, mostly located in intergenic regions, which showed low similarity to B. gedangense (Table S2). These regions are related to large deletions, which become obvious in a multiple alignment of related whole cp genome sequences. The largest deletion, which is in trnE-trnT linker, is of a size of about 298 bp (at position 32, 368 bp of the B. gedangense sequence). The additional large deletions were in different intergenic regions, e.g., trnV-rps12 (237 bp), rps16-trnQ (155 bp), and ycf4-cemA intergenic regions (189 bp) (Table S2). Several shared deletions were in introns, e.g., the intron of atpF, the intron of trnL-UAA. Only a very few deletions are in the coding regions, e.g. ycf2 gene (115 bp) (Table S2).

3.3.3. Border contraction and extension

We analyzed the IR/single copy (SC) region border positions and their adjacent genes in the four Bulbophyllum cp genomes (Fig. 4). In all cp genomes of the four Bulbophyllum species, the genes rpl22, ndhF, ycf1, rps19, and psbA were located at the junction of the LSC/IRb, IRb/SSC, SSC/IRa, and IRa/LSC borders (Fig. 4). The expansion and contraction of the border regions of cp genomes were similar in the four Bulbophyllum species. There were 10–78 bp extension of rpl22 gene into IRb region. Moreover, the IRb/SSC borders were well conserved among the four cp genomes, of which the IRb region expanded into the gene ndhF with 70 bp. The ycf1 gene crossed the SSC/IRa junction expanding 1059 bp to IRb regions in three species of section Macrocaulia, while 1074 bp in B. gedangense. There were 206–267 bp between rps19 and the LSC/IRa border, meanwhile, the psbA generated a distance of 81–136 bp to another LSC/IRa junction.

3.3.4. Nucleotide diversity

The nucleotide diversity value (Pi) of the coding regions and intergenic regions were calculated using DnaSP for understanding the sequence divergence among these four Bulbophyllum species. The coding regions were more conserved with the Pi value ranging from 0 to 0.06329, compared to the intergenic spacer with the Pi value ranging from 0 to 0.12543 (Fig. 5). Most of highly variable sites were located in the LSC regions, followed SSC regions. The sequences in IR regions exhibited lower nucleotide diversity value. Among the coding regions, the highest nucleotide diversity (Pi = 0.06329) was in trnL-UAA gene. A total of 12 genes were detected as highly variable coding regions, e.g., trnL-UAA, rpl16, trnK-UUU (Pi > 0.03), and ycf1, psbM, matK, ndhK, petD, clpP, rps16, rps15, ndhA (Pi > 0.02) (Fig. 5A). Among the non-coding regions, the highest nucleotide diversity (Pi = 0.12543) was in the psbB-psbT region. A total of 20 intergenic regions, psbB-psbT, rps16-trnQ, psaJ-rpl33, trnS-trnG, trnR-atpA, ndhF-rpl32, psbI-trnS, rpl32-trnL, psaC-ndhE, trnG-trnR, trnE-trnT, atpI-rps2, trnS-rps4, psbM-trnD, rpl33-rps18, ndhC-trnV, atpB-rbcL, trnF-ndhJ, petN-psbM, and psbK-psbI were recognized as hotspot regions with nucleotide diversity > 0.04 (Fig. 5B).

3.4. Selection pressure analysis of evolution

Synonymous and non-synonymous nucleotide substitution patterns are very important markers for gene evolution studies. A ratio of dN/dS < 1 indicates purifying selection, dN/dS > 1 denotes probable positive selection, and dN/dS values close to one indicates neutral evolution. We estimated the selective pressure among four Bulbophyllum species. 79 protein-coding genes were used to calculate the ratio of synonymous and non-synonymous substitution, and thereby to determine whether the selection pressure has acted on individual sequences. Thus, we found all 79 genes of the dN/dS average ratio were lower than 1.0 might be under purifying selection, and no gene under positive selection in four Bulbophyllum species (Fig. 6). We found that ten genes (matK, rps15, rps16, ndhK, petA, rpoA, rpl16, rpl22, yfc1, and ycf2) of the dN/dS average ratio above 0.5, indicating relaxed selection (Fig. 6). Fourteen most conserved genes (psaC, psaI, psbF, psbI, psbJ, psbL, psbN, psbT, psbZ, petL, petG, atpH, rpl23 and rps7) with the dN/dS average values of 0.0, suggesting very strong purifying selective pressure. The remaining genes showed the dN/dS average value were ranged from 0.0 to 0.5, suggesting strong purifying selective pressure. In addition, according to the models analysis, we found a total of 6 genes 13 sites exhibiting site-specific selection (Table 7). Of these genes, the ndhK gene harbored 6 sites under positive selection, with three in ycf2, and the other four genes (ndhF, rbcL, rps4, and rps16) each had only one selection site (Table 7).

3.5. Phylogenetic analysis

The protein coding gene sequences and all the complete chloroplast genome sequences of 15 Malaxideae members were used to reconstruct the phylogenetic relationships, including the 13 Bulbophyllum cpDNA sequences, and two Dendrobium species (Table 1, Fig. 7). In our analysis, the topologies obtained based on the maximum likelihood (ML), the parsimony (MP) and Bayesian (BI) analyses were largely consistent with the different two datasets, and the bootstrap values and posterior probabilities were very high for each lineage (Fig. 7). These results suggested that there was no conflict among the entire genome data set and protein coding gene sequences. All the sampled species of Bulbophyllum were clustered into one clade with 100% bootstrap value or the Bayesian posterior probability. Within the Bulbophyllum clade, the Asian tropical and neotropical Bulbophyllum form a monophyletic clade, respectively. Within the Asian tropical clade, B. pentaneurum was sister to B. lingii with B. menghaiense as their closest relative, which strongly supported that section Macrocauliais was monophyletic; In another clade, B. gedangense and B. inconspicuum (native to Japan) grouped together (Fig. 7).

4. Discussion 4.1. Conservation and divergence of chloroplast genome in the section Macrocaulia

Based on our comprehensive comparative genomics analysis, the chloroplast genomes of section Macrocaulia showed both conservative and divergent. In terms of genome size, all of four cp genomes fell within the size of Orchidaceae species (Kim et al., 2015; Lin et al., 2015). The three species of the section Macrocaulia cp genomes were similar and smaller than that of B. gedangense. In addition, although loss of ndh genes has been previously in some orchids (Feng et al., 2016; Niu et al., 2017; Yang et al., 2013; Zavala-Páez et al., 2020), all the cp genomes of the section Macrocaulia had no evidence of ndh genes loss. The total length of angiosperm cp genomes is influenced by contractions and expansions of the IRs (Green, 2011; Weng et al., 2017). However, the decrease of cp genomes size in three species of section Macrocaulia might have been due to the large deletions. Several large deletions, more than 100 bp, were detected in the cp genomes of these miniature Bulbophyllum species, mostly in intergenic regions (Table S2). Compared with B. gedangense, 84 indels were shared in the three species of section Macrocaulia. The similarity between cp DNA sequences indicated the cp genome are highly conserved in the section Macrocaulia, and contributed to their monophyletic relationships. However, several exclusive large indels were also identified for the three species of section Macrocaulia respectively, e.g., an 87 bp-indel in rpoB-trnC intergenic region for B. menghaiense, a 65 bp-indel in the intron of rps16 gene for B. lingii, and a 43 bp indels in petN-psbM intergenic region for B. pentaneurum (Table S2). The differences between cp DNA sequences indicated interspecific differentiation in the section Macrocaulia, and might contributed to the species divergence and adaptive evolution. Large indels easily identified in whole cp genome alignments might be used for species-specific molecular markers as recently shown for a Populus tremula-specific marker that has been developed based on a 96 bp-indel (Kersten et al., 2016). The indels identified among the four Bulbophyllum cp genomes in this study might also serve as species-specific molecular marker for species delimitation and identification.

4.2. Simple sequence repeats diversity

Simple sequence repeats (SSRs) are the most frequent in genome sequences with high mutation rates, and play important roles in plants (Qin et al., 2015). The most abundant SSR type was the mononucleotide repeat, and the majority of SSRs in four Bulbophyllum species was composed of polyT and polyA repeats, corresponding with the angiosperm cp genome composition (Kim et al., 2015; Qin et al., 2015). In our study, the number and types of SSRs identified from the cp genomes were variable in four Bulbophyllum species (Fig. 2). Of di- to hexa-nucleotide SSRs among four Bulbophyllum species, most of SSRs were specific to each species, and only five SSRs shared in the section Macrocaulia. These findings indicated that SSRs exhibiting highly variable even in closely related species, might have few phylogenetic significance. Most of the previous studies revealed that the richness of SSR types is various in different species, which may contribute to the genetic variations differently among species (Kuang et al., 2011). SSRs are also effective molecular markers for interspecific and intraspecific polymorphisms and extensively used for species identification, speciation, population genetic diversity and phytogeography studies (Powell et al., 1995). In four Bulbophyllum species, SSR loci in the cp genome exhibited high level of interspecific polymorphisms. Several exclusive SSR loci presented in each species, e.g., B. menghaiense with 10 exclusive SSRs, B. lingii with six, and B. pentanurum with five. Thus, these specific SSRs are useful molecular markers for species identification and detecting genetic diversity.

4.3. Hotspot regions

We found the nucleotide sequence diversity was higher in the non-coding regions than the coding regions, which is generally consistent with most previous studies of the chloroplast genomes of angiosperms (Clegg et al., 1994). A number of cpDNA, e.g., matK, rbcL, psbA-trnH, trnS-trnG, trnL-trnF, trnF-ndhJ, trnE-trnD are previously used within Bulbophyllum to reconstruct phylogenetic relationships, but have not been shown useful for detecting infragenic or interspecific relationships (Fischer et al., 2007; Hosseini et al., 2012, 2016; Smidt et al., 2011). Smidt et al. (2011), for example, used two non-coding regions of cpDNA, psbA-trnH and trnS-trnG to examine phylogenetic relationships among Neotropical sections of Bulbophyllum, but did not support the monophyly of the neotropical group, offering very little phylogenetic value. Fischer et al. (2007), who studied Madagascan Bulbophyllum using four plastid markers, did not resolved some terminal nodes. In recent years, numerous orchid plastid markers based on the comparative cp genome have been proposed for Orchidaceae phylogeny at different taxonomical level (Dong et al., 2018; Li et al., 2019; Niu et al., 2017; Yang et al., 2013; Zavala-Páez et al., 2020; Zhu et al., 2018).

Based on comprehensive comparative analysis, a total of 20 sequence divergence hotspot regions located in intergenic regions will be helpful to analyze phylogenetic relationships especially between closely related taxa, or where recent divergence, rapid radiation in Bulbophyllum (Daniell et al., 2016; Tonti-Filippini et al., 2017). Compared with the previous study of cp genomes of Neotropical Bulbophyllum (Zavala-Páez et al., 2020), where 10 potential markers proposed for barcoding and phylogenetic studies in Bulbophyllum, there are four regions identical (psbB-psbT, trnS-trnG, trnR-atpA, rpl32-trnL) with our work. In addition, three protein coding genes (trnL-UAA, rpl16, trnK-UUU) also showed remarkably higher Pi values (> 0.03). These genes might be undergoing rapid nucleotide substitution at the species level, and might be applied for molecular evolution and adaptive evolution analysis.

4.4. Adaptive evolution

The dN/dS ratios have been widely used to infer the evolutionary dynamics and identify adaptive signatures among species (Yang and Bielawski, 2000). The dN/dS value in all of 79 protein coding genes in four Bulbophyllum species were lower than 1. The dN/dS values are higher in protein-coding genes as a result of more frequent non-synonymous nucleotide substitutions than the synonymous substitutions, and the adaptation evolution may have contributed to the elevated dN/dS ratios (Makałowski and Boguski, 1998). Thus, the elevated dN/dS ratios indicated some unknown selective forces on Bulbophyllum species and results in species divergence (Hurst, 2002), which may increase our understanding of the evolutionary speciation and environment adaptation in Bulbophyllum species. We found that ten genes under relaxed selection (Fig. 6) and six genes under site-specific selection (Table 7). Of these genes, petA encoding cytochrome f, ndhF and ndhK encoding subunits of NADH dehydrogenase, rbcL encoding the large subunit of ribulose-bisphosphate carboxylase, all functioned in photosynthesis (Green, 2011; Table 3). A previous study showed that rbcL is often under positive selection in land plants, and evolved under strong positive selection after the C3–C4 photosynthetic transition in grasses (Kapralov and Filatov, 2007; Piot et al., 2017). The genes rpoA, rps15, rps16, rpl16, rpl22 all functioned in self-replication. The other genes are matK acting as a maturase, ycf1 and ycf2 with unknown function (Daniell et al., 2016; Green, 2011). These genes play important role in improving energy conversion efficiency, therefore they might have some evolutionary and adaptive significance for Bulbophyllum species. Different from our results, selection analysis in eight Neotropical Bulbophyllum species showed that four genes (ycf1, ycf2, accD, and rps3) were under site-specific selection (Zavala-Páez et al., 2020). The differences of site-specific selection pattern between Asian and neotropical Bulbophyllum indicated that they might each had various adaptive strategies under unpredictable environmental conditions. Our findings from selective pressure analysis could lead to a better understanding of the evolution of Bulbophyllum species.

4.5. Phylogenetic relationship

Complete chloroplast DNA sequences were valuable for analyzing phylogenetic relationships (Daniell et al., 2016; Tonti-Filippini et al., 2017). Kim et al. (2020a), for example, utilized the chloroplast genomes of 124 species of orchids to evaluate the phylogenomics and relationship within the Orchidaceae family. Moreover, the whole cp genome phylogeny that was generated at different taxonomical level in several groups of Orchidaceae, such as subtribe Aeridinae (Kim et al., 2020b), subtribe Calypsoinae (Li et al., 2020), and Holcoglossum (Li et al., 2019). Considering that most species-level diversity of Orchidaceae, a whole cp genome phylogeny is highly desired for Bulbophyllum. We examined phylogenetic relationships among 13 Bulbophyllum species from two continental groups in this study. Both the complete cp genome sequences and 65 protein coding genes datasets produced the same topology with similar support. In addition, all the ML, BI and MP phylogenetic results strongly supported that B. pentaneurum is sister to B. lingii with B. menghaiense as their closest relative, and all 13 Bulbophyllum species were closely clustered into one lineage, which strongly support that section Macrocaulia and the genus of Bulbophyllum are monophyletic, respectively (Fig. 7). The monophyly of section Macrocaulia was also supported by Gravendeel et al. (2014) using nuclear ITS sequence. However, our phylogenetic analysis did not show that section Macrocaulia was the basal clade in Asian Bulbophyllum group as their work, might due to sampling limited. Furthermore, the phylogenetic relationships among the Neotropical clade were also in accordance with recent phylogeny research (Zavala-Páez et al., 2020). In the future, more complete cpDNA sequences of member individuals of Bulbophyllum are needed to exploit the power of the whole cp genome phylogenies, especially for differentiation at lower taxonomic levels.

Author contributions

YL and LL designed the study and revised the manuscript. HQT, LT drafted the manuscript and prepared tables and figures. HQT, LT and SCS performed comparative genome analysis. YLP performed phylogenetic analysis. All authors discussed the results and commented on the manuscript.

Declaration of competing interest

The authors declared no conflict of interest.

Acknowledgements

We are grateful to Wenbin Yu, Yu Song, from Xishuangbanna Tropical Botanical Garden, Chinese Academy of Sciences, Fei Zhao from Kunming Institute of Botany, Chinese Academy of Sciences, for their help in genomic data analysis. We also thank Lin Li, from South China Botanical Garden, Chinese Academy of Sciences, Xiaoyun Wang and Zhifeng Xu from the Orchid Conservation Center of Yunnan Fengchunfang Biotechnology Co., Ltd, for kindly providing photos and sample materials. This work was supported by the National Natural Science Foundation of China (No. 31870183, No.U1702235) and Southeast Asia Biodiversity Research Institute, Chinese Academy of Sciences (Y4ZK111B01).

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2021.01.003.