Next Article in Journal
The Reproductive Journey in the Genomic Era: From Preconception to Childhood
Next Article in Special Issue
Rooting in the Desert: A Developmental Overview on Desert Plants
Previous Article in Journal
The First Case of Congenital Myasthenic Syndrome Caused by a Large Homozygous Deletion in the C-Terminal Region of COLQ (Collagen Like Tail Subunit of Asymmetric Acetylcholinesterase) Protein
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Analysis of Whole Transcriptome RNA-seq Data Reveals Many Alternative Splicing Events in Soybean Roots under Drought Stress Conditions

1
Joint International Research Laboratory of Agriculture and Agri-Product Safety, Jiangsu Key Laboratory of Crop Genomics and Molecular Breeding, Co-Innovation Center for Modern Production Technology of Grain Crops, Yangzhou University, Yangzhou 225009, China
2
National Center for Soybean Biotechnology and Division of Plant Sciences, University of Missouri, Columbia, MO 65211, USA
3
National Key Laboratory of Crop Genetics and Germplasm Enhancement, National Center for Soybean Improvement, Nanjing Agricultural University, Nanjing 210095, China
4
Institute of Industrial Crops, Jiangsu Academy of Agricultural Sciences, Nanjing 210014, China
*
Authors to whom correspondence should be addressed.
Submission received: 9 October 2020 / Revised: 10 December 2020 / Accepted: 17 December 2020 / Published: 19 December 2020

Abstract

:
Alternative splicing (AS) is a common post-transcriptional regulatory mechanism that modulates gene expression to increase proteome diversity. Increasing evidence indicates that AS plays an important role in regulating plant stress responses. However, the mechanism by which AS coordinates with transcriptional regulation to regulate drought responses in soybean remains poorly understood. In this study, we performed a genome-wide analysis of AS events in soybean (Glycine max) roots grown under various drought conditions using the high-throughput RNA-sequencing method, identifying 385, 989, 1429, and 465 AS events that were significantly differentially spliced under very mild drought stress, mild drought stress, severe drought stress, and recovery after severe drought conditions, respectively. Among them, alternative 3′ splice sites and skipped exons were the major types of AS. Overall, 2120 genes that experienced significant AS regulation were identified from these drought-treated root samples. Gene Ontology term analysis indicated that the AS regulation of binding activity has vital roles in the drought response of soybean root. Notably, the genes encoding splicing regulatory factors in the spliceosome pathway and mRNA surveillance pathway were enriched according to the Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis. Splicing regulatory factor-related genes in soybean root also responded to drought stress and were alternatively spliced under drought conditions. Taken together, our data suggest that drought-responsive AS acts as a direct or indirect mode to regulate drought response of soybean roots. With further in-depth research of the function and mechanism of AS in the process of abiotic stress, these results will provide a new strategy for enhancing stress tolerance of plants.

1. Introduction

Alternative splicing (AS) is an important posttranscriptional regulatory process for altering the structure of pre-mRNA that modulates the function, structure, and cell location of genes. AS has been found in many tissues and different development stages of eukaryotes and is essential in coordinating the growth and development of plants [1,2,3,4]. It has been reported that developmental stage-dependent AS events dramatically affect rice grain yield [5]. Additionally, many genes have been selected for increased AS complexity during the domestication of maize [6]. Therefore, AS has been assumed to have an important role in the developmental plasticity of plants. Additionally, the alternate protein variants derived from AS events may compete for interactions with the same factors and therefore function differently to cope with environmental changes [7].
With the advent of next-generation sequencing-based approaches, a large number of plant transcripts have been shown to be alternatively spliced under abiotic stress conditions in plants [8,9,10,11]. For example, the number of AS events increased significantly under cold or salt stress in Arabidopsis and cassava [12,13,14]. AS analysis in maize suggested that drought stress induced large developmental splicing changes in a tissue-dependent pattern [9]. AS regulation was shown to coordinate with transcriptional regulation to fine-tune the heat or heat/drought stress response in wheat [15]. In rice, the conserved AS of dehydration responsive element binding protein 2B (DREB2B) introduced a premature translation-termination codon (PTC) into the splicing isoform OsDREB2B1 under non-stress conditions, which resulted in the production of a nonfunctional isoform. However, under high temperature or drought stress, the functional isomer OsDREB2B2 was significantly induced to improve the stress resistance of rice [16]. Wheat WDREB2, an orthologue of OsDREB2B, showed similar AS patterns under drought conditions, which highlighted the conservation of AS modulation during response to drought stress among different species [17,18]. Increased AS complexity in the rice maintainer HuHan2B also contributed to additional drought resistance during the breeding process [19]. Therefore, an in-depth analysis of the regulation mechanism of AS at the post-transcriptional level will help to understand the complex regulation of the plant response to environmental changes [20,21].
Soybean is a major crop that provides oil and protein for humans [22]. A recent study has shown that the majority of RNA splicing events occur in developing soybean embryos [23]. In soybean, more than 63% of multiexonic genes undergo AS [24]. There are five AS types in plants: intron retention (IR), skipped exon (SE), alternative 5′ splicing site (A5SS), alternative 3′ splicing site (A3SS), and mutually exclusive exon (MXE) [20]. IR is the dominant type of AS, even in response to abiotic stress [14,25,26,27,28].
Soybean is a drought-sensitive crop during the seedling stage, and drought stress is a major factor limiting global soybean growth and production [29]. The developmental plasticity of the soybean root system provides an opportunity to improve water and nutrient uptake and usage during environmental changes [30]. However, it is barely known how AS and transcriptional regulation participate in drought response in soybean root, although more recent research in plants provides a new perspective for AS as a promising regulator of root development. For example, the different patterns of AS impact root distribution and are related to the deep root system in rice [31]. Splicing factor 3b involves pre-mRNA splicing associated with root hair development in response to light signals in Arabidopsis [32]. In addition, different AS informs may respond differently to stress. The GATA gene family is involved in nitrate assimilation, and one of the members, OsGATA23a, showed high expression levels under salinity and drought conditions, whereas another alternative splice variant, OsGATA23b, did not respond to the above conditions [33]. Therefore, investigating the AS regulation and expression patterns of different AS isoforms under stress conditions may shed some light on the mechanism of abiotic stress response in soybean. Understanding the impact of AS on the plasticity of root development under abiotic stress conditions in soybean could also offer promising guidance to soybean breeders to develop drought-resistant cultivars [34].
In the present study, soybean (Williams 82 variety) roots under very mild, mild, and severe drought stress as well as water recovery after severe drought stress were used as materials in an RNA-sequencing (RNA-seq) experiment to identify the alternatively spliced isoforms and quantify the differential drought-responsive AS events. All five AS types were found in soybean roots in response to drought stress. The differentially spliced genes were determined and used to identify the enriched Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Our results indicate that the skipped exon (SE) type was the most abundant AS event type. We also identified some novel AS events that might function in the drought response. These findings provide insight into AS regulation during soybean response to drought in the post-transcriptional process.

2. Materials and Methods

2.1. Datasets and Treatments

The RNA-seq dataset (http://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/sra/SRP067593, [35]) was used to test AS with replicate multivariate analysis of transcript splicing (rMATs) program in this study. Soybean variety Williams 82 (G. max) was used in this study. Four treatments—very mild stress (VMS), mild stress (MS), severe stress (SS), and water recovery after severe stress (SR)—were applied and monitored on the basis of the leaf water potential [35]. All plants were well-watered until V3 (three unfolded trifoliate leaves) stage before treatment. Each stress treatment corresponded to a well-watered control at the same growth stage. The VMS treatment involved 5 days of withholding water while other plants were well-watered as a corresponding control (VMS-WW). For MS treatment, water was withheld for 12 days while well-watered plants were used as corresponding control (MS-WW). For SS treatment, water was withheld for 19 days while other well-watered plants were used as a corresponding control (SS-WW). For SR treatment, plants were re-watered for 2 days after SS treatment while well-watered plants were used as a corresponding control (SR-WW). Therefore, the transcriptome datasets from a total of 24 libraries [(4 treatments + 4 controls) × 3 reps] were used in this study.

2.2. Reference-Based Transcriptome Structure Assembly

RNA isolation, library construction, and RNA sequencing data were consistent with the previous description in [35]. More mapping processes of RNA-seq reads were performed in this study to match the rMATs analysis. Except for fastx_toolkit_0.0.14 (http://hannonlab.cshl.edu/fastx_toolkit/), additional filtering for poor-quality bases was performed through SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle). Furthermore, in order to identify novel splice sites with direct mapping to known transcripts and produce sensitive and accurate alignments, we used TopHat2 (http://0-ccb-jhu-edu.brum.beds.ac.uk/software/tophat/index.shtml, [36]) instead of TopHat for alignment with the G. max reference genome (G. max Wm82.a2.v1) and Phytozome v13 gene model in this study. The software StringTie (http://0-ccb-jhu-edu.brum.beds.ac.uk/software/stringtie/) was used to assemble the mapped reads, which were used to compare with the known transcripts in order to obtain novel transcripts. Novel transcripts were named as MSTRG plus numbers.

2.3. AS Detection and Identification of Drought Stress-Responsive AS Events

The software rMATs (version 4.0.1) [24] (http://rnaseq-mats.sourceforge.net/index.html) was used to identify AS events and analyze the differential AS events between samples. The junction count-only method was used in this study. We identified AS events with a false discovery rate (FDR) < 0.05 in comparison to AS events with a significant drought response. The classification of AS was as follows: SE: skipped exon, MXE: mutually exclusive exon, A5SS: alternative 5′ splice site, A3SS: alternative 3′ splice site, RI: retained intron.

2.4. Differential Gene Expression Analysis

The differential expression analysis of RNA-seq data was performed between 2 different groups by using DESeq2 software (and by edgeR between 2 samples) [37,38]. The genes/transcripts with a false discovery rate (FDR) of <0.05 and absolute fold change ≥ 2 were considered as differentially expressed genes/transcripts. The expression values of the resulting genes or transcripts were presented as transcripts per million (TPM) values.

2.5. Gene Ontology and Pathway Enrichment Analysis

All DSGs were mapped to GO terms in the Gene Ontology database (http://www.geneontology.org/) and gene numbers were calculated for every term. Significantly enriched GO terms in DSGs compared to the genome background were defined by hypergeometric test. The calculation of the p-value was performed with FDR correction, taking an FDR ≤ 0.05 as a threshold. The significantly overrepresented KEGG pathways were determined by Fisher’s exact test (p-value < 0.05) compared with the entire genome background.

3. Results

3.1. Overview of the Sequencing Quality of RNA-seq Data for Identifying AS Events

Twenty-four RNA-seq libraries were used to examine the regulation of AS in soybean roots under various drought stress conditions. These libraries included very mild drought stress (VMS), mild drought stress (MS), severe drought stress (SS), and recovery after severe drought (SR), with each stress treatment corresponding to a well-watered control at the same growth stage (WW). Each sample had three biological replicates and was sequenced separately [35]. Our previous study showed that thousands of genes were differentially expressed in soybean roots in response to various levels of drought stress [35]. Here, we used the same RNA-seq dataset to further examine the comprehensive profiles of AS events under those stress conditions. Instead of TopHat, TopHat2 was used in this study to identify novel splice sites with direct mapping to known transcripts, since the software can produce sensitive and accurate alignments, even for highly repetitive genomes or for those with pseudogenes [36]. In total, over 1.54 billion clean reads were utilized in downstream analysis for each sample, and the Q20 (base ratio > 20), Q30 (base ratio > 30), and GC content of the obtained clean data were calculated. The high percentage of Q30 (over 86.2%) in each sample suggested that a high accuracy of sequencing was achieved (Table S1). As shown in Figure S1A, the total filtered reads were uniformly distributed across all 20 soybean chromosomes. In addition, the coverage of sequencing results indicated that there was no significant bias in sequencing (Figure S1B). Therefore, the quality of assembly and sequencing data was high, which is consistent with previous reports [35], and it could be used for the detection and analysis of AS events.

3.2. Identification of AS Events in Soybean Roots under Various Drought Treatments

High-quality reads were first mapped to the soybean reference genome (Williams 82, Ensembl version Glycine_max_v2.1) and then AS events were identified and quantified by the junction counts only method. All five AS types were identified in soybean roots with or without stress treatments. There was no significant difference in the number of these five AS types between well-watered control of very mild drought stress (VMS_WW) and VMS conditions or between well-watered control of mild drought stress (MS_WW) and MS conditions. However, the numbers of AS events increased significantly under the SS condition and decreased significantly under the SR condition (Figure 1). The differences in terms of AS numbers between replicates were smaller than that of treatments in SS- and SR-treated samples. Specifically, 10,312 SE-type AS events and 5448 A3SS-type AS events were detected under the SS condition, which was significantly greater than that of AS events under the SS_WW condition (p < 0.01). On the contrary, the number of AS events under the SR condition significantly decreased compared with that of the SR_WW and SS conditions. For example, only 7862 SE-type AS events were identified under the SR condition (Figure 1). In addition, the distributions of AS types under each condition were comparable—SE was the most abundant (41–45%) of the AS events out of all five AS types, followed by A3SS (24–26%); the number of A5SS and retained intron (RI) (14–16%) events were similar, but the number of MXE events (2%) was far smaller than that of the other AS types (Figure S2). These results suggest that the number of AS events increased significantly under the severe drought condition and SE was the major splicing type in soybean roots under various drought conditions. The number of splicing events increased with the deepening of the degree of drought while it decreased significantly after re-watering.

3.3. Identification of Drought Stress-Responsive AS Events in Soybean Root

The splicing isoforms of the genes that showed differential expression after stress treatments were called stress-responsive (drought) AS (DAS) events. This kind of event may be involved in regulating the plant response to drought stress. In this study, the “junction reads-based FDR (false discovery rate calculated from p-value) < 0.05” criteria were used to determine the DAS events. In total, 385, 989, 1429, and 465 stress-responsive DAS events were identified under VMS, MS, SS, and SR conditions, respectively. This result indicates that AS responses under the MS and SS treatments were enhanced compared with the VMS and SR treatments (Figure 2). Among all types of AS events, A3SS was the most abundant AS type under the VMS (45.2%), MS (33.0%), and SR (41.5%) conditions, but the number of SE-type (especially exclusion-type) events significantly increased under the MS and SS conditions (39.8%) compared to all the other AS types. In addition, the numbers of RI events under the MS and SS treatments were more than those of the VMS and SR treatments. Moreover, the number of alternative exon exclusions was more than the number of exon inclusions in A3SS and A5SS types under MS and SS conditions but was less than those of the VMS and SR conditions (Figure 2). These results suggest that the distribution of AS types and ratios of different isoforms under various drought conditions were different.
The IEP (isoform expression percentage) changes have been used to classify the groups of stress-responsive AS events in wheat [15]. In the present study, the IEPs for all stress-responsive AS events were further clustered to reveal the magnitude of AS pattern changes (Figure 3). Seven groups were identified on the basis of AS drought-response profiles. Group 1, representing a total of 240 AS events, indicated responses to VMS treatment. Group 2 showed AS pattern changes under both VMS and MS treatments. Group 3 showed slight AS pattern changes under MS, SS, and SR treatments. There were 569 AS events in group 4, which included response to MS and SS treatments. Group 5 only showed significant changes of AS pattern under SR treatment, suggesting that SR treatment created a specific AS response. Groups 6 and 7 only exhibited significant changes of AS pattern in response to the VMS or SS treatments, respectively. Therefore, our present study revealed seven groups of AS events that responded differentially to different drought conditions.

3.4. Identification of Potential Drought Regulators via Analyzying Novel Drought-Responsive DSGs and DEGs

The genes with significantly differential drought stress-responsive AS events were named differentially spliced genes (DSGs). In this study, a total of 2120 DSGs (345, 801, 1130, and 398 under VMS, MS, SS, and SR treatments, respectively) were identified in soybean roots (Table S2). First, the DSGs and differentially expressed genes (DEGs) were compared and the results indicated that 183 genes showed significant drought-responsive expression patterns and AS regulation patterns (Figure 4A). Second, the DSGs were compared under various drought conditions, and the results showed that MS and SS shared the largest overlap—12% (255) of the DSGs overlapped between MS and SS treatments, while only 2.2–5.3% (46-111) of the DSGs overlapped between VMS/MS, SS/SR, VMS/SS, VMS/SR, or MS/SR (Figure 4B). Similarly, more stress-responsive AS gene overlap was found between MS and SS conditions (Figure S3). These results suggested that the mild drought-induced AS response was similar to those induced by severe drought stress, but distinct from VMS and SR. Notably, many genes encoding hormone-related proteins, drought-related proteins, splicing-related proteins, and heat shock-related proteins exhibited significant AS pattern changes (Figure S4). For example, drought-induced protein 19 exhibited significant AS pattern changes in response to the VMS, SS, and SR treatments (Figure S4A), and auxin efflux carrier family protein (Glyma.19G128800) underwent different AS events under different drought stresses (Figure S4D,E). These results indicate that many pathways enriched in DEGs were also found in DSGs.
To further investigate whether any novel AS events were identified and involved in soybean drought response in this study, we studied three types of AS (MXE, RI, and SE) events that significantly responded to drought stress (Table S3). First, more novel AS events were identified under the MS and SS conditions than under the VMS and SR conditions. Second, a total of 232 novel SE–type AS events were identified under all treatments (Table S3), suggesting the SE-type AS may be more involved in drought response of soybean roots. These results further indicate that AS may play an important role in regulating the drought response in soybean roots.

3.5. Comparative Analysis of the Biological Functions of All DSGs Regulated at AS and Transcription Levels

To investigate the biological functions that changed in response to drought and AS modulation, we performed Gene Ontology (GO) enrichment analysis was performed for all DSGs. As shown in Figure 5A, many genes involved in metabolic and cellular biological processes, as well as in the binding and catalytic activity involved in molecular function, were subject to AS regulation. Of particular interest were the many DSGs enriched in the following GO terms: nucleic acid binding, organic cyclic compound binding, and heterocyclic compound binding, suggesting that the AS regulation of binding activity had vital roles in the drought response of soybean root (Figure 5B). In addition, we found that the biological process related to the substance metabolic process and the primary metabolic process were significantly overrepresented (Figure S5). Those enriched biological processes, metabolic pathways, and biochemical activities provided an overview of the possible regulation of drought response in soybean roots through AS.

3.6. Comparative Analysis of the Selected Differentially Expressed and Alternatively Spliced Genes

Splicing regulatory factors (SPFs), such as Ser/Arg-rich proteins, pre-mRNA splicing factors, and heterogeneous nuclear ribonucleoproteins (hnRNPs), are not only responsive to environmental stresses but are also alternatively spliced themselves under various environmental stresses [12,26,39,40]. KEGG pathway analysis was applied to evaluate whether such genes and pathways were enriched in the DSGs. We found that the spliceosome pathway (p value < 4.12 × 10−10) and mRNA surveillance pathway (p value < 3.28 × 10−5) were enriched in the DSGs (Figure 6A). These results indicate that the SPF-related genes in soybean also responded to drought stress and were alternatively spliced under drought conditions.
In the present study, a total of 56 SPFs-related genes were also found to be significantly differentially spliced under various drought treatments. Among them, A3SS and SE types comprised the majority of AS types (Figure 6B). In addition, the number of upregulated SPF-related genes was greater than that of the downregulated genes (Figure 6C). The response of 49 SPF-related genes to various drought treatments was significantly different. Most SPF-related genes were upregulated in response to the MS and SS treatments compared with the VMS and SR treatments (Figure S6). These results indicate that the genes encoding the splicing factors were severely spliced in response to the drought stress of soybean roots.
To further reveal how the different stress-responsive transcripts responded to the various drought stresses, we investigated the expression patterns of four SPFs genes. As indicated in Figure 7A, the MS and SS treatments slightly induced the expression of Glyma.13G241800 (encoding the splicing factor 3b subunit), but greatly induced the expression level of a novel transcript (named MSTRG.29150.1). Conversely, the expression pattern of the novel transcript (named MSTRG.13253.3) of Glyma.06G324500 gene was downregulated under the MS condition but upregulated under the SS treatment (Figure 7B). In addition, the expression level of MSTRG.44538.1 was significantly higher than that of Glyma.20G236800 (Figure 7C). On the contrary, the expression level of MSTRG.19304.3 was significantly lower than that of Glyma.09G104200 (Figure 7D). These results indicate that different transcripts responded differently to drought and therefore may have different roles in regulating soybean response to drought.

4. Discussion

4.1. Many AS Events Occurred in Soybean Roots under Various Drought Conditions

Abiotic stress-regulated AS events have not been systematically analyzed and reported at the whole transcriptome level in soybean, although AS has been studied in other plants and implicated in regulating plant stress responses [9,12,13,14,15]. In the present study, 385, 989, 1429, and 465 drought stress-responsive AS events were identified in soybean roots under VMS, MS, SS, and SR conditions, respectively, using whole transcriptome RNA-seq data. More AS events occurred under the MS and SS treatments than under the VMS and SR treatments, indicating that the number of AS events increased under the severe drought condition and decreased after re-watering. Therefore, our findings not only revealed a large number of AS events in soybean roots in response to various drought treatments, but also support the important roles of AS events in regulating soybean response to drought stress.

4.2. A3SS and SE Were the Main AS Types in Soybean Roots in Response to Drought Stress

Although all five AS types were identified in the present study, SE and A3SS were the most abundant AS types and represented 67.1–69.0% of the total AS events under all the drought stress conditions (Figure S2). Further investigation showed that those two AS types accounted for 63.7–68.8% of all the drought stress-responsive AS events (Figure 2). Additionally, an approximately equal number of the observed AS events was found between A5SS and RI types under all conditions, but the number of A5SS type (85-294) was more than that of RI type (23-127) in all the stress-responsive AS events (Figure 2). On the contrary, RI was the most abundant AS event in different developmental stages of soybean [41,42,43]. RI was the major AS type that responded to drought, heat, and their combination in wheat and to nitrogen in maize [15,44]. The difference in dominant AS types between different developmental stages in soybean and between different plant species may be related to the regulation of SPF-related genes (both gene expression level and splicing events). Further work is needed to understand the observed difference.

4.3. A Number of Key GO Terms and Pathways Were Enriched in DSGs

DEGs were reported in a previous study and the significantly enriched KEGG pathways were related to the biosynthesis of secondary metabolites and phenylpropanoid biosynthesis, as well as metabolic pathways [35]. In the present study, the most enriched pathways in the DSGs were related to splicing and mRNA surveillance (Figure 6). In total, 85 (4%) DSGs were related to the spliceosome and mRNA surveillance pathways. The spliceosome processes pre-RNA to mature mRNA and is extensively involved in the splicing process [20,21]. Recently, it was reported that the modification of key spliceosome components is consistent with the transcriptional and proteome changes under drought stress in Arabidopsis [45]. In the present study, RNA-binding proteins (Glyma.03G214400, Glyma.05G003400) were also identified to be regulated at the splicing level, although a detailed composition of the spliceosome is yet to be characterized in soybean under drought stress. Therefore, it appears that specific spliceosome changes also occur in soybean roots under drought stress.
The mRNA surveillance pathway is involved in the quality control of the accuracy of gene expression, as it detects and degrades abnormal mRNAs. Nonsense-mediated decay (NMD) is a major pathway of mRNA surveillance, which eliminates mRNAs containing premature translation termination codons (PTCs) [46]. It was reported that at least 17.4% of all multi-exon, protein-coding genes produce splicing variants that are targeted by NMD in Arabidopsis. UP FRAMESHIFT1 (UPF3), an NMD factor, could result in the de-capping and rapid exonucleolytic digestion of the mRNA [47]. UPF3 expression is feedback-regulated at multiple levels, and its balanced expression is essential for coping with salt stress in Arabidopsis [48]. Furthermore, the AS-coupled NMD was modulated by NaCl stress in Arabidopsis [47]. In our present study, one gene (Glyma.06G324500) encoding UPF3 regulator of nonsense transcripts homolog was found to be regulated by drought stress and AS (Figure 7D). Therefore, the mRNA surveillance pathway may also be involved in regulating soybean response to drought stress.
Plant splicing regulatory factors are known to be essential in regulating the AS of drought-related genes [12,26,39,40]. The interactions between Ser/Arg-rich protein and pre-mRNA are essential to constitutive and alternative splicing and have crucial contributions to maintaining cell and tissue homeostasis [12,49,50]. Interestingly, plant SPF-related genes themselves also undergo AS in a developmental and tissue-specific manner as well as in response to various hormones and abiotic stresses [10,12,51]. Overexpression of BrSR45a in Arabidopsis not only increased the abundance of the drought stress-inducible genes but also influenced the splicing patterns of target genes [52]. The transcript abundance and the ratio of exon-skipped transcripts of AtSR45a were promoted by heat and drought stress in Arabidopsis [53]. Notably, a number of genes encoding the Ser/Arg-rich proteins involved in pre-mRNA splicing were found to generate AS events in response to drought in the present study (Figure 6), for example, Glyma.16G121300, a homologue of SR45a. Likely the stress-dependent AS of SPF-related genes may also function in regulating soybean response to drought stress.
The most enriched GO molecular functional terms in DSGs were related to nucleic acid binding, organic cyclic compound binding, and heterocyclic compound binding (Figure 5B). It was reported that nucleic acid binding and organic cyclic compound binding were significantly enriched in the OsMYBs regulatory network in drought response in rice [54]. The genes related to “organic cyclic compound binding and heterocyclic compound binding” were also over-represented under water deficit conditions in Chilean quinoa [55]. Therefore, in addition to drought-induced changes in gene expression levels, the drought-induced AS events may also play important roles in soybean root development. The enriched pathways may be conserved in different plant species in response to drought stress.

4.4. The Complexity of the Steady State Transcriptome in Soybean Roots under Drought Conditions

Our present and previous studies have demonstrated that the soybean transcriptome was not only regulated at the expression level, but also through the AS regulation. Previously, we showed that 6609 genes were differentially responsive to various drought conditions in soybean root [35]. Here, we revealed not only a large number of soybean genes that produced splicing variants that responded to drought conditions at transcript levels, but also a large number of genes that produced differentially spliced variants in response to different drought levels. Furthermore, 183 genes appeared to experience both changes in expression level and AS regulation in response to drought stresses. These results indicate the complexity of the transcriptome in soybean roots under drought conditions.

4.5. AS Has an Important Regulatory Role during Plant Development and Response to Abiotic Stress

AS has been studied extensively in several plant species. The concept of deep learning and gene editing technology (Clusters of Regularly Interspaced Short Palindromic Repeats, CRISPR) have been used to identify AS events and gene function [52,56,57]. Moreover, differences in the recognition sequences of splice sites have been reported between plants and animals [26,58,59]. This might occur because plants have to endure unfavorable environments due to a sessile lifestyle. In particular, the conserved AS regulation model of key abiotic stress regulators among different plant species have been found, such as heat stress transcription factors (HsfA2) and dehydration-responsive element binding protein 2B (DREB2B) [17,18,60]. In addition, environmental stress regulates the changes of those AS regulators, and further regulates the AS and expression of other stress-related genes in plants, which ensures that plants can quickly respond and adapt to environmental changes [61,62,63,64,65].
With the development of high-throughput sequencing technology, more and more plant genome sequences and transcriptome data can provide technical support for the identification of species, tissues, developmental stages, and environmental specific regulation of AS, which help researchers to determine the dynamic changes of AS in different development stages and environmental conditions. Moreover, the research on the function and mechanism of AS in the process of abiotic stress can provide new strategies for enhancing stress tolerance of plants.

5. Conclusions

In summary, the present study revealed a large number of AS events that occurred in soybean roots in response to drought stress through analyzing the whole transcriptome RNA-seq data, providing a comprehensive view of AS in soybean roots under drought stress conditions. Therefore, in addition to the regular differential gene expression, our data suggest that drought-responsive AS events and/or their expression levels both act in a direct or indirect regulation mode to modulate soybean drought responses. Systemic approaches are further needed to dissect and understand soybean responses to drought stress and the roles of DSGs as well as DEGs in regulating such responses. The obtained knowledge will help develop future drought-resistant soybean cultivars.

Supplementary Materials

The following are available online at https://0-www-mdpi-com.brum.beds.ac.uk/2073-4425/11/12/1520/s1: Figure S1: Quality and feature of the RNA-seq datasets used in this study. (A) Chromosome-wise read count distribution of all filtered reads in root samples. Reads were obtained throughout the soybean genome and distributed on different chromosomes using Williams 82 sequence as a reference. Raw reads were filtered using the criteria of MAPQ > 20 and insert size 100–1000 bp. (B) Normalized gene-body coverage of RNA-seq libraries generated from soybean root. For each library, the average coverage is shown at each relative position along the transcripts’ length. Figure S2: The number of different AS events in genes expressed in seedlings based on RNA-seq data compared with the annotated gene models under all drought conditions. Blue, exclusion; red, inclusion. Figure S3: Venn diagrams of the differentially spliced AS events derived from five AS types under various drought conditions: (A) SE type, (B) A3′SS type, (C) A5′SS type, (D) RI type, and (E) MXE type. Figure S4: AS profiles of significantly stress-responsive genes under various drought conditions. The bar charts show relative expression levels of alternatively spliced isoform that was well-watered (ww) (blue) and isoform stress (yellow) of AS genes under various drought conditions revealed by RNA-seq. The numbers in each bar indicate the expression percentages of different isoforms. Figure S5: The top 20 GO enrichment terms (biological process) of DEGs under all conditions. Notes: the blue bar chart represents the gene percent (%) of the GO terms. Figure S6: Clustering analysis of differentially expressed SPF-related genes under various drought treatments. The transcripts per million (TPM) values of expression levels of those genes were used. Upregulated (red) and downregulated (blue) SPF-related genes are presented by different color scales in the heatmap. Table S1: Raw data and clean data statistics of RNA sequencing. Table S2: The ID list of DSGs. Table S3: DSGs encoding novel transcripts under various drought conditions.

Author Contributions

Conceptualization, L.S. and H.C.; methodology, L.S.; software, Y.D.; validation, Z.P., L.C., and L.S.; formal analysis, L.S.; investigation, L.S.; resources, L.S. and H.T.N.; data curation, H.T.N.; writing—original draft preparation, L.S. and H.C.; writing—review and editing, J.W., H.Y., H.T.N., and G.Z.; visualization, Z.P. and L.C.; supervision, L.S.; project administration, L.S. and H.C.; funding acquisition, L.S. and H.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Key R&D project of Jiangsu Province (BE2019376); the National Natural Science Foundation of China (31871540); the Natural Science Foundation of Jiangsu Province (BK20191438); Jiangsu Agriculture Science and Technology Innovation Fund (CX(20)2007); the Open Project Program of Joint International Research Laboratory of Agriculture and Agri-Product Safety; the Ministry of Education of China, Yangzhou University(JILAR-KF202005); and Graduate Practice Innovation Program of Yangzhou University (no. XSJCX19_096).

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviation

ASalternative splicing
A5SSalternative 5′ splicing site
A3SSalternative 3′ splicing site
DASstress-responsive AS events
DREB2Bdehydration-responsive element binding protein 2B
DSGsdrought spliced genes
PTCpremature translation termination codon
FDRfalse discovery rate calculated from p-value
GOGene Ontology
hnRNPsheterogeneous nuclear ribonucleoproteins
IEPsisoform expression percentage changes
IRintron retention
KEGGKyoto Encyclopedia of Genes and Genomes
MXEmutually exclusive exon
MSmild drought stress
NMDnonsense-mediated decay
rMATsreplicate multivariate analysis of transcript splicing
SEskipped exon
SSsevere drought stress
SRrecovery after severe drought
SPFssplicing regulatory factors
UPF3UP FRAMESHIFT1
VMSvery mild drought stress
WWwell-watered control

References

  1. Dugas, D.V.; Monaco, M.; Olson, A.; Klein, R.R.; Kumari, S.; Ware, D.; Klein, P.E. Functional annotation of the transcriptome of Sorghum bicolor in response to osmotic stress and abscisic acid. BMC Genom. 2011, 12, 514. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Gracheva, E.O.; Cordero-Morales, J.F.; González-Carcacía, J.A.; Ingolia, N.T.; Manno, C.; Aranguren, C.I.; Weissman, J.S.; Julius, D. Ganglion-specific splicing of TRPV1 underlies infrared sensation in vampire bats. Nature 2011, 476, 88–91. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Kriechbaumer, V.; Wang, P.; Hawes, C.; AbellRobert, B.M. Alternative splicing of the auxin biosynthesis gene YUCCA4 determines its subcellular compartmentation. Plant J. 2012, 70, 292–302. [Google Scholar] [CrossRef] [PubMed]
  4. Pikaard, C.S.; Scheid, O.M. Epigenetic Regulation in Plants. Cold Spring Harb. Perspect. Biol. 2014, 6, a019315. [Google Scholar] [CrossRef] [PubMed]
  5. Xing, Y.; Zhang, Q. Genetic and Molecular Bases of Rice Yield. Annu. Rev. Plant Biol. 2010, 61, 421–442. [Google Scholar] [CrossRef]
  6. Huang, J.; Gao, Y.; Jia, H.; Liu, L.; Zhang, D.; Zhang, Z. Comparative transcriptomics uncovers alternative splicing changes and signatures of selection from maize improvement. BMC Genom. 2015, 16, 363. [Google Scholar] [CrossRef] [Green Version]
  7. Posé, D.; Verhage, L.; Ott, F.; Yant, L.; Mathieu, J.; Angenent, G.C.; Immink, R.G.H.; Schmid, M. Temperature-dependent regulation of flowering by antagonistic FLM variants. Nature 2013, 503, 414–417. [Google Scholar] [CrossRef]
  8. Feng, J.; Li, J.; Gao, Z.; Lu, Y.; Yu, J.; Zheng, Q.; Yan, S.; Zhang, W.; He, H.; Ma, L.; et al. SKIP Confers Osmotic Tolerance during Salt Stress by Controlling Alternative Gene Splicing in Arabidopsis. Mol. Plant 2015, 8, 1038–1052. [Google Scholar] [CrossRef] [Green Version]
  9. Thatcher, S.R.; Danilevskaya, O.N.; Meng, X.; Beatty, M.; Zastrow-Hayes, G.; Harris, C.; Van Allen, B.; E Habben, J.; Li, B. Genome-Wide Analysis of Alternative Splicing during Development and Drought Stress in Maize. Plant Physiol. 2016, 170, 586–599. [Google Scholar] [CrossRef] [Green Version]
  10. Zhu, J.; Wang, X.; Guo, L.; Xu, Q.; Zhao, S.; Li, F.; Yan, X.; Liu, S.; Wei, C. Characterization and Alternative Splicing Profiles of the Lipoxygenase Gene Family in Tea Plant (Camellia sinensis). Plant Cell Physiol. 2018, 59, 1765–1781. [Google Scholar] [CrossRef] [Green Version]
  11. Zhu, F.-Y.; Chen, M.-X.; Ye, N.-H.; Shi, L.; Ma, K.-L.; Yang, J.-F.; Cao, Y.-Y.; Zhang, Y.; Yoshida, T.; Fernie, A.R.; et al. Proteogenomic analysis reveals alternative splicing and translation as part of the abscisic acid response in Arabidopsis seedlings. Plant J. 2017, 91, 518–533. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  12. Ding, F.; Cui, P.; Wang, Z.; Zhang, S.; Ali, S.; Xiong, L. Genome-wide analysis of alternative splicing of pre-mRNA under salt stress in Arabidopsis. BMC Genom. 2014, 15, 431. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  13. Calixto, C.P.; Guo, W.; James, A.B.; Tzioutziou, N.A.; Entinze, J.C.; Panter, P.E.; Knight, H.; Nimmo, H.G.; Zhang, R.; Brown, J.W. Rapid and Dynamic Alternative Splicing Impacts the Arabidopsis Cold Response Transcriptome. Plant Cell 2018, 30, 1424–1444. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Li, S.; Yu, X.; Cheng, Z.; Zeng, C.; Li, W.; Zhang, L.; Peng, M. Large-scale analysis of the cassava transcriptome reveals the impact of cold stress on alternative splicing. J. Exp. Bot. 2020, 71, 422–434. [Google Scholar] [CrossRef]
  15. Liu, Z.; Qin, J.; Tian, X.; Xu, S.; Wang, Y.; Li, H.; Wang, X.; Peng, H.; Yao, Y.; Hu, Z.; et al. Global profiling of alternative splicing landscape responsive to drought, heat and their combination in wheat (Triticum aestivum L.). Plant Biotechnol. J. 2017, 16, 714–726. [Google Scholar] [CrossRef] [Green Version]
  16. Matsukura, S.; Mizoi, J.; Yoshida, T.; Todaka, D.; Ito, Y.; Maruyama, K.; Shinozaki, K.; Yamaguchi-Shinozaki, K. Comprehensive analysis of rice DREB2-type genes that encode transcription factors involved in the expression of abiotic stress-responsive genes. Mol. Genet. Genom. 2010, 283, 185–196. [Google Scholar] [CrossRef]
  17. Egawa, C.; Kobayashi, F.; Ishibashi, M.; Nakamura, T.; Nakamura, C.; Takumi, S. Differential regulation of transcript accumulation and alternative splicing of a DREB2 homolog under abiotic stress conditions in common wheat. Genes Genet. Syst. 2006, 81, 77–91. [Google Scholar] [CrossRef] [Green Version]
  18. Terashima, A.; Takumi, S. Allopolyploidization reduces alternative splicing efficiency for transcripts of the wheat DREB2 homolog, WDREB2. Genome 2009, 52, 100–105. [Google Scholar] [CrossRef]
  19. Wei, H.; Lou, Q.; Xu, K.; Yan, M.; Xiaosong, M.; Ma, X.; Yu, X.; Luo, L. Alternative splicing complexity contributes to genetic improvement of drought resistance in the rice maintainer HuHan2B. Sci. Rep. 2017, 7, 1–13. [Google Scholar] [CrossRef]
  20. Syed, N.H.; Kalyna, M.; Marquez, Y.; Barta, A.; Brown, J.W.S. Alternative splicing in plants—Coming of age. Trends Plant Sci. 2012, 17, 616–623. [Google Scholar] [CrossRef] [Green Version]
  21. Reddy, A.S.; Marquez, Y.; Kalyna, M.; Barta, A. Complexity of the Alternative Splicing Landscape in Plants. Plant Cell 2013, 25, 3657–3683. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  22. 2019 Soystats. Available online: http://soystats.com/wp-content/uploads/2019-SoyStats-Web.pdf (accessed on 10 August 2019).
  23. Aghamirzaie, D.; Nabiyouni, M.; Fang, Y.; Klumas, C.; Heath, L.S.; Grene, R.; Collakova, E. Changes in RNA Splicing in Developing Soybean (Glycine max) Embryos. Biology 2013, 2, 1311–1337. [Google Scholar] [CrossRef] [PubMed]
  24. Shen, Y.; Zhou, Z.; Wang, Z.; Li, W.; Fang, C.; Wu, M.; Ma, Y.; Liu, T.; Kong, L.-A.; Peng, D.-L.; et al. Global Dissection of Alternative Splicing in Paleopolyploid Soybean. Plant Cell 2014, 26, 996–1008. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. Modrek, B.; Lee, C. A genomic view of alternative splicing. Nat. Genet. 2002, 30, 13–19. [Google Scholar] [CrossRef]
  26. Filichkin, S.A.; Priest, H.D.; Givan, S.A.; Shen, R.; Bryant, D.W.; Fox, S.E.; Wong, W.-K.; Mockler, T.C. Genome-wide mapping of alternative splicing in Arabidopsis thaliana. Genome Res. 2010, 20, 45–58. [Google Scholar] [CrossRef] [Green Version]
  27. Marquez, Y.; Brown, J.W.S.; Simpson, C.; Barta, A.; Kalyna, M. Transcriptome survey reveals increased complexity of the alternative splicing landscape inArabidopsis. Genome Res. 2012, 22, 1184–1195. [Google Scholar] [CrossRef] [Green Version]
  28. Walters, B.; Lum, G.; Sablok, G.; Min, X.J. Genome-Wide Landscape of Alternative Splicing Events in Brachypodium distachyon. DNA Res. 2013, 20, 163–171. [Google Scholar] [CrossRef]
  29. Farooq, M.; Gogoi, N.; Barthakur, S.; Baroowa, B.; Bharadwaj, N.; Alghamdi, S.S.; Siddique, K.H.M. Drought Stress in Grain Legumes during Reproduction and Grain Filling. J. Agron. Crop Sci. 2017, 203, 81–102. [Google Scholar] [CrossRef]
  30. Dastidar, M.G.; Jouannet, V.; Maizel, A. Root branching: Mechanisms, robustness, and plasticity. Wiley Interdiscip. Rev. Dev. Biol. 2012, 1, 329–343. [Google Scholar] [CrossRef]
  31. Wei, H.; Lou, Q.; Xu, K.; Zhou, L.; Chen, S.; Chen, L.; Luo, L. Pattern of alternative splicing different associated with difference in rooting depth in rice. Plant Soil 2020, 449, 233–248. [Google Scholar] [CrossRef]
  32. Ishizawa, M.; Hashimoto, K.; Ohtani, M.; Sano, R.; Kurihara, Y.; Kusano, H.; Demura, T.; Matsui, M.; Sato-Nara, K. Inhibition of Pre-mRNA Splicing Promotes Root Hair Development in Arabidopsis thaliana. Plant Cell Physiol. 2019, 60, 1974–1985. [Google Scholar] [CrossRef] [PubMed]
  33. Gupta, P.; Nutan, K.K.; Singla-Pareek, S.L.; Pareek, A. Abiotic Stresses Cause Differential Regulation of Alternative Splice Forms of GATA Transcription Factor in Rice. Front. Plant Sci. 2017, 8, 1944. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Stauffer, E.; Maizel, A. Post-transcriptional regulation in root development. Wiley Interdiscip. Rev. RNA 2014, 5, 679–696. [Google Scholar] [CrossRef] [PubMed]
  35. Song, L.; Prince, S.; Valliyodan, B.; Joshi, T.; Dos Santos, J.V.M.; Wang, J.; Lin, L.; Wan, J.; Wang, Y.; Xu, D.; et al. Genome-wide transcriptome analysis of soybean primary root under varying water-deficit conditions. BMC Genom. 2016, 17, 1–17. [Google Scholar] [CrossRef] [Green Version]
  36. Kim, D.; Pertea, G.; Trapnell, C.; Pimentel, H.; Kelley, R.; Salzberg, S.L. TopHat2: Accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013, 14, R36. [Google Scholar] [CrossRef] [Green Version]
  37. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef] [Green Version]
  38. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2009, 26, 139–140. [Google Scholar] [CrossRef] [Green Version]
  39. Palusa, S.G.; Ali, G.S.; Reddy, A.S.N. Alternative splicing of pre-mRNAs of Arabidopsis serine/arginine-rich proteins: Regulation by hormones and stresses. Plant J. 2007, 49, 1091–1107. [Google Scholar] [CrossRef]
  40. Tanabe, N.; Yoshimura, K.; Kimura, A.; Yabuta, Y.; Shigeoka, S. Differential expression of alternatively spliced mRNAs of Arabidopsis SR protein homologs, atSR30 and atSR45a, in response to environmental stress. Plant Cell Physiol. 2007, 48, 1036–1049. [Google Scholar]
  41. Shen, S.; Park, J.W.; Lu, Z.-X.; Lin, L.; Henry, M.D.; Wu, Y.N.; Zhou, Q.; Xing, Y. rMATS: Robust and flexible detection of differential alternative splicing from replicate RNA-Seq data. Proc. Natl. Acad. Sci. USA 2014, 111, E5593–E5601. [Google Scholar] [CrossRef] [Green Version]
  42. Wang, L.; Cao, C.; Ma, Q.; Zeng, Q.; Ewang, H.; Cheng, Z.; Zhu, G.; Qi, J.; Zhang, X.; Nian, H.; et al. RNA-seq analyses of multiple meristems of soybean: Novel and alternative transcripts, evolutionary and functional implications. BMC Plant Biol. 2014, 14, 169. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  43. Íñiguez, L.P.; Ramírez, M.; Barbazuk, W.B.; Hernández, G. Identification and analysis of alternative splicing events in Phaseolus vulgaris and Glycine max. BMC Genom. 2017, 18, 1–17. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Wang, Y.; Xu, J.; Ge, M.; Ning, L.; Hu, M.; Zhao, H. High-resolution profile of transcriptomes reveals a role of alternative splicing for modulating response to nitrogen in maize. BMC Genom. 2020, 21, 1–19. [Google Scholar] [CrossRef] [PubMed]
  45. Marondedze, C.; Thomas, L.; Lilley, K.S.; Gehring, C. Drought Stress Causes Specific Changes to the Spliceosome and Stress Granule Components. Front. Mol. Biosci. 2020, 6, 163. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  46. Lareau, L.F.; Brooks, A.N.; Soergel, D.A.W.; Meng, Q.; Brenner, S.E. The Coupling of Alternative Splicing and Nonsense-Mediated mRNA Decay. Adv. Exp. Med. Biol. 2007, 623, 190–211. [Google Scholar] [CrossRef] [Green Version]
  47. Drechsel, G.; Kahles, A.; Kesarwani, A.K.; Stauffer, E.; Behr, J.; Drewe, P.; Ratsch, G.; Wachter, A. Nonsense-Mediated Decay of Alternative Precursor mRNA Splicing Variants Is a Major Determinant of the Arabidopsis Steady State Transcriptome. Plant Cell 2013, 25, 3726–3742. [Google Scholar] [CrossRef] [Green Version]
  48. Vexler, K.; Cymerman, M.A.; Berezin, I.; Fridman, A.; Golani, L.; Lasnoy, M.; Saul, H.; Shaul, O. The Arabidopsis NMD Factor UPF3 Is Feedback-Regulated at Multiple Levels and Plays a Role in Plant Response to Salt Stress. Front. Plant Sci. 2016, 7, 1376. [Google Scholar] [CrossRef] [Green Version]
  49. Reddy, A.S.; Ali, G.S. Plant serine/arginine-rich proteins: Roles in precursor messenger RNA splicing, plant development, and stress responses. Wiley Interdiscip. Rev. RNA 2011, 2, 875–889. [Google Scholar] [CrossRef]
  50. Laloum, T.; Martín, G.; Duque, P. Alternative Splicing Control of Abiotic Stress Responses. Trends Plant Sci. 2018, 23, 140–150. [Google Scholar] [CrossRef] [Green Version]
  51. Zhang, X.-N.; Mount, S.M. Two Alternatively Spliced Isoforms of the Arabidopsis SR45 Protein Have Distinct Roles during Normal Plant Development. Plant Physiol. 2009, 150, 1450–1458. [Google Scholar] [CrossRef] [Green Version]
  52. Muthusamy, M.; Yoon, E.K.; A Kim, J.; Jeong, M.-J.; Lee, S.I. Brassica Rapa SR45a Regulates Drought Tolerance via the Alternative Splicing of Target Genes. Genes 2020, 11, 182. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  53. Gulledge, A.A.; Roberts, A.D.; Vora, H.; Patel, K.; Loraine, A.E. Mining Arabidopsis thaliana RNA-seq data with Integrated Genome Browser reveals stress-induced alternative splicing of the putative splicing regulator SR45a. Am. J. Bot. 2012, 99, 219–231. [Google Scholar] [CrossRef] [PubMed]
  54. Esmita, S.; Ekatiyar, A.; Echinnusamy, V.; Epandey, D.M.; Bansal, K.C. Transcriptional Regulatory Network Analysis of MYB Transcription Factor Family Genes in Rice. Front. Plant Sci. 2015, 6, 1157. [Google Scholar] [CrossRef] [Green Version]
  55. Morales, A.; Zurita-Silva, A.; Maldonado, J.; Silva, H. Transcriptional Responses of Chilean Quinoa (Chenopodium quinoa Willd.) Under Water Deficit Conditions Uncovers ABA-Independent Expression Patterns. Front. Plant Sci. 2017, 8, 216. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  56. Louadi, Z.; Oubounyt, M.; Tayara, H.; Chong, K.T. Deep Splicing Code: Classifying Alternative Splicing Events Using Deep Learning. Genes 2019, 10, 587. [Google Scholar] [CrossRef] [Green Version]
  57. Butt, H.; Piatek, A.; Li, L.; Reddy, A.S.N.; Mahfouz, M. Multiplex CRISPR Mutagenesis of the Serine/Arginine-Rich (SR) Gene Family in Rice. Genes 2019, 10, 596. [Google Scholar] [CrossRef] [Green Version]
  58. Moore, M.J.; Proudfoot, N.J. Pre-mRNA processing reaches back to transcription and ahead to translation. Cell 2009, 136, 688–700. [Google Scholar] [CrossRef] [Green Version]
  59. Reddy, A.S. Plant serine/arginine-rich proteins and their role in pre-mRNA splicing. Trends Plant Sci. 2004, 9, 541–547. [Google Scholar] [CrossRef]
  60. Hu, Y.; Mesihovic, A.; Jiménez-Gómez, J.M.; Röth, S.; Gebhardt, P.; Bublak, D.; Bovy, A.; Scharf, K.D.; Schleiff, E.; Fragkostefanakis, S. Natural variation in HsfA2 pre-mRNA splicing is associated with changes in thermostolerance during tomato domestication. New Phytol. 2020, 225, 1297–1310. [Google Scholar] [CrossRef] [Green Version]
  61. Achard, P.; Chen, D.; Steele, A.D.; Lindquist, S.; Guarente, L. Integration of Plant Responses to Environmentally Activated Phytohormonal Signals. Science 2006, 311, 91–94. [Google Scholar] [CrossRef]
  62. Baena-Gonzalez, E.; Rolland, F.; Thevelein, J.M.; Sheen, J. A central integrator of transcription networks in plant stress and energy signalling. Nature 2007, 448, 938–942. [Google Scholar] [CrossRef] [PubMed]
  63. Barta, A.; Kalyna, M.; Reddy, A.S. Implementing a Rational and Consistent Nomenclature for Serine/Arginine-Rich Protein Splicing Factors (SR Proteins) in Plants. Plant Cell 2010, 22, 2926–2929. [Google Scholar] [CrossRef] [PubMed]
  64. Jia, Y.; Mu, J.C.; Ackerman, S.L. Mutation of a U2 snRNA Gene Causes Global Disruption of Alternative Splicing and Neurodegeneration. Cell 2012, 148, 296–308. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  65. Scharf, K.-D.; Berberich, T.; Ebersberger, I.; Nover, L. The plant heat stress transcription factor (Hsf) family: Structure, function and evolution. Biochim. Biophys. Acta 2012, 1819, 104–119. [Google Scholar] [CrossRef]
Figure 1. Categories and numbers of different types of alternative splicing (AS) events identified in the transcriptomes of soybean Williams 82 under various drought stress conditions. VMS: very mild drought stress, MS: mild drought stress, SS: severe drought stress, SR: water recovery after severe drought stress, SE: skipped exon, MXE: mutually exclusive exon, A5SS: alternative 5′ splice site, A3SS: alternative 3′ splice site, RI: retained intron. The data are expressed as mean ± standard deviation of three biological replicates. **, p < 0.01; *, p < 0.05.
Figure 1. Categories and numbers of different types of alternative splicing (AS) events identified in the transcriptomes of soybean Williams 82 under various drought stress conditions. VMS: very mild drought stress, MS: mild drought stress, SS: severe drought stress, SR: water recovery after severe drought stress, SE: skipped exon, MXE: mutually exclusive exon, A5SS: alternative 5′ splice site, A3SS: alternative 3′ splice site, RI: retained intron. The data are expressed as mean ± standard deviation of three biological replicates. **, p < 0.01; *, p < 0.05.
Genes 11 01520 g001
Figure 2. Comparison of the number and proportions of differential AS events under various drought conditions. (A) Under VMS condition, (B) under MS condition, (C) under SS condition, and (D) under SR condition. Exclusion is shown in blue, and inclusion is shown in red. The y-axis shows the number of AS events. The x-axis shows the types of AS.
Figure 2. Comparison of the number and proportions of differential AS events under various drought conditions. (A) Under VMS condition, (B) under MS condition, (C) under SS condition, and (D) under SR condition. Exclusion is shown in blue, and inclusion is shown in red. The y-axis shows the number of AS events. The x-axis shows the types of AS.
Genes 11 01520 g002
Figure 3. Cluster analysis of stress-responsive AS events in soybean root displayed on the basis of their isoform expression percentage (IEP) variations under various drought conditions. (IEP = Average ∆PSI (ww)/(Average ∆PSI (ww) + Average ∆PSI (stress)). The number of AS events in each cluster is listed. The x-axis shows various stress conditions, and the y-axis shows IEP value. The black lines indicate the trends of average IEP value for all AS events in each cluster.
Figure 3. Cluster analysis of stress-responsive AS events in soybean root displayed on the basis of their isoform expression percentage (IEP) variations under various drought conditions. (IEP = Average ∆PSI (ww)/(Average ∆PSI (ww) + Average ∆PSI (stress)). The number of AS events in each cluster is listed. The x-axis shows various stress conditions, and the y-axis shows IEP value. The black lines indicate the trends of average IEP value for all AS events in each cluster.
Genes 11 01520 g003
Figure 4. Venn diagram of differentially expressed genes (DEGs) and differentially spliced genes (DSGs). (A) Comparison analysis of overlapping between DEGs and DSGs under all drought conditions. (B) Comparison of DSGs under VMS, MS, SS, and SR conditions.
Figure 4. Venn diagram of differentially expressed genes (DEGs) and differentially spliced genes (DSGs). (A) Comparison analysis of overlapping between DEGs and DSGs under all drought conditions. (B) Comparison of DSGs under VMS, MS, SS, and SR conditions.
Genes 11 01520 g004
Figure 5. Representation of the Gene Ontology (GO) classification for biological process, molecular function, and cellular component. (A) GO terms assigned to DSGs under various drought conditions (false discovery rate (FDR) value < 0.05). (B) The top 20 GO enrichment terms (molecular function) of DSGs under various drought conditions. Notes: the blue bar chart represents the gene percentage (%) of the GO terms.
Figure 5. Representation of the Gene Ontology (GO) classification for biological process, molecular function, and cellular component. (A) GO terms assigned to DSGs under various drought conditions (false discovery rate (FDR) value < 0.05). (B) The top 20 GO enrichment terms (molecular function) of DSGs under various drought conditions. Notes: the blue bar chart represents the gene percentage (%) of the GO terms.
Genes 11 01520 g005
Figure 6. Expression and AS analysis of splicing regulatory factors (SPF)-related genes in response to various drought treatments. (A) The 20 most enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (G. max. L.) of DSGs. The −Log10 (p value) was used to indicate the ratio of the DSG number and the number of genes annotated in this pathway. The greater the −Log10 (p value), the greater the degree of enrichment. **, p < 0.01; *, p < 0.05. (B) The number of significantly differentially spliced SPF-related genes identified under various drought retreatments. (C) The number of significantly differentially expressed SPF-related genes identified under various drought retreatments.
Figure 6. Expression and AS analysis of splicing regulatory factors (SPF)-related genes in response to various drought treatments. (A) The 20 most enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (G. max. L.) of DSGs. The −Log10 (p value) was used to indicate the ratio of the DSG number and the number of genes annotated in this pathway. The greater the −Log10 (p value), the greater the degree of enrichment. **, p < 0.01; *, p < 0.05. (B) The number of significantly differentially spliced SPF-related genes identified under various drought retreatments. (C) The number of significantly differentially expressed SPF-related genes identified under various drought retreatments.
Genes 11 01520 g006
Figure 7. Visualization of representative transcript isoforms of four SPF-related genes to indicate expression changes under various drought stresses. MSTRG indicates the novel transcripts in this study. (A) Splicing factor 3b subunit (Glyma.13G241800) switched between MSTRG.29150.1 and Glyma.13G241800.1 under MS and SS conditions. Novel transcripts are illustrated as MXE events. (B) UPF3 regulator (Glyma.06G324500) switched between MSTRG.13253.3 and Glyma.06G324500 under MS and SS conditions. Novel transcripts are illustrated as A3SS events. (C) Serine arginine-rich splicing factor (Glyma.09G104200) switched between MSTRG.19304.3 and Glyma.09G104200.1, Glyma.09G104200.2, and Glyma.09G104200.3 under VMS, MS, and SS conditions, respectively. Novel transcripts are illustrated as A5SS events. (D) Serine/arginine-rich splicing factor SR30 (Glyma.20G236800) switched between MSTRG.44538.1 and Glyma.20G236800.1 under MS and SS conditions. Novel transcripts are illustrated as MXE events. The color of the functional gene ID is blue, while the color of nonfunctional gene ID is black.
Figure 7. Visualization of representative transcript isoforms of four SPF-related genes to indicate expression changes under various drought stresses. MSTRG indicates the novel transcripts in this study. (A) Splicing factor 3b subunit (Glyma.13G241800) switched between MSTRG.29150.1 and Glyma.13G241800.1 under MS and SS conditions. Novel transcripts are illustrated as MXE events. (B) UPF3 regulator (Glyma.06G324500) switched between MSTRG.13253.3 and Glyma.06G324500 under MS and SS conditions. Novel transcripts are illustrated as A3SS events. (C) Serine arginine-rich splicing factor (Glyma.09G104200) switched between MSTRG.19304.3 and Glyma.09G104200.1, Glyma.09G104200.2, and Glyma.09G104200.3 under VMS, MS, and SS conditions, respectively. Novel transcripts are illustrated as A5SS events. (D) Serine/arginine-rich splicing factor SR30 (Glyma.20G236800) switched between MSTRG.44538.1 and Glyma.20G236800.1 under MS and SS conditions. Novel transcripts are illustrated as MXE events. The color of the functional gene ID is blue, while the color of nonfunctional gene ID is black.
Genes 11 01520 g007
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Song, L.; Pan, Z.; Chen, L.; Dai, Y.; Wan, J.; Ye, H.; Nguyen, H.T.; Zhang, G.; Chen, H. Analysis of Whole Transcriptome RNA-seq Data Reveals Many Alternative Splicing Events in Soybean Roots under Drought Stress Conditions. Genes 2020, 11, 1520. https://0-doi-org.brum.beds.ac.uk/10.3390/genes11121520

AMA Style

Song L, Pan Z, Chen L, Dai Y, Wan J, Ye H, Nguyen HT, Zhang G, Chen H. Analysis of Whole Transcriptome RNA-seq Data Reveals Many Alternative Splicing Events in Soybean Roots under Drought Stress Conditions. Genes. 2020; 11(12):1520. https://0-doi-org.brum.beds.ac.uk/10.3390/genes11121520

Chicago/Turabian Style

Song, Li, Zhenzhi Pan, Lin Chen, Yi Dai, Jinrong Wan, Heng Ye, Henry T. Nguyen, Guozheng Zhang, and Huatao Chen. 2020. "Analysis of Whole Transcriptome RNA-seq Data Reveals Many Alternative Splicing Events in Soybean Roots under Drought Stress Conditions" Genes 11, no. 12: 1520. https://0-doi-org.brum.beds.ac.uk/10.3390/genes11121520

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop