Jasmonate induced alternative splicing responses in Arabidopsis.

Jasmonate is an essential phytohormone regulating plant growth, development, and defense. Alternative splicing (AS) in jasmonate ZIM-domain (JAZ) repressors is well-characterized and plays an important role in jasmonate signaling regulation. However, it is unknown whether other genes in the jasmonate signaling pathway are regulated by AS. We explore the potential for AS regulation in three Arabidopsis genotypes (WT, jaz2, jaz7) in response to methyl jasmonate (MeJA) treatment with respect to: (a) differential AS, (b) differential miRNA targeted AS, and (c) AS isoforms with novel functions. AS events identified from transcriptomic data were validated with proteomic data. Protein interaction networks identified two genes, SKIP and ALY4 whose products have both DNA- and RNA-binding affinities, as potential key regulators mediating jasmonate signaling and AS regulation. We observed cases where AS alone, or AS and transcriptional regulation together, can influence gene expression in response to MeJA. Twenty-one genes contain predicted miRNA target sites subjected to AS, which implies that AS is coupled to miRNA regulation. We identified 30 cases where alternatively spliced isoforms may have novel functions. For example, AS of bHLH160 generates an isoform without a basic domain, which may convert it from an activator to a repressor. Our study identified potential key regulators in AS regulation of jasmonate signaling pathway. These findings highlight the importance of AS regulation in the jasmonate signaling pathway, both alone and in collaboration with other regulators. SIGNIFICANCE STATEMENT: By exploring alternative splicing, we demonstrate its regulation in the jasmonate signaling pathway alone or in collaboration with other posttranscriptional regulations such as nonsense and microRNA-mediated decay. A signal transduction network model for alternative splicing in jasmonate signaling pathway was generated, contributing to our understanding for this important, prevalent, but relatively unexplored regulatory mechanism in plants.

INTRODUCTION

Plants have developed inducible defense mechanisms and complex signaling networks for efficient, precise and fast response to ever‐changing environmental stimuli and stresses. Jasmonate is a phytohormone induced by herbivore or pathogen attack, as well as abiotic stresses such as cold, salt, UV light, and ozone (Browse & Howe, 2008; Howe & Jander, 2008; Hu, Jiang, Wang, & Yu, 2013; Katsir, Chung, Koo, & Howe, 2008; Valenzuela et al., 2016). Jasmonate regulates a variety of plant events including photosynthesis (Attaran et al., 2014), root and shoot morphogenesis (Gasperini et al., 2015; Zheng et al., 2017), flowering time (Zhai et al., 2015; Thatcher, Danilevskaya, et al., 2016), stamen development (Qi, Huang, Song, & Xie, 2015a; Song et al., 2011), seed germination (Linkies & Leubner‐Metzger, 2012), and senescence (Jiang, Liang, Yang, & Yu, 2014; Miao & Zentgraf, 2007; Qi, Wang, et al., 2015; Yu et al., 2016). In general, jasmonate inhibits plant growth, promotes growth to defense transition and triggers early reproduction.

Biologically active (+)‐7‐iso‐Jasmonoyl‐L‐isoleucine (JA‐Ile), Jasmonate ZIM‐domain (JAZ) proteins and the F‐box protein of the E3 ubiquitin‐ligase complex Coronatine Insensitive 1 (COI1) form the key three component regulatory module mediating jasmonate signal transduction (Chini et al., 2007; Fonseca et al., 2009; Thines et al., 2007). Upon JA‐Ile accumulation, the three‐component module (JA‐Ile, COI1 and JAZ) is formed, and the JAZs are polyubiquitinated and degraded through the 26S proteasome, which releases suppression of the JAZ‐interacting transcription factors that regulate early‐jasmonate‐responsive genes (Chini et al., 2007; Chini, Gimenez‐Ibanez, Goossens, & Solano, 2016; Thines et al., 2007). Thirteen JAZ genes have been identified in the Arabidopsis genome (Bai, Meng, Huang, Qi, & Chen, 2011; Thireault et al., 2015), and they have both functional redundancy and specificity (Chini et al., 2016). Members of the JAZ family share a conserved N‐terminal domain (Moreno et al., 2013), a TIFY (previously known as ZIM) domain (except JAZ13), and a C‐terminal Jas domain (Chini et al., 2007; Thines et al., 2007; Yan et al., 2007). The TIFY domain is responsible for homo‐ and hetero‐dimerization (Chini, Fonseca, Chico, Fernández‐Calvo, & Solano, 2009) and interaction with NOVEL INTERACTOR OF JAZ (NINJA), which further recruits the co‐repressor TOPLESS (TPL) through its conserved ERF‐associated amphiphilic repression (EAR) domain (Pauwels et al., 2010). However, JAZ5, JAZ6, JAZ7, JAZ8, and JAZ13 could directly recruit TPL by EAR in their N‐terminal domain (Causier, Ashworth, Guo, & Davies, 2012; Shyu et al., 2012; Thireault et al., 2015; Thatcher, Cevik, et al., 2016). The Jas domain is the major domain that is responsible for interaction with transcription factors (Chini et al., 2016), and COI1, except in JAZ7 and JAZ8 where the Jas domain has diverged (Chini et al., 2016; Shyu et al., 2012). Interestingly, almost all the JAZs (except for JAZ1, JAZ7, and JAZ8) share a homologous intron (Jas intron) that divides the Jas domain into a 20‐amino acid N‐terminal motif and a 7‐amino acid C‐terminal motif (X5PY; Chung et al., 2010). AS around this conserved intron can lead to a truncated protein lacking the X5PY motif or lacking the entire Jas domain (Chung et al., 2010; Chung & Howe, 2009; Monte et al., 2019; Moreno et al., 2013; Yan et al., 2007). These truncated proteins retain their ability to interact with transcription factors through the 20‐amino acid N‐terminal motif, or a similar sequence in their N‐terminal domain, but have reduced (or abolished) ability of being recognized by COI1 (Chung & Howe, 2009; Chung et al., 2010; Moreno et al., 2013; Zhang, Ke, et al., 2017). As a result, these AS isoforms avoid degradation and function as dominant repressors in the presence of jasmonate. AS events around the Jas intron are conserved among monocots, eudicots, and bryohpytes, underscoring their functional importance (Chung et al., 2010; Monte et al., 2019).

It was predicted that ~60% of Arabidopsis protein‐coding genes are alternatively spliced (Zhang, Calixto, et al., 2017). However, other than JAZ repressors, little is known about AS regulation in the jasmonate signaling pathway. AS generates multiple transcripts from a single gene by regulated differential selection of splice sites. This process is regulated in part by interactions between cis‐elements (intronic/exonic splicing enhancer/silencer) and trans‐factors (e.g. Serine/Arginine‐rich (SR) splicing activators, heterogeneous nuclear ribonucleoprotein (hnRNP) splicing repressors, spliceosome components, and additional RNA binding factors; Kornblihtt et al., 2013; Meyer, Koester, & Staiger, 2015). Regulation of AS influences the relative proportion of different isoforms for a given gene. Moreover, AS interacts with other regulatory mechanisms such as non‐sense mediated decay (NMD; Kalyna et al., 2012; Kervestin & Jacobson, 2012) and microRNA (miRNA) regulation (Reddy, Marquez, Kalyna, & Barta, 2013). AS frequently results in a pre‐mature termination codon (PTC; Kalyna et al., 2012) that can mark the transcript for degradation through the NMD pathway (Kervestin & Jacobson, 2012). Thus, NMD coupled with AS can function as an important posttranscriptional mechanism to regulate protein levels (Kervestin & Jacobson, 2012). Inclusion/exclusion of the alternative region may introduce/remove miRNA targeting sites within one isoform relative to another (Reddy et al., 2013), thus genes producing AS isoforms may be subject to both AS and miRNA regulation. Finally, if the AS isoform can make a protein product, it may generate the same protein or a different protein depending on whether AS occurs within the untranslated region (UTR) or coding sequence (CDS). Generation of different protein products may result in nonfunctional, partially‐functional, redundantly‐functional, or neo‐functional alterations (Reddy et al., 2013; Staiger & Brown, 2013). Isoforms with partial‐ or neo‐functions are especially interesting as they may lead to different gene functions, as is the case for JAZs (Chung et al., 2010; Chung & Howe, 2009; Moreno et al., 2013; Yan et al., 2007).

Here we integrated transcriptomics and proteomics of WT, jaz2 and jaz7 mutant Arabidopsis in response to methyl jasmonate (MeJA) to identify three aspects of AS‐dependent regulation potentially impacting the jasmonate signaling pathway: (a) differential AS in response to MeJA treatment; (b) miRNA regulation edited by differential use of AS isoforms; and (c) AS splice variants with novel functions. The jaz2 and jaz7 mutants are within genes that are in different subgroups of JAZ family (Chung et al., 2010). In addition, jaz2 is a knock‐down mutant and jaz7 is an over‐expression mutant, both of which exhibit obvious phenotypic changes compared to WT (Yan et al., 2014). RNA‐Seq analysis of the jaz2 and jaz7 mutants have been described previously (Yan et al., 2014) and the underlying sequence data are available in the NCBI short‐read archive (SRP026541). We re‐analyzed these data to discover AS isoforms and examine these for differential expression under MeJA treatment—many of which we further explored and describe herein. The AS events identified from the RNA‐seq data were further validated with previously unpublished proteomics data that was generated from the same samples used by Yan et al. (2014) for the RNA‐Seq data acquisition.

EXPERIMENTAL PROCEDURES

Experimental design and transcriptome/proteome library preparation and sequencing

Thirty‐six RNA‐seq samples corresponding to three replicates each from shoots and roots obtained from the Arabidopsis thaliana mutants jaz2 (SALK079895C) and jaz7 (SALK040835C) and WT (Col‐0) with (10 μM) or without (0 μM) MeJA treatment were used for transcriptome and proteome profiling. The transcriptome data were published (Yan et al., 2014) and available in the NCBI Sequence Read Archive under SRP026541. The proteome data that were newly generated have been deposited in the ProteomeXchange MassIVE partner repository with the data set identifiers PXD014563 and MSV000084071.

Transcriptome assembly and differential AS detection

Raw single end reads (75 bp) obtained from an Illumina Hiseq2000 instrument were filtered with FASTX‐Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/index.html) and mapped to the Arabidopsis reference sequence and annotation (ver. TAIR10 http://www.arabidopsis.org) with GSNAP v2013‐07‐20 (Wu & Nacu, 2010) allowing a maximum intron size of 8,000, a minimum intron size of 20, and a maximum of 5% mismatch. The uniquely mapped reads from each library were assembled into transcripts with Cufflinks v2.2.1 (Trapnell et al., 2013) using the reference‐guided method and the same intron/exon size constraints used by GSNAP (above). The 36 assemblies were merged to a single transcript reference set by Cuffmerge v2.2.1 (Trapnell et al., 2013) with the minimum isoform abundance parameter set to five and further filtered with two criteria: (a) each junction of the transcript is supported by>= 3 reads; (b) The minimum average mapping depth within the alternative region of intron retention events was required to be four if 100% covered, five if 90% covered, and six if 80% covered. The AS event types were classified by AStalavista v3.2 (Foissac & Sammeth, 2007). The percentage of novel junction reads (“complete novel” indicates both the 5’ and 3’ splicing sites are novel, “partial novel” indicates only one splice site is novel) were calculated by RSeQC v2.6.2 (Wang, Wang, & Li, 2012). Differential AS and expression were identified by Cuffdiff v2.2.1 (Trapnell et al., 2013) with default parameters. Supplementary‐dataset1‐differential‐expression‐splicing.xlxs contains the cuffdiff results of all significant differential genes (expression and splicing) for all tests. Heatmaps illustrating the expression of transcription factors and splicing factors were produced using data from cuffdiff results. The ratio of FPKM values from two treatments was used as input for the heatmaps, which were generated using the function heatmap.2 from the gplots R package. Supplementary‐dataset2‐heatmapData contains the data used to generate the heatmaps.

Open reading frame prediction

TransDecoder v3.0.0 (https://transdecoder.github.io/) was used to integrate BLAST and Pfam searches for protein prediction and ORF searches on the filtered transcript assemblies.

Protein interaction network analysis

Genes with significant differential expression and/or AS in response to MeJA treatment were used for protein interaction network construction using STRING v10.0 (Szklarczyk et al., 2015). Genes or their homologs with experimental evidence of protein interaction were connected with lines. Protein interaction networks were displayed using Cytoscape v3.4.0 (Shannon et al., 2003) and annotated with gene expression heatmaps.

miRNA target predication

The collection of 427 Arabidopsis miRNAs from the miRNA database (miRBase Release 21, Kozomara & Griffiths‐Jones, 2014) were used as queries to psRNATarget (Dai & Zhao, 2011) to search for potential target sites within the filtered transcript assemblies. Default parameters in scoring schema V2 were used.

Protein extraction, digestion, iTRAQ labeling, and LC‐MS/MS

Proteins were quantified as previously described (Koh et al., 2012), and dissolved in denaturant buffer (0.1% SDS [w/v]) and dissolution buffer (0.5 M triethylammonium bicarbonate, pH 8.5) in the iTRAQ 8‐plex kit (AB sciex Inc.). For each sample, a total of 100 μg of protein were reduced, alkylated, trypsin‐digested, and labeled according to the manufacturer’s instructions (AB Sciex Inc.). The control samples from wild type, jaz2, and jaz7 were labeled with iTRAQ tags 113, 114, and 115, respectively, and the corresponding treated samples were labeled with iTRAQ tags 116, 117, and 118 respectively. In addition, all six samples were mixed and labeled with iTRAQ tag 119 as an internal control. Labeled peptides were desalted with C18‐solid phase extraction and dissolved in strong cation exchange (SCX) solvent A (25% (v/v) acetonitrile, 10 mM ammonium formate, and 0.1% (v/v) formic acid, pH 2.8). The peptides were fractionated using an Agilent HPLC 1260 with a polysulfoethyl A column (2.1 × 100 mm, 5 µm, 300 Å; PolyLC). Peptides were eluted with a linear gradient of 0–20% solvent B (25% (v/v) acetonitrile and 500 mM ammonium formate, pH 6.8) over 50 min followed by ramping up to 100% solvent B in 5 min. The absorbance at 280 nm was monitored and a total of 12 fractions were collected. The fractions were lyophilized and resuspended in LC solvent A (0.1% formic acid in 97% water (v/v), 3% acetonitrile [v/v]). A hybrid quadrupole Orbitrap (Q Exactive Plus) MS system (Thermo Fisher Scientific) was used with high energy collision dissociation (HCD) in each MS and MS/MS cycle. The MS system was interfaced with an automated Easy‐nLC 1000 system (Thermo Fisher Scientific). Each sample fraction was loaded onto an Acclaim Pepmap 100 pre‐column (20 mm × 75 μm; 3 μm‐C18) and separated on a PepMap RSLC analytical column (250 mm × 75 μm; 2 μm‐C18) at a flow rate at 350 nl/min during a linear gradient from solvent A (0.1% formic acid [v/v]) to 30% solvent B (0.1% formic acid (v/v) and 99.9% acetonitrile [v/v]) for 95 min, to 98% solvent B for 15 min, and hold 98% solvent B for additional 30 min. Full MS scans were acquired in the Orbitrap mass analyzer over m/z 400–2000 range with resolution 70,000 at 200 m/z. The top ten most intense peaks with charge state ≥2 were isolated (with 2 m/z isolation window) and fragmented in the high energy collision cell using a normalized collision energy of 28%. The maximum ion injection time for the survey scan and the MS/MS scans were 250 ms, and the ion target values were set to 3e6 and 1e6 respectively. The selected sequenced ions were dynamically excluded for 60 s.

Proteomics data analysis

The raw MS/MS data files were processed by a thorough database searching approach considering biological modification and amino acid substitution against customized Arabidopsis database using the ProteinPilot v4.5 with the Fraglet and Taglet searches under ParagonTM algorithm (Shilov et al., 2007). The following parameters were considered for all the searching: fixed modification of methylmethane thiosulfonate‐labeled cysteine, fixed iTRAQ modification of amine groups in the N‐terminus, lysine, and variable iTRAQ modifications of tyrosine. The false discovery rate at the peptide level was estimated with the integrated PSPEP tool in the ProteinPilot Software to be 1.0%. The identified peptide reads were screened for confidence no less than 95%. The screened peptide reads were mapped to TAIR10 primary protein database and peptides which failed to map to the database were regarded as candidates for supporting the AS isoforms. These candidates were manually validated and only reads spanning the AS junction were regarded as evidence for that AS isoform.

RESULTS

Transcriptomic analysis and genome‐guided assembly

452 million reads from 36 RNA‐seq Arabidopsis libraries (three each from WT, JAZ2 and JAZ7 roots or shoots w/wo 10 μM methyl jasmonate) were uniquely mapped to the Arabidopsis TAIR10 genome assembly (Table S1). These assembled into 20,524 transcripts from 13,647 genes. 15,947 transcripts were previously known (TAIR10) while the remaining 4,577 were novel (Figure 1a). Among all the splice junctions identified in this project, 26% of the junctions have at least one boundary novel to the TAIR10 annotation (Figure 1b). A total of 4446 genes (32.58% of the identified 13,647 genes) demonstrate AS with the majority (98%) of them producing two to five splice variants (Table 1). The most abundant AS events are alternative acceptor (38%) and intron retention (30%; Figure 1c).

FIGURE 1

Characterization of assembled transcripts. (a) Comparison of the assembled transcripts from the MeJA RNA‐seq data and TAIR10 annotation. (b) Comparison of the mapped junction reads from the MeJA RNA‐seq data to the TAIR10 annotated junctions. (c) Identified AS patterns detected in the MeJA RNA‐seq assemblies. (d) ORF prediction of the 20,524 transcripts assembled from the MeJA RNA‐seq data with the TransDecoder program (Haas et al., 2013)

TABLE 1

Gene isoform number in TAIR10 annotation and in RNA‐seq assembly from this project

Isoform no.	TAIR10		MeJA project
	Gene no.	Percentage	Gene no.	Percentage
1	21402	78.67	9201	67.42
2	4251	15.63	2998	21.97
3	1133	4.17	881	6.46
4	291	1.07	345	2.53
5	89	0.33	115	0.84
6	26	0.10	65	0.48
7	7	0.03	21	0.15
8	5	0.02	11	0.08
9	1	0.00	6	0.04
10	1	0.00	0	0.00
11	0	0.00	1	0.01
12	0	0.00	2	0.02
13	0	0.00	0	0.00
14	0	0.00	0	0.00
15	0	0.00	1	0.01
Total	27206	100.00	13647	100.00

Jasmonate‐Related Protein Interaction Network

Significant differential gene expression and differential AS were identified with Cuffdiff (Trapnell et al., 2013) in response to MeJA treatment, between tissues, and between genotypes (Table 2). In each case there were always more genes with significant expression changes than genes with significant AS changes. The largest number of differences was observed between tissues. Genes with significant differential expression and/or AS in response to jasmonate treatment were used to generate a protein interaction network with STRING (Szklarczyk et al., 2015). A subset of the network including splicing related proteins, jasmonate key regulatory proteins (JAZs, COI1, and NINJA) as well as genes interacting with them were identified and examined further (Figure 2). We arbitrarily divided the network into four modules. One module contains splicing related genes, and this is connected to three other modules (named Module 1, 2, 3; Figure 2). Module 1 contains the key jasmonate regulatory factors, JAZ3, JAZ10, TIFY7, COI1, and NINJA, as well as transcription factors such as bHLHs and R2R3‐MYBs. Based on the network, the communication between jasmonate signaling and splicing signaling was shown to be mediated by bHLHs, R2R3‐MYBs and three splicing‐related proteins: CBF1‐interacting co‐repressor (CIR), pre‐mRNA‐splicing factor ISY1‐like protein (LSY1) and SNW/SKI‐interacting protein (SKIP). Module 2 mainly contains kinases and transcription factors, and its interaction with the splicing‐related proteins is through kinases and SKIP. Module 3 is centered by a topoisomerase (TOPII) and a ubiquitin (UBQ11) protein. Interaction between the third module and the splicing‐related proteins is through Embryo Defective 2816 (emb2816) and Always Early 4 (ALY4).

TABLE 2

Differential AS or expression of genes undergoing AS in treatment comparisons, tissue comparisons and between mutant backgrounds

Genotype	Tissue	Differential alternative splicing				Differential expression
		Splicing diff.	Promoter diff.	CDS diff.	Total
Untreated vs. 10 μM MeJA
WT	shoot	44	3	6	48	537
jaz2	shoot	32	5	8	37	586
jaz7	shoot	234	90	72	326	1400
WT	root	49	25	15	75	886
jaz2	root	32	14	9	48	557
jaz7	root	43	22	13	67	1131

FIGURE 2

Protein interaction network of genes that undergo AS which show differential expression or differential AS in response to MeJA treatment. Edges indicate experimentally validated interactions of the two connected proteins or their homologs. Yellow boxes indicate proteins related with splicing; blue boxes indicate transcription factors; pink boxes indicate kinase proteins. Heatmaps under the gene name indicate gene expression patterns in response to MeJA treatment

Regulation of Transcription Factors (bHLHs and MYBs) and Splicing Factors (SRs and hnRNPs)

The protein interaction network indicates the inclusion of several members of the bHLH and MYB gene families, suggesting their importance in the jasmonate pathway. Moreover, the SR and hnRNP gene families are key trans‐factors involved in the selection of splice junctions and play essential roles in AS regulation (Barbazuk, Fu, McGinnis, 2008). We further analysed the expression patterns of these four important gene families in response to MeJA treatment as well as between tissues and genotypes (Figure 3). Genes without any significant expression changes in any of the three comparisons were not included.

FIGURE 3

Heatmap of differentially expressed transcription factor (bHLH and MYB) and splicing factor (SR and hnRNP) gene family members that undergo AS. Differentially expressed bHLH and MYB genes under MeJA treatment (a), between tissues (shoot relative to root; b), and between mutant backgrounds (mutant relative to WT; c). Differentially expressed SR and hnRNP genes under MeJA treatment (d), between shoots and roots; e), and between genotypes (mutants relative to WT; f)

The expression profiles of the MYB and bHLH gene family reveal the impacts of expression in response to MeJA treatment (up or down), as well as whether they act downstream of JAZ2 and JAZ7. For example, while the expression of bHLH116 is strongly downregulated in WT root, it is slightly upregulated in jaz2 root, and higher still in jaz7 root in response to MeJA (Figure 3a). This expression pattern indicates that bHLH116 might be regulated by JAZ2 and JAZ7 in the root and is also supported by the expression change between genotypes (Figure 3c). Similarly, bHLH77 might be regulated by JAZ2 and JAZ7 in the root tissue; bHLH153, bHLH155, and bHLH137 might be regulated by JAZ2 and JAZ7 in the shoot tissue; MYB48, and bHLH36 are likely regulated by JAZ7 in both the shoot and root tissues (Figure 3a).

The majority of the splicing factors were downregulated in response to MeJA treatment (Figure 3d). The expression changes were less dramatic in the SR and hnRNP splicing factors compared with the bHLH and MYB transcription factors. However, we observed several splicing factors that may be involved in jasmonate responses: RBP45a, RS2Z33, RNPA/B_3, and RSZ21. The expression of RBP45a was downregulated in response to MeJA in the WT shoots but exhibits increased expression in the jaz7 mutant and even higher expression in the jaz2 mutant. Another interesting case is RS2Z33, which is dramatically downregulated in response to MeJA in jaz7 mutant root tissue but upregulated in the jaz2 mutant and WT in both shoot and root tissues (Figure 3d). Moreover, the expression of RS2Z33 in shoot tissue was greatly induced in jaz2 mutants with respect to WT in the shoot tissue (Figure 3f). The expression of RNPA/B_3 did not change much in response to MeJA in WT shoot. However, its expression is greatly downregulated in the shoot of jaz2 and greatly upregulated in the shoot of jaz7. These altered expression patterns suggest possible involvement of these splicing factors in the jasmonate pathway regulated by JAZ2 and JAZ7.

Differential alternative splicing in response to MeJA treatment

The highest number of significantly differential AS events in response to MeJA treatment was observed in shoot from jaz7 mutants (Table 2). Across each of the six treatments (two tissues sampled from each of three genotypes), only a fraction of the differential AS events are shared in response to MeJA (Figure S1). Genes with significant differential AS in response to MeJA in more than one genotype of a specific tissue (shoot or root) were analyzed further (Figure S2). Among the 16 genes, seven genes have AS events in the coding region, while five have events that could lead to a PTC that could target these transcripts for NMD (Kervestin & Jacobson, 2012). Based on the AS pattern and expression profiles, we observed cases where AS regulation alone may regulate gene expression (through NMD), such as NUDIX HYDROLASE HOMOLOG 9 (NUDX9; Figure 4a), and cases where AS and transcriptional regulation together regulate gene expression, such as NITRATE TRANSPORTER 1.8 (NRT1.8; Figure 4d).

FIGURE 4

Two genes under AS regulation in response to MeJA treatment. (a) Expression profile of NUDX9. Error bars indicate standard deviation. Pie charts under each condition indicate the proportion of each AS isoform relative to the total expression from the locus. Red lines under pie charts indicate significantly differential AS in the shoot tissue. (b) Gene structure of the two isoforms of NUDX9. Angled lines indicate introns; thin black boxes indicate UTRs; blue boxes indicate CDS; pink boxes indicate regions which were converted to noncoding regions as a result of an AS induced PTC. Gray line connected red boxes indicate mapped reads supporting the novel junction in the alternative isoform. (c) Regulation of NUDX9 in response to MeJA treatment by AS. (d) Expression profile of NRT1.8. Blue lines under pie charts indicate significantly differential AS patterns in the root tissue. Red and blue lines above the barplot indicate significantly differential gene expression in shoot and root respectively. (e) Gene structure of the two isoforms of NRT1.8. (f) Regulation of NRT1.8 in response to MeJA treatment by transcription and AS

Two splice variants were identified for gene NUDX9 with an alternative acceptor event in the fifth exon leading to a frame shift that results in a PTC in one of the AS transcripts. There are 17 mapped reads supporting the AS junction (Figure 4b). While no significant changes in the gene expression level were identified in response to MeJA, we observed significant differential AS in the shoot tissue of WT and jaz7 upon MeJA treatment (Figure 4a). In shoot of WT and jaz7, When there is no MeJA treatment the two splice variants are present at similar levels, with one isoform (the primary isoform annotated by TAIR10) generating a protein product, while the alternative acceptor AS isoform is potentially subject to NMD. However, upon MeJA treatment almost all isoforms are the primary isoform providing the potential to generate twice as much productive translation of the protein product under MeJA treatment relative to untreated tissue with no required increase in the rate of transcription (Figure 4c).

NRT1.8 exhibited both differential expression and AS in response to MeJA treatment in the root of the three genotypes (Figure 4d). Two splice variants were identified for NRT1.8. The AS isoform contains an intron retention and an alternative polyadenylation event compared with the primary annotated isoform. The removal of the retained intron in the AS isoform, which is supported by 10 mapped reads (Figure 4e), leads to PTC and thus the AS isoform is unproductive. Upon MeJA treatment, the expression level of NRT1.8 was up‐regulated and AS generated more of the productive isoform. In the case of NRT1.8, differential transcription and AS both play a role to potentially generate more protein products in response to increased jasmonate (Figure 4f).

Alternative splicing variants differentially targeted by miRNA

We used psRNATarget (Dai & Zhao, 2011) to predict miRNA targets on the filtered transcript assemblies. A total of 508 genes were predicted to be targeted by miRNAs, among which 171 have evidence of AS (Figure S3). Sixty‐four genes contain miRNA target sites subjected to AS regulation, suggesting the potential of being regulated by a combination of AS and miRNA. Among them, 21 showed significantly changed isoform proportions (significantly differential AS identified by Cuffdiff) in MeJA treatment/tissue/genotype comparisons (Figure S4). Most of the changes in isoform ratio occur between tissues, with only five genes—SCHLAFMUTZE (SMZ), ALDEHYDE OXIDASE 2 (AAO2), PLANT U‐BOX 4 (PUB4), LON PROTEASE 2 (LON2), At3g02740—exhibiting differences in relation to MeJA treatment and three (AAO2, LON2, and CEST) exhibiting differences between genotypes (Figure S4).

Alternative splicing variants with novel functions

In order to identify AS isoforms with novel functions, the ORF and protein products of the 20,524 transcripts were predicted with TransDecoder (Figure 1d; https://transdecoder.github.io/). For transcripts shared with TAIR10 database, 99.6% predicted products from TransDecoder match the TAIR10 annotation, supporting the efficacy of TransDecoder. We used the ORF annotation from TAIR10 for the 15,947 known transcripts and the ORF prediction from TransDecoder for the remaining 4,577 novel transcripts. Transcripts with no ORF prediction or multiple ORF predictions were excluded from further analysis, thus 3,220 genes with evidence of undergoing AS (a total of 7,638 transcripts) with ORF predictions from either TAIR10 or TransDecoder were used for further analysis.

Not all AS leads to multiple protein products of a gene. For example, AS events that affect UTR regions will not generate a novel protein, and those that result in a PTC leading to NMD are not translated. In this study, 37.6% of the AS events occurred in UTRs (Figure S5). Long 3’‐UTRs and introns in the UTR region are cis‐elements that trigger NMD (Kalyna et al., 2012; Kervestin & Jacobson, 2012). Two experimentally validated criteria for NMD in Arabidopsis were applied to identify potential cases of NMD (Kalyna et al., 2012): (a) 3’ UTR length >350 nt; and (b) the distance of stop codon to last exon junction >55 nt (Figure S5). As a result, a total of 1,745 transcripts from 1,106 genes (34.3% genes with evidence of AS) were predicted to be targets of NMD. After eliminating genes undergoing AS but producing only a single protein product, there were 1,464 genes with AS transcripts that can potentially generate multiple protein products.

Among these 1,464 genes, 171 were identified to be splicing‐related genes, jasmonate‐related genes, or transcription factors, and also show significant differential expression and/or AS in response to jasmonate treatment/tissues/genotypes comparison. Predicted proteins from the 372 transcripts of these 171 genes were analyzed with the NCBI Conserved Domain Database (CDD, Marchler‐Bauer et al., 2015) and Simple Modular Architecture Research Tool (SMART, Letunic, Doerks, & Bork, 2015) to determine their domain structure. We identified potential novel AS isoforms by two criteria: (a) the AS protein has conserved domain(s) which suggests functional importance, and (b) the domain arrangement of the AS variant is different from the primary protein, indicating functional divergence. A total of 30 genes were identified that satisfy both criteria (Figure 5). We observed conservation between members of the same gene family for domain pattern changes. For example, the MADS‐box genes FOREVER YOUNG FLOWER (FYF) and FLOWERING LOCUS M (FLM) both have AS isoforms lacking the MADS domain; R2R3‐MYB genes MYB59, MYB48, MYB28, MYB15, and MYB47 all have AS isoforms lacking the R2 domain; 3R‐MYB genes MYB3R1 and MYB3R4 both have AS isoforms lacking the C‐terminal repression motif (Figure 5).

FIGURE 5

AS Genes with different domain structures predicted between transcript isoforms. Heatmaps indicate up‐ or down‐regulated transcript expression in response to MeJA treatment in six conditions (SW: shoot of WT; S2: shoot of jaz2; S7: shoot of jaz7; RW: root of WT; R2: root of jaz2; R7: root of jaz7)

Alternative splicing variant of bHLH160 with a potential novel function

Among the 30 genes that were predicted to produce AS isoform(s) with novel functions, we identified bHLH160 as an interesting case. Like all bHLH proteins, the bHLH160 contains a basic‐Helix‐Loop‐Helix domain, where the basic region is responsible for DNA binding, and the HLH region is responsible for dimerization (Figure 6a, Carretero‐Paulet et al., 2010). When the bHLH dimer stably binds with the DNA recognition sequence, it serves as a transcriptional activator/repressor. Three transcript isoforms of bHLH160 were assembled with two of them (TCONS_00034502 and TCONS_00034498) predicted to generate the primary protein product and the third isoform (TCONS_00034499) predicted to produce a product lacking part of the N‐terminal region, including 13bp basic region and its upstream sequences (Figure 6a). There are 45 mapped reads supporting the AS junction of the third isoform (Figure S6). Given the function of the basic and HLH regions, a bHLH protein that lacks the basic region is able to dimerize but unable to bind DNA. In addition, by dimerizing with normal bHLH proteins, bHLH160b‐ is expected to decrease the active bHLH dimers (Figure 6b). Thus bHLH160b‐ functions in a manner opposite to that of the primary isoform. We observed significantly upregulated gene expression of bHLH160 in response to MeJA treatment in WT and jaz2 mutant roots (Figure 6c).

FIGURE 6

Arabidopsis bHLH160 AS patterns and the proposed regulatory function of the splice variant. (a) Gene structure of the transcript isoforms of Arabidopsis bHLH160. (b) Protein 3D structure of the Arabidopsis bHLH160 proteins adapted from Ma et al., 1994 and their regulatory roles. (c) Expression profile of Arabidopsis bHLH160. Error bars indicate the standard deviation. Pies charts under each condition indicate the proportion of each AS isoform present relative to the total expression from the locus. Red lines above the barplot indicate significantly changed gene expression. (d) Gene structure of the transcript isoforms of Camelina sativa LOC104712692, LOC104701783, and LOC104750970. Angled lines indicate introns; black thin boxes indicate UTRs; thick boxes indicate CDS; orange, green, and purple indicate basic, helix, and loop, respectively, in the bHLH domain. (e) Multiple sequence alignments of the protein isoforms from Arabidopsis thaliana bHLH160, Camelina sativa LOC104712692, LOC104701783, and LOC104750970. The bHLH domain is indicated within black boxes

To investigate whether the AS regulatory pattern in bHLH160 is conserved among other bHLH genes, all other bHLH genes expressed in these data with evidence of AS were examined and no similar AS pattern was observed. Phylogenetic analysis of the bHLH gene family indicates bHLH160 is a newly evolved gene in Brassicaceae (Carretero‐Paulet et al., 2010). We further checked the bHLH160 orthologs in other Brassicaceae species and three genes with similar AS patterns based on the gene annotation (NCBI database) of Camelina sativa were identified (Figure 6d). Protein sequence alignments of the primary and AS isoforms of the Arabidopsis bHLH160 and the three Camelina genes (LOC104712692, LOC104701783, LOC104750971) suggest a conservative alternative start codon at the first amino acid of Helix1 (Figure 6e). A similar AS pattern observed in multiple species suggests that AS generation of both an activator and a repressor might be a conserved and important regulatory mechanism in bHLH160 and its orthologs.

Proteomics validation of alternatively spliced isoforms

An attempt was made to validate AS events identified at the transcript level with proteomics data. We applied three criteria to select peptide sequences supporting identified AS events: (a) The peptide sequence should have at least a 95% confidence level; (b) The peptide sequence is uniquely mapped to a single gene; (c) The peptide sequence is mapped to the AS junction of the AS isoform. With the above criteria, AS events of nine genes were identified to have evidence supporting the protein level expression (Figure 7). Of the nine AS events, two (PAPP5, AKINBETA1) are exon skipping; two (RCA, GRXC2) are intron retention; two (SYP43, GSTZ1) are alternative acceptor, two (MORF1, PYR6) are alternative polyadenylation, and one (RPAC14) is alternative promoter (Figure 7). Among them, GSTZ1 has support for both the AS junction and the primary annotated splicing junction. Two (MORF1 and RPAC14) of the nine genes showed differential AS in response to MeJA and/or genotype comparison (Figure S7).

FIGURE 7

Proteomics validated AS isoform expression. Gene structure was displayed by Gene Structure Display Sever 2.0 (http://gsds.cbi.pku.edu.cn, Hu et al., 2013): thick black line—UTR; blue box—exon; angled line—intron; red angled line—intron splitting the junction of the transcript isoform supported by proteomics data; black/gray line under blue box—peptide mapped region. Shaded isoform name indicates primary annotated isoform from TAIR10 annotation

DISCUSSION

Alternative splicing profiles and networks

Focusing on jasmonate signaling core regulators and splicing related genes, a protein network was identified (Figure 2). The bHLH and MYB transcription factors likely mediate transmission of the jasmonate signal to the splicing‐related genes through splicing factors CIR (At2g44200), LSY1 (At3g18790), and SKIP (Ag1g77180; Figure 2). SKIP is an SNW domain‐containing protein, which has conserved functions as both a transcriptional co‐regulator (Li et al., 2016; Lim et al., 2009) and a splicing factor (Cui et al., 2017; Feng et al., 2015; Li et al., 2016; Wang, Wu, et al., 2012). In Arabidopsis, SKIP is a key regulator of AS across a broad collection of regulated processes, such as abiotic stress response (Feng et al., 2015; Lim et al., 2009), circadian clock (Wang, Wu, et al., 2012), light signaling (Zhang et al., 2014), and flowering time (Cui et al., 2017). SKIP is a spliceosome component that physically interacts with SR family splicing factors, such as SR45 (Wang, Wu, et al., 2012), and a variety of its target pre‐mRNAs (Cui et al., 2017; Feng et al., 2015; Wang, Wu, et al., 2012). SKIP regulates AS by controlling recognition and cleavage of the 3’ and 5’ splicing sites (Feng et al., 2015) and the skip‐1 mutant in Arabidopsis causes genome‐wide splicing defects (Wang, Wu, et al., 2012). In the regulatory network (Figure 2), SKIP is a hub protein connecting with the most proteins including splicing‐related proteins, transcription factors (module 1), and kinases (module 2). Considering its function and placement in the network, perhaps SKIP may serve as a key regulator mediating jasmonate signaling through alternative splicing regulation (Figure 8). The third module is centered by a polyubiquitin UBQ11 (At4g05050) and a topoisomerase TOPII (At3g23890), which are connected to the splicing‐related proteins mainly through two genes: emb2816 (At2g03870) and ALY4 (At5g37720; Figure 2). emb2816 is a SM‐like gene that encodes a component for the U6 snRNP (Perez‐Santángelo et al., 2014). Homologs of emb2816 in Arabidopsis and the human ortholog of emb2816 regulate circadian rhythms by regulating RNA splicing (Perez‐Santángelo et al., 2014). Jasmonate is an important phytohormone involved in plant circadian regulation that demonstrates peak expression at midday in preparation for potential insect feeding, which is most probable at that time (Lu, McChung, & Zhang, 2017). The presence of emb2816 in the protein network may be indicative of a possible regulatory pathway for jasmonate signaling in circadian rhythm regulation although further experiment validation is needed. ALY4 is a key regulator involved in nucleo‐cytosolic mRNA transport (Pfaff et al., 2018) and transcriptional regulation (Storozhenko, Inzé, Montagu, & Kushnir, 2001). In Arabidopsis, mutation in ALY genes leads to multiple vegetative and reproductive defects (Pfaff et al., 2018). Based on its function, ALY4 could potentially regulate TOPII and UBQ11 through transcriptional regulation and/or mRNA transportation (Figure 2). It is interesting that SKIP and ALY4, which are key regulators connecting splicing signaling with other genes, have both RNA‐binding and DNA‐binding affinity. The broad binding affinity may confer these genes with greater impact on the gene cascade by serving as a hub which can integrate multiple signals.

FIGURE 8

Model of alternative splicing regulation in the jasmonate signaling pathway mediated by JAZ2 and JAZ7. JA‐Ile (JA) negatively regulates JAZ2 by SCF COI1 complex and JAZ7 by other unknown pathways. JAZ2/7 represses the expression of transcription factor MYBs and bHLHs. These two types of transcription factors further regulate downstream pathways mediated by transcription factors, kinases, or ubiquitin. SKIP could interact with the transcription factor and kinases mediated pathways and ALY4 could interact with the ubiquitin mediated pathways. The SKIP and ALY4 are directly or indirectly involved in splicing regulation, which regulates alternative splicing of targeted genes (e.g. NUDX9, NRT1.8—Figure 4; bHLH16—Figure 6) by splicing factors SR or hnRNP family proteins

Taken together, the three modules interacting with the splicing‐related proteins suggest three possible signal transduction pathways between jasmonate and mRNA splicing regulation: transcription factors (module 1), kinases (module 2), and ubiquitin pathways (module 3; Figures 2 and 8). The first two modules are likely mediated by a key regulatory SKIP and the third module by another key regulator ALY4. Further functional validation of the two proteins will be important for our understanding of AS in the jasmonate signaling pathway.

Alternative splicing coupled miRNA regulation

In this study, we identified 64 genes which contain miRNA binding sites potentially regulated by AS, including the known SPL4 (Figure S8; Wu & Poethig, 2006). In Arabidopsis, alternative polyadenylation on the 3’UTR of SPL3/4/5 affects the miR156 binding site and generates miR156‐sensitive and miR156‐insensitive isoforms. Increased abundance of miR156‐insensitive isoforms causes juvenile‐to‐adult transition (Wu & Poethig, 2006). Although jasmonate is known to be an important negative regulator in juvenile‐to‐adult transition (Hibara et al., 2016), we did not observe a significantly changed level of the miR156‐insensitive form of SPL4 (TCONS_00010799) in response to jasmonate treatment (Figure S8). Whether AS‐coupled miRNA regulation in SPL4 is regulated by jasmonate remains to be explored.

We identified a small proportion of genes (32 genes, 6.3% of miRNA regulated genes) that were predicted to contain multiple miRNA target sites. Most plant miRNA regulated genes were regulated by only one miRNA with a single target site (Axtell, Westholm, & Lai, 2011; Voinnet, 2009). Plant genes containing multiple miRNA binding sites were rarely reported (Campo et al., 2013). In rice, the dominant splicing variant of Nramp6 (Natural resistance‐associated macrophage protein 6, Os01g31870), Nramp6.8, contains two miR7695 binding sites in its 3’UTR region. The expression level of Nramp6.8 is negatively correlated with miR7695, and the over‐expression of miR7695 confers pathogen resistance in rice (Campo et al., 2013). Gene AAO2 identified in this work showed significantly changed isoform proportion in response to MeJA treatment in the shoot of jaz7, and it was predicted to contain two miRNA target sites recognized by different miRNAs, miR414 and miR5021 (Figure S3). The miRNA binding sites are predicted to be on the 3’UTR of one of the two transcript isoforms. AAO2 is one of the four Arabidopsis aldehyde oxidase genes (AAO1, AAO2, AAO3, and AAO4) potentially involved in the last step of biosynthesis of indole‐3‐acetic acid (IAA) and/or abscisic acid (ABA; Akaba et al., 1999; Seo et al., 2004). One possible explanation is that AAO2 is under AS and miRNA regulation and may be involved in the interactions among jasmonate, IAA, and ABA signals.

In animals, AS‐associated miRNA regulation is known to play important roles in a wide variety of pathways, such as immune cell activation, cell proliferation, muscle stem cell function, and heart development (Boutet et al., 2012; Kalsotra, Wang, Li, & Cooper, 2010; Sandberg, Neilson, Sarma, Sharp, & Burge, 2008). In Arabidopsis, more than 12.4% of miRNA binding sites are predicted to be regulated by AS (Yang, Zhang, & Li, 2012). However, evidence identified for AS‐coupled miRNA regulation in response to MeJA treatment in this study is limited (Figure S4). Moreover, many of the predicted miRNA target sites were located at the long 3’UTRs of the AS isoforms (e.g. AAO2, GRF3, STR16, Figure S3). These long 3’UTRs themselves can trigger NMD (Kalyna et al., 2012). It is possible that miRNAs may be important regulators, working alongside the NMD pathway, to reduce the abundance of transcripts with long 3’UTRs. It has been shown that miRNAs target introns in Arabidopsis and rice (Meng, Shao, Ma, & Wang, 2013). Taken together, these results suggest that miRNA may specifically repress or degrade erroneous mRNA with intron or long 3’UTR generated by AS.

Functional AS regulation

Several cases were identified where AS may play functional roles in response to MeJA. In terms of AS regulation—determined by a changed proportion of splice variants under different conditions—we identified genes NUDX9 and NRT1.8 (Figure 4). Both genes exhibited significantly altered isoform proportions upon increased jasmonate, suggesting existence of cis‐elements (intronic/exonic splicing enhancer/silencer) and different involvement of trans‐factors (e.g. SR splicing activators and hnRNPs splicing repressors). Specifically, we observed significantly changed expression of SR genes RS2Z33, RSZ21, and hnRNP genes RBP45a, RNPA/B_3 in response to MeJA treatment (Figure 3). These splicing factors are candidate genes responsible for differential AS induced by jasmonate. Moreover, we explored a protein interaction network of jasmonate‐responsive genes (Figure 2). In this scenario, increased jasmonate would trigger the jasmonate signaling pathway. The signal could then be mediated and transferred via transcription factors, kinases and ubiquitin signaling pathways and ultimately transmitted to splicing‐related factors that would regulate AS of the responsive genes, such as NUDX9 and NRT1.8 (Figures 4 and 8). Interestingly, the AS isoforms of both genes are potentially subjected to NMD due to a PTC. NMD‐coupled AS is an efficient mechanism to downregulate gene expression (Kalyna et al., 2012; Kervestin & Jacobson, 2012) and it was predicted to occur in 34.3% of the genes demonstrating AS in this analysis (Figure S5). These AS isoforms would escape proteomics validation due to lack of protein products. NRT1.8 is a nitrate transporter which enhances nitrate uptake by mediating nitrate unloading from the xylem vessels (Li et al., 2010). Specifically, NRT1.8 not only shows differential AS but also increased transcription in response to MeJA treatment. The two regulation points (transcription and AS) function synergistically to increase NRT1.8 protein product (Figure 4). Increased NRT1.8 protein level indicates its possible involvement in jasmonate‐triggered plant nitrate uptake. NUDX9 is a GDP‐D‐Man pyrophosphohydrolase, which indirectly modulates ammonium responses by hydrolysis of GDP‐D‐Man in the root (Tanaka et al., 2015). Here, we reported a possible involvement of NUDX9 in the jasmonate signaling pathway by AS regulation in the shoots of Arabidopsis.

In addition to regulating the relative abundance of isoforms, AS also plays important roles by generating functionally novel splice variants, as is the case for the well documented JAZ repressors (Chung et al., 2010; Chung & Howe, 2009; Moreno et al., 2013; Yan et al., 2007). Notably, in response to wounding, the ratio of the expressed JAZ10 splice variants did not change significantly (Yan et al., 2007). The AS regulation of the JAZ repressors lies in generating a stable repressor (lacking the domain recognized by COI1) in addition to the primary degradable repressor rather than modulating the proportion of isoforms. AS has the potential to generate isoforms with different functions by modulating domain structures (S. R. Thatcher, Danilevskaya, et al., 2016). By identifying AS genes with different domain arrangements between splice variants, an attempt was made to identify cases similar to JAZ repressors in the jasmonate signaling pathway. We identified a copper‐responsive transcription factor bHLH160 (Bernal et al., 2012; Yamasaki, Hayashi, Fukazawa, Kobayashi, & Shikanai, 2009), which has the potential to generate an activator and a repressor through AS by inclusion/exclusion of sequences coding for the DNA‐binding domain (Figure 6). The splice variant bHLH160 contains the dimerization domain but lacks the DNA‐binding domain, thus it would repress the primary isoform function by dimerizing with primary bHLH160 transcripts but not able to bind with their cis‐regulatory elements. Human Id genes, which only contain the HLH domain like Arabidopsis bHLH160, function in a similar manner to repress gene expression (Sikder, Devlin, Dunlap, Ryu, & Alani, 2003). Interestingly, similar AS regulation was observed in MADS‐box genes, FYF and FLM. The AS isoforms of FYF and FLM have lost the MADS‐box domain, which is the DNA‐binding domain, but retained the K‐box domain, which facilitates dimerization of MADS‐box proteins (Par̆enicová et al., 2003; Figure 5). Thus, the splice variants of FYF and FLM may negatively regulate the original function of the gene in a manner similar to bHLH160. AS could also result in isoforms that have the same outcome but different strengths, which is the case for MYB3R1 and MYB3R4 (Figure 5). MYB3R1 and MYB3R4 are transcriptional activators which upregulate G2/M transition in cell cycle (Haga et al., 2007). AS of MYB3R1 and MYB3R4 leads to truncated proteins lacking the repression motif, which would serve as hyper‐activators compared with the primary proteins (Feng, Burleigh, Braun, Mei, & Barbazuk, 2017; Kato et al., 2009). In most of the identified cases, the primary protein product of the genes has multiple domains and the AS isoform has lost one or more domain(s; Figure 5). It is not surprising that many of the cases we identified are transcription factors and splicing factors, as they usually contain at least two domains: the DNA/RNA‐binding domain and the regulatory domain. Based on these observations, AS could serve as an important regulator on multi‐domain proteins that could generate splice variants with different functions by altering domain arrangement. Further functional validation of these cases would help to elucidate regulatory roles of different isoforms and their involvement in the jasmonate signaling pathway.

AUTHOR CONTRIBUTIONS

GF, SC, WBB designed the work. GF, MJY, JK, RD generated data. QF, MJY, JK, RD, LB, SC, WBB analyzed and reviewed the data. GF, SC, and WBB wrote the manuscript.

Supplementary Material

Click here for additional data file.

Supplementary Material

Click here for additional data file.

Supplementary Material

Click here for additional data file.