Identification of RNA editing sites
A series of species that represent distant evolutionary clades were selected judging by two criteria, one is enough transcriptomic data of leaf tissue in the SRA database, and another is the availability of sequenced chloroplast genome. Hence, 21 species, consisting of 6 ferns, 4 gymnosperms, and 11 angiosperms (see Table 1), and corresponding 317 SRA accessions were chosen finally. Detailed information of SRA data and chloroplast genome accessions were listed in Table S1. Mapping results (Figure S1) showed that RNA-seq data volume varied among different species, compared with angiosperms, ferns and gymnosperms have lower mapping depths, especially for Adiantum aleuticum and Histiopteris incisa, indicating that the actual number of RNA editing in ferns and gymnosperms might be underestimated. We also observed that the read density also varied widely among different genomic regions, ranging from less than 30 to more than 800 in a few species, demonstrating varied expression levels of genes or sequence biases. Based on results of RNA-seq data mapping and SNP calling, an automated bioinformatics pipeline implemented in REDO tool (Wu, Liu et al. 2018) was conducted under default thresholds. Consequently, there were a total of 6,011 raw editing sites located chloroplast genome was detected in leaf. However, despite multiple stringent filters, sequence mismatches that accord with RNA editing occasionally appeared, hence, we manually examined all mismatches to eliminate false positive, only kept C-to-U and U-to-C editing types, editing sites with editing efficiency close to 100% that may be derived from genomic variation were also eliminated. Finally, a total of 5,205 RNA editing sites, consisting of 237 U-to-C and 4,968 C-to-U edit sites, were screened out for further analysis. To evaluate the reliability of results, PREPACT webserver (Lenz, Hein et al. 2018) was also used to verify the editing sites. We found that the distribution of the predicted number of RNA editing based on PREPACT agreed with that of our prediction based on RNA-seq data basically (see Figure 1), reflecting our pipeline offered high performance with reliable results. To describe the attributes of RNA editing sites, we illustrated by one example of samples of Adiantum aleuticum (see Figure S2), which depicts the reliability of RNA editing sites by REDO tools statistically.
Table 1
The summary of RNA editing sites across three clades
|
Species
|
chloroplast genome
|
Number of raw editing sites
|
Number of filtered
editing sites
|
Percent of edited genes
|
MD1
|
Ferns
|
Pteris vittata
|
MH500228
|
349
|
319
|
74%
|
102
|
Adiantum aleuticum
|
MH173079
|
422
|
362
|
82%
|
23
|
Selaginella moellendorffii
|
HM173080
|
1,316
|
1,316
|
96%
|
73
|
Histiopteris incisa
|
NC_040220
|
588
|
414
|
85%
|
27
|
Cibotium barometz
|
NC_037893
|
562
|
540
|
85%
|
32
|
Cyrtomium fortunei
|
NC_037510
|
778
|
551
|
91%
|
40
|
Gymnosperms
|
Ginkgo biloba
|
NC_016986
|
301
|
291
|
81%
|
228
|
Picea abies
|
HF937082
|
123
|
121
|
61%
|
450
|
Cycas revoluta
|
NC_020319
|
173
|
132
|
44%
|
27
|
Pinus massoniana
|
MF564195
|
96
|
93
|
59%
|
182
|
Angiosperms
|
Liriodendron tulipifera
|
NC_008326
|
320
|
260
|
75%
|
113
|
Nelumbo nucifera
|
NC_025339
|
154
|
121
|
63%
|
295
|
Nicotiana tabacum
|
Z00044
|
107
|
74
|
23%
|
134
|
Glycine max
|
NC_007942
|
64
|
59
|
34%
|
195
|
Populus tremula
|
KP861984
|
93
|
89
|
43%
|
244
|
Arabidopsis thaliana
|
KX551970
|
50
|
45
|
24%
|
389
|
Gossypium hirsutum
|
NC_007944
|
127
|
97
|
37%
|
171
|
Helianthus annuus
|
NC_007977
|
59
|
57
|
29%
|
148
|
Phoenix dactylifera
|
NC_013991
|
121
|
118
|
42%
|
781
|
Zea mays
|
NC_001666
|
105
|
79
|
37%
|
814
|
Oryza sativa
|
NC_001320
|
103
|
65
|
26%
|
550
|
Note: 1 mapping depth.
Distribution of RNA editing across three clades
After manual inspection and elimination of mismatches, filtered RNA editing sites across 21 species were screened out, the summary was listed in Table 1, Table S2, detailed corresponding annotation information was produced simultaneously, as listed in Table S3. All the filtered editing sites were located in 1,038 genes across all species, the average percent of edited genes for the three clades (ferns, gymnosperms, and angiosperms) are 0.85, 0.61, and 0.39 respectively. Compared to the latter two clades, more sites and genes were edited in chloroplast transcripts from ferns remarkably (Figure 2a). In ferns, Selaginella moellendorffii has the largest number of editing sites, up to 1,316, which is nearly 100-fold more abundant than that of flowering plants, nearly all the chloroplast genes (77 genes, ~96%) suffered effective editing. Differ from Selaginella moellendorffii that exclusively belongs to Lycopsida, the other five fern plants are members of Leptosporangiopsida, have relatively smaller numbers of editing sites, represented by Cyrtomium fortune, which owned the second-largest number of editing sites, with 551 editing sites and 79 edited genes (~91%). Whereas for gymnosperms, the average number of editing sites and percent of edited genes were all less than that of ferns. Compared with the other three gymnosperms, Ginkgo biloba has the most editing sites, with 291 editing sites and 68 edited genes (~81%). On the opposite end, angiosperms have the lowest average numbers of editing sites, only a part of genes were effectively edited, lower than 50%. It was noticeable that Liriodendron tulipifera and Nelumbo nucifera distinguished them from other angiosperms with more editing events. A total of 260 editing sites were detected in Liriodendron tulipifera, which is nearly 3-fold more abundant than that of other angiosperms and gymnosperms except for Ginkgo biloba, the percent of edited genes was up to 75%, well above the average. The numbers of editing sites among other 9 angiosperms showed no significant differences. The above results illustrated the diversity of RNA editing distribution across the three clades, early-diverging clades/species showed higher numbers of editing sites compared with that of later-diverging. The statistics of RNA editing sites showed that 21.97%, 67.16%, and 10.58% of sites were edited at first, second, and third codon positions respectively (Figure 2b). In terms of editing types, C-to-U was the dominant editing type (nearly ~95.1%), the next is U-to-C type, and other mismatch types were rare (Figure 2c). The statistics of editing types showed that the majority (~95%) of the editing events resulted in non-synonymous codon changes, and the changed amino acids tend to be hydrophobic, the frequency of changes from hydrophilic to hydrophobic was the highest, followed by changes from hydrophobic to hydrophobic (Figure 2d). The most common amino acid change types were Ser-to-Leu and Pro-to-Leu, serine is hydrophilic, whereas leucine and proline are both hydrophobic. The above results demonstrated that the RNA editing exhibited a selective advantage in the overall increase in hydrophobicity of the resulting proteins, which was also in good agreement with our previous studies (Zhang, Jiang et al. 2020).
Variability of RNA editing among species and genes
To further explore the evolutionary trajectory of RNA editing for the three clades, for each gene, we summed the number of its editing sites across 21 species (see Table S4), and picked out the top 30 genes with most editing sites for cluster analysis. Base on the matrix of numbers of editing sites across 21 species, a hierarchically clustered heatmap was plotted in Figure 3. The clustering relations showed that the 21 species were divided into three subgroups. The remarkable abundant of editing events make two exceptional species (Selaginella moellendorffii and Liriodendron tulipifera) cluster far away from their own clades. We found that not every gene were edited in all the species, by grouping the genes based on their function, genes encoding membrane subunits of the chloroplast NDH complex (ndh) exhibited larger average numbers of editing sites, this is consistent with previous studies that RNA editing occurred preferentially in genes encoding membrane-bound proteins under strong selection (Mower and Palmer 2006). Whereas ribosomal subunit genes showed lower numbers of editing sites. Due to the well-studied background and abundant editing sites in the plant, ndhB gene is assumed to be a good case for the study of RNA editing evolution, in our study, ndhB was confirmed to possess the most editing sites, with 333 editing sites spread across 17 species. Despite this, there is a biased distribution of RNA editing sites in ndhB among three clades. In ferns, 50 sites were detected in Selaginella moellendorffii, only 4~17 sites in the other five ferns; in angiosperms, there were about 20 editing sites for each of the 11 angiosperm species; in gymnosperms, RNA editing in ndhB was only detected in Ginkgo biloba, for Picea abies and Pinus massoniana, no ndhB gene annotated in their chloroplast genome, whereas for Cycas revolute, no RNA editing events were detected in its ndhB gene, which may result from loss of editing or too low depth around genomic regions of its ndhB gene. Based on alignments of merged protein sequence for top 30 genes across 21 species, phylogenetic tree by Maximum Likelihood method was constructed, as shown in Figure 3b, which showed that the clustering relations agreed with that of the matrix of numbers of editing sites roughly. The above observation implied the conservation of RNA editing evolution across different clades. Hierarchically clustered heatmaps of numbers of RNA editing sites from all/each clade were shown in Figure S3-6 respectively.
Variability of RNA editing level
RNA editing level was used to measure to what extent the edited transcripts among all transcriptomes for one gene. In this study, we explored the distribution of editing levels among three aspects: codon positions, species, and edited genes. We observed that the distribution of RNA editing level for three codons did not comply with the normal distribution, featuring a peak around ~0.2 and fat tails, as shown in Figure 4. The comparison between codon positions showed that the average editing level in the second codon position (~0.72) is higher than that of the first (~0.63) and third codon positions (~0.61) significantly, suggesting non-synonymous substitution occurred in second codon position tend to be effectively edited, it was higher editing level that dominated the landscape of RNA editing. In addition, we also found that the average editing level also varied widely across 21 species, ranging from 0.42 to 0.85 (Figure 5a), the trend of plotting of filtered sites was consistent with that of raw sites roughly. Selaginella moellendorffii has the lowest editing level (~0.42), far below that of other species. In gymnosperms and angiosperms, Ginkgo biloba and Nelumbo nucifera had the lowest editing efficiencies (~0.6) respectively. It seemed that abundant editing sites detected in those early-branching plants might have a negative impact on their editing level. However, as an ancient plant in angiosperms, Liriodendron tulipifera was one exception, with an editing level up to 0.81. The editing level was also analyzed in each gene individually, we averaged the RNA editing level among the top 30 genes across 21 species, the result also demonstrated the diversity of distribution, as shown in Figure 5b, Maturase K gene (matK) has the lowest editing level (~0.5), oppositely, editing level of atpA gene was the highest, up to 0.88.
Reduced cytosines content with the evolution of plants
Considering the large differences in the scale of RNA editing events along with evolution, we analyzed the nucleic acid base composition of edited genes shared by the 21 species. There were a total of 51 genes shared in all the chloroplast genomes of 21 species. The percent of cytosines for each gene was calculated by the formula: the number of cytosines (C)/total number of bases (A/T/C/G), as listed in Table S5. Afterward, the cytosines content of genes was averaged for each clade, remarkable significances (p-value < 0.05) were detected for pairwise comparisons between every two clades except for gymnosperms-angiosperms, as shown in Figure 6a. The percent of cytosines in ferns was far below that of angiosperms, followed by gymnosperms. The highest average of cytosines content was found in Selaginella moellendorffii (~0.26), the lowest average was found in Glycine max (~0.17). Ginkgo biloba and Liriodendron tulipifera have the highest averages (~0.19) in their own clade, showing a positive correlation with their high numbers of RNA editing sites. We further compared the cytosines content for each gene across the 21 species, as shown in Figure 6b, which demonstrated that the percent of cytosine dramatically declined roughly along with evolution with a few exceptions, such as 50S ribosomal protein L23 (rpl23) gene. One striking example was Photosystem II reaction center protein I (psbI) gene, which has the highest percent of cytosines in Selaginella moellendorffii (~0.37), and dropped to about 0.18 in other species. The above results indicated that the declination of cytosine content with evolution was the significant driver of loss of editing sites for later-branching plants.
Illustration of atpA RNA editing
We illustrated the RNA editing of atpA gene to help understand the evolution trajectory vividly. All editing sites identified in atpA gene were marked by a yellow color, as shown in Figure 7, it demonstrated that numbers of editing sites, as well as cytosine content, declined from ferns to angiosperms. We found that editing sites at the third codon position only occurred in certain species, were poorly conserved, for example, RNA editing in locus-a (Figure 6) only occurred in Adiantum aleuticum despite the existence of cytosine in other fern members, indicating that synonymous substitutions at third codon position were not conserved. In contrast, RNA editing at the second codon position was relatively conserved, RNA editing in locus-b and -c occurred in all three clades. Whereas in locus-d, RNA editing occurred in two members of gymnosperms, and the base types of other 19 species were all corrected to thymine in the genome, furthermore, we found many sites occurred in certain species for the three clades, implying the independent origins of RNA editing across different clades. In extreme cases, when all the cytosines with potential editing capacity of atpA gene were completely edited in fern plants, the overall uracil levels in edited transcripts will increase to the comparable levels of angiosperms and gymnosperms on the whole. However, in ferns’ transcripts, most of the sites were partly edited, hence, their uracil level is still lower than that of angiosperms and gymnosperms.