Genome-wide transcriptomic analysis identifies candidate genes involved in jasmonic acid-mediated salt tolerance of alfalfa

Plant materials

Cultivar XinJiangDaYe alfalfa used in this study was provided by professor Zhizhong Cao, Gansu Agricultural University. Alfalfa seeds were sowed in pots with a mixture of peat and vermiculite (v:v, 2:1) and placed in the artificial climate chamber. The growth condition was as follows: 30 °C at day and 25 °C at night, 16 h light and 8 h darkness, a relative humidity at 70%. 3-weeks old plants with identical growth characteristics were transplanted into 1/2 Hoagland nutrient solution and treated with salt (S, 150 mM), Jasmonic acid (J, 5 μM), or salt+JA (JS) after 1-week adaptional growth. The same nutrient solution without salt or JA was used as the control (C). To perform transcriptomic analysis, root samples were collected from three biological replicates after 24 h of treatment.

Isolation of RNA and transcription in reverse

The total RNA was isolated using the Plant total RNA purification kit (Gmbiolab. Co., Ltd, Dali City, Taiwan). RNA degradation and contamination were assessed using 1% agarose gels. To assess the purity and quantity of RNA, a NanoPhotometer spectrophotometer was used (IMPLEN, Calabasas, CA, USA). The input material for RNA sample preparations was 5 g of RNA per sample. Total RNA (500 ng) was reverse transcribed to cDNA using a ReverTra Ace qPCR RT Kit with genome-DNA-removing enzyme (Toyobo, Osaka, Japan).

Construction of a cDNA library and sequencing

After separating the total RNA, ribosomal RNA was removed to get messenger RNA (mRNA). The RNA was then fragmented using a fragmentation buffer. cDNA was generated using random hexamers containing dNTPs (dUTP, dATP, dGTP, and dCTP) and DNA polymerase I, and then purified using AMPure XP beads. The purified double-stranded DNA was then subjected to end repair. To create the final cDNA library, a tail well was introduced to connect the sequencing adaptor, AMPure XP beads were used for fragment size selection, and PCR enrichment was performed. The sequence information was uploaded in figshare (DOI 10.6084/m9.figshare.21628505).

Library quality assessment

After the development of the library, its quality was assessed to confirm that it met sequencing standards. To complete the library inspection, the detection techniques include: (1) the use of Qubit2 for further integrity testing; and (2) the -PCR method to appropriately estimate the effective concentration of the library (library effective concentration >2 M). Following validation of the onboard sequencing library, several libraries were pooled and sequenced on the Illumina HiSeq platform based on the desired off-board data volume. The approach was then proceeded by a bioinformatics analysis in which off-board data was filtered using FASTP to produce Clean Data (Chen et al., 2018). (1) Removal of amplicons with adapters; (2) removal of paired reads when the N content in any sequencing read surpasses 10% of the total number of bases in the read; and (3) removal of paired reads when any sequencing read contained low-quality (Q-20) and the total number of bases in the read exceeded 50% of the total number of bases in the read. Additional structural level analyses were performed on the resulting clean mapped data.

Alignment of data, genes expression profiles, and differentially expressed genes

The clean reads were aligned to the gene set and compared to the reference genome using HISAT2 software to calculate the mapping event on both the sense and antisense strands. We used the following method for calculating FPKM (Fragments Per Kilobase of Transcript per Million Mapped Fragments) as a measure of transcript or gene expression levels: FPKM = mapped fragments of transcript/Total Count of mapped fragments (Millions) × Length of transcript (kb). To identify differentially expressed genes (DEGs), the raw counts were used as input and analyzed with the DESeq R package (Varet et al., 2016). The thresholds for identifying DEGs were false discovery rate (FDR) <0.05 and |log2Fold Change| ≥ 1. Gene expression comparisons were conducted among different treatments.

Gene expression cluster analysis

Cluster analysis is used to determine the expression patterns of DEGs under different treatments, and to identify the functions of unknown genes or unknown functions of known genes with the same or similar expression patterns.

Cluster analysis was conducted to investigate gene expression patterns under different treatments. The FPKM values were first centralized and normalized, and then K-means clustering analysis was done. The same class of genes had similar trends under different experimental treatments and may have similar functions.

Functional annotation and enrichment analysis of differentially expressed genes

The GO function classification was evaluated and visualized, hypergeometry was assessed, and validation and screening of genes enriched with GO significant terms were performed. Finally, DEGs were subjected to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using the KOBAS online software and statistical analysis. The DEGs were further functionally annotated by comparing the protein or cDNA sequence to the KOG database and then obtaining the annotation from the KOG database.

Variants of a single nucleotide and alternative splicing

We assessed the expression of alternative splicing events in biologically replicated samples using rMATS, estimated the P value using the likelihood-ratio test to indicate the difference between the two sets of samples, and then applied Benjamini Hochberg Methods. FDR was obtained by utilizing multiple hypothesis testing to adjust the P value. Differential events were defined as alternative splicing events with an FDR less than 0.05.

Validation of genes and expression analysis

To confirm the results, selected DEGs were subjected to real-time quantitative PCR (q-PCR) using primers identified by Primer Premier software (Premier Biosoft International, Palo Alto, CA, USA) (Table S1). RT-qPCR was used to amplify cDNAs in a final volume of 20 μL, which included 2 μL of cDNA, 10 μL of SYBR Green master mix with low Rox, 0.5 μL forward and reverse primers, and 7 l nuclease-free water. On an ABI Quantstudio 6 Flex real-time PCR apparatus, 50 °C for 20 s, 95 °C for 10 min, followed by 45 cycles of 3 min at 94 °C, 15 s at 94 °C, 15 s at 58 °C, and 20 s at 72 °C were used to standardize amplification (Applied Biosystems, Foster City, CA, USA). The housekeeping gene, MsACTIN2 was used as an internal control for normalizing gene expression. To ensure reproducibility and reliability, q-PCR analysis was performed on three independent biological replicates.

Analytical statistics

A pairwise comparison of phenotypic and relative expression levels between samples was performed using Duncan’s multiple range test with a significance threshold of 0.05. All tests were carried out with IBM SPSS Statistics for Windows 19.0 (IBM Corp, Armonk, NY, USA).

Transcript profile of DEGs in the root of alfalfa under different treatment

In order to analyze the effects of JA on the salt stress response, the differentially expressed genes (DEGs) were compared between C_vs_S and S_vs_JS, as well as between C_vs_J and J_vs_JS. There were 4,286 and 2,896 genes that were significantly regulated by salt and JA treatment when compared to control, respectively. C_vs_S and S_vs_JS shared 151 genes and C_vs_J and J_vs_JS shared 185 genes. Importantly, 507 and 1,365 genes were detected exclusively in S_vs_JS and J_vs_JS analysis respectively (Figs.1A and 1B). Heatmap of the FPKM of all DEGs obtained by pairwise comparison further revealed that a series of genes were specifically regulated by JA and salt combined treatment (Fig. 1C, Table S2). We focused on the DEGs regulated by salt alone relative to control and JS relative to S. Volcano plots revealed that most DEGs showed in C_vs_S and S_vs_JS comparisions (Figs. 1D and 1E). Remarkably, a set of genes were further altered when JA was added to salt stress compared to salt stress alone (Fig. 1E, Table S2).

Transcription pattern analysis of all DEGs was analyzed by kmeans and six major expression patterns were identified. Total of 1,129 genes in cluster 1 showed upregulation under salt stress, whereas downregulation under JS conditions compared to S (Fig. 2A). Cluster 2 contains 1,249 genes which were upregulated by salt treatment alone and further induced by JS treatment (Fig. 2B). Expression levels of genes in other four clusters were equal under S and JS conditions (Figs. 2E and 2F). In addition, the DEGs regulated by both S and JS treatment but in opposite directions were analyzed. As shown in Table S3, there were 74 genes were upregulated by salt stress, which were compromised by JS treatment. Some of them were related to plant growth response, such as two Expansin genes (MS.gene011579, MS.gene051061, marked in Figs.1D and 1E). Meanwhile, 28 genes down-regulated by salt stress were upregulated in S_vs_JS analysis, such Fe-Mn superoxide dismutase (Fe-Mn SOD, Gene ID: MS.gene013412, marked in Figs.1D and 1E) and peroxidase family protein (POD, Gene ID: MS.gene02080, marked in Figs.1D and 1E). Meanwhile, 42 were commonly upregulated in both C_vs_S and S_vs_JS analysis, such as some salt-responsive TFs (bZIP, MS.gene07054; ERF, MS.gene50632) and an auxin efflux carrier (PIN-LIKES 7, MS.gene86524). Another seven genes were commonly downregulated.

Functional annotation for DEGs

To identify critical gene functions and pathways in large collections of DEGs, we performed a gene set enrichment analysis and retrieved individual genes from overrepresented gene ontology (GO) categories. For each treatment, we looked at the top 50 enriched ontologies. Intriguingly, most ontologies enriched in C_vs_S and S_vs_JS were different, excluding glutathione transferase activity, hexitol dehydrogenase activity, and mannitol dehydrogenase activity in molecular functions (Fig. 3). In detail, the most significant ontology in salt-only regime was a biological process with 28 enrichments. Response to hydrogen peroxide, cellular amino acid catabolic process, and pathogenesis had the most significant number of genes (30), followed by the alpha amino acid catabolic process, which had 28 genes. The molecular function was the second most annotated ontology, with 22 classifications. The two enrichments were oxidoreductase activity, acting on single donors with the incorporation of molecular oxygen and two atoms of oxygen, with 46 and 37 genes, respectively. Nutrient reservoir activity and protein kinase binding were the following enrichments with 35 and 32 genes, respectively (Fig. 3A). In S_vs_JS analysis, biological process and molecular function ontologies showed an equal number of enrichments (23). Remarkably, the enriched biological process in S_vs_JS analysis was completely different from that in C_vs_S analysis. Genes related to oxylipin biosynthetic and metabolic process, and cellular oxidant detoxification were the most significantly enriched by JS treatment. Glutathione transferase activity was the most enriched molecular function (13 genes), followed by FMN binding, NADH dehydrogenase (quinone) activity, NADH dehydrogenase activity, and oxidoreductase activity, acting on NAD(P)H, quinone or similar compound as acceptor (Fig. 3B).

DEGs from the KOG database were concurrently classified into 26 functional classes. The most common category was general function prediction, with 406 and 84 genes in the C_vs_S and S_vs_JS analysis, respectively (Figs. 4A and 4B). Control of the cell cycle, cell division, chromosomal partitioning, and glucose transport and metabolism were all significantly enriched in C_vs_S and S_vs_JS analysis (Figs. 4A and 4B). A total of 13 genes controlling defense mechanisms, and three genes related to the nuclear structure were only enriched in C_vs_S analysis (Fig. 4A).

A KEGG pathway analysis of frequent DEGs was undertaken to better understand their functional relevance in salt response. Results showed that the most significantly enriched pathways were metabolic pathways (737 in C_vs_S and 119 in S_vs_JS) and secondary metabolite biosynthesis (509 in C_vs_S and 79 in S_vs_JS) in both C_vs_S and S_vs_JS comparisons. In addition to the shared enriched-pathways, there were some specific pathways enriched exclusively in C_vs_S or S_vs_JS comparisons. For example, genes involved in starch and sucrose metabolism (106) and ribosome (100) were enriched when comparing S to C, rather than comparing JS to S. Whereas 10 genes related to ascorbate and aldarate metabolism were regulated only by JS treatment (Fig. 5). We then focused on the genes regulated by JS relative to S. DEGs in metabolic pathways and secondary metabolite biosynthesis, the common-enriched pathways, were analyzed in S_vs_JS comparison. As shown in Figs. 6A and 6B, DEGs in the same pathway showed different regulatory patterns, roughly half upregulated and half downregulated. Ascorbate and aldarate are essential antioxidants, and genes related to their metabolism were specifically detected by comparing JS to S. Intriguingly, seven genes in this pathway were upregulated by JA addition to salt treatment (Fig. 6c).

Discovery of alternative splicing events

Following the differential analysis of alternative splicing events, rMATS was used to analyze the categories of alternative splicing events, compute their number, and determine the expression level of each type of alternative splicing event individually for each differential grouping. When detected by junction counts (JC), salt stress resulted in 5,568 alternative splicing (AS) events when compared to control, 95.6% of which were skipped exon (SE). Mutually exclusive exons (MXE) events accounted for 4.4% (Fig. 7A). Among all AS events, 238 events were differentially expressed, with 79.0% SE events (23 upregulated and 165 downregulated) and 21.0% MXE events (21 upregulated and 29 downregulated) (Fig. 7B). Similarly, SE events (95.4%) were predominant when compared JS to S treatment. The rest were MXE events (Fig. 7C). There were 114 differentially expressed AS events caused by JS treatment, with 72 SE and 42 MXE events. Up- and down-regulated MXE events were equal, whereas most SE events were upregulated (Fig. 7D). There was no alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS) and retained intron (RI) events were detected in our test. Results obtained by JCEC were similar to JC.

We further discovered SE events in various essential genes involved in salt stress response, including transcription factors (bZIPs and AP2), Gluthatione S-transferases, cyclin-dependent kinase, auxin efflux carrier family protein (Table S2). JS treatment resulted in some unique differentially expressed SE events, such as genes related to Heat shock proteins (HSP), ubiquitin-conjugating enzyme, ATP-binding cassette, and transcription factors (MYBs and MYC2) (Table S3).

Gene expression validation

qRT-PCR was used to verify the fold change of DEGs. Among the 20 DEGs picked were enzyme genes from the secondary metabolic pathway (MS.gene015472), a new gene (novel.17466), and transcription factors (MS.gene074091, MS.gene94790, MS.gene30384, MS.gene026605, and MS.gene20637). All of the randomly picked genes were identified using PCR. The consistency between the fold change of these DEGs and the RNA-seq results indicates the data’s repeatability (Fig. 8, raw data is shown in Table S4).