Dynamic landscape of mitochondrial Cytidine-to-Uridine RNA editing in tobacco (Nicotiana tabacum) shows its tissue specificity

RNA editing is a prevalent nucleotide modification at the RNA level in higher plants. However, little is known about the dynamic distribution of RNA editing among tissues. In this study, we explored the tissue specificity of mitochondrial RNA editing in tobacco (Nicotiana tabacum) based on publicly available RNA-seq data from four tobacco tissues: root, stem, leaf, and flower. As a result, 473 RNA editing sites involved in 60 mitochondrial genes were identified. The results showed an uneven distribution of editing sites among tobacco tissues, a total of 106 sites and 11 genes were identified as tissue-specific editing in the four tissues, and a total of 11 sites located in six genes were detected differentially edited statistically (p-value < 0.01). The expression of RNA edited genes and RNA editing factors was analyzed, and we observed that most tissue-specific edited genes were expressed at a low level. There were about ~  20 RNA editing factors that were differentially expressed between different tissues, indicating that the heterogeneity of RNA editing in different tissues might result from the expression regulation of RNA editing factors. Our analyses provide insights into the understanding of landscape, regulation, and function of RNA editing events in higher plants. Dynamic landscape of conserved editing sites revealed the tissue specificity of mitochondrial RNA editing in tobacco.


Introduction
RNA editing is defined as the insertion, deletion, and replacement of nucleotide bases that occurs after transcription, which usually results in differences between the RNA genetic information and the genome template (Liscovitch-Brauer et al. 2017;Walkley and Li 2017;Zahn 2017;Peng et al. 2018). The specific RNA positions affected by RNA editing and the corresponding DNA positions are called as editing sites (Edera et al. 2018a). For instance, the transcript of maize chloroplast gene 50S ribosomal protein L2 (rpl2) can only create the initiation codons after uridine replacement of cytidine (Hoch et al. 1991). In 1986, Benne et al. Communicated by Joyce Van Eck. biologists discovered an RNA editing event in the mitochondrial gene cytochrome c oxidase subunit II (cox2) of trypanosomes firstly (Benne et al. 1986), they found that the insertion of four uridines caused the cox2 gene to form a continuous open reading frame, resulting in changes in genetic information. In 1989, Hiesel et al. revealed that RNA editing also exists in primrose mitochondria (Hiesel et al. 1989). Subsequently, increasing RNA editing events were confirmed, to date, RNA editing has been found in primitive eukaryotes, vertebrates, plants, fungi and viruses (Palladino et al. 2000;Bahn et al. 2012;Alon et al. 2015;Guo et al. 2015;Riemondy et al. 2018;Zaidan et al. 2018).
Numerous studies have shown that RNA editing alters the genetic information of the genome and enriches the expression products of genes; on the other hand, it provides new genetic structures and functions for the evolution of organisms (Small et al. 2020;Lukes et al. 2021). For mammals, the regulation of RNA editing events has been widely studied. For example, RNA editing in mammals has been proved to be a dynamic landscape across tissues and additional potential editing factors have been discovered to be involved in the editing events (Tan et al. 2017;Blanc et al. 2019). In the case of higher plants, RNA editing occurs mainly in the coding regions of mitochondrial and chloroplast genes, where the most common types of editing is cytidines substituting uridines (C-to-U) (Covello and Gray 1989;Tillich et al. 2006;Takenaka et al. 2013). So far, in higher plants, the total number of RNA editing sites in the chloroplast genome (20-100) is much less than that in the mitochondrial genome (300-600) (Chen et al. 2017;Edera et al. 2018a;Brenner et al. 2019;Nawae et al. 2020).
RNA editing in the plant is mainly mediated by editing complex involving editing factors such as Pentapeptide repeat (PPR) proteins and multiple organellar RNA editing factors (MORF) (Haag et al. 2017;Tian et al. 2019). Editing factors recognize a 20-25 nucleotide sequence upstream of the editing site, and the editing factors required for different RNA editing sites may differ (Yagi et al. 2014;Yan et al. 2018). The PPR protein is the largest class of RNA editing factors in plants, which is encoded by nuclear genes and located in mitochondrion and/or chloroplast, playing a role of site-specific recognition in editing (Yagi et al. 2013(Yagi et al. , 2014Ichinose and Sugita 2017). PPR protein has two subfamilies, P-type and PLS-type, the P-type consists of the 35-amino acid classic PPR (P) motif, while the PLS-type consists of the classic P-motif and its variants L (35 or 36 amino acids) and S (31 amino acids) (Manna 2015). A plantspecific conserved E domain often exists at the C-terminus of the PLS-type PPR protein. Generally, in plant organelles, it is the PLS-type PPR protein that recognizes the specific editing sites (Shikanai 2015;Yan et al. 2018). MORF is a small protein family, with 10 members in Arabidopsis, and MORF mutants exhibit reduced efficiency at multiple sites (Takenaka et al. 2012). However, the target relationship of editing factors catalyzing RNA editing sites is still unclear.
Previous studies have revealed that RNA editing is specific in terms of tissues and developmental stages, suggesting that RNA editing may function as a regulatory device in plastid gene expression (Bock et al. 1993;Zeltz et al. 1993;Miyata and Sugita 2004). Tissue-specific and developmentspecific RNA editing were detected in the bryophyte ribosomal protein S14 (rps14) gene in previous studies (Miyata and Sugita 2004). Tseng also observed that the RNA editing efficiency of Arabidopsis plastid mRNA was variable among tissues (Tseng et al. 2013). Actually, earlier researches on RNA editing mainly focused on selected genes based on experimental methods with a disadvantage of low throughput. With the rapid improvement of genomic and transcriptome sequencing technology, a large number of plant genomes have been sequenced and numerous RNA sequencing (RNA-seq) data have also been generated, offering an opportunity to test the function and regulatory mechanism of RNA editing in plant growth and development. In the future, the regulation and function of RNA editing in plants and their effects on traits, especially some essential agronomic traits, will attract more people's attention (Small et al. 2020).
As a model plant, tobacco (Nicotiana tabacum) plays a key role in plant molecular research and is also an economically important plant worldwide. In this study, to illustrate the tissue distribution characteristics of RNA editing in the higher plant, we identified and analyzed RNA editing sites in four tobacco tissues (root, stem, leaf and flower) based on transcriptome data. And then, the dynamic landscape of RNA editing sites was analyzed, we discovered that editing sites were distributed differently in the four tissues, among which the root had the least editing sites, which might be related to the physiological function requirements for the root. Simultaneously, we also analyzed the expression levels of RNA editing factors in different tissues, its heterogeneity of expression levels might partly explain the varied RNA editing events happening in different tissues. Our findings showed that RNA editing was differentially regulated in various types of tissues, and may contribute to the functional differentiation of tobacco tissues.

Data collection and pre-treatment
The transcriptome data of four tissues (root, stem, leaf, and flower) and whole genome re-sequencing (WGS) data from Nicotiana tabacum cultivar TN90 (common tobacco) were downloaded from NCBI (https:// www. ncbi. nlm. nih. gov/ biopr oject/ PRJNA 208209), each sample consists of three replicates (Sierro et al. 2014), detailed information was listed in Table S1. In addition, the genome (assembly: GCF_000715135.1) with mitochondrial sequences (accession number: NC_006581.1) of Nicotiana tabacum cultivar Ntab-TN90, as well as their annotation files were also downloaded from NCBI data repository (https:// www. ncbi. nlm. nih. gov) (Sugiyama et al. 2005). We utilized the FastQC v0.11.8 tool to check the quality of the transcriptome data firstly (Brown et al. 2017). To increase sequencing depth, we merged the three duplicates from each tissue into one sample. The transcriptome data were mapped to mitochondrial genome reference using HISAT2 v2.1.0 software with default parameters (Kim et al. 2015). Afterward, each SAM file was converted into a BAM file, sorted and duplicatesremoved with SAMtools v1.9 (Li et al. 2009;Danecek et al. 2021). The variant calling process was conducted by SAMtools 'mpileup' command, the single nucleotide polymorphisms (SNPs) were identified by BCFtools 'call' command (Danecek and McCarthy 2017).

Identification of RNA editing sites
Based on the SNP results and genome annotation files, RNA editing sites were identified by using the REDO tool (Wu et al. 2018). REDO identifies RNA editing events in plant organelles and adopts a series of stringent rule-dependent filters to evaluate the reliability of editing sites (Wu et al. 2018), such as (1) quality control filter: the low-quality sites are filtered based on reads quality; (2) mapping depth filter: sites with total reads depth more than 4 are kept; (3) alt proportion filter: alt_proportion-error_proportion] <0.1; (4) multiple alt filter: only the variant with one alt allele is retained for RNA editing detection; (5) distance filter: given the possible positional interference for RNA editing, the variant sites in short distance (< 3 bp) are filtered out; (6) spliced junction filter: variants within short spliced anchor (< 2) are removed; (7) indel filter: the indel variants are excluded by default; (8) likelihood ratio (LLR) test filter (LLR < 10): LLR test is a probabilistic test incorporating error probability of bases (error probability is obtained using adjacent nonvariant sites in specific window) for detecting RNA editing sites (Chepelev 2012;Sun et al. 2016); (9) Fisher's exact test filter (p-value < 0.01): the significance of the RNA editing site is evaluated by using Fisher exact test to compare the expected level of itself, and the p value is used for filter; (10) complicated filter model: according to the statistical results of the RNA editing sites that have been experimentally verified and the attributes of the codon table, the complicated filter model consists of five characteristics of RNA editing sites, which are RNA editing types (C→T, A→G, T→C, etc), alt proportion, amino acids change, codon phase, and hydrophobic/hydrophilic change (Wu et al. 2018). Finally, the corresponding annotation information files for all identified RNA editing sites were generated. To improve the accuracy of RNA editing sites filtering, genome SNP-calling result was further used to eliminate genomic variation. Prediction of RNA editing events using by PREPACT (Lenz et al. 2018) was further compared with the above results to test the reliability of the final editing sites.

Characteristic statistics and comparison of RNA editing sites
For each tissue, the resulted editing sites were used for further characterization, including statistics on the number of editing, editing types, editing efficiency, codon positions, amino acid changes, and involved genes. The nucleotide substitution types of total editing sites were counted. For the triplet codon, we calculated the substitutions of base at each position and counted the editing events separately. In addition, according to codon functional traits, we counted the number of start codons and termination codons and analyzed their changes at the amino acid level. For one site, the value of RNA editing efficiency was expressed as the proportion between edited transcripts and total transcripts, if the site was edited, the C/G base (wild type) should be altered to the T/A base (edited type), since one editing site could be detected hundreds of times via RNA sequencing, the number of wild type (C/G) or edited type (T/A) of bases could then be counted at this particular site, hence, the editing efficiency at one site could then be calculated by the formula: depth of edited bases (T and A)/total read depth of bases. Furthermore, we compared the RNA editing efficiency between every two tissues, RNA editing sites with statistical significance (p-value < 0.01) were identified. Aiming to decipher the tendency of RNA editing efficiency for each tissue, cluster analysis was performed based on the matrix of RNA editing efficiency, which was normalized by subtracting the row-wise mean from the values in each row of data and multiplying all values in each row of data by standard deviation value documented in our previous study (Zhang et al. 2020). 'Pheatmap' function in R (http:// www.R-proje ct. org/) was used to plot the heatmap of editing efficiency, 'dist' function was used to calculate the distance matrix of different samples with the default 'euclidean' method, 'hclust' function was used to compute the hierarchical clustering.

Expression analysis
Furthermore, to measure and compare the gene expression level of RNA edited genes and RNA editing factors, we also performed transcriptome analysis as described in the previous study for the RNA-seq data (Pertea et al. 2016). Two main types of RNA editing factors, including PPR proteins and MORF proteins, were chosen to explore their expression 1 3 in RNA editing activity. To obtain all tobacco PPR and MORF genes, two representative protein sequences (Uni-ProtKB ID: Q9SAD9, O49429) from Arabidopsis thaliana were used as queries to search against the tobacco protein databases using BLASTP v2.10.1 (http:// blast. ncbi. nlm. nih. gov). All positive hits were retrieved for gene function annotation to blast against the Swiss-Prot protein database. The clean reads of RNA-seq data from each sample were mapped against the tobacco genome reference with HISAT2 v2.1.0 (Kim et al. 2015), each SAM file was converted into a BAM file, and sorted, duplicates-removed with SAMtools v1.9 (Li et al. 2009;Danecek et al. 2021). Further transcript assembly and quantification of the read alignments were performed using Stringtie v1.3.3b (Pertea et al. 2015). Gene expression levels were measured using FPKM (fragments per kilobase of transcript per million mapped reads), expression values were also normalized by the method mentioned above. EdgeR package (Robinson et al. 2010) was used to determine the differentially expressed genes between every two tissues. A heatmap with all samples was plotted using 'pheatmap' function in R (v.3.3.1) environment.

Statistical analysis
A two-tailed Wilcoxon rank-sum test was used to perform the pairwise comparison of RNA editing efficiency between every two tissues. As for the pairwise comparison for each editing site, Fisher's exact test was used.

Alignment of transcriptome data
A total of 12 RNA-seq extracted from four tobacco tissues (root, stem, leaf, and flower) and one WGS samples from PRJNA208209 were used in our study. There are three replicates for each tissue. We merged the replicates of each tissue and aligned the transcriptome data against the tobacco mitochondrial genome as reference. The tobacco mitochondrial genome used in our study is 430,597 bp in size and has 116 Open Reading Frames (ORFs). To avoid data bias in sequencing depth that might affect the comparability among four tissues, we unified the number of reads in four tissues RNA-seq data, then a total of 130,415,792 reads were selected randomly for transcriptome alignment. At last, there were 56,774, 81,378, 125,616 and 92,165 reads in the root, stem, leaf, and flower tissues mapped onto the tobacco mitochondrial genome, with corresponding overall alignment rate was about 0.04%, 0.06%, 0.09% and 0.07% respectively ( Table 1). The alignment results of transcriptome data revealed that roots have a lower mapping rate, one reason is that genes were expressed at a lower level in root, the other reason is that variation in tissue composition, like the mitochondrial density, may affect the RNA source of different tissue sample. We observed that the average mapping depths were around 20× that meets the minimum requirement of identification of editing sites.

Characteristic statistics of tobacco mitochondrial RNA editing
We used REDO tool that is an automated approach based on multiple layers of filtering to detect RNA editing sites, therefore, for merged RNA-seq data from four tobacco tissues, a total of 487 raw editing sites were detected. To improve the accuracy of RNA editing sites filtering, genome SNP-calling result was further used to eliminate genomic variation, 14 sites were both detected in RNA and genome level with the same substitution type, we excluded them as genomic variations. After that, we also inspected the remaining editing sites manually to ensure the reliability of results. Finally, a total 473 editing mitochondrial RNA editing sites were predicted that located in 60 mitochondrial genes in four tobacco tissues, However, we only detected RNA editing sites that occurring within mitochondria genes, and excluded the sites occurring within non-coding regions, hence, the number of RNA overall editing sites is lower than that of a previous work (Grimes et al. 2014). We also compared our results with that predicted by PREPACT method (http:// www. prepa ct. de/ prepa ct-main. php) (Lenz and Knoop 2013), the resulting number and distribution agrees with that of our results, see Table S2, hence, we deemed that the identification results from the REDO tool in our study were reliable.
The statistics of RNA editing types revealed that C-to-U is the most common editing type, accounting for more than  Figure 1A). There were also other substitution types, such as U-to-C, which was recently assumed to be RNA editing type in Arabidopsis thaliana (Ruchika et al. 2021), other mismatches might be the result of sequencing errors and double transcription of nuclear genes and mitochondrial genes (Edera et al. 2018a). Detailed information on editing sites is listed in Table S3. The average numbers of editing sites varied widely among genes. On average, NADH dehydrogenase subunit 7 (nad7) gene had the largest number of editing sites, approximately 17, followed by 16 editing sites found in cytochrome oxidase subunit 1 (cox1) gene. Out of the 473 identified editing sites, 143 (30.23%), 281 (59.41%), and 48 (10.15%) were located in first, second, and third codon positions respectively. It means that RNA editing events occurred preferentially on the second codon and followed by the first codon ( Figure 1B). The statistics of editing type showed that there were far more editing events of non-synonymous mutations than synonymous mutations, with non-synonymous editing events up to 88%, as shown in Table S3. In addition, the editing level of non-synonymous mutations was significantly higher than that of synonymous mutations ( Figure 1C), indicating that non-synonymous C-to-U editing was favored by natural selection. In addition, two sites located in cox1 and ribosomal protein S10 gene (rps10) were detected to produce functional start codon, another five sites were found to produce premature termination codon, see Table 2. Edits at cox1_2 and rps10_2 were found in all tissues, and the average editing efficiency was 87.2% and 94.2%, respectively. Sugiyama et al. similarly reported that the genomic ACG codon in the rps10 gene was changed to the start ATG codon by RNA editing, but the same change was not reported in cox1 gene (Sugiyama et al. 2005). The edited positions of five premature termination codons were all located at the first position of codon, the editing efficiency ranges from 17.3% to 97%, four out of five codons were translated into glutamine before editing.
Most of the amino acid changes tended to be hydrophobic, the most amino acid changes were Ser-to-Leu (21.35%) and Pro-to-Leu (18%), see Figure 1D, these changes of hydrophobicity could greatly restore the conservation of proteins and were essential for maintaining normal protein function. The proportion of hydrophobic amino acids in the root, stem, leaf, and flower issues were 40.98%, 37.38%, 41.67%, and 42.24%, after editing, the proportions were increased to 81.97 %, 80.84%, 81.67%, and 81.9% respectively (Supplementary Figure S1). For RNA editing efficiency, the average is 80%, ranging from 10 to 100% (see Figure S2).

Tissue-specific analysis of tobacco mitochondrial RNA editing
For the four tissues (root, stem, leaf, and flower), the number of RNA editing sites identified in mitochondrial corresponded to 204, 294, 234, and 343, involving 43, 45, 45 and 54 genes respectively, as shown in Table S3, more editing sites were detected in flower, whereas root had the least number. As for RNA editing efficiency, the RNA editing efficiency of root, stem, leaf, and flower corresponds to 85.6%, 88.5%, 89% and 86.6% respectively (Table S1), leaf had the highest editing efficiency, followed by stem, flower, and root. Clustering and heatmap plotting was also performed based on the normalized matrix of RNA editing efficiency, as is shown in Figure 2, root clustered far away from other tissues. These results demonstrated an uneven distribution and dynamic landscape among different tissues, for example, compared with other tissues, root had not only the least number of editing sites but also the lowest editing efficiency.
The venn diagram analysis (Figure 3) showed that there were 131 common editing sites and 34 edited genes shared in all four tissues, considerable editing sites were only detected in certain tissues, the tissue-specific editing sites in root, stem, leaf, and flower correspond to 16, 18, 13, 59 respectively, more sites were specifically edited in flower. Detailed information on tissue-specific editing sites is listed in Table S5. For example (Figure 4), RNA editing in cob_114 (genome position: 41934) occurred in all the tissues except root; matR_326 (genome position: 5687) was only edited in leaf tissue. As for tissue-specific edited genes, a total of 11 genes, containing 14 editing sites, were specifically edited in certain one tissue, including one gene (orf105) in root, two genes (rrn26 (rRNA), orf306) in stem, one gene (succinate dehydrogenase subunit 3 (sdh3)) in leaf, and seven genes (30S ribosomal gene (rps14), orf121a, orf158, orf111, orf101a, orf274, cytochrome c maturation genes (ccmB)) in flower, which owned more tissue-specific edited genes (Table 3). We found that all these genes contained only one or two editing sites, nearly all the sites were located in the second position of codons except two sites (orf306_242, orf121a_327).
In the root, both editing sites of orf105 genes did not change the amino acid type, however, there were still more sites that altered the amino acid type through RNA editing, such as orf306_242, rps14_305, ccmB_204, and so on. In addition, these editing sites demonstrated a high editing efficiency except for two sites (orf306_242, rrn26_295). The above results revealed that the editing events are tissue-specific, it seems editing is more active in flower and weaker in the root, however, the vast majority of genes were conservative and edited in at least two tissues, tissue-specific edited genes were mostly orf genes whose functions remain unclear.

Changes in RNA editing efficiency among different tobacco tissues
Despite tissue-specific editing, we also observed that a batch of sites demonstrated varied editing efficiency among different tissues, some of them were remarkably reduced or increased in certain tissues. Pairwise comparisons of editing efficiency between any two tissues were conducted, a total of 11 sites located in six genes were detected differentially edited statistically (p-value < 0.01), see Table 4. The six genes consisted of three ribosomal genes (rps10, rpl16, rps4), one hypothetical gene (orf125e), and two other functional genes (atp9, nad5). One site of ribosomal protein S10 (rps10_2) was nearly completely edited in stem, leaf and flower, but reduced to 69.2% in root, in addition, another site in ribosomal protein L16 (rpl16_488) was completely edited in the stem, but reduced to 54.5%, 77.3%, and 82% in leaf, root, flower respectively; another one site of ATP synthase F0 subunit 9 (atp9_223) in root and flower corresponds to 97.5% and 100%, but reduced to 78.9% and 69.2% in leaf and stem respectively. Remarkably, two sites in orf125e (orf125e_153, orf125e_192) demonstrated an opposite trend of change, compared with the other three tissues, orf125e_153 was up-edited in leaf, whereas orf125e_192 was down-edited, indicating that the heterogeneity and complexity of editing not only occurred in different tissues, but also in different sites of the same gene. Detailed information on these editing sites can be found in Table S6.

Expression analysis of RNA edited genes in different tissues
To investigate the reason for RNA editing variation in different tissues, we examined the expression of RNA edited genes and RNA editing factors. For the 60 RNA edited genes, we observed that most of the RNA edited genes (54 out of 60) were expressed in all tissues with no bias, however, the overall expression level was highest in stem, but lowest in root ( Figure 5A). There were only six genes that were not expressed in certain tissues, including four orf genes (orf132, orf158, orf190, orf111), one tRNA gene (trnD(guc)), and cytochrome c maturation gene (ccmFc). In terms of the 11 tissue-specific edited genes, we observed that they were nearly expressed in all tissues at a low level except rrn26 and rps14 ( Figure 5B), although some genes were expressed at a low level compared with that of other tissues, they were still only edited in this tissue, such as  orf105 gene, which was highly expressed in stem, but only edited in root. Pairwise differential expression analysis was conducted for these RNA edited genes, a total of 16 genes were differentially expressed (Table S7), several of them were tissue-specific edited genes or contained differentially editing sites, such as orf306, rps14, rps4, nad5, orf125e, however, no associations between them was detected. One striking example is rps14 gene, which showed vastly different expression levels between tissues, highest in leaf, lowest in root, but only edited in flower. The above observation suggested that although both expression and editing of transcripts are tissue specifically regulated in different tissues, however, there was no causality between these two events, different types of tissues may generate different kinds of signals to affect RNA editing and expression of organelles genes separately. Our observations also agree with our previous study about the dynamic response of RNA editing to heat stress in grape, no expression difference of RNA edited genes were detected under different temperatures.

Expression analysis of RNA editing factors in different tissues
To investigate the correlation between RNA editing pattern and editing factors or their interacting proteins, we also evaluated the expression of RNA editing factors, including PPR and MORF proteins, which function as catalytic machines and play key roles in RNA editing process (Yan et al. 2018). After blast searching, a total of 17 MORF and 170 PPR proteins were identified in tobacco. We found that approximately 20 RNA editing factors were differentially expressed between leaf and other three tissues, whereas only 10 RNA editing factors were differentially expressed between tissues of flower, stem and root, see Table S8. Based on the normalized expression matrix of all RNA editing factors and differentially expressed factors, heatmap plottings were performed, as is shown in Figure 6, Figure S3, which demonstrated an uneven distribution among different tissues. A batch of RNA editing factors was highly expressed in leaf significantly, we found a large proportion (8 out of 17) of MORF proteins were differentially expressed, five of them were highly expressed in leaf, the other two morf genes were highly expressed in leaf and flower, including multiple organellar RNA editing factor 2 (morf2, gene symbol: LOC107777961) and multiple organellar RNA editing factor 3 (morf3, gene symbol: LOC107760068). However, there were still five genes that were expressed in leaf at a low level, but highly expressed in root, like gene LOC107787019, LOC107782030. MORF protein family were once considered to be essential components of plant editosomes, a recent study has proved that RNA editing events in rosette leaves and flowers were reduced by morf9 mutation (Tian et al. 2019). Hence, the above results suggested that these RNA editing factors also exhibited tissue-specific expression, especially morf genes, which may partly explain the reason for the varied RNA editing events among different tissues.

Discussion
As an important epigenetic mechanism that modified genome-encoded transcripts, RNA editing diversifies genomically encoded information to expand the complexity of the transcriptome (Takenaka et al. 2013). Previous studies revealed that RNA editing has various biological functions, including promoting RNA splicing by affecting the intron structures, and playing a central role in plant development and evolutionary adaptation (Ichinose and Sugita 2017;Tang and Luo 2018). If the proteins produced by edited and unedited transcripts had different functions, it would be beneficial for plants to have regulated RNA editing. Hence, it was conceivable that dynamic RNA editing, not changes in the genomic DNA sequence, might provide a flexible way to increase plant fitness and adaptation to the constantly changing environment. Occasionally, transcripts that contained a partially edited site might generate different forms of a protein with distinct or overlapping functions and further increase the diversity of the encoded proteins, but the underlying mechanism of functional differentiation remains unclear. Until now, the discrepancy in RNA editing sites among tissues has only been detected in certain genes, such as photosystem II protein VI (psbF) and photosystem II protein L (psbL) genes (Bock et al. 1993). It is still difficult to make general conclusions to what extent that discrepancy exists between them. Hence, to obtain an overview of the distribution of RNA editing in the model plant tobacco, we used the transcriptome data and performed bioinformatics analysis to examine the editing profiles of RNA editing and gene expression in root, stem, leaf, and flower tissue. A total of 473 RNA editing sites in 60 mitochondrial genes were identified. The tissue specificity of RNA editing was examined, we found 106 tissue-specific RNA editing sites and 11 tissue-specific edited genes in the four tissues. Despite tissue-specific editing, we also observed that a batch of sites demonstrated varied editing efficiency among different tissues. Our results demonstrated the fact of heterogeneity of RNA editing patterns among different tissues. However, the tissue-specific edited genes were mostly expressed at a low level, no associations were detected between their expression and editing. Oppositely, considering that the differentially expressed RNA editing factors between tissues, we proposed that the tissue-specific RNA editing may be mediated by them via expression regulation.
Compared to the other three tissues, we also found that root had not only the least number of editing sites (204), but also the lowest average editing efficiency, although we observed that nearly all the RNA edited genes were expressed in four tissues and the overall expression level was the lowest in root ( Figure 5A). Hence, transcript abundance might be an important limiting factor in identifying editing sites and the increase of sequencing depth of RNAseq will effectively avoid loss of editing sites. However, given the observation that a few sites occurred in low-level expressed genes, we still consider that they were tissuespecific edited and irrelevant to expression. Although the overall average editing efficiency was presented at a high level (~ 80%), certain sites still had varied editing efficiency among different tissues. For example, the editing efficiency of rps10_2 in root was 0.69, whereas it was nearly fully edited in the other tissues (≥0.95), as shown in Table S4. Compared with root, the other three tissues had more editing events and higher editing efficiency on the main respiratory chain (consisting of complexes I, III and IV), indicating that root tissues may differ from other tissues in terms of energy requirements. Root is an underground organ, functions in fixation of plant to the soil and absorption of water, it does not perform photosynthesis, and simultaneously lose the necessity of correction of normal functions of certain genes. Moreover, compared with root, transcripts of ribosomal genes and maturase-related genes owned more editing sites in stem, leaf, and flower, indicating that RNA editing plays a key role in rebuilding functions of protein biosynthesis and RNA splicing for these tissues. Our results were consistent with the results of previous studies (Tseng et al. 2013), the editing of respiratory complex I transcripts in root was significantly downregulated compared with other tissues. It is worth noting that more editing sites and genes were detected in flower, as a propagative organ, more energy-relevant proteins may be required in mature flowers to ensure normal physiological functions. The functional discrepancy between tissues makes different selection pressure on RNA editing in respiratory genes, which may partly explain the tissue specificity of RNA editing. Regardless of the dramatic discrepancy between tissues, there are still many conservative editing sites that are highly edited, up to 0.9 in all tissues, such as cytochrome oxidase subunit genes (cox1, cox2, cox3), they play key roles in the respiratory electron transport chain, their stable and high editing ensures their edited proteins are homogeneous in assembling into functional complexes that are equally important for any tissues.
The current hypothesis is that several types of RNA editing factors constitute an editosome to catalyze the C-to-U conversion (Schmitz-Linneweber and Small 2008). MORF proteins were demonstrated to affect RNA editing event, and they were necessary for seedling survival in rice and seed development in maize (Liu et al. 2020;Zhang et al. 2021). Therefore, it is probable that RNA editing is regulated by the expression of editing factors and their interacting proteins, such as PPR, MORF proteins. In our study, we also observed that RNA editing factors demonstrated tissuespecific expression pattern, especially morf genes, which may partly explain the reason for tissue specificity of RNA editing. However, more ppr genes were highly expressed in leaf, but not flower, indicating that the regulation relationships between RNA editing factors and editing sites were not one-to-one, a striking feature of the plant RNA editosome is its diversity and complexity in the composition without consistency between individual RNA targets, hence, further studies are still needed to support this points.