Analysis of Whole-Transcriptome RNA-Seq Data Reveals the Involvement of Alternative Splicing in the Drought Response of Glycyrrhiza uralensis

Abstract

Alternative splicing (AS) is a post-transcriptional regulatory mechanism that increases protein diversity. There is growing evidence that AS plays an important role in regulating plant stress responses. However, the mechanism by which AS coordinates with transcriptional regulation to regulate the drought response in Glycyrrhiza uralensis remains unclear. In this study, we performed a genome-wide analysis of AS events in G. uralensis at different time points under drought stress using a high-throughput RNA sequencing approach. We detected 2,479 and 2,764 AS events in the aerial parts (AP) and underground parts (UP), respectively, of drought-stressed G. uralensis. Of these, last exon AS and exon skipping were the main types of AS. Overall, 2,653 genes undergoing significant AS regulation were identified from the AP and UP of G. uralensis exposed to drought for 2, 6, 12, and 24 h. Gene Ontology analyses indicated that AS plays an important role in the regulation of nitrogen and protein metabolism in the drought response of G. uralensis. Notably, the spliceosomal pathway and basal transcription factor pathway were significantly enriched with differentially spliced genes under drought stress. Genes related to splicing regulators in the AP and UP of G. uralensis responded to drought stress and underwent AS under drought conditions. In summary, our data suggest that drought-responsive AS directly and indirectly regulates the drought response of G. uralensis. Further in-depth studies on the functions and mechanisms of AS during abiotic stresses will provide new strategies for improving plant stress resistance.

Keywords: Glycyrrhiza uralensis; alternative splicing; drought regulators; drought response; splicing regulatory factor; transcriptome analysis.

FIGURE 1

Evaluation of RNA sequence sample dataset used to identify alternative splicing (AS) events. (A) Principal component analysis of RNA sequence sample data. Samples from aboveground and underground parts formed different clusters at different time points (APDS_2h, APDS_6h, APDS_12h, APDS_24h, UPDS_2h, UPDS_6h, UPDS_12h and UPDS_24h). (B) Correlation analysis between samples. Correlation matrices were calculated by comparing the mean values of entire transcriptome in samples from aboveground and underground parts at different time points during drought treatment. Pearson’s correlation coefficients between samples were analyzed using an R script.

FIGURE 2

Classes and numbers of different types of AS events detected in the transcriptome of Glycyrrhiza uralensis at different time points during drought treatment. (A) Schematic diagram of the AS: SE, A3SS, A5SS, RI, and MXE. Each event produces two types of isomers: isoforms 1 and 2. (B) Number of five different types of AS events in AP of G. uralensis at different time points under drought stress. (C) Number of five different types of AS events in UP of G. uralensis at different time points under drought stress. Data are mean ± standard deviation of three biological replicates. Mean values for different time points under drought stress were tested using the Bonferroni test. Different letters indicate highly significant differences among treatments.

FIGURE 3

Identification and comparative analysis of stress-responsive AS events in G. uralensis. (A) Comparison of number and proportion of different AS events in different parts of G. uralensis. Exclusion is shown in blue, inclusion in red. y-axis shows number of AS events, x-axis shows types of AS events. (B) Cluster analysis of stress-responsive AS events in G. uralensis, as indicated by changes in isoform expression percentage (IEP) under various drought conditions. [IEP = mean PSI (ck)/(mean PSI (ck) + mean PSI (stress)]. Number of AS events in each cluster is listed. x-axis shows various stress conditions and y-axis shows IEP values. Red line shows the trend of average IEP value for all AS events in each cluster.

FIGURE 4

Comparative analysis of differentially spliced genes (DSGs) and differentially expressed genes (DEGs) in G. uralensis at different time points under drought stress. (A) Comparative analysis of overlapping DEGs and DSGs under all drought conditions. (B) Number of DSGs and DEGs in AP and UP, and those overlapping between AP and UP, at different time points during the drought treatment.

FIGURE 5

Gene Ontology (GO) analyses of DSGs. (A) GO terms in the biological process and cellular component categories assigned to DSGs under various drought conditions [false discovery rate (FDR) values <0.05]. (B) Top 20 GO terms (molecular function) enriched with DSGs under different drought conditions. Note: Histograms represent the percentage of genes (%) for GO terms.

FIGURE 6

Expression and AS analysis of SPF-related genes in different drought treatments. (A) Top 20 pathways enriched with DSGs. −Log10 (p-value) expresses the ratio of the number of DSGs to the number of genes annotated in the pathway. The larger the −Log10 (p-value), the greater the enrichment. **, p < 0.01; *, p < 0.05. (B) Number of significantly differentially spliced SPF-associated genes at different time points under drought stress. (C) Number of significantly differentially expressed SPF-associated genes at different times points under drought stress.

FIGURE 7

Experimental validation of stress response events in licorice AP and UP by RT-PCR. (A) Gene encoding 22 kDa heat shock protein; (B) Gene encoding probable WRKY transcription factor 19; (C) Gene encoding probable transcription factor 20; (D) Gene encoding probably protein phosphatase 2C 52. In the left panels, bars show relative expression levels of alternatively spliced isoforms 1 (blue) and 2 (red) in AS at 0, 2, 6, 12, and 24 h of drought stress as determined from RNA-seq data. Panels on the right show results of RT-PCR analyses to detect splice variants using specific primer pairs.