Journal of Horticulture
Open Access

Genome-Wide Analysis of mRNA Splicing vaRiants in Higher Plants

Ying Li, Qianhuan Guo, Chengchao Zheng, Shizhong Zhang^* and Kang Yan^{* *}Correspondence: Shizhong Zhang, College of Life Sciences, Shandong Agricultural University, Taian, Shandong 271018, P.R. China, Tel: 86-538-8242894, Fax: +86-538-8226399, Email: Kang Yan, College of Life Sciences, Shandong Agricultural University, Taian, handong 271018, P.R. China, Tel: +86-538-8242894, Fax: +86-538-8226399, Email:

Author info »

Introduction

Alternative splicing (AS) is a common and fundamental process that contributes to both transcriptome and proteome diversities [1]. AS generates two or more mRNA variants from a single pre-mRNA with multiple introns through different splice sites [2,3]. Most plant genes (80% to 85%) are interrupted by introns [4]. However, the intron and exon organization of higher plant genes is not similar to that of animal genes [5]. Introns in plant genes are generally much shorter in length and fewer in number compared with those in animal genes [6]. In contrast to reports on the differences between plant and animal species, studies that compare the gene structures of different plant species such as food crops and horticultural crops are scarce.

Moreover, AS has been found to be an important regulatory process in different cell types, different developmental stages, and environmental responses. Therefore, individual mRNA variants may play specific spatial or temporal roles [7,8]. The crucial role of AS has been extensively investigated in animals and plants [9,10].

AS is associated with the human genetic diseases, and also involved in a range of functions in plants, such as seed germination, stress response [11,12]. In eukaryotes, AS of pre-mRNAs significantly contributes to the proper expression of the genome and results in new protein products (13). However, not all of the products are functional [13,14]. Several AS variants contain premature termination codons (PTCs) that potentially lead to unproductive transcripts, truncated proteins or mRNA decay [15,16].

Conversely, recent studies indicated AS act as a regulatory mechanism influencing the expression of non-coding RNA [17]. Primary transcripts of miRNAs (pri-miRNAs) contain introns, and their AS events has been detected [18,19]. In Arabidopsis, AS event disrupts the secondary structures of pri-miR162a and provides a mechanism for miR400 expression in response to environmental cues [20,21]. Tomato as a kind of very important horticultural crops, research has found that the AS of miR4376 through regulating the expression of an autoinhibited Ca²⁺- ATPase, ACA10, affect the reproductive growth of tomato [22,23]. Thus, growing evidence has indicated that AS in plants is much more prevalent than previously thought and plays crucial roles in the diversification of gene regulation.

Genome-wide studies of AS events in various organisms have relied on using the traditional expressed sequence tag and cDNA libraries [24]. However, the occurrence of AS events might be underestimated because data from these libraries do not provide a complete coverage of the transcriptome. High-throughput RNA sequencing (RNA-seq) technologies allow transcriptome analyses at unprecedented levels of sensitivity and precision. Such technologies also provide large data sets for comprehensive analyses [8].

Studies on plant RNA-seq data sets have identified thousands of AS events and confirmed most annotated introns in several vertebrate species. A previous RNA-seq analysis revealed that approximately 95% of human multi-exon genes express multiple splicing variants [25]. Recently, it was shown that 61% of Arabidopsis intron-containing genes and 56% of maize multiexonic genes are subject to AS [26,27]. A large number of variable splicing events of genes have been found in horticultural crops. For example, in strawberries, 66.43% detected multi-exon genes undergo AS [28,29]. And in cucumber, recent studies have found that 58% of the multi-exon genes underwent AS [30]. Also, a large number of genes were alternatively spliced in the soybean genome [31]. To our knowledge, the extent and frequency of plant mRNA variants remain poorly characterized from a global perspective.

According to different AS types, AS events can be mainly classified as exon skipping (ES), intron retention (IR), alternative donor sites (AD), and alternative acceptor sites (AA). IR is the most common splicing form in plants [32]. However, one AS gene variant might simultaneously have multiple different AS types. In this study, to achieve a nonbiased and complete analysis of the plant transcriptome, we analyzed the mRNA splicing variants of 12 plant species (8 dicotyledons and 4 monocotyledons) including food crops and horticultural crops which subjected to high-throughput sequencing database, including Poplar (Populus trichocarpa), Medicago (Medicago truncatula), soybean (Glycine max), Arabidopsis (Arabidopsis thaliana), cotton (Gossypium raimondii), cacao (Theobroma cacao), rose gum (Eucalyptus grandis), blue columbine (Aquilegia coerulea), sorghum (Sorghum bicolor), maize (Zea mays ), rice (Oryza sativa), and Brachypodium (Brachypodium distachyon). All mRNA variants from the annotated genes were confirmed by genome alignment and investigated the relationship between the number of AS mRNA variants and their effects on gene features at the genome-wide level.

Materials and Methods

Database set-up

A total of 12 recently published RNA-seq data sets from higher plants were selected and downloaded at Phytozome, which is the Plant Comparative Genomics portal of the Department of Energy's Joint Genome Institute (https://phytozome.jgi.doe.gov/ pz/portal.html). Detail information about mRNA variants of 12 plant species is provided in the Supplemental section (Table 1 and Supplementary Table 1).

Table 1: The mRNA variants of 12 plant species used in this study.

mRNA variants alignment to the reference plant genome

A total of 589,034 mRNA variants were aligned to perfectly match gene locations in the reference genome with a maximum of two mismatches. Among the mRNA splicing products, only the annotated genes were selected and classified into groups according to the number of mRNA variants. The correlative gene features (introns, coding sequence, exons, or UTRs) were generated, including their distribution, location, length, and numbers. Detail information about the gene features are given in the Supplemental section (Supplementary Table 1).

Clustering of AS gene features

The gene features in the different groups were calculated.

Results

Identification and classification of mRNA variants in higher plants

In general, RNA-seq provides broad and deep sequencing of the transcriptome with low false discovery rates. Therefore, to facilitate further investigation of plant mRNA variants from a global perspective, we mined the RNA-seq data sets of 12 plant species to identify thousands of mRNA variants and all the gene sets from 12 annotation genomes downloaded at https://phytozome.jgi.doe. gov/pz/portal.html (Figure 1 and Table 1). All mRNA variants were aligned to perfectly match the corresponding gene locations in the reference genome (Supplementary Table 1). To ensure the authenticity of each mRNA product, only the annotated mRNA were pick up and classified according to the number of mRNA variants into groups (Table 2). The correlative gene features of introns, coding sequence, exons, or untranslated regions (UTRs) were recorded. These features include distribution, location, length and number (Tables 1- 3 and Supplementary Table 1).

Figure 1: Characteristics of plant annotated gene features and AS variants.
Statistics of the annotated gene features and AS events in 12 plant species. The line represents the evolutionary relationships among various biological species in the phylogenetic tree.

Table 2: Classification of annotated mRNA according to the number of mRNA variants into groups.

Table 3: The correlative gene features of 12 plant species used in this study.

Generally, AS can be classified as IR, ES, AD and AA by the type of AS events. However, one AS gene variant might have multiple different AS types, simultaneously. This means that one AS gene can generate two, three, four, or more variants, all AS variants were classified into four groups on the basis of the number of AS variants: 2V (two variants), 3V (three variants), 4V (four variants), and 5V+ (five or more AS variants) (Table 2).

In total, 442,541 annotated genes from 589,034 mRNA products were obtained in the genome annotation of the 12 plant species (Figure 1 and Table 1). Our analysis identified 224,593 annotated mRNA variants. The events were distributed in 78,190 AS genes, which accounted for 17.67% of the total annotated genes (Table 1 and Supplementary Table 1). To eliminate the bias effect of the different species, we calculated the AS gene number and ratio (AS gene number/total gene number) for each species. We generated 3,063 genes in Populus (7.5%), 1,797 in Medicago (3.5%), 6,611 in soybean (14.3%), 5,804 in Arabidopsis (21.2%), 15,403 in cotton (41.1%), 7,630 in Theobroma (25.9%), 6,110 in Eucalyptus (16.8%),7,487 in Aquilegia (30.2%), 1,489 in sorghum (5.4%), 13,137 in maize (33.8%), 6,458 in rice (11.6%), and 3,201 in Brachypodium (12.6%) from the pooled data (Figure 1). A comparison of the AS variants among these species revealed that the AS gene ratios were higher in cotton and maize but lower in Populus, Medicago, and sorghum than in the other species (Figure 1 and Table 1).

Global analysis of gene structure

Gene structure analysis from a global perspective differentiates the architecture of plant and animal genes. In general, the composition of introns differs in plants and animals. Therefore, we evaluated the characteristics of gene structure to determine the putative differences in splicing among these plant species (Table 1).

Our analysis revealed that the composition of introns differed in dicots and monocots, especially in terms of the number of introns. The average number of introns in the different dicotyledonous plants ranged from approximately 2.5 to 5.8 (Table 1). This wide range indicated their divergence despite the close phylogenetic relationship. Interestingly, each gene contained an average of four introns in the monocotyledonous plants. These numbers were relatively stable (Table 1). However, the average length of each intron ranged from 0.16 kb to 0.72 kb and greatly differed across the different species. We also found that more than 85% of the plant introns were located in the coding regions (Table 1 and Table 3).

Unlike the composition of introns, the average number and length of exons were highly similar in plants with subtle differences. These differences determined the sizes of the coding regions and proteins (Table 1). The above mentioned exon analysis was consistent with the similar average size of proteins in multiple species (Table 3). This finding explained why the different average gene sizes of various plant species generated similar protein sizes. In addition, the average sizes of UTRs (both 5′UTR and 3′UTR) differed among plant species. The length of the 3′UTR was much longer than that of the 5′UTR, except in Medicago (Table 1 and Table 3).

To gain clues about the features of multiple AS variants in 4 variant groups, we checked the assembled raw numbers and frequencies of AS among the different plant species. Group 2V had the most abundant AS variants, followed by groups 3V, 4V, and 5V+. A high proportion of AS variants in group 2V were generated from each species: Populus, 76.4%; Medicago, 70.8%; soybean, 71.7%; Arabidopsis, 73.2%; cotton, 39.3%; Theobroma, 51.9%; Eucalyptus, 62.5%; Aquilegia, 47.9%; sorghum, 81.8%; maize, 55.0%; rice, 63.8%; and Brachypodium, 72.0% (Figure 2 and Table 2). It is noteworthy that almost more than 50% pre-mRNAs which underwent alternative splicing events produced two variants in plants, whereas the proportions decrease with the increasing numbers of variants in other groups (Figure 2).

Figure 2: Features of multiple AS variants in plants.
From outside to inside, each circle represents the gene model of multiple AS variants, the distribution of AS variants along annotated genes features, and the proportion of multiple AS variants in plants. See also the Supplemental section (for details including the gene structure, location, number, length, and number of AS events).

Examination of the characteristics of gene structure from different AS groups in plants showed that the average number of exons (ranged from 6.25 to 9.74) and introns (ranged from 4.65 to 9.64) in AS genes were more than the average number of the total annotated genes in each genome, but the average lengths of these exons and introns were slightly shorter (Figure 2 and Table 2). Additionally, the average sizes of UTRs were longer in the genes that generated multiple AS variants. These results implied the presence of a relationship between gene structure and AS. However, the total length of annotated genes did not significantly change (Table 2).

Distribution and frequency of AS variants

The above-mentioned analysis suggested that gene features such as the number and length of introns can influence AS variants. To further explore other genomic features, we performed correlation analysis to investigate the relationship between the AS variants and gene density of higher plants. The distributions of the AS variants in the four groups were pooled in 1 Mb contiguous regions across the chromosome where the genes are located (Figure 3). The analysis revealed extensive transcriptional activity in the chromosomes of the different plant genomes. The AS variants were principally located in the chromosome arms and were associated with gene distribution. As expected, the AS variants in group 2V that matched with multiple locations were predominantly mapped to the most transcribed chromosomal regions. These results were consistent with the high number of group 2V variants in all AS genes (Figure 3).

Figure 3: Statistics of multiple AS variants and AS distribution with gene density along chromosome length in plants.
Each vertical bar represents the frequency of an annotated gene in the plant genome. A schematic of each chromosome and its features is shown below. Different colors represent different groups on the basis of the number of AS variants. Red, two variants; deep blue, three variants; yellow, four variants; green, five or more variants.

A comparison of the AS variants to the annotated plant genome showed that most of the AS variants that matched the genome can be mapped to annotated genic regions in plant species. However, approximately 50% of the AS variants cannot be mapped to the right chromosome in Aquilegia (Figure 3). Two reasons may explain this discrepancy. First, the depth of coverage of the Aquilegia genome was lower than the gene location coverage. Second, information regarding the annotated genes in each particular chromosome was limited. Overall, the distribution of multiple AS variants significantly positively correlated with gene density. Similar correlation patterns were observed in all plant species.

Multiple AS variants exert different effects toward improving proteome diversity

AS variants improve gene representation at the transcriptional and proteomic levels. The overall proteome from all the pooled data sets contained 211,558 protein isoforms, whereas each AS gene generated 2.87 mRNA variants and 2.71 protein variants. The number of variants per AS gene was not significantly different among the different plant species (Figure 1 and Table 1). These data indicated that AS enhanced transcriptome plasticity and proteome diversity in higher plants.

To explore the effects of multiple AS variants on proteome diversity, we determined the number of protein isoforms that were generated by the AS variants in the four groups. Intriguingly, the proportion of the effective AS events to generate various proteins was the highest in 2V group, but not in 4V or 5V+. Within group 2V, approximately 30% AS genes generated only one protein product, and similar proportions were observed for all plant species, except in Medicago (Figure 4). By contrast, approximately half of the genes from group 3V generated three different proteins, whereas the remaining genes produced one or two proteins (except in Medicago and maize) (Figure 4). The proportion of effective AS events noticeably increased with the number of mRNA variants. The observed proportion was even higher in groups 4V and 5V+ (Figure 4). In brief, a high proportion of the annotated genes which generated the isofroms numbers less than mRNA variants were found, indicating that the effects of AS on the improvement of proteome diversity were lower than previously thought, especially in groups 3V, 4V, and 5V+. The AS variants in the different AS groups exerted distinct influences on proteome diversity (Figure 4).

Figure 4: Comparison of studies on proteome diversity with multiple AS variants. From outside to inside, each circle represents the classification of AS as defined by the number of AS variants.

Discussion

In this study, the total number and ratio of AS were below those of AS in previous studies on some species. Among the multiexonic genes, approximately 42% to 61% of those in Arabidopsis, approximately 56% of those in maize, and approximately 63% of those in soybean underwent AS. Two possible reasons may explain this result. First, several mRNA variants occurred in specific developmental states might have triggered in response to environmental conditions. Second, the basic information of the annotated gene in some species would be limited. Therefore, alignment quality, low transcript abundance, and transcriptome coverage directly influenced the number of annotated AS genes.

We noticed that the features of multiple AS variants in 4 variant groups differed in most plants food and horticultural. With the increasing numbers of variants, the AS genes contained more but shorter exons and introns, implying AS might improve gene representation by changing the extent and complexity of gene structure (Figure 2 and Table 2). We speculated that these trends may be related to the evolution of various species from a long-term perspective.

Our analysis suggested that a high proportion of AS events might occur in the UTRs. We hypothesized that a large proportion of the AS events in the UTRs plays crucial roles in gene regulation [33,34]. Growing evidence has indicated that splicing is a regulatory mechanism for noncoding RNA [35]. In a recent study, we found a stress-induced AS event in the 5′UTR that influences intronic miRNA expression in Arabidopsis [21,36]. However, we also found a few introns in the UTRs (both 5′UTR and 3′UTR). Further investigation should be conducted to determine whether or not the AS variants are governed by the same regulators in the UTRs. In summary, we investigated the relationship between the number of AS mRNA variants and their effects on gene features at the genome-wide level through a comprehensive comparison of 12 plant species, as a valuable resource for functional research on food crops and horticultural crops. The present reanalysis and reclassification of previously reported RNA-seq data sets evaluated the global existence of AS variants in plants and suggested the need for future analyses on mRNA variants [37,38].

Acknowledgments

This work was supported by the National Natural Science Foundation (grant numbers 31570271 and 31870234).

Authors Contributions

Ying Li and Kang Yan performed the research and wrote the manuscript. Shizhong Zhang, Qianhuan Guo and Chengchao Zheng conceived of the study, and participated in its design and contributed to revisions of the manuscript. All authors approved the final version of the manuscript.

References

1. Filichkin SA, Priest HD, Givan SA, Shen R, Bryant DW, Fox SE, et al. Genome-wide mapping of alternative splicing in Arabidopsis thaliana. Genome Res. 2010;20(1):45-58.
2. Wang Z, Burge CB. Splicing regulation: from a parts list of regulatory elements to an integrated splicing code. Rna. 2008;14(5):802-813.
3. Lareau LF, Green RE, Bhatnagar RS, Brenner SE. The evolving roles of alternative splicing. Curr Opin Struct Biol. 2004;14(3):273-282.
4. Brown JW, Simpson CG. Splice Site Selection in Plant PRE-mRNA Splicing. Annu Rev Plant Biol. 1998;49:77-95.
5. Lorkovic ZJ, Wieczorek Kirk DA, Lambermon MH, Filipowicz W. Pre-mRNA splicing in higher plants. Trends in plant science. 2000;5(4):160-167.
6. Reddy AS. Alternative splicing of pre-messenger RNAs in plants in the genomic era. Annu Rev plant Biol. 2007;58:267-294.
7. Moore MJ, Proudfoot NJ. Pre-mRNA processing reaches back to transcription and ahead to translation. Cell. 2009;136(4):688-700.
8. Ramani AK, Calarco JA, Pan Q, Mavandadi S, Wang Y, Nelson AC, et al. Genome-wide analysis of alternative splicing in Caenorhabditis elegans. Genome Res. 2011;21(2):342-328.
9. Braunschweig U, Gueroussov S, Plocik AM, Graveley BR, Blencowe BJ. Dynamic integration of splicing within gene regulatory pathways. Cell. 2013;152(6):1252-1269.
10. Luco RF, Allo M, Schor IE, Kornblihtt AR, Misteli T. Epigenetics in alternative pre-mRNA splicing. Cell. 2011;144(1):16-26.
11. Sugliani M, Brambilla V, Clerkx EJ, Koornneef M, Soppe WJ. The conserved splicing factor SUA controls alternative splicing of the developmental regulator ABI3 in Arabidopsis. The Plant cell. 2010;22(6):1936-1946.
12. Mach J. Alternative splicing produces a JAZ protein that is not broken down in response to jasmonic acid. The Plant cell. 2009;21(1):14.
13. Lewis BP, Green RE, Brenner SE. Evidence for the widespread coupling of alternative splicing and nonsense-mediated mRNA decay in humans. Proceedings of the National Academy of Sciences of the United States of America. 2003;100(1):189-92. 
14. Modrek B, Lee C. A genomic view of alternative splicing. Nature genetics. 2002;30(1):13-9. 
15. Drechsel G, Kahles A, Kesarwani AK, Stauffer E, Behr J, Drewe P, et al. Nonsense-mediated decay of alternative precursor mRNA splicing variants is a major determinant of the Arabidopsis steady state transcriptome. The Plant cell. 2013;25(10):3726-3742.
16. Kalyna M, Simpson CG, Syed NH, Lewandowska D, Marquez Y, Kusenda B, et al. Alternative splicing and nonsense-mediated decay modulate expression of important regulatory genes in Arabidopsis. Nucleic Acids Res. 2012;40(6):2454-2569.
17. Rogers K, Chen X. Biogenesis, turnover, and mode of action of plant microRNAs. The Plant cell. 2013;25(7):2383-2399.
18. Jia F, Rock CD. MIR846 and MIR842 comprise a cistronic MIRNA pair that is regulated by abscisic acid by alternative splicing in roots of Arabidopsis. Plant Mol Biol. 2013;81(4-5):447-460.
19. Szarzynska B, Sobkowiak L, Pant BD, Balazadeh S, Scheible WR, Mueller-Roeber B, et al. Gene structures and processing of Arabidopsis thaliana HYL1-dependent pri-miRNAs. Nucleic Acids Res. 2009;37(9):3083-93.
20. Hirsch J, Lefort V, Vankersschaver M, Boualem A, Lucas A, Thermes C, et al. Characterization of 43 non-protein-coding mRNA genes in Arabidopsis, including the MIR162a-derived transcripts. Plant Physiol. 2006;140(4):1192-204.
21. Yan K, Liu P, Wu CA, Yang GD, Xu R, Guo QH, et al. Stress-induced alternative splicing provides a mechanism for the regulation of microRNA processing in Arabidopsis thaliana. Mol Cell. 2012;48(4):521-531.
22. Wang Y, Itaya A, Zhong XH, Wu Y, Zhang JF, van der Knaap E, et al. Function and Evolution of a MicroRNA That Regulates a Ca2+-ATPase and Triggers the Formation of Phased Small Interfering RNAs in Tomato Reproductive Growth. The Plant cell. 2011;23(9):3185-3203.
23. Ben Abdallah RA, Jabnoun-Khiareddine H, Nefzi A, Daami-Remadi M. Evaluation of the Growth-Promoting Potential of Endophytic Bacteria Recovered from Healthy Tomato Plants.  J Horticul. 2018;05(02).
24. Chen FC, Wang SS, Chaw SM, Huang YT, Chuang TJ. Plant Gene and Alternatively Spliced Variant Annotator. A plant genome annotation pipeline for rice gene and alternatively spliced variant identification with cross-species expressed sequence tag conservation from seven plant species. Plant Physio. 2007;143(3):1086-1095.
25. Pan Q, Shai O, Lee LJ, Frey BJ, Blencowe BJ. Deep surveying of alternative splicing complexity in the human transcriptome by high-throughput sequencing. Nature genetics. 2008;40(12):1413-1415.
26. Marquez Y, Brown JW, Simpson C, Barta A, Kalyna M. Transcriptome survey reveals increased complexity of the alternative splicing landscape in Arabidopsis. Genome research. 2012;22(6):1184-1195.
27. Li P, Ponnala L, Gandotra N, Wang L, Si Y, Tausta SL, et al. The developmental dynamics of the maize leaf transcriptome. Nature Gen. 2010;42(12):1060-1067.
28. Li Y, Dai C, Hu C, Liu Z, Kang C. Global identification of alternative splicing via comparative analysis of SMRT- and Illumina-based RNA-seq in strawberry. Plant J. 2017;90(1):164-176.
29. Verma S, Whitaker V. A New Technology Enabling New Advances in Strawberry Genetics. J Horticul. 2016;3(2).
30. Sun Y, Hou H, Song H, Lin K, Zhang Z, Hu J, et al. The comparison of alternative splicing among the multiple tissues in cucumber. BMC Plant Biol. 2018;18(1):5.
31. Shen Y, Zhou Z, Wang Z, Li W, Fang C, Wu M, et al. Global dissection of alternative splicing in paleopolyploid soybean. The Plant cell. 2014;26(3):996-1008.
32. Chang CY, Lin WD, Tu SL. Genome-wide analysis of heat-sensitive alternative splicing in Physcomitrella patens. Plant Physiol. 2014.
33. Grozeva S. Effect of copper levels in the culture medium on shoot regeneration in pepper. Banats J Biotechnol. 2015;6(12):86-91.
34. Barazesh F, Oloumi H, Nasibi F, Kalantari KM. Effect of spermine, epibrassinolid and their interaction on inflorescence buds and fruits abscission of pistachio tree (Pistacia vera L.), "Ahmad-Aghai" cultivar, Banats J Biotechnol. 2017;8(16):105-115.
35. Idris A. Comparative analysis of 16SrRNA genes of Klebsiella isolated from groundnut and some american type culture collections, Banats J Biotechnol. 2016;7(13):34-40.
36. Hariri Moghadam F, Khalghani J, Moharramipour S, Gharali B, Mostashari Mohasses. Investigation of the induced antibiosis resistance by zinc element in different cultivars of sugar beet to long snout weevil, Lixus incanescens (Col: Curculionidae), Banats J Biotechnol. 2018;9(17):5-12.
37. Ayadi Hassan S, Belbasi Z. Improvemnet of hairy root induction in Artemisia annua by various strains of agrobacterium rhizogenes, Banats J Biotechnol. 2017;8(15):25-33.
38. Kosev VI. Multivariate analysis of spring field pea genotypes. Banats J Biotechnol. 2015;6(11):23-29.