Genome-Wide Identification of U-To-C RNA Editing Events for Nuclear Genes in Arabidopsis thaliana

Cytosine-to-Uridine (C-to-U) RNA editing involves the deamination phenomenon, which is observed in animal nucleus and plant organelles; however, it has been considered the U-to-C is confined to the organelles of limited non-angiosperm plant species. Although previous RNA-seq-based analysis implied U-to-C RNA editing events in plant nuclear genes, it has not been broadly accepted due to inadequate confirmatory analyses. Here we examined the U-to-C RNA editing in Arabidopsis tissues at different developmental stages of growth. In this study, the high-throughput RNA sequencing (RNA-seq) of 12-day-old and 20-day-old Arabidopsis seedlings was performed, which enabled transcriptome-wide identification of RNA editing sites to analyze differentially expressed genes (DEGs) and nucleotide base conversions. The results showed that DEGs were expressed to higher levels in 12-day-old seedlings than in 20-day-old seedlings. Additionally, pentatricopeptide repeat (PPR) genes were also expressed at higher levels, as indicated by the log2FC values. RNA-seq analysis of 12-day- and 20-day-old Arabidopsis seedlings revealed candidates of U-to-C RNA editing events. Sanger sequencing of both DNA and cDNA for all candidate nucleotide conversions confirmed the seven U-to-C RNA editing sites. This work clearly demonstrated presence of U-to-C RNA editing for nuclear genes in Arabidopsis, which provides the basis to study the mechanism as well as the functions of the unique post-transcriptional modification.


Introduction
RNA editing, one of the most promising means of post-transcriptional gene regulation, has been widely investigated in various protozoa [1], mammalian apolipoprotein-B [2], animals [3], fungi [4], bacteria [5,6], and viruses [7,8] as well as in plants [9][10][11]. A-to-I (Inosine) RNA editing is observed in animal nuclear genes, while C-to-U RNA editing is not limited to animals but is also spreading in plant organelles. The mechanism of cytidine-touridine (C-to-U) RNA editing in plant organelles is completely different from that in animal nucleus but also reasonably well understood, mainly owing to the characterization of many RNA editing factors in model systems such as Arabidopsis thaliana and Physcomitrella patens [12]. In flowering plants, the RNA editing machinery, collectively described as the editosome, consists of at least four proteins including pentatricopeptide repeat (PPR) protein, Multiple Organellar RNA editing factor (MORF)/RNA editing factor interacting protein (RIP), Organelle RNA Recognition Motif (ORRM) proteins, and organelle zincfinger protein (OZ1).
PPR proteins constitute a large family of proteins, with more than 400 members [13] and are either directly or indirectly responsible for RNA editing [14][15][16]. Direct selection of target sites is governed by PLS subclass PPR proteins with additional E1 and E2 of contaminating genomic DNA. After DNase treatment, the samples were purified by phenol-chloroform extraction, followed by ethanol precipitation. The purified RNA was quantified using a NanoDrop Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Subsequently, cDNA was synthesized using reverse transcriptase (Superscript III; Invitrogen, Carlsbad, CA, USA) and a random hexamer (oligo dT) primer. The sequences of forward and reverse primers are given in Table S3.

Library Preparation for Transcriptome Sequencing
The mRNA from 12-d-and 20-day-old samples were enriched using oligo (dT) beads. A total amount of 3 µg RNA per sample was used as input material for the RNA sample preparations. Then, total RNA was extracted and was sent to the company, Novogene Co., Ltd., for Next Generation Sequencing analysis. Sequencing libraries were generated using NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligoattached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer 5X. First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H -). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase. After adenylation of 3 ends of DNA fragments, NEBNext Adaptor with hairpin loop structure was ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 150~200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, MA, USA). Then 3 µL USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 • C for 15 min followed by 5 min at 95 • C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index (X) Primer. At last, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. The workflow for library preparation and transcriptome sequencing is shown in supporting Figure S1.

Clustering and Sequencing
The clustering of the index-coded samples was performed on a cBot Cluster Generation System using HiSeq PE Cluster Kit cBot-HS (Illumina) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform and 125-bp/150-bp paired-end reads were generated.

Data Analysis
Quality control. Raw data (raw reads) of fastq format were firstly processed through in-house perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing poly-N, and low-quality reads from raw data. At the same time, Q20, Q30, and GC content from the clean data were calculated, as given in supporting Table S1. All the downstream analyses were based on the clean data with high quality.
Reads mapping to the reference genome. Reference genome (TAIR 10) and gene model annotation files were downloaded from genome website directly. Index of the reference genome was built using Bowtie v2.2.3 and paired-end clean reads were aligned to the reference genome using TopHat v2.0.12. We selected TopHat as the mapping tool, as it can generate a database of splice junctions based on the gene model annotation file and, thus, a better mapping result than other non-splice mapping tools.
Quantification of gene expression level. High-throughput sequencing (HTSeq v0.6.1, University of Heidelberg, Heidelberg, Germany) was used to count the reads' numbers mapped to each gene. Then the FPKM of each gene was calculated based on the length of the gene and reads count mapped to this gene. FPKM, expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced, considers the effect of sequencing depth and gene length for the reads' count at the same time and is currently the most commonly used method for estimating gene expression levels [25]. HTSeq software was used to analyze FPKM, indicating the gene expression levels in this experiment, using the union mode. The resulting files presented the number of genes with different expression levels and the expression level of single genes.
Differential expression analysis (For DESeq with biological replicates). Differential expression analysis of two conditions/groups (two biological replicates per condition) was performed using the DESeq R package (1.18.0). DESeq provide statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting p-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted p-value < 0.05 found by DESeq were assigned as differentially expressed. (For DEGSeq without biological replicates.) Prior to differential gene expression analysis, for each sequenced library, the read counts were adjusted by edgeR program package through one scaling normalized factor. Differential expression analysis of two conditions was performed using the DEGSeq R package (1.20.0). The p values were adjusted using the Benjamini and Hochberg method. Corrected p-value of 0.005 and log 2 (Fold change) of 1 were set as the threshold for significantly differential expression. SNP analysis. Picard-tools v1.96 and samtools v0.1.18 were used to sort, mark duplicated reads, and reorder the bam alignment results of each sample. Genome Analysis Toolkit, GATK2 (v3.2) software was used to perform SNP calling. The mapping status of reads was provided in BAM files, which were visualized using the Integrative Genomics Viewer (IGV) software.

Sanger Sequencing
After doing PCR with equal amounts of cDNA (100 ng) in each reaction of 20 µL volume, the PCR products were purified by 1% agarose gels and the bands were cut out and frozen. DNA was purified using the QIAquick Gel Extraction kit, and concentration was measured by Nano-Drop. Sequencing of the purified DNA was performed using the Big Dye Terminator v3.1 Cycle Sequencing Kit (Thermo Fisher Technologies, Waltham, MA, USA) using the forward and reverse primers (Table S3). The raw sequencing data were analyzed using the Sequence Scanner software version 2 (Applied Biosystems) and DNADynamo software.

Identification and Analysis of DEGs by RNA-Seq
The level of gene expression was measured by determining transcript abundance; the greater the transcript abundance, the higher the gene expression level [26]. In RNA-seq analysis, gene expression level is estimated by counting the number of reads mapped onto genes or exons. The lists of descriptions for all expressed genes are given in Supporting data S1. Read count was proportional not only to the actual gene expression level but also to gene length and sequencing depth. Transcriptome data indicated that a total of 33,641 genes were expressed, of which 2140 were specifically expressed in 12-day-old seedlings' genes and 1485 in 20-day-old seedlings' ( Figure 1A). The correlation coefficient is an important indicator of the reliability of the experiment: the closer the value of the correlation coefficient to 1, the greater the similarity between samples. The square of the Pearson's correlation coefficient (R2) should be greater than 0.92 under ideal experimental conditions. In this study, R2 was greater than 0.8, indicating a slight difference in gene expression between 12-and 20-d-old seedlings ( Figure 1B). Volcano plots were used to infer the overall distribution of differentially expressed genes (DEGs). In experiments without biological replicates, the threshold is normally set as |log 2 (Fold Change)| > 1 and q-value < 0.005. By contrast, in experiments with biological replicates, DESeq eliminates Cells 2021, 10, 635 5 of 15 biological variation; therefore, we set our threshold as adjusted p-value (padj) < 0.05. ( Figure 1C). in 12-d-old seedlings. A total of 2712 genes were differentially expressed, of which 1642 were upregulated and 1070 were downregulated ( Figure 1E), further indicating higher expression in 12-d-old seedlings. All DEGs are listed in Supporting data S2.
To compare gene expression levels under different conditions, an FPKM distribution diagram was used. The final FPKM value represents the mean of biological replicates. In general, an FPKM value of 0.1 or 1 was used as a threshold to determine whether a gene is expressed or not. The number of genes with different expression levels is shown in Fig

Comparison of Nucleotide Differences between Genomic DNA in Database and RNA-Seq of 12-or 20-D-Old Seedlings
Comparison of RNA-seq data of 12-day-or 20-day-old Arabidopsis Col-0 plants to genomic DNA sequence in the database identified 12 types of possible nucleotide conversion patterns in transcripts: G-to-A, C-to-U, U-to-C, U-to-A, A-to-G, C-to-A, A-to-T, G-to- The FPKM is the most well-known method of gene expression estimation in RNA-seq, as it takes into account the effects of both sequencing depth and gene length on read counts. Figure 1D shows that read counts and FPKM values were higher in 12-d-old seedlings than in the control sample (20-d-old seedlings), indicating higher expression of genes in 12-d-old seedlings. A total of 2712 genes were differentially expressed, of which 1642 were upregulated and 1070 were downregulated ( Figure 1E), further indicating higher expression in 12-d-old seedlings. All DEGs are listed in Supporting data S2.
To compare gene expression levels under different conditions, an FPKM distribution diagram was used. The final FPKM value represents the mean of biological replicates. In general, an FPKM value of 0.1 or 1 was used as a threshold to determine whether a gene is expressed or not. The number of genes with different expression levels is shown in Figure 1F.

Comparison of Nucleotide Differences between Genomic DNA in Database and RNA-Seq of 12or 20-D-Old Seedlings
Comparison of RNA-seq data of 12-day-or 20-day-old Arabidopsis Col-0 plants to genomic DNA sequence in the database identified 12 types of possible nucleotide conversion patterns in transcripts: G-to-A, C-to-U, U-to-C, U-to-A, A-to-G, C-to-A, A-to-T, G-to-T, C-to-G, A-to-C, G-to-C, and U-to-G. Among these patterns, U-to-C was the third most predominant after G-to-A and C-to-U. Single-nucleotide base differences are listed in Supporting data S3. RNA-seq analysis revealed 590 different sites, of which 79 sites (13%) represented possible U-to-C conversion. Out of 253 genes showing nucleotide differences, 50 contained possible U-to-C conversion (Figure 2A,B). A list of candidate U-to-C RNA editing sites detected in Arabidopsis seedlings is given in Table 1.
T, C-to-G, A-to-C, G-to-C, and U-to-G. Among these patterns, U-to-C was the third most predominant after G-to-A and C-to-U. Single-nucleotide base differences are listed in Supporting data S3. RNA-seq analysis revealed 590 different sites, of which 79 sites (13%) represented possible U-to-C conversion. Out of 253 genes showing nucleotide differences, 50 contained possible U-to-C conversion (Figure 2A,B). A list of candidate U-to-C RNA editing sites detected in Arabidopsis seedlings is given in Table 1. f candidate U-to-C RNA editing sites detected in Arabidopsis seedlings at different developmental stages.

Reads
Gene  Next, we analyzed the next-generation sequencing (NGS) data of Arabidopsis for expressed PPR genes using the Bioinformatics & Evolutionary Genomics website. The expressed PPR genes are listed in descending order of expression in Supporting data S4. The right column contains the genes containing all nucleotide base conversions. Out of 465 expressed PPR genes, 10 genes including AT3G62470, AT1G50270, AT1G16830, AT1G63080, AT1G06580, AT3G56550, AT1G09820, AT3G53360, AT2G22410, and AT4G32430 showed nucleotide conversion ( Figure 3A). Out of 54 U-to-C variant genes, one gene, AT4G32430, was found as PPR gene ( Figure 3B). The list of expressed genes, PPR genes that differed in base nucleotide conversions, the genes that differed in U-to-C base conversion, and the PPR gene that differed in U-to-C base conversion are shown in Figure 3C.    The list of expressed genes, PPR genes that differed in base nucleotide conversions, the genes that differed in U-to-C base conversion, and the PPR gene that differed in U-to-C base conversion are shown in (C).

Identification of Genes Harboring U-To-C RNA Editing Site
We selected the genes of both samples that had a minimum number of reads to be able to infer an editing event. This minimum number should be reasonably high to minimize the impact of sequencing artifacts. For example, the T-to-C change at position 14,198,871 in AT3G41768 was supported by 647(29%) and 240 (19%) in 12-d-old and 20-d-old seedlings, respectively. In addition, there were some variants that were supported by 100% of the reads in both samples (12-and 20-d-old). Therefore, these are several editing events that seem to be polymorphisms. For the same gene, we found many reads in the same U-to-C conversions. Genes with higher read coverage were further examined for the confirmation of U-to-C RNA editing sites. Genes, such as AT2G16586, AT5G02670, AT5G42320, AT4G16380, and AT5G08740, showed 249, 13, 55, 268, and 146 reads at the converted sites, respectively. Genes showing extremely low reads (0, 2) were also analyzed by RT-PCR. However, very few sites were confirmed as editing events. Because many reads mapped to each U-to-C conversion site, we considered that these nucleotide conversions were caused by RNA editing [27]. The flowchart for methodology for identification of U-to-C RNA editing sites is shown in Figure 4A. Furthermore, we validated the RNA-seq-based candidates experimentally by Sanger sequencing of both genomic, gDNA, and cDNA for all candidate genes. We extracted DNA and mRNA from the same aliquots of seedling samples. By sequencing the paired DNA and cDNA samples and analyzing each chromatogram by two individuals independently we confirmed the U-to-C RNA edited sites. The cDNA showed a double peak, representing T and edited C nucleotides, while no double peak was observed in gDNA sequencing. The sequencing was performed using sense primer targeting at the editing sites. Validation using PCR and Sanger sequencing verified seven genes, including AT2G16586, AT5G42320, AT5G02670, AT3G41768, AT4G32430, AT3G47965, and AT5G52530, containing U-to-C RNA editing sites. The Sanger sequence chromatograms for all seven edited genes are showed in Figure 5. The raw sequencing data were analyzed using the software, DNADyanamo and Sequence Scanner version 2 (Applied Biosystems). When the edited and unedited products were presented together, a dual peak (T (unedited) and C (edited)) was observed at the target site. The editing efficiency was calculated from peak area and a list of genes showing U-to-C RNA editing in 12-d-and 20-day old Arabidopsis seedlings is given in Table 2. Furthermore, we also investigated the editing efficiency at different developmental stages of Arabidopsis, such as four days and eight days. It was found that no editing occurred at early stages of development, like in four days, while a few editing could be identified in 8-day-old seedlings (Table S3). The U-to-C RNA editing sites were majorly located within the UTRs of mature mRNAs.
Cells 2021, 10, x FOR PEER REVIEW 10 of 1 found as PPR gene (B). The list of expressed genes, PPR genes that differed in base nucleotide conversions, the genes that differed in U-to-C base conversion, and the PPR gene that differed in U-to C base conversion are shown in (C).    ther read (N > 10%). (3) Discard reads when low-quality nucleotides (base quality less than 20) constitute more than 50% of the read. For mapping sequences, TopHat2 was chosen for plant genomes. The mismatch parameter was set to 2 and other parameters were set to default. Appropriate parameters were also set, such as the longest intron length. Only filtered reads were used to analyze the mapping status of RNA-seq data to the reference genome. Edited sites were further validated and confirmed by RT-PCR. (B) Clean reads for day 12 and day 20. (C) Percentage of reads mapped to genome regions for day 12 and day 20. Figure 5. The Sanger sequence chromatogram depicting the U-to-C types of RNA editing events in 12-d-and 20-day-old seedlings from the same tissues of Arabidopsis via cDNA and genomic, gDNA using forward primers. Arrows indicate the position of RNA editing.

2.
List of genes identified with U-to-C RNA editing sites in 12-day-and 20-day-old Arabidopsis seedlings.

Discussion
In our knowledge, this is the first report of U-to-C RNA editing for nuclear genes confirmed by both RNA-seq and Sanger sequencing in flowering plants. In this study, total RNA extracted from 12-d-and 20-d-old seedlings was examined by high-throughputing RNA-seq.
Total RNA isolated from 12-d-old seedlings was examined by NGS, and DEGs were identified based on FPKM values and read counts. The results showed that DEGs were expressed to higher levels in 12-d-old seedlings than in 20-d-old seedlings. This was confirmed by higher FPKM values and read counts and more upregulated genes in 12-d-old seedlings than in 20-d-old seedlings. The ANOVA test was performed for comparing the gene expression levels. The summary for regression analysis of differentially expressed genes among the replicates of 12-d-and 20-d-old seedlings is given in Table S5. Additionally, PPR genes were also expressed to higher levels in 12-d-old than in 20-d-old seedlings, as indicated by the log 2 FC values. These data suggest that DEGs are more likely to be expressed in young Arabidopsis seedlings than in older seedlings. Therefore, more mutations could occur at this stage of development because RNA editing events are more frequent in seedlings than in any other plant tissue.
While investigating for RNA editing events to create a global map of high-quality candidates, an appropriate balance between sensitivity (identifying a highly inclusive set of possible edits) and specificity (being more confident that a call is, in fact, a true RNA edit) is required. We considered it better to have a fewer number of candidate RNA editing events that are more likely to be true than to have a larger number with an increased percentage of false positives. We undoubtedly did not score a substantial number of true, low-level, U-to-C RNA editing events in the process. Up to 90% of nucleotide variants that are not SNPs (either in dbSNP or private genomic SNVs) are U-to-C calls; this suggests they are likely to be U-to-C editing candidates. Furthermore, more than 85% of these candidates are located in UTRs. Our candidate U-to-C RNA editing sites had a different variant frequency from known SNPs. They tended to cluster predominantly in the untranslated regions.
We investigated single-nucleotide base changes and the percentage of read coverage was calculated (Table 1). We predicted 12 types of nucleotide differences, including possible U-to-C conversions. RT-PCR products of the genes including the candidate U-to-C conversions were subjected to Sanger sequencing. A total of seven genes, AT2G16586, AT5G42320, AT1G05670, AT3G41768, AT4G32430, AT3G47965, and AT5G52530, were identified as targets for U-to-C RNA editing ( Table 2). The UTRs of genes encoding proteins involved in RNA metabolism and RNA binding, including PPR proteins, Zn-finger (ZnF)-related proteins, ribosomal protein L2, transmembrane proteins, and two hypothetical proteins, were identified as target of U-to-C editing. Interestingly, the ribosomal RNA, AT3G41768, was identified for 45.65% of U-to-C RNA editing efficiency. Since about 50% of genes are affected with editing, it might had had significant effect on their functions. Similarly, the transmembrane protein, AT2G16586, was identified with 77.3% of U-to-C RNA editing efficiency, which may affect its general physiology. In addition, the PPR gene, AT4G32430, was also identified with 20.43% U-to-C RNA editing.
While RNA editing in introns or UTR regions can affect mRNA stability, translation, or splicing activity because of the modification of its secondary structure, those in coding region can also affect the translated polypeptide sequence [28][29][30]. In this study, we demonstrated that most U-to-C RNA editing events are located in UTRs, which may affect the secondary structure and, consequently, the stability of mRNA.
In Arabidopsis, C-to-U and U-to-C RNA editing have been reported at the translation borders of nuclear transcripts, AT1G29930.1 and AT1G52400.1 [31]. These deamination(C-to-U) and amination (U-to-C) events accumulated at adjacent sites; therefore, the possibility that the deamination reaction serves as the amino group donor for the amination reaction was proposed, although the frequency of amination was higher than that of deamination [31]. Although this hypothesis is attractive, we could not detect the same RNA editing events in our RNA-seq data. Thus, the amino group donor of the U-to-C amination in plants is unclear. However, in cDNA AT3G47965 there is also a small T superposing with the C just upstream the edited T, showing the possible immediate donor of amino group. Previously, an extensive research on editing sites in nuclear transcripts for mRNA by Parallel Analysis of RNA Ends (PARE) and Massively Parallel Signature Sequencing (MPSS)data was reported. It showed that all 12 RNA editing patterns may exist in the nuclear genes and that perhaps the numbers of editing sites in a specific pattern may vary. The study suggested that RNA editing is an essential RNA-based regulatory layer not only for mitochondrial and chloroplast genes but also for nuclear genes. However, a global vision of RNA editing in plant nuclear protein-coding transcripts has not been realized. Therefore, this work intended to uncover the occurrence of RNA editing events in the nuclear genes of Arabidopsis.
We further compared the gene expression levels for seven identified U-to-C RNA editing target genes among different tissues ( Figure S2). The green bar shows the genes expressed in seedling stage of development of Arabidopsis. The day-specific characteristic of the U-to-C RNA editing events implied that these were post-transcriptional modifications, not genomic mutations. These editings were identified as a growth-dependent RNA editing efficiency alteration. Day 4 seedlings did not have RNA editing, at least (Table S4). It indicates that the enzyme important for this editing events might have been expressed at defined stages of seedling development. Next, to validate whether the identified RNA editing sites were true positive, we searched for evidence of the identified RNA editing sites in Arabidopsis RNA-seq data generated by public laboratories, using online software http://signal.salk.edu/atg1001/3.0/gebrowser.php, accessed on 5 January 2021. All seven identified U-to-C RNA editing sites AT2G16586, AT5G42320, AT5G02670, AT3G41768, AT4G32430, AT3G47965, and AT5G52530 were aligned against the publicly available RNA-seq databases (Supporting data S6) and confirmed our findings. The target T sites were identified as edited C sites in various databases. The comparative analysis of Arabidopsis RNA-seq is shown in Figure S3. The edited sites are indicated within red boxes. Further studies are needed to better understand the processes involved in U-to-C RNA editing, including the identification of cis or trans regulatory elements, isolation of editing enzymes, and validation of editing sites.

Conclusions
Our findings confirm the uridine-to-cytidine RNA editing sites in some nuclear genes in Arabidopsis thaliana. A comprehensive analysis of RNA-seq data to detect nucleotide base conversions was performed. In this study, we examined U-to-C RNA editing in Arabidopsis seedlings at different developmental stages. Sanger sequencing identified the sites and efficiency of seven U-to-C editing events. Most U-to-C RNA editing here identified occurred in the UTR of mature mRNAs. Thus, we confirmed the presence of U-to-C RNA editing in nuclear genes of plants. We provided the experimental basis to explore the mechanism involved in the amination of U-to-C editing and functions and effects of U-to C RNA editing on mRNA stability, other RNA modifications, and translation.
Supplementary Materials: The following are available online at https://www.mdpi.com/2073-4 409/10/3/635/s1. Figure S1: The workflow for library preparation and transcriptome sequencing. Figure S2: Comparative analysis of gene expression levels for seven identified U-to-C RNA editing target genes among different tissues. Green bar shows the genes expressed in seedling stage of development of Arabidopsis. Figure S3: Validation of target U-to-C RNA editing sites on Arabidopsis RNA-seq database. Table S1: Data table for quality control. Table S2: List of candidate U-to-C RNA editing sites detected in Arabidopsis seedlings showing the percentage of read coverage. Table S3: List of candidate genes for U-to-C RNA editing sites in Arabidopsis seedlings at different developmental stages showing the forward and reverse primer sequences. Table S4: List of genes identified with U-to-C RNA editing in Arabidopsis seedlings at different developmental stages. Table  S5: Summary for regression analysis of differentially expressed genes among the replicates of 12-dayand 20-day-old seedlings. Supporting data S1: Descriptions for all genes. Supporting data S2: Lists of DEGs. Supporting data S3: Lists of single-nucleotide conversions. Supporting data S4: Expressed PPR gene lists. Supporting data S5: Description for U-to-C conversion. Supporting data S6: RNA-seq database table.