1. Introduction

RNA editing is a post-transcriptional modification of RNA that occurs in nuclear and/or organellar genomes of some organisms (Maier et al., 1996; Gott and Emeson, 2000). In land plants, the most prevalent type of RNA editing in organellar genomes is cytidine-to-uridine (C-to-U) editing, while U-to-C editing only occurs in some plants, e.g., hornworts, lycophytes, and ferns (Chateigner-Boutin and Small, 2011; Guo et al., 2015; Knie et al., 2016; Ichinose and Sugita, 2017). The frequency of organellar RNA editing in plants varies from zero to thousands of sites (Ichinose and Sugita, 2017; Lenz et al., 2018). Most RNA editing events serve to restore the evolutionarily conserved amino acid residues in mRNAs or to create start and stop codons (Chateigner-Boutin and Small, 2011; Takenaka et al., 2013). Therefore, RNA editing is an essential process to correct defective genes and maintain genetic information (Knoop, 2011; Small et al., 2020). For plant molecular phylogeny, plastid DNA sequences are the primary data source (Gitzendanner et al., 2018), and over the last 30 years, protein-coding sequences or translated amino acid sequences have been the most commonly used sequences (Hasebe et al., 1995; Jansen et al., 2007; Ruhfel et al., 2014; Nie et al., 2020).

The abundance of RNA editing sites in DNA sequences raises the possibility that these sites bias phylogenetic analysis. Some early studies tested the effects of organellar RNA editing on phylogenetic reconstruction by comparing phylogenies generated from DNA and cDNA sequences. These studies generally suggested that RNA editing has no or little effect on DNA-based phylogenetic reconstruction (Bowe and dePamphilis, 1996; Varigerow et al., 1999; Szmidt et al., 2001; Petersen et., 2006, 2013). However, these studies only used a few organellar genes (mostly mitochondrial genes), and some sampled lineages had relatively low RNA editing frequency (Ichinose and Sugita, 2017; Lenz et al., 2018). Therefore, a reevaluation of the effect of RNA editing on phylogenetic reconstruction at a genomic level is essential, especially in highly RNA-edited lineages, such as hornworts, lycophytes, and ferns (Ichinose and Sugita, 2017; Lenz et al., 2018).

The frequency of RNA editing in plastid genomes (3415 editing sites) is highest in Selaginella uncinata (Desv. ex Poir.) Spring (Oldenkott et al., 2014); for mitochondrial genomes the highest frequency of RNA editing (2152 editing sites) is found in another species of Selaginella P. Beauv., Selaginella moellendorffii Hieron (Hecht et al., 2011). In both species, RNAs are exclusively modified from "C" to "U" and editing mainly occurs in protein-coding regions (Hecht et al., 2011; Oldenkott et al., 2014). Selaginella, the single genus in Selaginellaceae (Selaginellales or spike moss), is distributed worldwide and contains about 700–800 species (PPG, 2016; Weststrand and Korall, 2016). All sequenced Selaginella plastomes have a remarkably high GC content (50.7–56.5%) (Zhang et al., 2020), and because high GC content is thought to be correlated with organellar RNA editing, Selaginella species are supposed to have high levels of RNA editing (Malek et al., 1996; Tsuji et al., 2007; Smith, 2009; Hecht et al., 2011; Oldenkott et al., 2014).

Previous studies have shown that RNA editing sites can be verified by comparing genomic DNA sequences with cDNA sequences (Wolf et al., 2004; Oldenkott et al., 2014), or can be predicted by comparing genomic DNA sequences with verified DNA or amino acid sequences among closely related species (Lenz et al., 2018). Prediction of C-to-U RNA editing in coding regions is relatively reliable because most editing events tend to restore conserved codons (Oldenkott et al., 2014; Lenz et al., 2018; Small et al., 2020; Gerke et al., 2020). In this study, we used Selaginella as a research model to test the effect of extreme plastid RNA editing on phylogenetic reconstruction. We first predicted RNA editing sites in plastid protein-coding sequences and then compared phylogenetic reconstructions of Selaginella derived from plastome and RNA-edited DNA/amino acid sequences.

2. Materials and methods 2.1. Taxon sampling and predicting RNA editing sites

We downloaded 18 Selaginella plastome sequences from GenBank (Table 1), and then predicted the RNA editing sites in protein-coding regions of these plastomes using the online plant RNA editing prediction and analysis tool PREPACT 3 (Lenz et al., 2018). Predictions were conducted under BLASTX mode with default parameters, and the plastid protein database of S. uncinata (modified from AB197035.2, 76 protein sequences, 3488 RNA editing records) was selected as reference.

2.2. Data set construction

Predicted editing information of each plastome was transformed into feature file format using notepad++ (https://notepad-plus-plus.org/), and then was added to the original GenBank flat file using Sequin (NCBI). Six lycophyte plastomes (two species of Isoëtaceae or quillworts and four species of Lycopodiaceae or club mosses) were downloaded as outgroups (Table 1). Coding sequences of each gene were extracted using Geneious 8 (Kearse et al., 2012). Genes missing from more than half of sampled Selaginella were excluded in phylogenetic analyses. Coding sequences of each gene were aligned based on translated protein sequences, and obvious aligned gaps (occurred in > 50% samples) were trimmed by codon using Geneious 8 (Kearse et al., 2012). Seventy-five genes were finally aligned and concatenated to produce a coding sequences (CDS) matrix (Table 2). To test the effect of plastid RNA editing on phylogenetic reconstruction, we constructed three contrasting data set pairs, which consisted of the genomic DNA/amino acid sequences and RNA-edited (or corrected) sequences, respectively. Specifically, to produce a contrast CDSr matrix, predicted editing sites "C" were replaced by "T" using Geneious 8 (Kearse et al., 2012) (Table 2). To produce two new data sets named codon 1-2 and codon1r-2r, respectively, first and second codon positions of the CDS and CDSr matrices were extracted using Geneious 8 (Kearse et al., 2012) (Table 2). Finally, the matrices CDS and CDSr were translated into amino acids matrices named AA and AAr, respectively (Table 2), using MEGA 6 (Tamura et al., 2013).

2.3. Phylogenetic analyses

The best-fit substitution models of all the DNA matrices were determined by jModeltest 2 (Posada, 2008) (GTR + I + G), and all DNA matrices were partitioned by codon position. For the amino acid matrices, the best-fit substitution models were automatically determined by the phylogenetic inference software, i.e., substitution model set to PROTGAMMAAUTO. Maximum Likelihood (ML) method, implemented in RAxML 8.2.10 (Stamatakis, 2014), was used to infer the phylogenetic relationships within Selaginella. A thorough tree search for the best ML tree was performed, and bootstrap analyses were performed with 1000 replications; bipartition information from the bootstrap trees was drawn on the best ML tree. Possible correlations between number of plastid RNA editing sites, plastome GC content, CDS GC content, and the GC content of codon 1-2 data sets were tested using R ape package (https://www.r-project.org/). The non-independent effects of phylogenetic signal were corrected using the phylogenetically independent contrasts (PIC) method (Felsenstein, 1985).

3. Results 3.1. Predicted RNA editing sites

The numbers of predicted C-to-U RNA editing sites in protein-coding regions of 18 Selaginella plastomes are shown in Table 1. Detailed information of predicted RNA editing sites is provided in the Appendix. Of the 18 plastomes examined, S. lyallii had the least number of predicted editing sites (405); whereas S. hainanensis, Selaginella pennata, and S. uncinata had the highest numbers of editing sites, ranging from 3075 to 3100 (Table 1). Selaginella uncinata was predicted to have 3079 editing sites, and verified to have 3415 sites by transcriptome data (Oldenkott et al., 2014). Selaginella lepidophylla and S. kraussiana were predicted to have 761 and 1256 editing sites in our study, respectively, and were verified to have 720 and 1353 sites by transcriptomes data, respectively (Smith, 2020).

The characteristics of three data set pairs that were used in phylogenetic analyses are shown in Table 2. The CDS matrix has an aligned length of 64, 452 bp and contains 35, 285 editing sites. The codon 1-2 matrix has an aligned length of 42, 968 bp and contains 34, 889 editing sites.

3.2. Correlations between RNA editing and GC content

The GC content of 18 Selaginella plastomes varied from 50.0% to 56.5% (Table 1). The corrected data sets CDSr and codon1r-2r had relatively higher values of pairwise identity and lower GC content (Table 2). For the CDS and CDSr data set pairs, GC content of ingroups dropped from 51.9% to 48.7% (Table 2). For the codon 1-2 and codon1r-2r data set pairs, GC content of ingroups dropped from 53.9% to 49.1% (Table 2). Correlation analysis of RNA editing sites and GC content in different data sets is shown in Table 3.

The number of RNA editing sites in Selaginella plastomes was not correlated (P > 0.05) with the GC content of plastomes, but highly correlated (P < 0.01) with the GC content of the first and second codon positions (Table 3). For example, the plastomes of the related species Selaginella remotifolia and S. kraussiana were similar in size (131, 867 bp versus 129, 971 bp) and number of RNA editing sites (1474 versus 1256), but had very different GC content (56.5% versus 52.3%). In fact, S. remotifolia had the highest GC content among all 18 sampled Selaginella species, even higher than the most RNA-edited species examined, i.e., S. hainanensis (Table 1). The plastomes of S. lyallii, S. moellendorffii, and S. hainanensis have 405, 2357, and 3100 predicated editing sites, respectively, and the GC content of first and second codon positions are 50.1%, 53.3%, and 56.7%, respectively (Table 3). Interestingly, species of clade D in Selaginella (Fig. 1) had a relatively higher number of editing sites (2267–2866) and lower GC contents (50.2–50.8%), whereas species of clade A in Selaginella (Fig. 1) had a relatively lower number of editing sites (1776–1858) and higher GC contents (54.0–54.1%), although they share very similar GC contents at the first and second codon positions (52.7–54.4%) (Table 1). The present results also show that both the GC content of plastomes and the GC content of first and second codon positions are highly correlated (P < 0.01) with the GC content of CDS regions (Table 3).

3.3. Phylogenetic results

Phylogenetic trees resulting from three data set pairs are shown in Fig. 1a–f. Numbers of RNA editing sites in Selaginella plastomes were mapped on the phylogenetic tree of codon1r-2r (Fig. 1f). The main clades of Selaginella are indicated as clades A–E in Fig. 1a–f. The tree lengths of the phylogenetic trees derived from corrected data sets were relatively shorter than those derived from plastome data sets (Table 2). For example, compared with the tree lengths of 3.552255 derived from data set AA, that from data set AAr decreased to 2.872170 (Table 2, Fig. 1b, e). The phylogenetic results based on all three plastome data sets support the sister relationship between clade B and C (MLBS = 55/100/100; Fig. 1a–c), while the results based on all three corrected data sets support the sister relationship between clade B and A with high bootstrap values (all MLBS = 100; Fig. 1d–f). Five of six analyses strongly support the sister relationship of clade D and clade A + B + C (all MLBS = 100; Fig. 1a–c, e, f), but the analysis based on CDSr matrix supports the sister relationship of clade D and all other Selaginella species (MLBS = 100; Fig. 1d). The results of three analyses based on plastome data sets (CDS, codon 1-2, and AA) in the present analyses (Fig. 1a–c) were mostly consistent with the results of Zhang et al. (2020).

4. Discussion 4.1. Variation of predicted RNA editing in Selaginella plastomes

The predicted number of C-to-U RNA editing sites in plastomes of three Selaginella species (S. kraussiana, S. lepidophylla and S. uncinata) were basically equal to their actual numbers, as verified by previous studies (Oldenkott et al., 2014; Smith, 2020). Therefore, the predictions were largely accurate and reliable. Generally, the number of RNA editing sites were relatively uniform among species within each of the Selaginella clades A–E (e.g., 2876–3100 in the four species of clade B; Fig. 1f, Table 1), but varied greatly between different clades, e.g., the numbers of editing sites in the two species of clade A are 1776 and 1858, which are much lower than those of its sister group, clade B (2876–3100) (Fig. 1f, Table 1). The studies of Zhang et al., 2019, Zhang et al., 2020 suggested that species in clade E have lower substitution rates and simpler plastome structures, whereas species in clades A–D have higher substitution rates and more complex plastome structures. Our analyses indicate that the frequency of RNA editing in Selaginella plastomes may be positively correlated with the variation in plastome structure and the substitution rate. However, there was an exception in the two species of clades A. Compared with species in clade D, species in clades A had a relatively lower number of RNA editing sites (Fig. 1f, Table 1), but more complex plastome structures and an accelerated substitution rate (Zhang et al., 2019, Zhang et al., 2020). More comprehensive taxon sampling and transcriptomic evidence are needed to illuminate the evolutionary patterns of RNA editing in Selaginella.

4.2. Correlations between RNA editing and GC content

Correlation analysis indicated that the number of RNA editing sites in Selaginella plastomes was highly correlated with the GC content of the first and second codon positions, but not correlated with the GC content of plastomes (Table 3). The GC content of plastomes was only moderately correlated with the GC content of first and second codon positions (Table 3). A similar pattern has been observed in the mitochondrial genomes of ferns and gymnosperms (Guo et al., 2016, 2017). The GC content of the plastid genomes of lycophytes vary extensively, with ca. 35–36% in Lycopodiales, ca. 38% in Isoëtales, and ca. 51–57% in Selaginellales (Mower et al., 2019; Zhang et al., 2019, Zhang et al., 2020). The present study showed that even for the corrected data set consisting of the most frequently edited first and second codon positions (codon1r-2r), the average GC content of the sampled Selaginella species (49.1%) was still relatively higher than those of club mosses (41.6–42.7%) and quillworts (44.5%) (Table 1, Table 2). The GC content of non-coding sequences, which are usually about half the size of the plastomes, are generally higher than 50%; however, RNA editing is rarely found in non-coding sequences (Oldenkott et al., 2014). Our study implies that there may be other mechanism(s) associated with the GC content of both coding and non-coding regions of plastome sequences.

4.3. Impact of RNA editing on phylogeny

The prediction of RNA editing sites in plastomes of Selaginella showed that an average of about 10% of the translated amino acids had changed due to an extreme editing event (Table 1). Previous studies have indicated that genomic DNA sequences may carry considerable homoplasy but have a greater number of variable sites; in contrast, RNA-edited DNA sequences may be more phylogenetically informative, although they have fewer variable sites because C-to-U editing tends to restore conserved codons (Bowe and dePamphilis, 1996; Varigerow et al., 1999; Peterson et al., 2006). Our results show that all the corrected data sets had relatively higher values of pairwise identity and thus contained less phylogenetic information (Table 2); however, the bootstrap values did not decrease significantly (Fig. 1), although the tree length of the phylogenetic trees derived from corrected data sets were relatively shorter (Table 2).

Contrast analyses showed that the extreme C-to-U RNA editing in coding regions of plastomes definitely altered the effect of the sequences in phylogenetic reconstructions. There were some substantial differences in the relationships among clades A–E in Selaginella between the results of the plastome and corrected data sets (Fig. 1). The three analyses based on the plastome data sets (CDS, codon 1-2, and AA) (Fig. 1a–c) supported the sister relationship between clade B and C, which is consistent with the results of Zhang et al. (2020), which were obtained using similar data sets and tree inference strategies. However, all three analyses based on corrected data sets (CDSr, codon1r-2r, and AAr) highly supported the sister relationship between clade A and B (Fig. 1d–f). Moreover, the results based on the corrected data set CDSr supported the sister relationship of clade D and all other Selaginella species (Fig. 1d), which is similar to the results of Zhou et al. (2016); however, the five other analyses supported the sister relationship between clade D and the aggregate of clade A + B + C. Overall, the infra-generic relationships within this big genus Selaginella still need further research with both adequate taxon sampling and molecular characters, as well as the correction of RNA editing sites.

5. Conclusions

In this study, we found that the extreme RNA editing in plastomes can substantially affect or even alter the result of phylogenetic reconstruction. Specifically, when we compared phylogenetic reconstructions of Selaginella based on plastome and RNA-edited data sets, we found that phylogenetic relationships among Selaginella differed. Therefore, we recommend that researchers correct RNA editing sites when using plastid or mitochondrial genes in phylogenetic analyses, especially in those lineages with abundant organellar RNA editing sites, such as hornworts, quillworts, spike mosses, and some seed plants.

Author contributions

J.-M.L., and D.-Z.L. designed research; X.-Y.D. performed research and analyzed data; and X.-Y.D., J.-M.L., and D.-Z.L. wrote the paper.

Declaration of Competing Interest

The authors declare no conflicts of interest.

Acknowledgements

We thank Dr. Peng-Fei Ma for improving the manuscript. We also thank the two anonymous reviewers for their constructive comments and suggestions. The study was supported by the Strategic Priority Research Program, Chinese Academy of Sciences, China (XDB 31000000); the National Natural Science Foundation of China (31970232); the Large-scale Scientific Facilities of the Chinese Academy of Sciences, China (2017-LSF-GBOWS-02); and the technological leading talent project of Yunnan, China (2017HA014).

Appendix A. Supplementary data

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