Fresh leaves were collected from Cautleya gracilis (99.70°E, 24.18°N), Rhynchanthus beesianus (99.50°E, 22.48°N), Pommereschea lackneri (101.23°E, 21.99°N), Hedychium coronarium (planted variety, 102.72°E, 25.05°N), and H. villosum (101.23°E, 21.99°N) in Yunnan, China, and 45G sequence data were generated for each species using the Illumina Hiseq 2500 platform (San Diego, CA, United States). A total of 277,483,161, 691,955,913, 631,731,352, 309,816,484, and 309,816,484 reads were generated for C. gracilis, R. beesianus, P. lackneri, H. coronarium, and H. villosum, respectively. GetOrganelle was used to execute the de novo assembly of the five chloroplast genomes (− R 15 − k 105,121; Jin et al., 2020), and several previously reported chloroplast genomes from the members of the Zingiberaceae were used as references for automatic annotation and manual adjustment, which were performed using GeSeq and DOGMA, respectively (Wyman et al., 2004; Michael et al., 2017). To ensure accuracy, the coding sequences were further confirmed by online BLAST searches in NCBI. Finally, a circular map of each annotated complete chloroplast genome was drawn using Organellar Genome DRAW (Lohse et al., 2007).

A total of 13 representative chloroplast genomes, including the five newly sequenced ones, were aligned using the Mauve plugin (Darling et al., 2004) in Geneious R8 (Biomatters Ltd., Auckland, New Zealand), with the default parameters to detect inversions and rearrangements. As the chloroplast genome borders of different species typically exhibit varying degrees of contraction and expansion, SC/IR boundary maps and sequence differences were plotted according to the length differences of the four regions and the distribution of related genes.

Even though chloroplast genomes are relatively conserved, structural differences and internal mutations exist between species. To determine the sequence variation of protein-coding genes, we aimed to identify potential DNA barcode genes that may be available in the future. Protein-coding sequences were aligned using MAFFT version 7.308 (Standley, 2013), and genome divergence and variation hotspots were identified using mVISTA (Frazer et al., 2004). Finally, nucleotide diversity (π) was calculated through sliding window analysis using DnaSP version 5 (Librado and Rozas, 2009), with a window length of 600 bp and step size of 50 bp.

Mean amino acid usage frequency was mapped using Circos version 0.69 (Krzywinski et al., 2009), and amino acids were calculated using Geneious R8 (Biomatters Ltd.). To calculate rates of synonymous (Ks) and non-synonymous (Ka) substitution and their ratio (Ka/Ks), the nucleotide sequences of protein-coding genes shared among the four species (C. gracilis, R. tibetica, P. lackneri, and R. beesianus) were extracted and aligned separately using MAFFT version 7.308. Before calculation, gaps and stop codons between the compared sequences were removed. As the YN model considers sequence evolution characteristics (e.g., transition/transversion ratio and codon usage frequency), it has been used increasingly in molecular evolution research (Yang and Nielsen, 2000; Zeng et al., 2017; Zhang R. T. et al., 2020). Thus, the YN algorithm was chosen in KaKs_calculator (Zhang et al., 2006) to illustrate the Ka/Ks value and perform selective pressure analysis. Genes with evidence of positive selection (Ka/Ks > 1) along each branch were identified using the improved branch-site model in PAML (Yang, 2007). The targeted branch(es) was assigned as the foreground branch and the remains served as background branches (Zhang et al., 2005). Finally, a likelihood ratio test (LRT) was used to compare a model (model = 2, NSsites = 2, omega = 1, fix_omega = 0) of positive selection on the foreground branch with a null model (model = 2, NSsites = 2, omega = 1, fix_omega = 1), where no positive selection occurred on the foreground branch. The LTR and corresponding P values were calculated using the chi-squared module in PAML.

Previous studies have suggested that chloroplast RNA editing can improve transcript stability, contribute to the regulation of chloroplast gene expression, and enable genes to produce multiple protein products, thereby expanding the original genetic information (Hanson et al., 1996). To investigate the role of RNA editing mechanisms in the evolution of the Zingiberaceae, PREP-cp (Mower, 2009) was used to predict RNA editing sites, with a parameter threshold (cutoff value) of 0.8 to ensure prediction accuracy.

The Zingiberaceae phylogeny was reconstructed using the chloroplast genome (whole genome or protein-coding only) and internal transcribed spacer (ITS) sequences. In addition to the five newly sequenced chloroplast genomes, other chloroplast genomes and all ITS sequences were downloaded from the NCBI database (Supplementary Table 1). In total, 47 chloroplast genomes and 54 ITS sequences, which each represented 20 genera were selected and aligned using MAFFT. Sequences from species in the Costaceae and Musaceae were also obtained for use as ingroups and outgroups, respectively. Modeltest version 3.7 (Posada and Crandall, 1998) was used to determine the best fitting model, based on Akaike Information Criterion (AIC) score (David and Buckley, 2004). Maximum-likelihood (ML) phylogenetic analysis was conducted using RAxML version 8 (Alexandros, 2014), with 1,000 bootstrap replicates, and Bayesian inference (BI) analysis was performed using the Markov Chain Monte Carlo (MCMC) algorithm in MrBayes version 3.2 (Ronquist and Huelsenbeck, 2003), with 1,000,000 generations and sampling once every 1,000 generations. The first 25% of trees from all runs were discarded as burn-in, and the remaining trees were used to construct a majority-rule consensus tree.

The five newly sequenced chloroplast genomes (162,878–163,831 bp, 36.0–36.1% GC content) possessed the typical quadripartite structure, including an LSC region (87,918–89,237 bp, 33.8–33.9% GC content), SSC region (15,707–16,720 bp, 29.3–29.6% GC content), and a pair of IR regions (IRa and IRb; 28,994–29838 bp, 41.0–41.4% GC content). Except for C. gracilis, which was missing the rps19 gene, the chloroplast genomes contained 133 genes, including 87 protein-coding genes, eight ribosomal RNA genes, and 38 tRNA genes (Supplementary Table 2). Of the 133 genes, 15 (atpF, petB, petD, ndhA, ndhB, rpoC1, rps16, rpl2, rpl16, trnA-UGC, trnI-GAU, trnV-UAC, trnL-UAA, trnG-UCC, and trnK-UUU) contained a single intron and 3 (rps12, clpP, and ycf3) contained two introns. The annotated complete chloroplast genome sequences were deposited in NCBI (GenBank accession numbers: MW769779–MW769783). Meanwhile, the lengths and GC contents of chloroplast genomes from all Zingiberaceae taxa (13 species and 12 genera) ranged from 161,920 bp (Alpinia pumila) to 164,068 bp (Wurfbainia longiligularis; Figure 1 and Supplementary Figure 1) and from 36.0 to 36.2%, respectively. More specifically, the lengths and GC contents of the LSC regions ranged from 86,982 bp (Curcuma amarissima) to 89,237 bp (C. gracilis) and from 33.7 to 34.0%, whereas those of the SSC regions ranged from 15,317 bp (A. pumila) to 16,720 bp (R. beesianus) and from 29.2 to 30.0%, and those of the IR regions ranged from 28,994 bp (R. beesianus) to 30,117 bp (Stahlianthus involucratus) and from 40.9 to 41.4% (Supplementary Table 3).

Moreover, variation at the SC-IR boundary and contraction and expansion were observed (Figure 1). The rpl22 and rps19 were located at the LSC-IRb junction, and ycf1 and ndhF were located at the SSC-IRb junction. In R. beesianus, P. lackneri, H. villosum, and H. coronarium, the rpl22 gene crossed the LSC-IRb boundary, with 52, 41, 53, and 53 bp located in the IRb region, respectively. Interestingly, in C. gracilis, the rps19 gene, which was represented by a copy in both the IRa and IRb regions of the other genomes, was only represented by a single copy in the LSC region. In P. lackneri and S. involucratus, the ndhF gene in crossed the SSC-IRb boundary, with 39 and 14 bp in the IRb region, respectively, the ycf1 gene crossed the SSC-IRa boundary in all 13 chloroplast genomes, with variable sequence lengths in the SSC region. The IRa-LSC boundary was relatively stable, except that the C. gracilis genome lacked an rps19 gene (Figure 1). No gene rearrangements or inversions were observed (Supplementary Figure 2).

Sequence mutates indicated that the chloroplast genomes of Zingiberaceae taxa were highly conserved (Figure 2). The coding regions were more conserved than the non-coding regions, and the IR regions were less variable than the single-copy regions. Four protein-coding regions (psbM, rps12, rpl22, and ycf1), which possessed > 25% variability (Supplementary Figure 3) could be used for DNA barcode research in the future.

Leucine (10.3%), isoleucine (8.8%), and serine (7.9%) were the most frequently used amino acids, whereas cysteine (1.1%) and tryptophan (1.7%) were the least frequently used amino acids (Supplementary Figure 4 and Supplementary Table 4). The nucleotide diversity of the four montane taxa was ∼0.01 (Supplementary Figure 5).

As some genes yielded Ks values of 0, which resulted in invalid Ka/Ks ratios, only 49 genes were included in the Ka/Ks analysis. KaKs_calculator suggested that four genes (atpF, rpoA, rps15, and ycf2) possessed Ka/Ks ratios of > 1 in at least one pairwise comparison among the four montane taxa (Figure 3). The genes atpF and rpoA were detected in P. lackneri and R. beesianus, respectively, whereas rps15 was detected in R. tibetica and P. lackneri. The gene ycf2 was detected in C. gracilis and R. beesianus. Further verification of the branch-site model revealed that the P-values of the targeted branches (rpoA and ycf2) were significant and retrieved sites under positive selection using the Bayes Empirical Bayes (BEB) method (Supplementary Table 5).

A total of 76–81 RNA editing sites were predicted in 25–27 genes (Supplementary Table 6). The ndhB gene contained the most predicted editing sites (9–11), which is consistent with findings in other plants, such as rice, maize, and tomato (Freyer et al., 1995). Meanwhile, ndhD contained 7–9 predicted editing sites, whereas ndhF contained 5–7 predicted editing sites, and the other genes contained between 0 and 7 predicted editing sites (ndhA, 4–7; rpoB, accD, 4–5; ycf3, 4; rpoC2, matK, 3–5; rpl20, rpoA, rps14, 3; ndhG, 2–3; petB, rpoC1, 2; atpB, atpI, psbB, rps16, 1–2; atpA, atpF, ccsA, psbF, rps8, 1; clpP, rpl2, rps2, 0–1). All predicted editing sites were C-to-U transitions, and most of the editing sites were predicted to greatly increase protein hydrophobicity but maintain the original function. While maintaining stability, it also provided a basis for adapting to different environments. More work is needed in this area in the future.

Pommereschea, Rhynchanthus, Cautleya, Roscoea, and Hedychium formed a monophyletic clade in the chloroplast genome tree, with BI support of 0.8, and the taxa were also closely related in the ITS tree (Figure 4). In both trees, the sister relationship of Pommereschea and Rhynchanthus was strongly supported (100% ML support and 1.0 BI support), and Roscoea was closely related to the Pommereschea-Rhynchanthus clade in the chloroplast genome tree, and the sister relationship of Cautleya and Roscoea was strongly supported in the ITS tree (100% ML support and 1.0 BI support).

Previous studies have reported that the chloroplast genomes of herbaceous plants have undergone rapid evolution, with certain structural changes, such as inversions (Doyle et al., 1992) and gene losses (Takayuki et al., 2004; Saski et al., 2005). No inversions or gene rearrangements were detected in the chloroplast genomes of the Zingiberaceae taxa included in this study. However, although most angiosperms, including most members of the Zingiberaceae, possess two copies of the rps19 gene at the boundaries of the LSC and IR regions (Xu et al., 2015), the Cautleya chloroplast genome only contained a single copy of the rps19 gene in the LSC region. Changes in rps19 genes have been reported in several other genera, including Dianthus (Caryophyllaceae; Raman and Park, 2015), Cardiocrinum (Liliaceae; Lu et al., 2016), Prunus (Rosaceae; Zhao et al., 2019), and Colobanthus (Caryophyllaceae; Androsiuk et al., 2020). However, the changes observed in the rps19 copies of Cautleya were different from those reported in other genera in two respects. First, the rps19 copy in the IRa region of Cautleya was completely lost, whereas those in the IRa regions of other genera were reportedly shortened and pseudogenized. Second, the rps19 gene in the IRb region of Cautleya was located in the LSC region, whereas in other taxa, the rps19 gene remained in the IRb region.

The rps19 protein is a component of the 40S small ribosomal subunit and is essential to both the maturation of the 3′-end of 18S rRNA and the assembly and maturation of pre-40S particles, which are related to chloroplast transcription and translation (Soulet et al., 2001; Matsson et al., 2004). The loss of rps19 has also been observed in a few other dicot taxa (e.g., Morus, Nicotiana, Vitis, and Tetrastigma) but is relatively rare in monocots (Ravi et al., 2006; Li et al., 2015), which suggests that rps19 is more likely to be lost or pseudogenized in dicots. The changes in rps19 could be due to (1) partial gene duplication (Lu et al., 2016; Zhao X. et al., 2019) or (2) the contraction and expansion of IR regions (Zhao X. et al., 2019). It was suggested that there are two evolutionary mechanisms of the IR region boundary: the small amplitude amplification of the boundary gene and the recombination repair of the boundary of the LSC region. The former is an important factor for maintaining the stability of IR regions (Goulding et al., 1996). The expansion and contraction of chloroplast IR regions are relatively common (Hansen et al., 2007). Except for Cautleya, other Zingiberaceae taxa included in this study possessed two complete rps19 copies, which suggests that the presence of two copies is the ancestral state within the Zingiberaceae. Cautleya also possesses the longest LSC region among the included taxa, which suggests that large changes in the rps19 of Cautleya should be the result of LSC region expansion and repair. Previous studies have suggested that rps19 cannot be completely removed from the IRa region through the expansion of LSC or IR regions (Raman and Park, 2015; Lu et al., 2016; Zhao X. et al., 2019; Androsiuk et al., 2020). Therefore, the complete loss of rps19 in Cautleya is more likely than the suppression of rps19 duplication by the LSC region expansion.

The rpoA and ycf2 genes are commonly associated with positive selection, which suggests that the chloroplast genomes of Cautleya, Roscoea, Rhynchanthus, and Pommereschea have undergone adaptive evolution. Notable adaptive divergence was noted for rpoA in the chloroplast genomes of the sister genera Rhynchanthus and Pommereschea. The rpoA gene encodes the α subunit of plastid-encoded RNA polymerase, which is responsible for the expression of most genes involved in photosynthesis and is essential for chloroplast gene expression and chloroplast development (Purton and Gray, 1989; Hajdukiewicz et al., 2014; Zhang et al., 2018). The evolution of rpoA is complicated in angiosperms. In the Annonaceae, Passifloraceae, and Geraniaceae, rpoA divergence was caused by structural rearrangement and purifying selection (Blazier et al., 2016). In Passiflora (Passifloraceae), rpoA is subject to either positive or purifying selection, depending on the specific clade (Shrestha et al., 2020). In Rehmannia (Orobanchaceae), rpoA is under positive selection (Zeng et al., 2017). In this study, Rhynchanthus, members of which are typically epiphytic on limestone or tree trunks in forest understories at lower elevations, when compared with Pommereschea. Habitat differentiation, in regard to sunlight exposure, suggests that these sister genera have experienced selection based on the utilization of different light intensities.

In angiosperms, ycf2 is the largest chloroplast gene (Huang et al., 2010) and is subjected to positive or purifying selection (Yan et al., 2019; Zhong et al., 2019). Even though previous studies have suggested that ycf2 has been lost from the chloroplast genomes of monocots (Drescher et al., 2000; Wang et al., 2018; Mishra et al., 2019), two ycf2 copies were present in the chloroplast genomes of the Zingiberaceae taxa included in this study. Furthermore, even though the specific function and role of ycf2 remain unclear, studies have suggested that the gene is not essential to either photosynthesis (Drescher et al., 2000; Zhang Y. et al., 2020) or leaf patterning and is, instead, related to cell survival and possibly ATPase metabolism (Kikuchi et al., 2018; Wang et al., 2018; Zhang et al., 2018; Zhang Y. et al., 2020). The ycf2 gene was also reported to contribute to encoding the 2-MD AAA-ATPase complex, which is a motor protein for generating ATP required for inner membrane translocation (Kikuchi et al., 2018), and to plant cell survival (Drescher et al., 2000). The positive selection of ycf2 suggests that the gene is involved in the adaptive evolution of the montane investigated here.

Even though the chloroplast-based Zingiberaceae phylogeny reconstruction was strongly supported and consistent with previous systematic studies (Kress et al., 2002), the phylogenetic positions of Cautleya and Roscoea in the chloroplast genome and ITS trees were inconsistent. Hybridization and incomplete lineage sorting are the most likely factors to underly phylogenetic conflict between nuclear and chloroplast genome signals (Degnan and Rosenberg, 2009; Joly et al., 2009; Petit and Excoffier, 2009). For example, Roscoea could be a hybrid descendant of Cautleya and the ancestor of Rhynchanthus and Pommereschea (Figure 4). However, incomplete lineage sorting is also possible because incomplete lineage sorting could be present at deeper-divergence lineages in angiosperms (Yang et al., 2020). Either way, this study confirmed the close phylogenetic relationships of the genera Pommereschea, Rhynchanthus, Cautleya, and Roscoea.