RNA editing is a post-transcriptional process that alters the genetic information of mRNAs, so that the translated proteins deviate from those predicted by the genomic DNA templates. RNA editing occurs in a wide range of eukaryotes, including protists, multicellular animals, and land plants (Gray, 2012; Ichinose and Sugita, 2016). RNA editing in animals, dominated by A-to-I synonymous editing with low editing efficiency (20–30%), might play important roles in the regulation of gene expression (Wang et al., 2013). In contrast, RNA editing in plant organelles is usually non-synonymous, affecting mostly the 2nd and 1st codon positions with very high editing efficiency (~80%). In plants, RNA editing is thought to affect plant phenotype and growth by restoring evolutionarily conserved amino acids and ensuring the correct folding of trans-membrane proteins on the respiratory chain complex (Sloan, 2017).

RNA editing occurs in the organellar genomes of all land plant lineages except marchantioid liverworts (Rudinger et al., 2008). Generally, mitochondrial genes have more RNA editing sites than their plastid counterparts. For example, in angiosperms, mitochondrial genomes contain between 300 and 500 RNA editing sites (Edera et al., 2018), whereas plastid genomes contain only between 30 and 50 RNA editing sites (Oldenkott et al., 2014). In addition, the number of RNA editing sites varies in different plant lineages. Lycophytes, which contain the largest numbers of RNA editing sites, possess 2139 RNA editing sites in 18 mitochondrial genes in Selaginella moellendorffii (Hecht et al., 2011), and 400–3100 in plastid genes in Selaginella (Du et al., 2020). Gymnosperms contain 11–418 (average: 245) RNA editing sites in mitochondrial genes, and 3–255 (average: 78) in plastid genes (Liu et al., 2022). Finally, early diverging lineages tend to have more RNA sites, e.g., note the high numbers of RNA editing sites in liverwort Haplomitrium (Rüdinger et al., 2012), lycophyte Selaginella (Hecht et al., 2011), ferns Ophioglossum and Psilotum, gymnosperm Ginkgo (Guo et al., 2016), angiosperms Amborella (Rice et al., 2013) and Liriodendron (Richardson et al., 2013).

RNA editing site abundance is correlated with the diversity of PPR PLS protein genes encoded by the nuclear genome (Rüdinger et al., 2012; Schallenberg-Rüdinger et al., 2014), the heterogeneity of GC content (Guo et al., 2017; Hecht et al., 2011; Malek et al., 1996; Smith, 2009), and substitution rates of genes (Cuenca et al., 2010). It has been well established that among-species heterogeneity of GC content and substitution rate result in systematic errors and affect phylogenetic analyses (Liu et al., 2014). Phylogenetic analyses of plant taxa are likely to be greatly affected as hundreds of RNA editing sites in plant organelles undergo biased conversions between C and T (Edera et al., 2018), introducing GC heterogeneity (Dong et al., 2022). Moreover, because there are 5–10 times more RNA editing sites in mitochondrial than in plastid genomes, the impact of RNA editing sites on phylogenetic analyses based on mitochondrial sequence data should be evaluated with more caution. However, studies on this issue have failed to yield a consensus (as reviewed in Dong et al., 2022). The prevailing opinions regarding the impact of RNA editing on phylogenetic reconstructions are that (1) RNA editing sites do not affect phylogenetic reconstructions (Boweand dePamphilis, 1996) and are reliable sources of phylogenetic information (Schmidt et al., 2002); (2) RNA editing sites in DNA sequences provide additional phylogenetic information, and phylogenomic inferences should be based on DNA sequences rather than on cDNA sequences (Petersen et al., 2006; Picardi and Quagliariello, 2008); and (3) RNA editing sites contain homoplasious signals (Groth-Malonek et al., 2005) and should be modified based on cDNA sequences (Du et al., 2020; Guo et al., 2016; Richardson et al., 2013) or be completely excluded in phylogenetic analyses (Bell et al., 2020).

Gymnosperms are mostly tall and woody trees, comprising about 1000 species in 12 families and 83 genera (Wang and Ran, 2014). Extant gymnosperms are distributed in all continents except Antarctica, and two-thirds of gymnosperms are conifers, constituting over 39% of the world's forests (Armenise et al., 2012). Gymnosperms represent the largest land-based biome, and accordingly are of high ecological and economic significance. Conservation and utilization of gymnosperm germplasm resources requires a well-resolved gymnosperm phylogeny. However, gymnosperm phylogenetic reconstruction has long been ambiguous, especially regarding the unstable placement of gnetophytes (Wang and Ran, 2014; Yang et al., 2022). The gnetophyte lineage has been proposed to be the sister to all seed plants (Soltis et al., 2002), sister to all other gymnosperms (Burleighand Mathews, 2004), sister to cupressophytes (Ran et al., 2010), and sister to Pinaceae (Zhong et al., 2010; Wu et al., 2011; Ran et al., 2018). Recent developments in sequencing techniques and genomic resources have greatly advanced our understanding of gymnosperm phylogeny (Gitzendanner et al., 2018; Ran et al., 2018; Leebens-Mack et al., 2019; Liu et al., 2022) and promise to resolve phylogenetic uncertainties. Indeed, phylogenomic studies have led to two hypotheses regarding the phylogenetic position of gnetophytes: nuclear phylogenomic analyses have indicated that gnetophytes are sister to Pinaceae (Ran et al., 2018), whereas organellar phylogenomic analyses have indicated that gnetophytes are sister to cupressophytes (Zhong et al., 2010). This cyto-nuclear incongruence remains unresolved (Liu et al., 2022).

In our recent study of the cycad genome (Liu et al., 2022), we generated empirical RNA editing data sets (by comparing the DNA and RNA sequences) for gymnosperm representatives, and analyzed both the original DNA data set and a data set with all RNA editing sites modified (based on the RNA sequences). The mitochondrial RNA editing-modified tree recovered the same topology as the original tree, only with slightly lower support for the sister relationship of gnetophytes and cupressophytes, i.e., BPP (Bootstrap Percentage Supports) and GCF (Gene Concordance Factor) decreased from 100% to 41% in original tree to 97% and 33% in the RNA editing-modified tree (Liu et al., 2022). The two data sets yielded the same topology, which might suggest that the impact of RNA editing was not completely mitigated. Here we revisit this empirical editing site data and further examine the impact of RNA editing on phylogenetic analysis. Specifically, we aim to re-evaluate the impact of mitochondrial RNA editing on the phylogenetic reconstruction of gymnosperms and construct a more reliable and robust mitochondrial genomic tree for gymnosperms. For this purpose, we completely excluded all RNA editing sites in our phylogenetic analysis of gymnosperms. Mitochondrial RNA editing data (Table S1), including all the annotated mitochondrial genomes for 46 gymnosperms, were retrieved from the data repository of the cycad genome study (Liu et al., 2022). The mitochondrial genome sequence of Liriodendron tulipifera was downloaded from NCBI organellar database.

Phylogenomic analyses of 46 gymnosperms were based on five mitochondrial genome data sets. The first data set consisted of original protein-coding genes these data included RNA editing sites. The second data set excluded RNA editing sites. The third data set consisted of sequences with degenerated codons (Regier et al., 2010). The fourth data set, called NT12, excluded 3rd codon positions. The final data set consisted of only RNA editing sites. For each of these data sets, we included 40 conserved protein-coding genes: atp1, atp4, atp6, atp8, atp9, ccmB, ccmC, ccmFC, ccmFN, cob, cox1, cox2, cox3, matR, nad1, nad2, nad3, nad4, nad4L, nad5, nad6, nad7, nad9, rpL2, rpL5, rpL16, rps1, rps2, rps3, rps4, rps7, rps10, rps11, rps12, rps13, rps14, rps19, sdh3, sdh4, and tatC.

Each protein-coding gene was extracted and aligned with MAFFT v5.0, with the '-localpair –maxiterate 1000' parameter, implemented in Geneious v.10.0.2 (http://www.geneious.com). The individual gene alignments and the corresponding RNA editing site annotation table were exported. Then an in-house python script was used to remove and retrieve all columns containing at least one RNA editing site annotation from the original data set. FASconCAT-G (Kück and Longo, 2014) was used to produce the codon degenerated and the NT12 data set, as well as to concatenate all the individual matrices into the concatenated super matrices for each data set.

Each concatenated data matrix was analyzed separately using PartitionFinder (Lanfear et al., 2012) for the optimal partitioning scheme and corresponding best-fit models, and subjected to maximum likelihood (ML) phylogenetic reconstruction using RAxML v.8.2.12 (Stamatakis, 2014) with 100 bootstrap replicates. The phylogenetic trees were visualized and modified with FigTree v.1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/).

In total, 11, 486 RNA editing sites were identified in the 40 conserved mitochondrial protein-coding genes of the 46 gymnosperms (Table S1). The number of mitochondrial RNA editing sites of gymnosperms is abundant, ranging from 11 to 418 sites, and averaged 245 sites per species. The abundance of RNA editing sites appears to be related to the function of genes, ranging from 7 sites in sdh4 to 147 sites in nad4, with an average of 63 sites per gene. Genes encoding the mitochondrial respiratory chain complex have the highest number of RNA editing sites, whereas genes encoding ribosomal proteins have the fewest RNA editing sites (Fig. 1), indicating that RNA editing may be important for plant respiration and energy metabolism.

To create a data set that excluded RNA editing sites, we removed 2537 sites (or 6.49% of the total length of the protein-coding sequence alignment) that contained at least one RNA editing site from the original protein-coding gene matrix of the 40 conserved protein-coding gymnosperm genes (Fig. 1). This alignment without RNA editing sites consisted of 36, 531 sites. The alignment of degenerated codons consisted of 39, 068 sites; the alignment that excluded third codon positions consisted of 26, 059 sites; and the alignment of RNA editing sites contained 2537 sites. GC content was highest in the data set that consisted only of RNA editing site sequences (41.29%), followed by the data set of original protein-coding gene sequences (28.97%). The ratio of parsimony-informative sites was highest in the RNA editing site sequences (77.57%), followed again by the original protein-coding gene sequences (23.78%). The data set with degenerated codons contained the lowest ratio of parsimony informative sites with only 13.34%.

The ML trees inferred from three data sets (i.e., the original sequence data set, NT12 data set, and codon degenerated data set) collectively supported the sister relationship of gnetophytes and cupressophytes with maximum supports (Fig. 2a, c, and d). This relationship was also strongly supported by the ML tree recovered from a data set that consisted of only RNA editing site sequences (Fig. S1). However, after RNA editing sites were excluded from the matrix, the mitochondrial ML tree (Fig. 2b) strongly supported the sister relationship of gnetophytes and Pinaceae with high bootstrap percentages (BPP = 90%). This result is consistent with recent nuclear phylogenomic results (Ran et al., 2018; Liu et al., 2022). Thus, the removal of RNA editing sites provide evidence from the mitochondrial genome for the sister relationship of the gnetophytes and Pinaceae. Nevertheless, excluding RNA editing sites did not change the sister relationship of cycads to the rest of gymnosperms, even though the sister relationship of cycads and Ginkgo has been recovered by nuclear and plastid data sets (Liu et al., 2022). This finding suggests that RNA editing may not be the sole cause of cyto-nuclear incongruences in gymnosperm phylogentic analyses.

Cyto-nuclear incongruences on the backbone phylogeny of gymnosperms mainly pertain to the placement of gnetophytes and cycads (Liu et al., 2022). Several approaches have been proposed to reconcile the incongruences related to the position of gnetophytes, including the curation of data biases (Burleighand Mathews, 2004; Wu et al., 2011; Zhong et al., 2010), the tree inference method (Huelsenbeck, 1995) and substitution models (Cox et al., 2014). However, the incongruence of the position of cycads might reflect ancient hybridization events. In extant cycads, Ginkgo, and gnetophytes both plastids and mitochondria are maternally inherited, whereas conifers have unstable maternal and paternal inheritance systems for mitochondria despite a consistent paternal inheritance system for plastids (Owens and Wilson, 1999). The complex cytoplasmic inheritance system in gymnosperms provides the possibility of cycads' ancestor hosting maternal plastid and paternal mitochondria, which may explain the mitochondrion-plastid conflict that obscures the phylogenetic position of cycads.

Nuclear and organellar RNA editing systems have opposite effects on protein diversity (Sloan, 2017). In animals and fungi, nuclear RNA editing increases protein diversity and creates novel protein types (Sloan, 2017). In plants, organellar RNA editing, which is functionally important, restores conserved amino acid sequences (Knoop, 2011; Bentolila et al., 2012; Sloan, 2017; Edera et al., 2018; Dong et al., 2019). Substitution of these edited sites, which is biased toward frequent changes from C and T (Edera et al., 2018), strongly deviate from neutral evolution and likely mislead phylogenetic inferences (Dong et al., 2022). An additional way in which RNA editing sites may affect phylogenetic reconstructions is that the abundance of these sites is positively correlated with GC content (Dong et al., 2019) and negatively correlated to substitution rates (Cuenca et al., 2010). Thus, RNA editing sites should be carefully considered in phylogenetic reconstructions.

In our study, GC content and the percentage of parsimony informative sites were both high in a data set that consisted of only RNA editing site sequences. The phylogenetic tree inferred from this data set recovered significantly long branches and aberrant relationships, indicating that RNA editing sites skew phylogenetic inferences. Although some studies have proposed revision of RNA editing sites rather than exclusion prior to phylogenetic reconstruction, this treatment might be too mild to thoroughly mitigate the impact of RNA editing, especially in plant groups containing large number of editing sites like gymnosperms and lycophytes (Richardson et al., 2013; Du et al., 2020; Dong et al., 2022; Liu et al., 2022). For example, in our study, revision of RNA editing sites generated a phylogenetic tree with the same topology as that based on RNA editing sites only. In fact, the tree topologies of phylogenies generated from several data sets (i.e., original sequences, NT12 data set, and sequences with degenerated codons), including revised RNA editing site sequences (Liu et al., 2022), did not differ. However, excluding RNA editing sites from the mitochondrial protein-coding gene data set restored the sister relationship of gnetophytes and Pinaceae.

Our findings are a reminder of the impact of RNA editing on plant phylogenetic reconstructions. Specifically, our work highlights how organellar genome sequences used as molecular markers in phylogenetic analyses must be carefully evaluated, taking into consideration the effect of RNA editing sites, especially in lineages with many RNA editing sites, like gymnosperms, lycophytes, and hornworts. We recommend that efforts to fully mitigate the effects of RNA editing in phylogenetic analysis exclude, instead of revise, RNA editing sites from sequence data.

Acknowledgements

This work is supported by the Shenzhen Urban Management Bureau Fund (No. 202106 and 202302) and Fairy Lake Science Fund (No. FLSF-2021-02).

Data availability

The mitochondrial genes and RNA editing site annotation for 46 gymnosperms were downloaded from the figshare repository, the matrices and trees newly generated for current study were deposited in figshare repository https://doi.org/10.6084/m9.figshare.20485698. Other supporting results are included within the article.

Author contributions

S.D. and Y.L. designed the study. S.D. and X.Z. carried out bioinformatics and phylogenetic analyses. S.D. drafted the manuscript. Y.L. and T.P. modified the final manuscript, and all authors reviewed it.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

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