1. Introduction

In photosynthetic plants, the plastid genome (plastome) is typically conserved in size, structure, and gene content (Wicke et al., 2011; Mower and Vickrey, 2018). However, in response to altered habitats and changes in external nutrient sources, the plastome may experience mild to severe functional disruption by progressive loss of genes in different functional categories due to altered selective regimes (Barrett and Davis, 2012; Graham et al., 2017; Kim et al., 2020). The early stage of genome degradation often includes the pseudogenization and/or loss of NADH dehydrogenase-like (ndh) genes, which are proposed to function as protection against photo-oxidative stress by balancing cellular redox levels (Martín and Sabater, 2010). The ndh genes act to relieve the effects of frequent episodes of environmental stresses but could be dispensable under favorable environments (Ruhlman et al., 2015).

Compared to the evolution of plastid genome architecture, relatively little is known about the evolutionary dynamics of posttranscriptional modifications. Plastid RNA editing, which mostly converts cytidine to uridine (C-to-U) at the posttranscriptional level, is widespread in land plants (Freyer et al., 1997; Small et al., 2020). Plastid transcripts in early land plants are generally heavily edited, yet most ancestral edits tend to be lost in higher plants (Tillich et al., 2006; Small et al., 2020). In angiosperms, RNA editing can affect dozens to > 180 sites in the plastome (Hein et al., 2016; Ishibashi et al., 2019; Smith, 2020). Comparisons of homologous editing sites have indicated that many sites are conserved in different angiosperm clades (Tsudzuki et al., 2001; Tillich et al., 2006). Slight differences in editing positions have also been observed in Arabidopsis ecotypes (Tillich et al., 2005). However, the evolutionary rates of plastid RNA editing sites may also be rapid. For example, plastid RNA editing is extremely high (up to ~ 3500 sites) and extraordinarily variable (> 6-fold change) in a lycophyte genus, Selaginella (Oldenkott et al., 2014; Smith, 2020).

How is the frequency of plastid RNA editing influenced by external and internal stimuli? Do sites of posttranscriptional modifications evolve in tight association with genetic variation? We answered these questions by examining the evolution of plastid genome architecture and RNA editing profiles for 11 orchids from Cymbidium, which comprises approximately 52 recognized species that are mainly distributed in Asian tropical and subtropical regions and northern Australia (Du et al., 2007; Yang et al., 2013). Cymbidium species have diverse life forms; ~70% are epiphytic, whereas a minority are terrestrial or lithophytic, and there is one partial mycoheterotroph, Cymbidium macrorhizon Lindl (Motomura et al., 2008; Kim et al., 2018; Suetsugu et al., 2018). These differences in life forms probably reflect changes in ecological adaptations to photosynthetic fluctuations and altered nutrient sources. Epiphytic plants are better able to compete for sufficient light when living on trees but need to balance transpirational water loss and photosynthetic rates. Terrestrial and mycoheterotrophic orchids compensate for their limited photosynthesis via increased fungal dependence. It has been reported that mycoheterotrophic orchids obtain 30–50% more carbon than terrestrial orchids to make up part of the nutrients needed for growth (Motomura et al., 2010). Thus, the diversity of Cymbidium plants provides an excellent opportunity to explore the genomic and posttranscriptional variation that may be associated with lifestyle-associated ecological adaptations.

2. Materials and methods 2.1. Taxon sampling

We sampled 11 Cymbidium species (Cymbidium dayanum Rchb. f., Cymbidium ensifolium (L.) Sw., Cymbidium faberi Rolfe, Cymbidium floribundum Lindl., Cymbidium lancifolium Hook., Cymbidium lowianum Rchb. f., C. macrorhizon, Cymbidium mannii Rchb. f., Cymbidium sinense (Jackson ex Andr.) Willd, Cymbidium tracyanum L. Castle, Cymbidium tortisepalum Fukuyama) with a diversity of lifestyle strategies (epiphytic, terrestrial and mycoheterotrophic). Stems and flower tissues were collected for the leafless mycoheterotroph C. macrorhizon, and fresh leaves were collected for other orchids. Sampling and sequencing details are included in Table S1.

2.2. Assembly and annotation of plastomes

Genome sequencing, assembly and annotation were performed as described previously (Fan et al., 2019; Liu et al., 2020). Briefly, total cellular DNA was isolated using the cetyltrimethylammonium bromide (CTAB) procedure (Doyle and Doyle, 1987), and sequenced on an Illumina HiSeq X Ten machine, generating 1.2–68.6 Gb of 150 bp paired-end reads per sample. Plastomes were assembled with NOVOPlasty v.3.2 (Dierckxsens et al., 2017) and annotated using the Geneious v.9.1.3 trial version (Kearse et al., 2012) with manual adjustment. The assembled plastomes are deposited in GenBank under accession numbers: MW582681–MW582691.

2.3. Strand-specific RNA-Seq library construction and transcriptome sequencing

To efficiently capture plant organelle-derived transcripts that generally lack poly(A) tails, we adopted the Ribo-minus RNA sequencing strategy. Briefly, total RNA was isolated using TRIzol reagent (Invitrogen) and treated with Ribo-Zero™ rRNA Removal Kit (Epicentre, Madison, WI, USA) for rRNA depletion and then with DNase I (TaKaRa Bio Inc., Japan) to ensure the removal of trace amounts of DNA. Strand-specific RNAseq libraries were constructed from purified mRNA using random hexamers, generating approximately 12 Gb of 150 bp paired-end reads for each species (Table S1). All raw sequencing data have been deposited in Sequence Read Archive under accession numbers SRR13674023–SRR13674033 and SRR13734467–SRR13734477.

2.4. Identification of editing sties

Raw RNAseq reads were filtered using Fastp v.0.20.1 (Chen et al., 2018) and then aligned to the respective plastomes with bowtie2 v.2.1.0 (Langmead and Salzberg, 2012). Expression signals were estimated via read coverages and calculated using Bedtools (Quinlan and Hall, 2010). RNA editing sites were detected using VarScan v.2.4.4 (Koboldt et al., 2012) and visually examined with Tablet v.1.19.09.03 (Milne et al., 2013) as described previously (Fan et al., 2019).

2.5. Phylogenetic construction and ancestral states

A 104.4-kb alignment of plastome sequences (including sequences of single copy regions and one copy of the inverted repeat) was performed using MAFFT v.7.475 (Katoh and Standley, 2013) that included all 11 Cymbidium species, along with Erycina pusilla (L.) N.H. Williams & M.W. Chase (NC_018114) and Oncidium sphacelatum Lindl. (NC_028148) as outgroups. A maximum likelihood-based plastid phylogeny was constructed using RAxML v.8.2.12 (Stamatakis, 2014) with a GTRGAMMAI substitution model and 100 bootstrap replicates. Ancestral states of genes and RNA editing were estimated using Dollo parsimony implemented in Count (Csűös, 2010).

3. Results 3.1. The trajectory of Cymbidium plastome evolution

To assess the evolutionary dynamics of plastid genome features, which may be associated with variation in lifestyles and other photosynthetically-related traits among Cymbidium orchids (Table S1), 11 plastomes were sequenced and assembled. All sequenced plastomes are largely syntenic in structure (Fig. S1). In general, plastomes are slightly larger in epiphytes compared with terrestrial and mycoheterotrophic orchids (Fig. 1), with genome sizes ranging from 149, 728 bp in C. macrorhizon (a mycoheterotrophyte) to 156, 616 bp in Cymbidium tacyanum (an epiphyte). Cymbidium plastomes encode 70–78 functional protein-coding genes. The differences in total gene numbers are due to different retention, loss and pseudogenization of the 11 ndh genes. All other typical plastid protein-coding genes are intact in Cymbidium.

The extent of ndh gene retention and loss, including pseudogenization via point mutations, fragmentation and gene deletions, varies among Cymbidium species (Fig. 1; Table S2). A previously documented 1-bp insertion at the 5′-end of ndhB (Kim et al., 2018) is shared with all 11 sampled Cymbidium species, and this frameshift might represent the initial degradation of the NDH complex in these genes. Other than the frameshifted ndhB, epiphytic Cymbidium orchids have generally retained nearly all plastid ndh genes, while most ndh genes have become progressively pseudogenized or lost in terrestrial and mycoheterotrophic taxa. Moreover, consistent with previous studies (Luo et al., 2014; Kim et al., 2018), ndh genes were also extensively pseudogenized in C. mannii, an epiphyte with strong crassulacean acid metabolism (CAM photosynthesis) (Zhang et al., 2016). Notably, the pseudogenized ndh genes are dispersed along the plastome (Fig. S1), indicating multiple mutational hotspots.

3.2. Extensive loss of plastid C-to-U RNA editing in Cymbidium plastomes

To further investigate plastome evolutionary dynamics at the posttranscriptional level, we performed non-polyA enriched transcriptome sequencing to efficiently capture nascent and mature organellar RNAs. Interestingly, RNAseq read mapping revealed that nearly all of the plastome is transcribed (mean depth = 2.0–5.9 × 10⁴), with the regions closest to ycf1 and ycf2 showing the lowest expression signals relative to others (Fig. 2). Moreover, intergenic regions and pseudogenes have relatively high expression levels, presumably due to co-transcription with neighboring genes.

The deep levels of expression obtained for most of the plastome provides ample data for the systematic identification of plastid RNA editing events. A total of 57–100 genomic sites are C-to-U edited in Cymbidium plastome transcripts. C. macrorhizon comprises the smallest set (57) of edits while epiphytic Cymbidium orchids (excluding C. mannii) generally have ≥90 edited sites. As many as 25–39 functional protein-coding genes are edited at one or more sites (Fig. 1; Table S3), mostly at their 1st and 2nd codon positions (Fig. S2). Interestingly, many pseudogenized ndh genes are still edited to some extent (Fig. 1). For example, of the seven pseudogenized ndh genes in C. macrorhizon, four were edited. A small number of intergenic and intronic edits (12.3–21.0%) were also identified (Table S3).

Comparisons of homologous editing sites revealed that the nearly two-fold differences in editing numbers were mainly due to progressive loss of ancestral edits, and to a much lesser extent, lineage-specific gains (Figs. 3 and S3). Not surprisingly, the ndh genes, both intact and pseudogenized, exhibit the most variation in their editing status in Cymbidium, requiring editing at 2–11 sites in ndhB, 0–3 sites in ndhC, 1–4 sites in ndhD, and 0–4 sites in ndhF (Table S3). Reconstruction of ancestral states indicates that editing of ndh transcripts is dramatically reduced in terrestrial and mycoheterotrophic orchids. Nearly 57–80% of the 30 edits that were probably present ancestrally in coding regions of ndh transcripts are lost in these orchids (Fig. S2). The epiphyte C. mannii, which also exhibits independent loss and pseudogenization of its ndh genes (Fig. 1), has only retained 9 ancestral editing sites in ndh transcripts. In contrast, the epiphytic orchids that have retained mostly functional ndh genes show only mild variation in ndh editing frequency.

Editing in non-ndh genes, such as other photosynthetic genes (pet, psa and psb) and ribosomal protein genes (rpl and rps), generally fluctuate in their editing profiles in a species-specific pattern (Fig. S2). Genes in other functional groups, such as those encoding atp, rpo and ycf proteins, are relatively conserved in RNA editing in all 11 Cymbidium species.

3.3. The extent of plastid RNA editing reduction varies among Cymbidium lineages and phylogenetic depths

Using the inferred ancestral editing status for each internal node, we compared the percentage of RNA editing reduction at different phylogenetic depths (Fig. 3A). The rampant RNA editing reduction of terrestrial and mycoheterotrophic orchids begins prior to the divergence of terrestrial orchids and the epiphytic Cymbidium floriundum, in which nearly 21.4% ndh and 2% non-ndh ancestral edits were lost (Fig. 3A). An additional 36.4% reduction of ndh edits occurred in the common ancestor of all terrestrial and mycoheterotrophic Cymbidium orchids. Subsequent RNA editing reduction is extensive for the mycoheterotrophic C. macrorhizon, while terrestrial orchids generally demonstrate lineage-specific editing reduction.

Among epiphytic orchids, the evolution of RNA editing was found to contrast between C. mannii and other taxa. Nearly 68% of RNA editing in ndh genes in C. mannii were lost, considerably higher than that in non-ndh genes. In contrast, only a small number of editing sites were affected in other epiphytic orchids, and these affected both ndh genes and non-ndh genes.

3.4. Loss of plastid editing is mainly due to reduced editing efficiency

The reduction of editing sites could not be simply explained by gene deletion. In fact, sequence deletion only eliminates the editing of six genomic sites (one in ndhA, one in ndhB and four in ndhF) (Fig. S2). To examine the cause of editing loss, we examined the nucleotide states at each homologous editing site (Fig. 3B). Nucleotide substitutions (i.e., converting an edited cytosine (C) to another nucleotide) make up only a small fraction (< 15%) of the editing loss at ancestral sites. Instead, the majority of ancestral edit sites that have lost editing have become unedited cytosines. These unedited cytosines account for 95% of loss events for terrestrial and mycoheterotrophic plants, 90% for C. mannii, and 85% for the other epiphytic orchids (Fig. 3B). This pattern indicates that losses of plastid editing sites in terrestrial and mycoheterotrophic Cymbidium orchids and C. mannii plastomes are mainly due to reduced editing efficiency rather than removal of the edited cytosine by deletion or nucleotide substitution. This is further supported by the extent of RNA editing in the transcript pool; ndh, petL, rpl genes tend to have low editing efficiency and reduced editing, whereas other genes (such as ATP synthase genes) have higher editing efficiency and less editing loss (Fig. S3).

4. Discussion

Collectively, these data suggest that throughout the diversification of Cymbidium species, both gene content and RNA editing profiles have changed to a certain extent, particularly in terrestrial and mycoheterotrophic species, which also have broader variation in genome architecture than epiphytic orchids (Fig. 1). This might be explained by increased opportunity for nutrient recovery from mycorrhizal partners when moving from trees to the ground, which ultimately may provide adaptive benefit to the switch to partial or full mycoheterotrophy (Motomura et al., 2010). One hallmark of orchids is their seeds generally lack endosperm and depend on fungi for nutrients during germination and the early stages of seedling growth (Du et al., 2007). Cymbidium plants with different lifestyle forms may have variable nutrient dependence on fungi. While epiphytic plants can absorb sufficient light, the terrestrial and mycoheterotrophic Cymbidium orchids generally grow in shade under the tropical forest all year round, which increases both their chances and their requirements for nutrients from fungal partners. This ecological adaptation could relax the functional constraints on photosynthetic components, typically on ndh genes, which function to protect against photo-oxidative stresses (Martín and Sabater, 2010).

However, we also found independent pseudogenization of most ndh genes in the plastome of the epiphytic C. mannii (Fig. 1). This may provide a good example of different ecological adaptation strategies leading to the convergent reduction of ndh genes. However, the underlying reasons for this reduction of ndh genes in C. mannii are still not well understood. One possibility is that C. mannii has evolved strong CAM photosynthesis, which could reduce transpirational water loss and improve photosynthetic rates by monitoring the opening and closing of stomatal aperture in a day–night fashion (Motomura et al., 2010; Zhang et al., 2016). This circadian regulation of photosynthesis could also relax selection on ndh genes by optimizing cellular oxidative redox. If this is true, we would expect such reduction of ndh genes would be repeatedly observed when more plastomes are sequenced in other strong CAM lineages.

How does the evolution of plastid RNA editing profiles associate with genetic variation in Cymbidium plastomes? To answer this question, we performed transcriptomic sequencing of all sampled Cymbidium taxa with a non-polyA enriched strategy, which efficiently captures early editing status in nascent RNAs (Germain et al., 2013). To our knowledge, our transcriptomic sequencing efforts represent the largest scale comparison of plastid RNA editing in a genus. Our findings provide strong evidence that non-random loss of ancestral edits in ndh genes in terrestrial and mycoheterotrophic species as well as C. mannii mainly accounts for the variation in total RNA editing sites (Fig. 1, Fig. 3), although several other sites may fluctuate in editing levels. Because of this correlated reduction of RNA editing sites in ndh genes accompanied with ndh pseudogenization, it seems that plastid RNA editing provides little or no advantage to counteract genetic variation at the DNA level.

How are RNA editing profiles in terrestrial and mycoheterotrophic orchids and C. mannii affected by relaxed selection? Among plastid genes, RNA editing is relatively abundant in ndh genes (Martín and Sabater, 2010). In our study, it appears that ndh genes continuously accumulated mutations in these lineages, suggesting increased mutation burden for the site recognition of editing machinery. If so, one would expect editing efficiency could be reduced at ancestral editing sites in the ndh genes for these particular lineages. By analyzing the proportion of edited reads within the total read coverage, we found relatively lower editing efficiency in ndh, petL, rpl genes compared with other genes (Fig. S3). Moreover, the contrasting evolution of ndh and non-ndh derived editing sites among Cymbidium plastomes further support this hypothesis (Fig. 3A). In summary, we conclude that the plastid RNA editing system could be fragile and the evolution of RNA editing profiles could provide additional evolutionary footprints left by altered external and internal environments.

Author contributions

A.Z. and J.Y. designed and conceived this study; M.Z., L. Z., F.L, Y.H and W.F performed this experiments and analyses; M.Z and A.Z. wrote the manuscript; funding acquisition, J.Y.

Declaration of competing interest

The authors declare that they have no conflict of interest.

Acknowledgements

We greatly thank Dr. Jia-Lin Huang and Mr. Ji-Dong Ya for sample collection and identification. We also appreciate Dr. Jeff Mower (University of Nebraska) for careful reading and language editing. This work was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (grant no. XDB31010000), by the Large-scale Scientific Facilities of the Chinese Academy of Sciences (grant no. 2017-LSF-GBOWS-02), by an open research project for "Cross-Cooperative Team" of the Germplasm Bank of Wild Species, Kunming Institute of Botany, Chinese Academy of Sciences, and by the CAS Pioneer Hundred Talents Program (to A.Z.).

Appendix A. Supplementary data

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