1. Introduction

Drought is an extreme abiotic stress condition for plants that results in a staggering 40.8% loss in crop production, based on an analysis of the last 40 years of agronomic data in the USA (Boyer, 1982). Therefore, enhancing drought tolerance in crops is a very desirable target for agriculture (Bailey-Serres et al., 2019) to help meet the demands of the worldwide human population (FAO et al., 2020).

Plants may exhibit a degree of drought memory in response to their exposure to recurring or chronic drought stress, which may increase their drought tolerance and survival under water scarcity conditions (Avramova, 2015; Li et al., 2019). Several studies have investigated drought stress memory by transcriptome deep sequencing (RNA-seq) in Arabidopsis (Arabidopsis thaliana) (Ding et al., 2012, 2013; Liu et al., 2016), maize (Zea mays) (Ding et al., 2014; Virlouvet et al., 2018), and rice (Oryza sativa) (Li et al., 2019). However, the molecular mechanisms underlying drought memory and how it is regulated remain poorly understood.

It is, however, clear that drought stress is associated with alternative splicing (AS) in plants, such as in maize (Thatcher et al., 2016), rice (Wei et al., 2017; Zhang and Xiao, 2018), and wheat (Triticum aestivum) (Liu et al., 2018). Alternative splicing, first discovered in 1977, is the process by which one pre-mRNA generates two or more distinct transcripts, and constitutes one of the most significant post-transcriptional mechanisms controlling gene expression (Lee and Rio, 2015; Zhang et al., 2015). Alternative splicing can be classified into five basic/simple types: intron retention (IR), alternative donor (A5SS), alternative acceptor (A3SS), exon skipping (ES), and mutual exon skipping (MXE) (Shen et al., 2014; Shi et al., 2019). Most previous reports have indicated that IR is the most frequent AS type in plants, based on genome-wide analyses (Filichkin et al., 2010, 2015; Reddy et al., 2013; Thatcher et al., 2016; Wei et al., 2017; Wang et al., 2019). By contrast, in animals the most frequent AS event is ES (Reddy et al., 2013; Zhang et al., 2017b; Wang et al., 2019; Jia et al., 2020). A handful of reports in plants have identified A3SS (Tang et al., 2016; Shi et al., 2019) or ES (Xu et al., 2019; Xia et al., 2020) as the most frequent AS event, when considering large datasets comprising over 1000 AS events. Among many software tools exist for AS detection, AStalavista and rMATS are widely accepted in plants (Song et al., 2019).

The canonical GT-AG dinucleotide sequence is highly over-represented at splicing sites, and non-canonical dinucleotides are rare (Sharp and Burge, 1997). While the U2-type spliceosome commonly uses GT-AG intron sites, AT-AC can be used by the U12-type spliceosome (Tarn and Steitz, 1996; Sharp and Burge, 1997). U12 and U2 spliceosomes are conserved between animals and plants (Lorkovic et al., 2005; Reddy et al., 2013), as are the U12- and U2-type intron positions (Basu et al., 2008). In Arabidopsis, U12-type introns are enriched in genes related to DNA/RNA metabolism (Marquez et al., 2012). The conserved motif at the 5′ splice site is AAG|GTATGTT, AAG|GTAAGTT, and CAG|GTAAGTA in budding yeast (Saccharomyces cerevisiae), fission yeast (Schizosaccharomyces pombe), human (Homo sapiens), and mouse (Mus musculus), respectively (Ast, 2004). Moso bamboo (Phyllostachys edulis) shows the conserved 5’ motif (A)_nG|GTAAG(T)_n (Wang et al., 2017).

Here, we reveal the connection between alternative splicing events and the molecular and functional characterization of drought memory genes in rice. We discovered that ES events may be induced by drought memory. Our results may provide new targets for breeding crops with enhanced drought tolerance in the future.

2. Material and methods 2.1. Plant materials and sampling

Seeds of the rice cultivar “Zhonghua 11” (Oryza sativa L. ssp. japonica) were surface sterilized with 2.5% NaClO for 15 min, rinsed with sterilized water, and then allowed to soak in sterilized water at 37 °C in an incubator. After germination, the seedlings were transferred to culture vessels containing one-quarter-strength modified Hoagland solution (Jones, 1982). Seedlings were grown in this solution for four weeks under a 12 h light: 12 h dark photoperiod at 28 °C. The Hoagland solution was changed every 2 days. We then selected four-week-old plants for the well-watered control samples (R0). We followed previously described procedures for drought memory experiments (Ding et al., 2012; Li et al., 2019) by exposing plants to air drying for 80 min at 28 °C; these samples were labeled as first drought stress treatment (S1). Plants were then re-watered after drought treatment for 22 h 40 min at 28 °C, constituting samples for first recovery treatment (R1). We repeated these stress/recovery cycles three times to enhance drought memory, annotating each sample stage as S2, R2, S3, R3, and S4 for the 2nd drought stress, 2nd recovery, 3rd drought stress, 3rd recovery, and 4th drought stress, respectively (Ding et al., 2012, 2013; Li et al., 2019).

2.2. Library preparations for RNA-seq

Total RNA from the leaves of R0, S1, R3, and S4 stage with three biological replicates was isolated with TRIzol reagent, according to the manufacturer's instructions (Invitrogen, Carlsbad, CA, United States), and then treated with RNase-free DNase I to remove trace DNA contamination. Total RNA was quantified on an Agilent Bioanalyzer 2100 instrument, and the libraries were constructed and sequenced on an Illumina HiSeq 2500 sequencing platform by Genedenovo Biotechnology, Co., Ltd. (Guangzhou, China) to generate 150-bp paired-end reads.

2.3. Data filtering, mapping of reads and differentially expressed genes (DEGs) analysis

The raw reads were filtered using Trimmomatic v.0.36 (Bolger et al., 2014) to remove adaptor sequences and low-quality reads. The filtered reads were then mapped with Hisat2 (Kim et al., 2015) to the rice japonica reference genome (ftp://ftp.ensemblgenomes.org/pub/plants/release-43/gtf/oryza_sativa/Oryza_sativa.IRGSP-1.0.43.gtf.gz, ftp://ftp.ensemblgenomes.org/pub/plants/release-43/fasta/oryza_sativa/dna/Oryza_sativa.IRGSP-1.0.dna.toplevel.fa.gz). The resulting output files were in SAM (Sequence Alignment/Map) format. We used SAMtools v.1.3.1 (Li et al., 2009) to convert SAM files into BAM (binary format of SAM file) format, followed by StringTie v.1.3.6 (Pertea et al., 2015) software to generate assembled GTF format files. We used DESeq2 v.1.32.0 (Love et al., 2014) to identify differentially expressed genes (DEGs) with log₂ (|fold change|)>1 and adjusted P value < 0.05.

2.4. Analysis of AS events

We used AStalavista v.4.0 software (Foissac and Sammeth, 2007) with the parameters (-t asta -i) to quantify the type of AS events from the merged GTF file for all assembled transcripts. We investigated the codes for AS events according to the notes from the software author. For simple AS events, AStalavista software defined AS code “0, 1–2ˆ” for exon skipping (ES), “1ˆ, 2ˆ” for alternative donor (A5SS), “1-, 2-” for alternative acceptor (A3SS), “'0, 1ˆ2-” for intron retention (IR), and “'1–2ˆ, 3–4ˆ” for mutual exon skipping (MXE). For transcripts with more than one AS event, AStalavista uses more complex codes (Sammeth et al., 2008; Foissac and Sammeth, 2015). The AStalavista webserver (v.2.2) uses “Others” for such complex events (http://astalavista.sammeth.net/). Thus, we extracted five basic AS events types (IR, A5SS, A3SS, MXE, and ES) and the remained complex AS events were grouped as “Others” type, as described previously (Foissac and Sammeth, 2007; Min et al., 2015; Li et al., 2017).

We also used rMATS v.4.0.2 software (Shen et al., 2014) to detect differential AS events with parameters --nthread 8 --cstat 0.0001 --tstat 6 --readLength 150. Splicing events were represented with rmats2sashimiplot (https://github.com/Xinglab/rmats2sashimiplot). Dinucleotides at splicing sites were counted by comparing our GTF files with the genome fasta file using in-house scripts. Flanking sequences of splicing sites were extracted with in-house scripts, then uploaded to WebLogo v.3 (http://weblogo.threeplusone.com/) to obtain conserved sequences.

We identified transcription factors (TFs) from genes showing differential alternative splicing (DAS) by querying the PlantTFDB 4.0 database (Jin et al., 2017).

2.5. Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) analysis

The agriGO (Du et al., 2010) software was used to identify enriched Gene Ontology (GO) terms. ClusterProfiler (Yu et al., 2012), and in-house scripts were employed to perform KEGG analysis. Protein–protein interaction (PPI) analysis was performed with the STRING database (Szklarczyk et al., 2019). Cytoscape 3.6.1 (Shannon et al., 2003) and TBtools v.0.66813 (Chen et al., 2018) software were used to generate figures.

2.6. RT-qPCR analysis

We initiated first-strand cDNA synthesis from 1 μg total RNA with oligo (dT) primer using the PrimeScript™ RT Reagent Kit (Takara, Japan). The relative expression level of each gene was measured with gene-specific primers (Supplementary Table S1). Real-time quantitative PCR (RT-qPCR) was performed with SYBR premix Ex Taq II kit (TaKaRa, Japan) in a total reaction volume of 20 μl on a CFX96™ system (Bio-Rad, CA, USA). The internal control was the rice ELONGATION FACTOR-1 (EF-1a, LOC_Os03g08020) gene (Li et al., 2019).

3. Results 3.1. Quality of the RNA-seq data

To identify alternative splicing (AS) events during drought stress, we collected rice samples from well-watered control plants (R0), plants exposed to drought stress once (S1) or four times (S4), and plants re-watered after the third drought stress cycle (R3). In total, we obtained 509 million raw reads corresponding to about 76.4 Gb data across all 12 samples. After quality-trimming, we retained 486 million reads (~65.3 Gb), with an average of 40 million reads per sample. We then aligned these clean reads to the rice reference genome (IRGSP-1.0) downloaded from the Ensembl website. On average, 96% of reads mapped to the rice genome, confirming the relative high quality of the RNA-seq data (Supplementary Table S2). Principal component analysis (PCA) profiling transcript patterns in rice leaf samples between well-watered and drought stress conditions revealed large contributions of repeated stress (50% of standing variance) and drought (12% of variance) to the total observed variance (Fig. 1A).

3.2. Analysis of differentially expressed genes (DEGs)

We identified 8679 DEGs when using well-watered samples (R0) as control (Fig. 1B and C). The number of DEGs almost doubled between the first and fourth drought-stress treatments (S1: 2549 DEGs; R3: 5781 DEGs; and S4: 4836 DEGs) (Supplementary Tables S3, S4 and S5). This large difference indicates that repeated cycles of drought and recovery influence gene expression in a more pronounced manner, presumably through drought memory.

To further investigate the functions associated with the three groups of DEGs mentioned above, we performed GO and KEGG enrichment analysis. We identified the GO terms GO:0019748 (“secondary metabolic process”), GO:0009875 (“pollen-pistil interaction”), GO:0007049 (“cell cycle”), GO:0009856 (“pollination”), GO:0051704 (“multi-organism process”), GO:0006259 (“DNA metabolic process”), GO:0006950 (“response to stress”), and GO:0006629 (“lipid metabolic process”) as being enriched in different drought conditions (Fig. 2A). After the 1st drought stress (R0_vs_S1), the enrichment levels of GO terms such as “secondary metabolic process” and “lipid metabolic process” significantly changed, suggesting that these two processes may be important for drought memory (Fig. 2A). KEGG enrichment analysis revealed a significant enrichment for “phenylpropanoid biosynthesis” (dosa00940), “plant hormone signal transduction” (dosa04075), “diterpenoid biosynthesis” (dosa00904), “MAPK signaling pathway – plant” (dosa04016), “DNA replication” (dosa03030), “zeatin biosynthesis” (dosa00908), “linoleic acid metabolism” (dosa00591), “starch and sucrose metabolism” (dosa00500), “galactose metabolism” (dosa00052), “brassinosteroid biosynthesis” (dosa00905), and “carotenoid biosynthesis” (dosa00906) (Fig. 2B, Supplementary Table S6). Of these KEGG terms, “phenylpropanoid biosynthesis” and “plant hormone signal transduction” showed the highest enrichment, suggesting that the associated pathways are highly affected by drought stress.

3.3. Alternative splicing events

We detected various types of AS events among our assembled transcripts, as determined by the AStalavista software. In this study, we extracted 24,558 AS events (Supplementary Table S7), which fell into six types: IR, A5SS, ES, A3SS, MXE, and Others. The numbers for each AS event type are listed in Table 1.

To investigate splicing events from another angle, we used the software rMATS. While AStalavista calculates the splicing events based on the GTF file of assembled transcripts at the whole-genome level, rMATS calculates differential AS events between two samples based on aligned BAM files. We thus subjected our samples to rMATS analysis to extract splicing events. Notably, the ES type was the predominant form of AS identified by rMATS. In addition, drought stress conditions (R0_vs_S1 and R0_vs_S4) resulted in more AS events (8456 and 9,586, respectively) than the stress recovery condition (R0_vs_R3, with 7265 AS events) (Fig. 3A, Table 2). This result indicates that drought stress induces differential alternative splicing.

To validate these results, we expanded our analysis of AS events to additional plant species exposed to drought stress with rMATS. We downloaded RNA-seq datasets from seven species from the Sequence Read Archive (SRA) database at the National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov/sra) (Supplementary Table S8), and processed the reads with the same pipeline used for our rice data to detect differential AS events. We discovered that ES was the largest AS type not only in rice (>70.73%), but also in sorghum (Sorghum bicolor, with 60.70%) and purple false brome (Brachypodium distachyon, with 45.76%). By contrast, IR was the largest class of AS events in Arabidopsis (47.31%) and maize (40.05%) (Fig. 3B). In wheat, ES (35.33%) and A3SS (35.61%) were the dominant AS types, while IR accounted for less than 13% of all AS events. We randomly selected one gene for each AS type to illustrate the expression and splicing events from AS genes (Fig. 4, Supplementary Fig. S1).

We investigated the dinucleotides at splicing sites: we observed that the canonical GT-AG dinucleotide had a frequency of 92% in our RNA-seq data, which was slightly lower than the genome-wide frequency of 95% in rice (Fig. 5B). The non-canonical dinucleotide sequences GC-AG and AT-AC both increased (Fig. 5A). AT-AC is recognized by U12-type splicing, which has been reported to participate in stress transduction (Bai et al., 2019). This observation indicates that drought stress may employ more non-canonical dinucleotide splicing sites. We further extracted flanking sequence (30 bp) around splicing sites and determined that the consensus for the 5′ splicing donor site was G|GTAAGT, while the 3’ acceptor site was (T)_nGCAG|G (Fig. 5C). Similar results were obtained in bamboo (Wang et al., 2017), suggesting a slight enrichment in conserved sequences at splicing sites for plants.

3.4. Drought memory genes with differential alternative splicing events

We previously identified 6885 drought memory genes (Li et al., 2019). In this study, 920 drought memory genes exhibited interesting patterns of differential alternative splicing events (DM-DAS) based on rMATS analysis (Fig. 6). GO and KEGG enrichment analysis showed different results compared to the initial set of 8679 DEGs mentioned above (Fig. 7).

Using these 920 DM-DAS genes as seeds, we performed a search in the STRING database, resulting in the identification of three well-known genes as hubs involved in drought memory: OS03T0738400-01, encoding serinehydroxymethyltransferase 1 (OsSHMT1, LOC_Os03g52840), OsJ_35887, encoding chloroplast stem-loop-binding protein of 41 kDa b (CSP41b, also called Light-Green Leaf 1 [LGL1], LOC_Os12g23180), and Phosphoinositide-dependent protein kinase 1 (PPDK1, LOC_Os05g33570) (Supplementary Fig. S2).

In order to better understand the transcriptional mechanisms underlying drought memory, we interrogated our DM-DAS genes for transcription factors (TFs), leading to the identification of 20 TFs. Two TFs belonged to the auxin response factor (ARF) family: LOC_Os05g48870 (OsETTIN1/OsARF15) and LOC_Os06g46410 (OsARF17). Another three genes encoded members of the plant-specific B3 superfamily (B3) of TFs: LOC_Os01g13300, LOC_Os07g37610 (GERMINATION DEFECTIVE 1, [GD1]), and LOC_Os10g17630. Four genes coded for members of the basic/helix-loop-helix (bHLH) superfamily: LOC_Os01g72370 (IRON-RELATED TRANSCRIPTION FACTOR 2, OsIRO2/OsbHLH056), LOC_Os03g43810 (PHYTOCHROME INTERACTING FACTOR 12, OsPIL12/OsbHLH103), LOC_Os07g43530 (OsbHLH1), and LOC_Os11g38870 (OsbHLH061). Three genes encoded bZIP TFs: LOC_Os01g64730 (ABSCISIC ACID [ABA] RESPONSIVE ELEMENT-BINDING FACTOR 1, ABF1), LOC_Os02g52780 (OsbZIP23), and LOC_Os11g05640 (OsbZIP80/OsZIP-2a). Three genes encoded members of the MYB TF family: LOC_Os12g13570, LOC_Os08g06110 (LATE ELONGATED HYPOCOTYL, OsLHY), and LOC_Os12g41920 (TELOMERE REPEAT-BINDING FACTOR 2, OsTRBF2). Finally, we identified one member each from the TF families G2-like, GRAS, NFYA, and YABBY, respectively: LOC_Os06g24070 (GOLDEN-LIKE 1, OsGLK1), LOC_Os03g48450 (OsGRAS17), LOC_Os02g53620 (NUCLEAR FACTOR-YA11, OsNF-YA11), and LOC_Os02g42950 (OsYABBY4) (Table 3).

We selected six genes for RT-qPCR analysis of their transcript levels. Our results suggested a good correlation between RNA-seq data and RT-qPCR data (Supplementary Fig. S3, Supplementary Table S1).

4. Discussion 4.1. DM-DAS genes might have very important function

Despite the identification of drought memory genes in multiple plant species such as rice (Li et al., 2019), Arabidopsis (Ding et al., 2012), and maize (Virlouvet et al., 2018), the molecular landscape of plant drought memory has remained unclear. Previous studies have suggested that drought memory may rely on alternative splicing of specific genes, for example CBL-INTERACTING PROTEIN KINASE (CIPK) (Wan et al., 2019), DEHYDRATION RESPONSE ELEMENT BINDING 2B (OsDREB2B) (Matsukura et al., 2010), and DROUGHT RESPONSIVE ANKYRIN 1 (DRA1) (Guerra et al., 2015). Here, we reveal a much larger complement of 920 DM-DAS genes in rice, providing a large candidate gene list for future analysis of drought resistance mechanisms (Fig. 6).

We also detected several hub DM-DAS genes which could be widely connected to other genes: for example, the protein encoded by LOC_Os03g52840 (OsSHMT1, OS03T0738400-01) was shown to enhance chilling tolerance and interacted with the heat shock protein Hsp70 (Fang et al., 2020). Transgenic rice plants overexpressing OsSHMT1 exhibited increased photosynthetic efficiency and improved plant productivity (Wu et al., 2015). The gene ZmPPDK1 encodes an essential photosynthetic enzyme in maize, a C4 plant. By contrast, the rice ortholog OsPPDK1 (LOC_Os05g33570, Os05g0405000) has no role in photosynthesis but rather plays an important role in starch biosynthesis and during seed development (Kang et al., 2005; Zhang et al., 2018; Wang et al., 2020).

It has been shown that TF genes also undergo AS events during stress response (Mastrangelo et al., 2012). In this study, we identified 20 TFs out of the set of 920 DM-DAS, some of which have previously been linked to drought stress response (Table 3). For example, OsbZIP23 is a key gene conferring ABA-dependent drought tolerance (Xiang et al., 2008). OsbZIP23 is also associated with histone H3K4me3 modification at the chromatin around dehydrin genes to regulate their expression (Zong et al., 2020). Likewise, OsABF1 regulates the expression of PROTEIN PHOSPHATASE 2C (PP2C) and bZIP23, which then feed back onto the drought/ABA signaling pathway (Zhang et al., 2017a). We therefore hypothesize that DM-DAS genes have great potential as drought-related genes, and provide a list of high-confidence candidates for future analysis.

4.2. Why is ES the most frequent AS events type under drought stress conditions in rice?

It is thought that intron retention (IR) constitutes the main type of AS event in plants, while exon skipping (ES) is the most prominent AS type in animals (Reddy et al., 2013; Zhang et al., 2015). However, our study suggests that ES is the dominant AS type among DAS events in rice plants under drought stress (Fig. 3A).

One possible reason to explain discrepancies between earlier analysis and our own was software bias, since rMATS and AStalavista were both designed to analyze animal data where ES is predominant. To eliminate this possibility, we investigated other plant datasets using rMATS. Notably, the reanalysis of these publicly available datasets showed the same high ES frequency in sorghum (60.70%), whereas ES frequency in maize, Arabidopsis and barley (Hordeum vulgare) were much lower (9.94%–16.67%), and the frequency in wheat and purple false brome was intermediate (35.33% and 45.76%, respectively) (Fig. 3B). These variable frequencies did not support a systematic bias caused by the software platform used.

A second potential reason for the differences in outcome called upon missing codes for complex AS events in AStalavista in a subset of earlier studies, as not all reports included complex/Others AS types when using AStalavista on plant samples. For example, when using only the simplest AStalavista code “0, 1–2ˆ” to extract ES events (Foissac and Sammeth, 2015), many ES events including under the “Others” category would be missed. A previous study in Oryza sativa ssp. japonica, Oryza sativa ssp. indica, sorghum and maize using AStalavista stated that IR constituted the “major common AS type” in these four species/subspecies, while the frequency of complex/Others types of AS events were 37.8%, 16.9%, 61.7%, and 20.4%, respectively (Min et al., 2015). In Arabidopsis, simple IR code: “0, 1ˆ2-” from AStalavista extracted just a little over half of all IR events, yielding a frequency of 24.21% when the total IR frequency is in fact ~40% (Marquez et al., 2012). In cabbage (Brassica oleracea), A3SS (10.31%) was the most frequent simple AS event, while complex AS (69.86%) was predominant (Xu et al., 2019). Thus, ES frequencies in plants may have been under-estimated. Since the rMATS software was released in 2014, few have paid attention to the relative frequencies of AS types, and our understandings of this area are under the ongoing development. When we applied AStalavista to all the data previously analyzed by rMATS, IR was the most prevalent type, while the frequency of complex/Others type varieties from 15.07% to 48.77% (Supplementary Table S9). These results therefore indicate that many ES events that belong to the complex/Others type may be under-represented in several earlier studies.

Last, another possible reason for the relative frequencies of AS events observed here was annotation bias in each genome. However, the majority of AS events corresponded to IR (83%) in 141 rice samples under different nutrient conditions that were analyzed by rMATS with two different annotations (MSU and RAPDB) of the rice genome (Dong et al., 2018). This study saw no annotation bias when using rMATS.

We thus hypothesize that 1) the high frequency of ES events during drought memory conditions in rice may be a direct consequence of drought stress on AS, and 2) the relative frequency of ES events is species specific. We propose the term “drought-stress related ES” for this hypothesis.

5. Conclusion

This study may be the first report to show that ES is the predominant form of alternative splicing in rice plants subjected to drought stress. This result differs greatly from those of most previous studies in plants, and may challenge currently accepted knowledge in the field. Based on our analysis and reasoning, we propose the hypothesis — “drought-stress related ES” — to explain why ES is the main AS event during drought stress in rice. This hypothesis should be further examined experimentally and computationally.

Declaration of competing interest

The author declares no conflict of interest.

Acknowledgements

This study was funded by National Natural Science Foundation of China (31571311 to CZ and 31971410 to LL) and the CAS Pioneer Hundred Talents Program (CZ and LL).

Appendix A Supplementary data

The following are the Supplementary data to this article: https://doi.org/10.1016/j.pld.2020.11.004.

Author contributions

CJ, LL, GJ and HY designed the experiments and analysis of data, drafting or revising the article. PL, DG and HY plant culture and performed the experiments. All authors reviewed and approved the manuscript.