Rice (Oryza sativa L.) is the second most widely grown cereal and one of the most important food sources for human nutrition. Most of rice production is concentrated in Asia. However, Brazil is the major producer outside this continent [1,2]. It is also important to notice that, apart from its social and economic importance, rice has been widely used as a model organism to investigate many aspects in plant biology. Thus, much of the knowledge acquired from this species can also be applied to other monocots. Several minerals are essential for plant growth and most of these are usually uptaken by roots, directly from the soil solution. Among these minerals, iron (Fe, from the Latin word "ferrum"), is a vital element for plant and animal life [3,4]. However, iron can be harmful when present at high concentrations [5], reason why its content is subjected to a precise control by the cells [6,7]. Iron excess is one of the most important constraints to rice production on soils irrigated by flooding. The low oxygen environment favors the reduction of this mineral, resulting in ferrous iron (Fe⁺²) [5], a quite soluble form present in higher concentrations in the soil solution. This excessive iron can be easily absorbed and transported along the xylem, carried by the transpiration stream, leading to the generation of reactive oxygen species (ROS) in leaves, where the first symptoms are observed. Different physiological disorders caused by this toxicity culminate with big losses in crop production [8,9]. To survive unfavorable conditions, plants actively employ pre-mRNA splicing as a mechanism to regulate expression of stress-responsive genes, reprogramming intracellular regulatory networks [10]. RNA splicing is a biological process responsible for removing introns from pre-mRNA, and joining exons. Splicing can lead to the generation of multiple protein isoforms from a single gene due to alternate use of 5’ and/or 3’ splice sites by the spliceosome, a phenomenon called alternative splicing (AS) [11,12]. In this sense, splicing sites, the nucleotide sequences surrounding exon-intron boundaries that determine the action of the spliceosome have an important role [13]. AS is an important post-transcriptional regulatory mechanism that modulates gene expression and, eventually, protein forms and functions [8,9]. This process consists in a posttranscriptional modification responsible for increasing the diversity of proteins generating two or more mRNA variants of a single gene causing the inclusion or exclusion of peptide sequences, and modulation of gene expression by producing mRNA variants [14]. Plant AS events have not been characterized in many different conditions, but there is strong evidence indicating that these can promote fast changes that can contribute to species-specific differences and adaptation, performing important roles in many events, including stress response [15-18]. Recent studies suggest that environmental stresses can induce AS or alter its efficiency and fidelity in a number of genes playing a role in plant stress response and tolerance [19-22]. The progress achieved in rice transcriptome analysis has enabled researchers to understand expression patterns and their relationship to function and regulation, as well as the relationship between global transcription and chromosome features [23]. Recent studies pointed out the extensive diffusion of these phenomena in plants and its importance in gene expression and stress response [24]. RNAseq technique has been used to achieve different objectives, including gene prediction [25,26], isoform identification and quantification [19,27-30], as well as discovery of non-coding transcripts [31-33]. Therefore, understanding transcriptome dynamics is essential for unveiling functional elements of the genome and interpreting phenotypic variation produced by combinations of genotypic and environmental factors [28]. Due to the little information available on this subject, in this study, we explore changes in splicing patterns that occur in response to iron excess, aiming to evaluate and quantify different kinds of junctions and alternative splicing events in the transcriptome of an iron-tolerant rice indica cultivar (BRS Querência) under high concentrations of iron.

Plant material

Seeds of rice (cv. BRS Querência) were germinated in a growth chamber for 7 days with a photoperiod of 16 hours and a temperature of 25 ± 2°C. After this period, seedlings were transferred to plastic trays (3 L) containing pre-washed sand, and kept in a greenhouse where irrigation, with water and nutrient solution [34], was maintained. Plants were subjected to iron treatment at the seedling stage [35], adding 300 mg L^-1 Fe⁺² to the nutrient solution [36], while untreated plants remained in normal nutrient solution (control). Plant material was collected 24 hours after iron treatment, frozen in liquid nitrogen and stored at -80°C for later RNA extraction.

Preparation of cDNA and sequencing

Total RNA was extracted from 100 mg of leaves, using Purelink® Plant RNA Reagent. A TruSeq RNA Sample Preparation v2 (Illumina®) kit was used for the preparation of libraries, following manufacturer's recommendations. RNA-Seq analysis was performed using Illumina HiSeq 2500 platform (Illumina®) with paired-end 2 × 100 reads.

Data analysis

In order to evaluate read quality, the software FastQC Ver. 0.11.2 [37] was used. The Trimmomatic Ver. 0.32 [38] software was used to remove the adapters and low quality base reads of each library. Then, reads were mapped against the reference genome of Oryza sativa cv. Nipponbare (IRGSP build 1.0) from The Rice Annotation Project Database (RAP-DB) [39], using the software TopHat Ver. 2.0.11 [40] with the argument “--report-secondary-alignments” to output additional or secondary alignments. TopHat uses software Bowtie Ver. 0.12.7 [41] to map the reads. To discover splice junction we mapped reads in the reference with the MapSplice program [42] using default parameters. This software detected novel canonical, semi-canonical and non-canonical splice junctions. After, a Python package Splice Grapher [43,44] was used to predict alternative splicing patterns by loci. For this analysis, first we reduced gene models and filtered alignments to reduce false-positive (with sam_filter.py), after we predicted splice graphs for all transcripts (with predict_graphs.py). Finally we compared the predicted with complete gene models to see how closely the predicted graphs match the complete gene models (with gene_model_to_ splicegraph.py). Finishing the analysis, MapMan Ver. 3.5.1.R2 [45] and the Predicted Rice Interactome Network - PRIN [46] were used to detect changes in metabolism, and in protein-protein interactions respectively. The identification of differentially expressed genes (DEGs) was performed using edge R Ver. 3.8.5 [46-48]. The expression levels were normalized by reads per kilobase per million reads (RPKM) method, with p-value < 0.01.

A total of 127,781 and 123,682 junctions were found in plants under normal and stressful conditions, respectively. When these junctions were distributed into the three main described categories (Figure 1), significant differences in junction site ratios were found. As expected, the majority of junctions belong to canonical (98.85% and 98.91%), followed by semi-canonical (0.73% and 0.70%) and non-canonical (0.42% and 0.40%) sites in control and stressed plants, respectively. The lower number of junction sites found in iron stressed rice could be explained by the high tolerant genotype used in this study. The AS is predominant in some gene families, while absent in others [49]. A larger number of genes expressed in plants under stressful conditions may belong to families that did not present AS. The analysis of 5’ and 3’ splice sites in all introns of Arabidopsis and rice indicates that these are very similar to human, with just some subtle differences in the frequencies of specific nucleotides at specific positions. In Arabidopsis, non-canonical splicing sites occur in only 0.7% of all splicing sites [13]. Other reports demonstrated that 97.5% of splicing sites are canonical GT-AG pairs, while 1% is GC-AG and 1.5% a combination of less frequent non-canonical sites [24]. Non-canonical AS comprised “variations affecting multiple exons”, which refer to deletions of large coding sequences spanning several exons, intraexonic deletions and generation of chimeric mRNAs. These join together portions of two separate mRNA molecules leading to frame shifting, which is caused by an intron excision, therefore a lower number of non-canonical sites junctions can be a defense strategy against the stress in this cultivar [10]. The four major kinds of alternative splicing (intron retention, exon skipping, alternative 5’ splicing site, and alternative 3’ splicing site), that cover the vast majority of the alternative splicing events [50], were investigated in control and Fe overload plants. A deeper alternative splicing analysis showed that the control plants presented a higher number of events (23,307) when compared to the treated ones (8,244). The AS events were organ-specific, indicating a strong association of AS events with organ- specific regulation and a major role for AS in the functional complexity of plants [28], being this one of the possible reasons that for having a low number of events in the treated samples. Furthermore, it is interesting to note that intron retention was the most frequent event in both situations, with 10,281 (44.1%) for control and 3,909 (47.4%) for stressed (Figure 2). Similar results were also detected in others studies, demonstrating that intron retention is the most frequent type of AS in plants [34,36,51]. The second major event is alternative 3’ splice site (5,261/22.6% and 1,802/21.9%), followed by exon skipping (4,413/18.9% and 1,429/17.3%) and alternative 5’ splice site (3,352/14.4% and 1,104/13.4%). The increase in intron retention percentage may be a reflection of the underlying mechanism to create novel function acting on many other genes in plants. The stress condition can affect the plant mechanisms, even with low effect, regulating splicing by multiple mechanisms, including alteration of the population or distribution of splicing factors, or induction of changes in phosphorylation status or expression of serine-arginine-rich protein [49]. The intron retention is the AS events observed with a presence of ~40% in Arabidopsis and rice, while only 9% is observed in human [52]. The most abundant splicing event in human is exon-skipping (42%). In plants, this event is considered relatively rare, as observed in experimental data in our study. This suggests that the splice regulatory components can influence the splicing mechanisms making difference between plant and animals. The obtained data fits the reported hypothesis [53] that organisms having typically small introns use an intron-definition splicing mechanism, performing most of AS events as intron retention, and that organisms that have large introns use an exon definition mechanism (exon skipping). Similar results were already found in which intron retention represents 47% of all splicing events in rice [28] and 77% in grape [24]. The results obtained are consistent with other studies [19,27,28,51] supporting the idea that intron retention is a common event in plants. Intron retention events and alternative 5’ and 3’ splice sites can produce more mRNAs with premature termination codons (PTCs). The presence of PTCs may result in mRNA degradation via the nonsense-mediated mRNA decay (NMD), and translation of truncated proteins [22,54,55]. Truncated proteins can have important roles in plant stress adaptation, once that truncated proteins derived from the PTC containing mRNAs are not necessarily functionless when compared to the full-length protein [37,56,57]. AS also affects a range of regulatory genes that are presumed to have important roles in abiotic stress adaptation [10]. Its influence is present in almost all aspects of protein functions, making it a central part of gene expression regulation. Recent studies suggest that epigenetic regulation not only determines which parts of the genome are expressed, but also how they are spliced. ROS and reactive nitrogen species (RNS) at an increased rate during stress, participating in signal transduction, modify cellular components and cause damage, being important groups of genes to analyse [58]. In this case, changes promoted by AS on the expression of proteins related to post-translational modifications were analysed. These proteins are responsible for diversifying the functions of other proteins and for the dynamic coordination of signaling networks [40]. A total of 25 differentially expressed genes were found associated to posttranslational modifications (Table 1). Of these, the first five genes (Os01g0631700, Os04g0660500, Os02g0777800, Os07g0693000, Os03g0125600) were upregulated, while the remaining 20 were downregulated. When these proteins were analysed, kinases and phosphatases were found to be the most common classes, where 14 assignments were found. These proteins are responsible for two important biological processes: phosphorylation with 14.29% of the assignments (with a total of three down-regulated genes) and dephosphorylation with 85.71% of the assignments (where 14 genes were down-regulated and four were upregulated). The most common molecular functions were transferase activity and ATP binding. In both cases, 16 downregulated and four upregulated assignments were observed. A change of scenario in genes responsible for posttranslational modifications can be seen when plants are subjected to stress. This change is perceived by decreasing the expression of most of these genes and by the activation of some other post-translational genes. These changes can be a distinct response of plants to stress. The activated post-translational modification genes are distributed into two main groups: protein kinases (Os01g0631700, Os04g0660500 and Os03g0125600) and receptor-like cytoplasmic kinases (Os02g0777800 and Os07g0693000). Protein kinases play key roles in most cellular activities, performing the phosphorylation in protein side chains, a process that usually results in a functional change of the target protein by changing enzyme activity, cellular location, or association with other proteins [59,60]. Receptor-like cytoplasmic kinases (RLCKs) have important roles in plant development and stress responses, where differential expression patterns suggest their involvement in diverse functions in rice [61]. These reports aid the understanding and confirmation of results obtained. Plants respond to adverse conditions through a series of signaling processes that often involves diverse protein kinases. We found calcineurin B-like protein-interacting protein kinases (CIPKs), mitogen-activated protein kinases (MAPKs), stress-activated protein kinase and receptor-like cytoplasmic kinases were downregulated. Differential expression of these genes had already been detected for stresses such as drought, salinity, cold, polyethylene glycol, or abscisic acid treatment [62-65]. Just one locus (LOC_ Os02g34600) showed interactions with other differentially expressed genes. Correlations with other six loci (Figure 3), five negative (coexpression values between -0.1879 and -0.0210), and one positive correlation (0.1756) were found. Metabolomic correlation networks, indicated negative co-expression between folding and post-translational modification proteins (Figure 3). Post-translational modifications had positive co-expression with proteins of RNA processing and also had negative co-expression with ribosomal protein synthesis and with some proteins with unknown metabolism (non-assigned).

Figure 1: Alternative splicing analysis in seedling stage rice leaves (Oryza sativa ssp. indica cv. BRS Querência) after 24 hours of iron excess exposure. (A) Types of alternative splicing junctions sites. (B) Pie chart showing the percentage distribution of alternative splicing events in control and stress condition (300 mg L^-1 Fe⁺²).

Figure 2: Gene ontology of differentially expressed genes related to posttranslational modifications in rice (Oryza sativa ssp. indica cv. BRS Querência) under iron excess (300 mg L^-1 Fe⁺²).

Figure 3: Interactions between genes involved in post-translational modifications in rice leaves at seedling stage (Oryza sativa ssp. indica cv. BRS Querência) after 24 hours of exposure to high concentrations of iron (300 mg L^-1 Fe⁺²). (A) Protein-protein interaction networks. Gray lines represent co-expression less or equal to 0.5 and red lines represent co-expression above 0.8. The circle of colors representing their possible subcellular localization: navy blue to mitochondria, dark blue for the nucleus, pink for cytoplasm and gray for unknown localization. (B) Co-expression metabolomics network.

Table 1: Genes related to post-translational modifications which are differentially expressed in rice (Oryza sativa ssp. indica cv. BRS Querência) under iron excess (300 mg L^-1 Fe⁺²).

Rice plants of cv. BRS Querência under stress with concentration 300 mg L^-1 Fe⁺² during 24 hours change the pattern of alternative splicing in leaves. Also, changes in the expression of post-translational modification genes occur. Changes of AS events may indicate a new tolerance strategy of plants. This work represents an important step for understanding how rice plants respond to stress and how this response can increase tolerance in this species. Further investigations are necessary in order to elucidate the influence of these up-regulated genes on plant tolerance responses under stress by iron overload. Alternative transcripts identified here constitute important targets for molecular breeding through genetic engineering and marker development.

Research supported by the following Brazilian funding agencies: Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and Fundação de Amparo à Pesquisa do Rio Grande do Sul (FAPERGS).

The authors have no potential competing interests to declare.