Integrated multi-omics framework of the plant response to jasmonic acid.

Understanding the systems-level actions of transcriptional responses to hormones provides insight into how the genome is reprogrammed in response to environmental stimuli. Here, we investigated the signalling pathway of the hormone jasmonic acid (JA), which controls a plethora of critically important processes in plants and is orchestrated by the transcription factor MYC2 and its closest relatives in Arabidopsis thaliana. We generated an integrated framework of the response to JA, which spans from the activity of master and secondary regulatory transcription factors, through gene expression outputs and alternative splicing, to protein abundance changes, protein phosphorylation and chromatin remodelling. We integrated time-series transcriptome analysis with (phospho)proteomic data to reconstruct gene regulatory network models. These enabled us to predict previously unknown points of crosstalk of JA to other signalling pathways and to identify new components of the JA regulatory mechanism, which we validated through targeted mutant analysis. These results provide a comprehensive understanding of how a plant hormone remodels cellular functions and plant behaviour, the general principles of which provide a framework for analyses of cross-regulation between other hormone and stress signalling pathways.

Introduction

Plant hormones are structurally unrelated small signaling molecules that play pivotal roles in a wide range of fundamental processes of plants spanning growth, development and responses to environmental stimuli [1]. Hormone perception by plants stimulates a cascade of transcriptional reprogramming that ultimately modifies cellular function and plant behavior [2-5]. This is initiated by one or a family of high-affinity receptors, followed by signal transduction through protein-protein interactions, post-translational modification events and regulation of TF activity that ultimately drives changes in gene expression [2,3,6].

One of the key plant hormones is jasmonic acid (JA), which regulates crucial processes including fertility, seedling emergence, the response to wounding and the growth-defense balance [7]. Jasmonates are perceived as jasmonoyl-isoleucine (JA-Ile) by the co-receptor COI1 (CORONATINE INSENSITIVE1)/JAZ (Jasmonate ZIM domain) complex [8-11]. COI1 is an F-box protein and part of a Skp-Cullin-F-box-E3 ubiquitin ligase complex (SCFCOI1) [12] that targets JAZ proteins for proteasomal degradation upon JA perception. JAZs are transcriptional repressor proteins that inhibit the activity of key TFs of the JA pathway such as the basic helix-loop-helix (bHLH) TF MYC2 and its closest homologues MYC3, MYC4 and MYC5 [13-15] in the absence of JA. The SCFCOI1-JAZ complex tightly controls the level of free non-repressed MYCs in a JA-dependent manner thereby determining the transcriptional output of the entire JA response [8,9,16]. The key regulatory step in the JA pathway is the hormone-triggered formation of a complex between the E3 ligase SCFCOI1 and JAZ repressors that are bound to the master TF MYC2. This results in degradation of JAZ repressors and permits the activity of a master regulatory TF MYC2, accompanied by MYC3, MYC4, MYC5 and numerous other TFs, all of which have distinct but overlapping roles in driving JA-responsive gene expression [13-20]. The result is a cascade of JA-induced genome reprogramming to modulate plant behavior such as plant immune responses [4,19,21] . However, our knowledge of the JA-responsive genome regulatory program and, more broadly, in the plants general response to environmental stimuli is limited currently by assessment of only one or a small number of components.

Here we aim to decipher the MYC2/MYC3-driven regulatory network using a multi-omics analysis that includes the direct targets of key TFs, chromatin modifications, global protein abundance and protein phosphorylation. Our analysis was conducted with etiolated seedlings in which the JA regulatory network is poorly characterized although MYC2 is active [21,22,23]. We discovered that MYC2/MYC3 directly target hundreds of TFs, resulting in a large gene regulatory network that not only amplifies the transcriptional JA response, but also facilitates extensive crosstalk with other signaling pathways. Furthermore, we found that MYC2 has a profound impact on the JA-dependent epigenome, proteome and phosphoproteome. We also generated a network model that predicted new components of the JA signaling pathway that we validated by targeted genetic analyses, demonstrating the power of our integrated multi-omic approach to yield fundamental biological insight into plant hormone responses.

Results

MYC2 and MYC3 target a large proportion of JA-responsive genes

To decipher the JA-governed regulatory network with its high degree of dynamic interconnectivity with other signaling pathways, we applied a multi-omic network approach that is comprised of five newly generated high-quality large-scale datasets (Fig. 1a, b, Extended Data Fig. 1a–i, 2a–d and Supplementary Table 1, 2). MYC2 is the master regulatory TF of JA responses and plants with a null mutation causing a clear decrease in JA sensitivity[15]. Thus, we included the myc2 (jin1–8 SALK_061267) mutant into our analyses (Fig. 1b) [15]. MYC2 is responsible for strong JA-responsive gene activation and acts additively with MYC3 and MYC4 [13,15-20]. myc3 and myc4 single mutants behave like wildtype with regards to JA-induced root growth inhibition. However, in combination with the myc2 mutant, myc2 myc3 double mutants exhibit an increased JA hyposensitivity, almost as pronounced as in myc2 myc3 myc4 triple mutants [13]. We consequently selected MYC3 for an in-depth analysis. In order to better understand how the master TFs MYC2 and MYC3 control the JA-induced transcriptional cascade, we determined their genome-wide binding sites using chromatin immunoprecipitation sequencing (ChIP-seq). Four biological replicates of JA-treated (2 hours) three-day-old etiolated Arabidopsis seedlings that express a native promoter-driven and epitope (YPet)-tagged version of MYC2 and three biological replicates of MYC3 (Col-0 MYC2::MYC2-Ypet, Col-0 MYC3::MYC3-Ypet) were used [24]. The rationale behind dissecting jasmonate signaling in etiolated seedlings is that although MYC2 is highly expressed in etiolated seedlings and regulates important processes such as photomorphogenesis and apical hook formation [21,22,23], a comprehensive characterization of this special developmental stage is still missing.

Figure 1.

Design of our study and key datasets utilized.

a, b, Overview of profiled regulatory layers (a) and detailed description of collected samples (b). c, AnnoJ genome browser screenshot visualizes the binding of MYC2 and MYC3 to the three example genes (IAR3/UR3, ACT1, JAZ9/TIFY7). MYC2/3 binding was determined with ChIP-seq using JA-treated (2 hours) Col-0 MYC2::MYC2-Ypet and Col-0 MYC3::MYC3-Ypet seedlings. Three independent biological ChIP-seq replicates are shown. In addition, mRNA expression of the three example genes WT seedlings (−/+ 2 hours JA) is shown as well. Expression data is derived from RNA-seq analysis. d, Venn diagram illustrates the overlap between MYC2, MYC3 target genes and genes that are differentially expressed after two hours of JA treatment (JA 2h DEGs). e, f, The top-ranked motif in MYC2 (e) and MYC3 (f) ChIP-seq data was the G-box (CAC/TGTG). Motifs were determined by MEME analysis using the top-ranked peaks that were identified by GEM.

Extended Data Fig. 1

Overview of quality metrics of generated ChIP-seq datasets

a-c, Correlation plot of the respective TF ChIP-seq samples is shown. The MYC2 and MYC3 ChIP-seq replicates are shown together in (a). Clustering is determined by the degree of correlation (Pearson correlation). ChIP-seq data is derived from at least three independent experiments: MYC2 (JA, n=4), MYC3 (JA, n=3), ZAT10 (air, n=3; JA, n=2), ANAC055 (JA, n=3). d-i, Cross-correlation (Pearson correlation) plot for the respective TF and histone ChIP-Seq sample is shown. NSC means normalized strand cross-correlation coefficient and RSC means relative strand cross-correlation coefficient. Qtag means quality tag based on thresholded RSC (codes = −2: very low, −1: low, 0: medium, 1: high, 2: very high). All shown TF ChIP-seq replicates are derived from independent experiments: MYC2 (JA, n=4), MYC3 (JA, n=3), ZAT10 (air, n=3; JA, n=2), ANAC055 (JA, n=3). Histone ChIP-seq data is derived from a single experiment (n=1).

Extended Data Fig. 2

Overview of quality metrics of generated RNA-seq and proteome data.

a, b, Multidimensional scaling (MDS) plots of replicate samples of the 24 h JA treatment RNA-seq time-series in WT (a) and the 4 h JA-treatment RNA-seq time-series in WT and myc2 seedlings (b) are shown. Both JA treatment time series consist of three independent samples (n=3) for each time point and genotype. c, d, Principal component analysis (PCA) plots of independent biological replicate samples analyzed by proteomics (c) and phosphoproteomics (d).

We identified 6,736 MYC2 and 3,982 MYC3 high confidence binding sites (p ≤ 1*10−25 and conserved in at least two independent biological replicates), equating to 6,178 MYC2 and 4,092 MYC3 target genes (within 500 nt of binding site center or nearest neighboring gene; Fig. 1c, d and Supplementary Table 1). Of the target genes identified, 3,847 were shared, meaning that almost all MYC3 target genes are also bound by MYC2 (Fig. 1c, d). Their target genes were enriched for JA-related gene ontology terms and for terms related to other hormones (Extended Data Fig. 3a). Target genes shared between MYC2 and MYC3 were significantly enriched (p < 0.05) for more JA-related gene ontology terms than target genes unique to either TF (Extended Data Fig. 3b). Proteins encoded by shared MYC2/MYC3 target genes were enriched for DNA-binding and transcriptional regulatory domains, but in contrast, proteins encoded by MYC2-only target genes were enriched for kinase domains (Supplementary Table 3). No significant protein domain or gene ontology term enrichment was detected amongst the small number of MYC3-only targets (Supplementary Table 3). Collectively, these data indicate that MYC2 and MYC3 have the potential to regulate 23.2% of genes in the Arabidopsis genome (27,655 coding genes). However, binding events are not necessarily regulatory [2,3,25]. We determined that 2,522 genes are differently expressed (false discovery rate, FDR < 0.05) after two hours of JA treatment using RNA-seq. A third (843 genes) of JA-modulated genes were directly bound by MYC2 or MYC3 (Fig. 1d and Supplementary Table 4). This is consistent with the important role of MYC2/3 in JA-responsive gene expression [13,15-17,19,20]. The majority of JA-responsive direct MYC2/3 target genes are transcriptionally upregulated after JA application indicating that MYC2 and MYC3 predominantly act as transcriptional activators (Extended Data Fig. 3c).

Extended Data Fig. 3

MYC2 and MYC3 act predominantly as activators for a functionally diverse range of target genes.

a, b, Gene ontology (GO) analyses using a hypergeometric distribution of all MYC2 and MYC3 targets (a) as well as MYC2 only and MYC2/MYC3 shared targets (b) are shown. Data is derived from four independent MYC2 (n=4) and three independent MYC3 (n=3) ChIP-seq samples. Analyses were conducted using clusterProfiler. c, Bar plots shows the portion of JA-induced and JA-repressed genes that are bound by MYC2 and MYC3. d, e, The CACG[A/C]G motif (286 sites, E = 2*10−52) (d) and the AT[A/T][A/T] [A/T]ATA motif (714 sites, E = 8.9*10−35) (e) were enriched in MYC2 high-confidence target regions that do not contain a G-box or the degenerate G-box motifs CATGTG or CACGTT.

The G-Box (CA[C/T]GT[G/T]) was the most common DNA sequence motif found at MYC2 or MYC3 binding sites, which is concordant with the observation that they shared a large proportion of their binding sites (Fig. 1e, f). This motif was also of similar sequence to a motif bound by MYC2 determined in vitro [26]. The majority of MYC2 and MYC3 binding sites contained the G-Box motif (MYC2: 4,240 of 6,736; MYC3: 3,072 of 3,982; Fig. 1e, f and Supplementary Table 5). However, the absence of the motif from a substantial number of MYC2 and MYC3 binding sites suggests the TFs may bind indirectly to some sites through partner protein(s). We identified putative partner TFs by determining DNA motifs enriched in MYC2 binding sites that did not contain a G-box motif. The most strongly enriched motifs were CACG[A/C]G (286 sites, E = 2*10−52), which may correspond to the TFs CAMTA1 (AT4G15090) or FAR1 (AT5G09410), and AT[A/T][A/T] [A/T]ATA (714 sites, E = 8.9*10−35), which may correspond to ARID family TFs (AT2G17410, AT1G04880) (Extended Data Fig. 3d, e). Molecular investigation of these TFs would be required to determine whether they do bind co-operatively with MYC2 to DNA.

Master TFs directly target the majority of signaling components in their respective pathway, a phenomenon which has already been observed for the ethylene, abscisic acid and cytokinin signaling pathways [2,3,27]. This pattern also holds true for the JA signaling pathway. Our MYC2/MYC3 ChIP-seq analyses determined that approximately two thirds of genes encoding for known JA pathway components (112 of 168 genes for MYC2 and 96 of 168 genes for MYC3) are bound by MYC2 and MYC3 (Extended Data Fig. 4a, b and Supplementary Table 6). Interestingly, the majority of all known JA genes that were differentially expressed following JA treatment were bound by MYC2 or MYC3 whereas fewer non-differentially expressed known JA genes were directly targeted (Extended Data Fig. 4b and Supplementary Table 6). MYCs initiate various feed forward loops that allow a rapid activation of the transcriptional JA response [19,28]. Our ChIP-seq approach revealed that beyond the autoregulation of MYC2 and MYC3, these TFs also regulate JA biosynthesis either directly by targeting the JA biosynthesis genes LOX2–4, LOX6 and AOS or indirectly through binding to the AP2-ERF TF gene ORA47 or (Supplementary Table 1, 6). In addition, MYCs simultaneously target various negative regulators enabling MYCs to efficiently dampen the JA response pattern (Extended Data Fig. 4c). Key negative regulators of JA signaling are the JAZ repressors, a gene family of 13 members in Arabidopsis [29], which can interact with the adaptor protein NINJA to confer TOPLESS-mediated gene repression [30]. Strikingly, all JAZs and also NINJA are directly bound by MYC2 and MYC3 (Extended Data Fig. 4c), with the probable effect of dampening the JA response thereby preventing excessive activation of JA signaling.

Extended Data Fig. 4

MYC2 and MYC3 regulate the majority of JA signaling pathway components.

a, Schematic overview of known MYC2/MYC3-targeted JA pathway components. Genes that are directly targeted by MYC2/MYC3 are highlighted in orange. b, Binding behavior of MYC2 and MYC3 at known JA genes (Supplementary Table 6) is shown. Known JA genes are grouped into non-differentially expressed and JA differentially expressed genes. c, AnnoJ genome browser screenshot visualizes MYC2 and MYC3 binding at all 13

MYC2 and MYC3 activate the JA response through a large TF network

To study the MYC2 and MYC3-governed transcriptional regulatory network in more detail, we investigated the relationship between MYC2/MYC3-bound TF-encoding genes and their transcriptional responsiveness to JA treatment. We conducted a JA time course experiment (time points 0, 0.25, 0.5, 1, 2, 4, 8, 12, 24 h post JA treatment) identifying a total of 7,377 differentially expressed genes at one or more time points within 24 h of JA treatment (Supplementary Table 4). Differentially expressed genes were categorized into clusters with similar expression trends over time to facilitate visualization of complex expression dynamics and enriched functional annotations (Extended Data Fig. 5a and Supplementary Table 7). The largest upregulated cluster was the “JA cluster” which was enriched for gene ontology (GO) terms associated with JA responses (Fig. 2a). In contrast, the “Cell wall cluster” was the largest cluster of downregulated genes and enriched for GO terms associated with cell wall organization, development and differentiation (Fig. 2b). These two main clusters illustrate the defense-growth trade-off when defense pathways are activated [31].

Extended Data Fig. 5

MYC2 and MYC3 target a large number of TFs.

a. Cluster analysis revealed the 5 other main clusters in the JA time course experiment. Clusters visualize the log2 fold change expression dynamics over the indicated 24 hours’ time period. The three strongest enriched gene ontology terms for each cluster are shown as well. b, Pie chart indicates the proportions of TFs that are transcriptionally induced by JA, bound by MYC2/MYC3, or both. c, d, Overview of MYC2/MYC3-bound plant hormone genes (c) and TFs (d) is shown. Plant hormones are abbreviated (ET (ethylene), BR (brassinosteroids), GA (gibberellic acid), ABA (abscisic acid), SA (salicylic acid), CK (cytokinin), AUX (Auxin), K (karrikin), SL (strigolactones)).

Figure 2.

MYC2 and MYC3 target a large proportion of JA-responsive genes that encode TFs.

a, b, Cluster analysis revealed the two main clusters in the JA time course experiment. The JA cluster (a) with 796 genes reflects the majority of JA-induced genes and the cell wall cluster (b) with 647 genes represents the largest cluster of JA repressed genes. Clusters visualize the log2 fold change expression dynamics over the indicated 24 hours’ time period. The three strongest enriched gene ontology terms for each cluster are shown as well. Clusters were identified by STEM clustering (Pearson correlation, minimum correlation 0.7, up to 50 permutations, significant clusters had Bonferroni corrected p < 0.05). For each of the indicated time points, the expression of three independent samples (n=3) was measured using RNA-seq. c, Bar plots illustrates the potential of MYC2 and/or MYC3 to bind to a portion of JA differentially expressed genes (JA DEGs) at the indicated time points. JA DEGs for all time points were identified with RNA-seq. MYC2/3 targets are derived from ChIP-seq analysis using Col-0 MYC2::MYC2-Ypet and Col-0 MYC3::MYC3-Ypet seedlings that were treated for two hours with JA. d, MYC2 and MYC3 target genes from a wide range of TF families. TF families are classified into four different groups; MYC2/MYC3 targets and differentially expressed after JA treatment (blue), MYC2/MYC3 targets and not differentially expressed (orange), not bound by MYC2/MYC3 but differentially expressed (grey) and not bound by MYC2/MYC3 but not differentially expressed (green). e, Nodes represent JA TFs for which direct binding data was generated. ChIP-seq data is indicated by presence of *, all other data was DAP-seq. Edges represent binding events and are directed. Self-loops indicate TF binds to its own locus, indicative of potential auto-regulation. Expression of the TF at 0.5 h after JA treatment is represented by color scale. f, Pearson correlation of TFs’ target sets of genes. Numerals in brackets indicate total number of target genes. ChIP-seq data is indicated by presence of *, all other data was generated by DAP-seq. ChIP-seq data is derived from at least three independent experiments: MYC2 (JA, n=4), MYC3 (JA, n=3), ZAT10 (air, n=3; JA, n=2), ANAC055 (JA, n=3). DAP-seq data is derived from a single experiment (n=1).

Based on our MYC2/MYC3 ChIP-seq data that is derived from a two hour-long JA treatment, up to 63% (0.5 h JA treatment) of differentially expressed genes at any given time point have the potential to be directly bound by MYC2 and/or MYC3 (Fig. 2c), highlighting the important role of MYCs in transcriptionally regulating JA responses. Our analysis also determined that 522 of 1,717 known or predicted TFs were differentially expressed within 24 h of JA treatment (Extended Data Fig. 5b). Half of these (268), representing 36 of 58 TF families, were also direct MYC2 or MYC3 targets (Fig. 2d and Extended Data Fig. 5b) indicating that MYC2 and MYC3 cooperatively control a massive TF network. The three most numerous families (ERFs, bHLHs and MYBs) in the Arabidopsis genome had the most JA-responsive MYC2 or MYC3 targeted members which is concordant with their previously annotated roles in JA responses (Fig. 2d) [32]. Plant hormone crosstalk is critical for an appropriate cellular response to environmental stimuli and numerous reports describe that MYC2 connects the JA pathway to other major plant hormone pathways [23,33]. To investigate this crosstalk function of MYC2 and MYC3 in more detail, we utilized our ChIP-seq data to determine the number of plant hormone TFs that are bound by MYC2 and MYC3. We found that 37 to 59% of annotated hormone pathway genes are bound by MYC2 and MYC3 and that their expression changes in response to 24 hours of JA treatment (Extended Data Fig. 5c). In addition, we discovered 122 annotated hormone TFs, with representatives from all hormone pathways, that are bound by MYC2 and MYC3 and 118 of these are differentially expressed (Extended Data Fig. 5d and Supplementary Table 1).

We next set out to better understand the target genes of the network of TFs downstream of MYC2 and MYC3. To do so we conducted ChIP-seq or DNA affinity purification sequencing (DAP-seq) on a subset of TFs that were direct MYC2/3 targets and rapidly upregulated (within 0.5 h) by JA treatment (DREB2B, ATAF2, HY5, RVE2, ZAT18; Fig. 2e) or were members of the upregulated “JA cluster” (TCP23; Fig. 2a). We also included TFs with known roles in JA signaling (ERF1, ORA59, ANAC055, WRKY51, ZAT10) [34-38]. These TFs formed a highly connected network, with all TFs except DREB2B targeting at least two TFs in the network and being themselves targeted by two TFs (Fig. 2e and Supplementary Table 8). Auto-regulation was common, with seven TFs targeting their own loci (Fig. 2e). The target genes of ZAT10, ANAC055 and ATAF2 were most similar to those of MYC2/3 (Fig. 2f). Consistent with this, their target genes shared several significantly enriched gene ontology terms (adjusted p < 0.05), suggesting related functions in jasmonate signaling (Extended Data Fig. 6a). ORA59 and ERF1, along with DREB2B, formed a distinct group that targeted a related set of genes (Extended Data Fig. 6a). Notably, ERF1 and ORA59 also shared significant enrichment of a separate set of gene ontology terms with one-another, but that were not enriched amongst MYC2/3 targets (Extended Data Fig. 6a). This is consistent with the joint role of ERF1 and ORA59 in controlling a pathogen defense arm of JA signaling [34,35]. No gene ontology terms were enriched amongst the targets of DREB2B. WRKY51 and RVE2 had relatively few enriched gene ontology terms but shared most of these with one-another (Extended Data Fig. 6a). Most of the terms related to anti-insect defense and were a subset of the enriched MYC2/3-ZAT10-ANAC055-ATAF2 gene ontology terms (Extended Data Fig. 6a). ZAT10 and ANAC055 are known regulators of anti-insect defense and our results suggest WRKY51 and RVE2 may also be involved in this component of jasmonate responses [39]. Interestingly, ZAT10 belongs to a group of genes that is inducible by the JA precursor 12-oxo-phytodienoic acid (OPDA) and not by JA [40]. We found that approximately one third of OPDA-specific response genes (45 genes) are targeted by MYC2 (Supplementary Table 3). Taken together, our analyses determine that MYC2 and MYC3 shape the dynamic JA response through the activation of a large TF network that includes various potentially coupled feedforward and feedback loops and that allows extensive cross-communication with other signaling pathways.

Extended Data Fig. 6

Overview of MYC-controlled TF network.

a. Significantly enriched (adjusted p < 0.05) gene ontology terms amongst the target of each TF. For each TF the 4 terms with the lowest p-value are shown, some of which are redundant between TFs. No enriched terms were detected for DREB2B targets. ChIP-seq data is indicated by presence of *, all other data was generated by DAP-seq. ChIP-seq data is derived from at least three independent experiments: MYC2 (JA, n=4), MYC3 (JA, n=3), ZAT10 (air, n=3; JA, n=2), ANAC055 (JA, n=3). DAP-seq data is derived from a single experiment (n=1).

We examined the effect of removing MYC2 activity on JA-responsive transcriptional regulation by generating transcriptomes from a myc2 null mutant (jin1–8) in an early JA response time series (0, 0.5, 1 and 4 h). The response of myc2 mutants to JA differed from that of WT plants. There were 2,905 unique genes differentially expressed between the myc2 and WT across the time-series (pairwise comparisons between genotypes at each time point; Supplementary Table 9). JA-responsive gene expression did occur in myc2 plants, consistent with the partially redundant function of MYC2, MYC3 and MYC4 [13]. However, JA-responsive genes were upregulated more highly in WT than myc2 (Supplementary Table 9). The JAZ/TIFY genes illustrate this, with 9 of the 12 genes being upregulated more highly in WT than myc2, as well as reaching peak expression at earlier time points in WT (0.5 or 1 h; (Extended Data Fig. 7a). Overall, a majority of the MYC2 target genes differentially expressed between myc2 and WT were more highly expressed in WT, indicating loss of MYC2 function reduces the JA responsiveness of these genes (Extended Data Fig. 7b). A total of 130 TFs targeted by MYC2 were differentially expressed in myc2 mutants compared with WT seedlings, including the TFs ATAF2, ERF1, ANAC055 and ZAT10 whose targets we had determined by DAP-Seq or ChIP-Seq (Supplementary Table 10). The myc2 mutation also affected expression of secondary, indirect MYC2 target genes (ie. genes targeted by MYC2-regulated TFs, but not by MYC2 itself). Between 23.6 and 26.3% of genes targeted by each of ATAF2, ERF1, ANAC055 and ZAT10, and not by MYC2, were differentially expressed in myc2 plants compared with WT (Extended Data Fig. 7c and Supplementary Table 11). Taken together, these data demonstrate that MYC2 regulates gene expression through a large network of downstream TFs during response to JA stimulus.

Extended Data Fig. 7

MYC2 partially controls expression of JAZ repressors.

a, Individual plots show expression of all JAZ/TIFYs and NINJA in WT (blue) and myc2 (orange) seedlings following JA treatment. log2 fold change (FC) was calculated relative to their respective 0 h (ie. non-treated) control samples. b, Bar chart shows the number of differentially expressed (DE) genes at each time point after JA treatment between WT and myc2 seedlings. The bar chart also indicates how many of these DE genes were direct binding targets of MYC2 (in ChIP-seq assays) and whether they were more highly expressed in WT (blue) or myc2 (orange) seedlings. c, Charts indicates of how MYC2 indirectly affects the expression of downstream genes through secondary TFs. The expression of genes in pairwise comparisons of WT and myc2 transcriptomes at 0, 0.5, 1 and 4 h was assessed. Only genes that were direct targets of the TFs ATAF2, ZAT10, ANACO55 and ERF1, and not direct targets of MYC2, were analyzed which are termed “non-MYC2 target genes”. ATAF2, ZAT10, ANACO55 and ERF1 are themselves direct targets of MYC2 and their expression levels were decreased in myc2 relative to WT, indicating they are directly regulated by MYC2. DE indicates differentially expressed genes.

MYC2 controls JA-induced epigenomic reprogramming

Reprogramming of the epigenome is an integral part of development and environmental stimulus-induced gene expression [41]. For example, activation of the transcriptional JA response requires the formation of MYC2/MED25-mediated chromatin looping [42]. To investigate the extent of JA-induced changes in chromatin architecture and the regulatory importance of MYC2 in this response, we conducted ChIP-seq assays to profile the genome wide occupancy of the histone modification H3K4me3 and the histone variant H2A.Z in untreated/JA-treated (4 h) WT and myc2 seedlings. Trimethylation of H3K4me3 marks active and poised genes and the histone variant H2A.Z confers gene responsiveness to environmental stimuli [43,44]. mRNA expression was monitored in parallel using RNA-seq. JA treatment led to a reprogrammed chromatin landscape with several thousand differentially enriched H3K4me3 and H2A.Z domains (Extended Data Fig. 8a–c and Supplementary Table 12). We identified 826 differentially expressed genes (675 induced, 151 repressed; WT control v. JA-treated) in that experiment. In line with the predominantly activating function of MYC2 (Extended Data Figure 3c), the JA-induced genes had a stronger promoter enrichment of MYC2 than the JA-repressed genes (Fig. 3a). The JA-induced genes had an increase of H3K4me3, whereas JA-repressed genes had no dynamic change in the level of H3K4me3 (Fig. 3b, d). Strikingly, myc2 mutants only display a compromised increase of H3K4me3 after JA treatment suggesting that the JA-induced trimethylation of H3K4me3 depends on functional MYC2 (Fig. 3b–d and Extended Data Fig. 8a). The impact of the myc2 mutation on JA-induced H3K4me3 changes can also be observed at JA-induced genes that are not directly targeted by MYC2 (Extended Data Fig. 8e, f), which is potentially caused by the decreased expression of MYC2-targeted TFs. The scenario of a direct MYC2 regulation is illustrated by two JA-induced genes, JAZ2 and GRX480, which are directly targeted by MYC2. Their expression depends on MYC2 and their JA-induced increase of gene body-localized H3K4me3 partially depends on MYC2 (Fig. 3d and Extended Data Fig. 8d). However, whether the MYC2-dependent changes in H3K4me3 levels precede transcription or rather reflect increased transcriptional activity cannot be addressed by these experiments. In contrast, JA-induced changes in H2A.Z occupancy are only slightly affected in myc2 mutants (Extended Data Fig. 8a, g, h) suggesting that JA-induced H2A.Z dynamics are either independent of MYC2, precede MYC2 binding or other MYCs such as MYC3, MYC4 and MYC5 are functionally redundant in regulating H2A.Z dynamics.

Extended Data Fig. 8

JA shapes the local chromatin architecture.

a, Bar plot shows the impact of two hours JA treatment on the genome-wide distribution of H3K4me3 and H2A.Z domains. Occupancy was determined in untreated/JA-treated WT and myc2 seedlings using ChIP-seq. SICER was used to identify the number of histone domains that show an increase (blue) or decrease (orange) of enrichment in response to JA. b, c, Heatmaps show the occupancy of H3K4me3 and H2A.Z from 1 kb upstream to 2 kb downstream of the transcriptional start site (TSS) at all Arabidopsis genes (TAIR10). Heatmaps are shown for H3K4me3 (b) and H2A.Z (c) in untreated and JA-treated (4 h) WT and myc2 seedlings. d, Quantification of H3K4me3 occupancy at JAZ2 and GRX480 is shown. It was calculated as the ratio between the respective ChIP-seq sample and the WT IgG control. e, f, Aggregated profiles show the log2 fold change enrichment of H3K4me3 at JA DEGs that are directly (e) and not directly targeted (f) by MYC2 from 2 kb upstream to 2 kb downstream of the transcriptional start site (TSS). g, h, Plot profiles show the log2 fold change enrichment of H2A.Z in WT (g) and myc2 mutants (h) from 2 kb upstream to 2 kb downstream of the transcriptional start site (TSS) at JA-induced and JA-repressed genes.

Figure 3.

The JA-responsive epigenome.

a, b, c, Aggregated profiles show the log2 fold change enrichment of MYC2 (a) and H3K4me3 (b, c) from 2 kb upstream to 2 kb downstream of the transcriptional start site (TSS) at JA-induced and JA-repressed genes. Profile of MYC2 is shown for Col-0 MYC2::MYC2-Ypet (a) seedlings and H3K4me3 profiles are shown are shown for WT (b) and myc2 (c) seedlings. d, AnnoJ genome browser screenshot visualizes MYC2 binding, mRNA expression and H3K4me3 occupancy at two example genes (JAZ2, GRX480) in WT and myc2 seedlings. All tracks were normalized to the respective sequencing depth.

JA extensively remodels the (phospho)proteome

We next explored how JA remodels the proteome and phosphoproteome of etiolated seedlings. Hormone signal transduction typically modifies phosphorylation of downstream proteins, changing their activity independent of transcript abundance [6]. Transcript abundances are also frequently weakly correlated with protein abundances [45,46]. Consequently, proteomic and phosphoproteomic analyses yield additional insight into gene regulatory networks. We determined that the loss of MYC2 caused substantial changes to the JA-responsive proteome and phosphoproteome; 1,432 proteins and 939 phosphopeptides (corresponding to 567 genes) were significantly differentially abundant in WT seedlings relative to myc2 seedlings after 2 h JA treatment (q < 0.1; Fig. 4a and Supplementary Table 13, 14). WT seedlings responded to JA (161 proteins, 443 phosphopeptides, WT JA v. WT air) and the response was smaller without functional MYC2 (79 proteins, 93 phosphopeptides, myc2 JA v. myc2 air) (Fig. 4a). These extensive changes in phosphopeptide abundance are consistent with the observation that 118 genes encoding protein kinases were differentially expressed between WT and myc2 seedlings in our transcriptome experiments (Supplementary Table 9).

Figure 4.

Loss of functional MYC2 impacts the global proteome and phosphoproteome.

a, Total significantly differentially abundant (False discovery rate q < 0.1 (estimated using a modified permutation plug-in method) proteins and phosphopeptides detected in comparisons between JA-treated (2 h) WT and myc2 seedlings and air controls. Three independent experiments (+/− 2h JA) (n=3) were conducted for WT and myc2 seedlings. For the third experiment, only the JA treatment was conducted. b, Venn diagram showing the overlap between significantly differentially abundant proteins, transcripts and differentially phosphorylated proteins in JA-treated WT seedlings compared to mock-treated WT controls. Also shown is the overlap with MYC2/3 target genes. c, Correlation between rank-normalized log2(FPKM)s of detected proteins and transcripts in WT seedlings treated with JA for 2 h (p-value cut-off was < 0.05 using paired t-tests). Scatter plot of log2 fold change in WT JA-regulated transcript levels versus log2 fold change in levels of corresponding proteins. Protein and transcript data is derived from three independent experiments (n=3) using WT and myc2 seedlings. d, Heatmap represents relative TPM of 137 isoform pairs exhibiting isoform switch events. Ratio calculated as logTPM (isoform 1/isoform 2). e, Plot shows an example of a transcript pair originating from AT2G43680 that had isoform switch events following JA treatment. Expression data is derived from a JA time course experiment. For each of the indicated time points, the expression of three independent samples (n=3) was measured using RNA-seq. Line indicates mean TPM of three independent samples. Shaded regions indicate standard error of these data.

Some direct overlap existed between proteins or phosphopeptides and transcripts responsive to JA treatment (Fig. 4b). Both transcripts and proteins encoded by 28 genes were differentially expressed in JA-treated WT seedlings relative to air controls (Fig. 4b). A further 33 differentially expressed proteins in JA-treated WT seedlings had no corresponding differentially expressed transcript but were encoded by genes that are targeted by MYC2 and MYC3 (Fig. 4b). Differentially abundant phosphopeptides were detected that corresponded with 15 differentially expressed transcripts (Fig. 4b). Transcript and protein abundances were weakly positively correlated (Pearson correlation 0.40341) in JA-treated WT seedlings (Fig. 4c), in agreement with prior studies [45,46]. The protein of only one known JA pathway component was differentially abundant in JA-treated WT seedlings relative to controls, and none were differentially phosphorylated. The fact that only a single JA-regulated protein and no phosphoproteins are annotated as JA pathway components may indicate that existing annotations are overly-dependent on transcriptome data and that consideration of (phospho)proteome data deepens our understanding of JA responses.

Alternative splicing can occur rapidly in response to environmental stimuli, contributing to transcriptome reprogramming and potentially fine-tuning physiological responses [47]. It is central to JA-mediated regulation of transcription, with an alternative isoform of the repressor JAZ10 creating a negative feedback loop that desensitizes cells to a JA stimulus [48,49]. However, the extent of alternative splicing in JA signaling beyond the JAZ repressors is poorly characterized. We observed that phosphorylation of proteins involved in RNA recognition and nucleotide binding was disrupted in JA-treated myc2 mutants compared with WT seedlings. The spliceosome was the only pathway significantly enriched amongst these differentially phosphorylated proteins (p < 0.05, 18 genes matched) suggesting that MYC2 may influence JA-responsive alternative splicing. Furthermore, 18 genes with splicing-related annotations were differentially expressed between myc2 and WT seedlings in our transcriptome experiments (Supplementary Table 9). None of the differentially phosphorylated spliceosome components was differentially expressed.

We examined isoform switching events across our JA transcriptome time-series, where the most abundant of two isoforms from a single gene changes, to determine the extent of JA-responsive alternative splicing (Fig. 4d, e, Supplementary Table 15). There were 151 switch events, corresponding to 137 isoform pairs from 120 genes, within 24 h of JA treatment. These were identified from 30,547 total individual transcripts detected (average TPM > 1; Supplementary Table 16). Two of the genes exhibiting isoform switches had prior JA annotations (RVE8/AT3G09600, SEN1/AT4G35770, Supplementary Table 15) and others were annotated to a variety of processes (including auxin, ABA, light signaling, disease response, amongst many others), but there was no significant enrichment of any gene ontology terms or pathways. This indicates that MYC2 influences alternative splicing that diversifies the transcriptome in response to a JA stimulus.

Multi-omic modelling of the JA response regulatory program

We wanted next to characterize the broader JA response genome regulatory program so that we could increase our understanding of the roles of known JA TFs within this and identify new potential regulatory interactions. To do so we generated a gene regulatory network model encompassing the (phospho)proteomic and time-series transcriptomic data (Extended Data Fig. 9a and Supplementary Table 17). Inclusion of the (phospho)proteomic data expanded the network by an additional 957 nodes (genes) compared with a transcript only network (3,409 v. 4,366 nodes, 28% larger) (Supplementary Table 17). The (phospho)proteomic and transcript data shared 217 nodes within the network, a relatively small proportion, indicating these datasets complement one-another when attempting to characterize the JA-response genome regulatory program.

Extended Data Fig. 9

The JA gene regulatory network.

a, Illustration of JA gene regulatory network for 1, 2 and 4 h time points. Edges were predicted using phosphoproteome (green), proteome (orange) and transcriptome (blue) data. Node sizes are scaled by normalized motif score, with larger nodes indicating greater scores and likely greater importance within the network. Edges predicted early in the time-series transcriptomic data are red (0.25 – 2 h), edges predicted late are blue (4 – 24 h). Proteome and phosphoproteome-data-predicted edges are grey and green, respectively.

Many known JA signaling components were present in the 100 most important predicted components of the network (MYC2, ERF1, JAZ1, JAZ2, JAZ5, JAZ10, ATAF2 and others; within top 100 of 4366 components by normalized motif score; Supplementary Table 17). MYC2 was predicted to regulate a subnetwork of 26 components, 23 of which were validated as directly bound by MYC2 in ChIP-seq assays (88.5%, Extended Data Fig. 10a and Supplementary Table 1, 17). We further validated the network by comparing the ChIP/DAP-seq data previously collected for the remaining 12 JA TFs to their targets in the gene regulatory network (Fig. 2e, f, Extended Data Fig. 10b and Supplementary Table 18). The gene regulatory network identified all of these TFs as components of the JA response, except MYC3 (Supplementary Table 17). It is likely that MYC3 was not part of the network due to it being only modestly differentially expressed following JA treatment and not being detected in the (phospho)proteome analyses (Supplementary Tables 4, 13, 14). The wider validation of targets was less strong than for MYC2, ranging from 0% to 33.3%. This could reflect the possibility that interactions predicted by the gene regulatory network may not identify all intermediate components. Lastly, we examined known genetic interactions. The MYC2 subnetwork included activation of JAZ10 within 0.5 h of a JA stimulus, with JAZ10 reciprocally repressing MYC2 (Extended Data Fig. 10a, b). This is consistent with the known role of JAZ10 in establishing negative feedback that attenuates JA signaling [49]. MYC2 was also predicted to activate AIB (ABA-INDUCIBLE BHLH-TYPE TRANSCRIPTION FACTOR, JAM1/bHLH017/AT2G46510) (Extended Data Fig. 10a, b), establishing a negative feedback loop in which AIB negatively regulates MYC2. This is in line with prior studies, which established AIB is dependent upon and antagonistic to MYC2, thereby repressing JA signaling [50,51]. Confirmation by both genetic data from the literature and our DAP/ChIP-seq experiments indicates that our gene regulatory network modelling approach is a useful tool to identify new regulatory interactions within JA signaling and to better understand known regulatory interactions.

Extended Data Fig. 10

Gene regulatory network validation against ChIP/DAP-seq data.

a, The MYC2 subnetwork is shown. Edges are directional and red edges exist at early time points (0.25 – 2 h), blue only at late time points (4 – 24 h). Thicker edges with chevrons indicate that MYC2 were directly bound to that gene in our ChIP-seq experiments. b, Validated edges are those between TFs and first neighbors in the JA gene regulatory network for which the first neighbor was also a direct target of the TF in ChIP/DAP-seq assays. These edges are indicated by chevrons. Early time-series transcriptome-predicted edges (0.25 – 2 h) are red and later edges (4 – 24 h) are blue. Edges detected in the proteomic data are grey and those detected in the phosphoproteomic data are green. c, Bar plot shows quantification of JA-induced root growth inhibition in the indicated T-DNA alleles. Seedlings were grown on LS media with or without 20µM MeJA. WT seedlings serve as independent controls for each tested T-DNA line. Sample size number n is shown within the respective bars. Samples are derived from three independent experiments. Asterisks represent significant differences between WT (−/+ JA) and indicated T-DNA lines (−/+ JA) (two-way ANOVA with Bonferroni post test, ns (not significant) p > 0.05, *p < 0.05, **p < 0.01, ***p < 0.001). d, Subnetwork of CYP708A2 is shown.

Crosstalk between hormone response pathways permits fine-tuning of plant growth and development in response to diverse environmental signals [1]. We examined the potential points at which MYC2 may interface directly with other hormone signaling pathways, since MYC2 is the master regulator of JA responses and one of the first TFs activated by JA. The MYC2 subnetwork identified a potential route for JA signaling to cross-regulate auxin hormone signaling. MYC2 activated ARF18 and ARF18 reciprocally activated MYC2 (Extended Data Fig. 10a and Supplementary Table 17). It also indicated that MYC2 may promote ethylene signaling by activating MAP kinase kinase 9 (MKK9) (Extended Data Fig. 10a). Prior genetic studies determined that MKK9 induces ethylene production, but had not examined a possible link with JA signaling [52]. Positive crosstalk is known to exist between JA and auxin signaling though the mechanism is not clearly determined [53]. RGL3, a regulator of gibberellic acid (GA) signaling previously associated with JA-GA crosstalk, was also present within the MYC2 subnetwork (Extended Data Fig. 10a), predicted to inhibit MYC2 but not to be reciprocally regulated by MYC2 [54]. These three interactions are potential points at which crosstalk can occur rapidly during a JA response with auxin, gibberellin and ethylene.

We next examined the broader gene regulatory network to identify additional predicted points of crosstalk between JA and other signaling pathways. The model predicted that STZ/ZAT10 is a key early hub through which JA signaling is prioritized over several other hormone and stress response pathways (Fig. 5a and Supplementary Table 17). STZ/ZAT10 is known to be a transcriptional repressor from genetic studies [55] and, consistent with this, our model predicted that it inhibited the majority of genes it regulates (25 of 34 genes). WRKY40, WRKY70, DDF and ERF6 were all predicted to be inhibited by STZ/ZAT10 within 0.25 h of a JA stimulus and GRX480 within 1 h. Direct binding of STZ/ZAT10 to ERF6 was detected in ChIP-seq assays (Supplementary Table 11). WRKY40 and WRKY70 are both early brassinosteroid response components that repress defense responses [56]. DDF1 promotes resistance to drought, cold, heat and salinity stress by reducing endogenous gibberellin abundance [57]. ERF6 similarly promotes drought resistance by reducing gibberellin abundance [58]. GRX480 regulates the negative crosstalk between salicylic acid and both JA/ethylene signaling through the direct interactions with TGA TFs [59,60]. The model also predicts that ERF6, WRKY70 and DDF1 exert negative feedback on STZ/ZAT10 by activating JAZ8 within 0.25 h of the JA stimulus (Fig. 5a and Supplementary Table 17). JAZ8 is a known repressor of JA signaling and is predicted to repress STZ/ZAT10 [61]. In sum the gene regulatory network predicts that STZ/ZAT10 is an important hub for JA signaling to be prioritized over other hormone and stress response pathways (Fig. 5a).

Figure 5.

JA response genome regulatory model positions known and new components.

a, b, Subnetworks of STZ (a) and ABO3 (b) are shown. Edges are directed. Red edges exist at early time points (0.25 – 2 h), blue only at late time points (4 – 24 h). Thicker edges with chevrons indicate a MYC2 directly bound that gene in ChIP-seq experiments. c, d, JA-induced root growth inhibition assay identified ABO3 as a positive JA regulator. Seedlings were grown on LS media with (d) or without 50µM MeJA (c). WT and myc2 seedlings served as controls. e, Quantification (± SEM) of JA-induced root growth inhibition in WT, myc2 and abo3 seedlings is shown. Sample size number n is shown within the respective bars. Samples are derived from three independent experiments. Asterisks represent significant differences between WT (−/+ JA) and abo3 mutants (−/+ JA) (two-way ANOVA with Bonferroni post test, ***p < 0.001). f, g, Root growth inhibition assay identified two cyp708A2 T-DNA alleles as JA hypersensitive. Seedlings were grown on LS media with (f) or without 50µM MeJA (g) and WT and myc2 seedlings serve as controls. h, Bar plot shows quantification (± SEM) of JA-induced root growth inhibition in WT and cyp708A2 seedlings. Sample size number n is shown within the respective bars. Samples are derived from three independent experiments. Asterisks represent significant differences between WT (−/+ JA) and cyp708A2 (−/+ JA) seedlings (two-way ANOVA with Bonferroni post test, **p < 0.01, ***p < 0.001).

Large-scale data-mediated identification of new JA regulators

We next utilized our regulatory network and large-scale datasets to identify novel regulators of the JA pathway using the JA root growth inhibition assay as our experimental readout. First, we focused on ABO3 (ABA overly sensitive 3), which is directly targeted by MYC2 and MYC3 (Supplementary Table 1) and whose subnetwork is comprised of 26 predicted regulated genes, the majority of which is positively regulated (22 of 26 genes) (Fig. 5b). ABO3 encodes the Arabidopsis WRKY TF gene WRKY63, which is involved in stress gene expression and drought tolerance [62]. To investigate the importance of the ABO3 subnetwork in JA signaling, we tested abo3 T-DNA mutant seedlings (SALK_075986C) in a JA-induced root growth inhibition assay. We found that abo3 mutants show a weak JA hypo-sensitive root growth inhibition phenotype (Fig. 5c–e) indicating that ABO3 is positive regulator of JA signaling and that our network approach is able to identify new pathway components.

Next, we expanded our phenotyping analysis to T-DNA lines of genes that display the strongest binding of MYC2 and MYC3 in their promoters (Supplementary Table 1, 18). The rationale behind this approach is that master TFs target the majority of key signaling components in their regulated respective pathways and that these are often the most strongly bound targets [2,3,27]. Of 99 genes tested (194 T-DNA lines in total, Supplementary Table 19), we discovered six genes that, when mutated, display mild JA root growth phenotypes (Extended Data Fig. 10c and Supplementary Table 19). Mild phenotypes as well as their low frequency were not surprising since gene redundancy is very common in the Arabidopsis genome and even the mutation of the master TF MYC2 only causes a mild JA-hyposensitive root growth phenotype (Fig. 5c–e)[15]. Among these genes was the cytochrome P450 enzyme CYP708A2 gene from which both tested T-DNA mutant alleles exhibit a JA hypersensitive root phenotype (Fig. 5f–h). Interestingly, our network analysis also discovered CYP708A2 as a regulatory hub (Extended Data Fig. 9a, 10d). CYP708A2 is involved in triterpene synthesis, which is known to be stimulated by JA [63]; future studies are however needed to further decipher the role of CYP708A2 in JA signaling. Another interesting uncharacterized gene that we discovered caused a JA phenotype is a Sec14p-like phosphatidylinositol transfer family protein (AT5G47730) (Extended Data Fig. 10c and Supplementary Table 19). Phosphatidylinositol transfer proteins (PITPs) are crucial for the phosphatidylinositol homeostasis in plants [64] and inositol polyphosphates have been implicated in COI1-mediated JA perception [65]. Taken together, these data show that our multi-omic approach goes beyond network description ultimately enabling the identification of novel JA pathway regulators.

Discussion

An important unanswered question in plant biology is how multiple signaling pathways interact to coordinate control of growth and development. In this study we have comprehensively characterized the cellular response to the plant hormone JA and generated a network-level understanding of the MYC2/MYC3-regulated JA signaling pathway. We used this to identify several new points at which JA signaling may have cross-regulation with other hormone and stress response pathways in order to prioritize itself. The results increase knowledge of how JA functions in the etiolated seedling, a less well characterized model for JA responses. Moreover, the general principles described here provide a framework for analysis of cross-regulation between hormone and stress signaling pathways. We provide our data in a web-based genome and network browsers to encourage deeper exploration (http://signal.salk.edu/interactome/JA.php, http://neomorph.salk.edu/MYC2).

The major insight provided by our study is that multiple points of crosstalk are likely to exist between JA signaling and other pathways. This was evident from the interactions within the genome regulatory network model and supported by our observation that many (37 to 59%) genes from other hormone signaling pathways are bound by MYC2/3 and JA-regulated. The WRKY family TF ABO3 was identified as a candidate JA response regulator and genetic analyses determined a mutant of the gene was JA hyposensitive. ABO3 is also a regulator of ABA responses [62] suggesting that ABO3 functions in the cross-communication between the JA and ABA pathway. The repressive zinc-finger family TF STZ/ZAT10, working with JAZ8, emerged as a potentially important point of contact with salt and drought stress, as well as the salicylic acid, brassinosteroid and gibberellin hormone signaling pathways. Combined these results illustrate the importance of transcriptional cross-regulation during a JA response in modulating the correct cellular output for the stimuli a plant perceives.

Our multi-omic analysis determined that the master TF MYC2 and its relative MYC3 directly target thousands of JA responsive genes including hundreds of JA-responsive TFs, thereby enabling a robust cascade of transcriptional reprogramming. Secondary TFs downstream of MYC2/3 directly targeted overlapping but distinct cohorts of genes, indicating they have distinct roles within the JA response. This illustrates the complexity of hormone-response genome regulatory programs; we have assayed only a fraction of the JA-responsive TFs and find that any individual JA-responsive gene may be bound by multiple TFs. How the final quantitative output of any individual gene is determined by combinatorial binding of TFs remains a major challenge to address. Achieving this will require analyses at cell-type resolution, resolving differences in TF activity between tissues that would be obscured by our bulk-tissue analyses. We further demonstrated the importance of MYC2/3 target genes in JA responses by analyzing JA root growth phenotypes in mutants of 99 genes strongly targeted by MYC2/3. Mutations in six genes caused clear disruption of JA responses, both hyper and hypo-sensitivity. It is probable that genetic redundancy accounts for a proportion of the mutants not causing phenotype changes. The structure of hormone-response genome regulatory programs will likely differ between cells and tissues and, while our findings can be translated between etiolated and light grown seedlings, exploration of other developmental stage-specific regulatory programs is needed to generalize these findings.

Our study also highlights that many different regulatory mechanisms are utilized by JA to exert its effects upon the cell. Expression of a large number of protein kinases was regulated by MYC2. Consistent with this, substantial MYC2-dependent changes in phosphopeptide abundance occurred in JA-treated seedlings. It is also likely that JA modulates alternative splicing through MYC2. Genes encoding splicing factors were differentially expressed between myc2 and WT and the spliceosome pathway was enriched amongst myc2-dependent JA-regulated phosphopeptides. Accordingly, isoform switch events occurred following JA-treatment. Collectively these findings indicate that investigation of post-transcriptional and post-translational layers of regulation are required to better understand the complexity of JA signaling. The targets of JA-regulated protein kinases are a notable prospect.

Another layer of regulatory complexity within the JA signaling pathway, and within signaling pathways in general, is the presence of multiple feedforward and feedback loops that are activated simultaneously. The interactions between these subnetworks through their kinetics and the strength of their regulatory impact on the broader network is not well understood. For example, we discovered that MYC2 and MYC3 stimulate JA biosynthesis but also target the entire JAZ repressor family from which the majority of members is also transcriptionally activated. Uncoupling these subnetworks would be an effective way to determine how they interact to drive very robust activation of the JA pathway. The combination of our multi-omic framework approach coupled with powerful genetic approaches such as the generation of the jaz decuple mutant [29] should significantly contribute to a better understanding of JA response pathways.

Material & Methods

Plant material and growth conditions

The myc2 mutant in this study is jin1–8 (SALK_061267) [15] and was obtained from the Arabidopsis Biological Resource Center (ABRC). Col-0 MYC2::MYC2-Ypet and Col-0 MYC3::MYC3-Ypet, generated by recombineering, have been described previously [66]. For the generation of all large scale datasets, three-day-old etiolated seedlings were used (Col-0 (WT), myc2, MYC2::MYC2-Ypet, MYC3::MYC3-Ypet). Seedlings were grown in the dark in closed lightproof containers. Gaseous methyl jasmonate treatments for the respective times were performed in these containers as previously described [17] with 1µl of methyl jasmonate (95% purity; Sigma-Aldrich) per 1 litre of container volume dropped on Whatman paper. For the JA-induced root growth inhibition assay, surface-sterilized WT, myc2 and T-DNA mutant seeds (Supplementary Table 19) were grown on agar plates containing LS medium supplemented with or without 20 µM MeJA (392707, Millipore Sigma) for 9 days. Plates were scanned afterwards and root length was measured using ImageJ.

ChIP-seq

Three-day-old etiolated Col-0 MYC2::MYC2-Ypet, Col-0 MYC3::MYC3-Ypet, Col-0 and myc2 seedlings were used for ChIP-seq experiments. ChIP assays were performed as previously described [67]. ChIP-seq assays were conducted with antibodies against H2A.Z (39647, Active Motif), H3K4me3 (04–745, Millipore Sigma) and GFP (11814460001, Millipore Sigma or goat anti-GFP supplied by David Dreschel, Max Planck Institute of Molecular Cell Biology and Genetics). As a negative control, mouse or goat IgG (015–000-003 or 005–000-003, Jackson ImmunoResearch) was used. The respective antibodies or IgG were coupled for 4–6 hour to Protein G Dynabeads (50µl, 10004D, Thermo Fisher Scientific) and subsequently incubated overnight with equal amounts of sonicated chromatin. Beads were washed twice with high salt buffer (50 mM Tris HCl pH 7.4, 150 mM NaCl, 2 mM EDTA, 0.5% Triton X-100) low salt (50 mM Tris HCl pH 7.4, 500 mM NaCl, 2 mM EDTA, 0.5% Triton X-100) and wash buffer (50 mM Tris HCl pH 7.4, 50 mM NaCl, 2 mM EDTA) before samples were de-crosslinked, digested with proteinase K and DNA was precipitated. Sequencing libraries were generated following the manufacturer’s instructions (Illumina). Libraries were sequenced on the Illumina HiSeq 2500 and HiSeq 4000 Sequencing system and sequencing reads were aligned to the TAIR10 genome assembly using Bowtie2 [68].

DAP-seq

DAP-seq assays were carried out as previously described [69,70] using recombinantly expressed ERF1 (AT3G23240, ERF1B, AtERF092), ORA59 (AT1G06160), ATAF1 (AT1G01720), DREB2B (AT3G11020), ZAT18 (AT3G53600), RVE2 (AT5G37260), WRKY51 (AT5G64810), HY5 (AT5G11260) and TCP23 (AT1G35560).

RNA-seq

Three-day-old etiolated seedlings were used for expression analyses. Total RNA was extracted with the RNeasy Plant Mini Kit (74903, Qiagen). cDNA library preparation and subsequent single read sequencing was carried as previously described [3].

RNA-seq analyses

Sequencing reads were quality trimmed using TrimGalore 0.4.5 https://github.com/FelixKrueger/TrimGalore) then aligned to the TAIR10 genome assembly using TopHat 2.1.1 [71]. Reads within gene models were counted using HTSeq [72]. Differentially expressed genes in time series RNA-seq were identified using EdgeR 3.6.2 with a likelihood ratio test (functions glmFit and glmLRT), batch correction Benjamini & Hochberg correction for multiple tests [73]. Differentially expressed genes in the Col-0 versus myc2 mutant RNA-seq were determined using EdgeR 3.18.1 and quasi-likelihood F-tests (function glmQLFit) [74]. Temporal co-regulation of transcripts was determined using the Short Time-Series Expression Miner (STEM) [75]. A minimum correlation coefficient of 0.7 was applied and up to 50 permutations permitted to identify correct cluster/gene matches. Significant clusters were those having a Bonferroni-corrected p-value < 0.05. Full STEM model parameters are given in Supplementary Table 7. Known A. thaliana TFs were identified by reference to PlantTFDB 4.0 [76].

ChIP-seq and DAP-seq analyses

ChIP-seq and DAP-seq sequence reads were mapped to the TAIR10 reference genome using Bowtie 2 v.2–2.0.5 with default parameters [77]. For TF and histone ChIP-seq, we first assessed the quality of the ChIP data by using PhantomPeakQualTools v.2.0 to calculate normalized strand correlation (NSC), relative strand correlation (RSC) and shift size [78]. Enriched binding sites were then identified using MACS2 v.2.1 (options -p 99e-2 --nomodel –shiftsize --down-sample --call-summits) against sequence reads from whole IgG control samples [79]. Subsequent analyses used summits only. Summit lists were filtered with a cut-off of p≤1*10−25 and remaining summits expanded from single nucleotides to 150 nt. Only summits with at least 10% nt overlap between at least two biological replicates were retained. These overlapping summits were merged between replicates using BEDtools v.2.17.0 to give the final set of high confidence binding sites, which were then annotated using ChIPpeakAnno v.2.2.0 to any gene within 500 nt of the center of the summit or, alternatively, the nearest neighbor if there was no gene within 500 nt [80,81]. Venn diagrams were drawn using Venny and Intervene (http://bioinfogp.cnb.csic.es/tools/venny/) [82]. Top-ranked MYC2/3 binding sites were identified by applying IDR to the summits from the two biological replicates that had the greatest number of summits above the MACS2 cut-off of p ≤ 1*10−25. TF binding motifs were determined using the MEME-ChIP webserver with default parameters on the sequences of the high-stringency MYC2 summits [83]. To identify potential partner TFs that may enable indirect MYC2 binding we removed all MYC2 high-stringency summits that contained the MYC2 motif (CACGTG, CATGTG or CACGTT). This was done by scanning them with FIMO set to default parameters (http://meme-suite.org/tools/fimo) against the position weight matrix for the MYC2 motif we previously identified by MEME-ChIP. We then conducted MEME-ChIP analyses on the remaining high-stringency summits as above. The Genome wide Event finding and Motif discovery (GEM) tool [84] was used to identify the target summits in DAP-seq data. Significant enrichments of histone modifications and histone variants were identified with the SICER software [85] using the TAIR10 genome assembly. The Intersect tool from BEDtools [80] was used to identify the genes in the histone ChIP-seq datasets most proximal to the binding sites. The fraction of reads in peak score (FRiP) score was calculated for DAP-Seq and histone ChIP-seq data using BEDtools and SAMtools [85,86]. For both ChIP-seq and DAP-seq gene ontology enrichment was assessed using clusterProfiler with default parameters[87]. Protein domain enrichment was assessed using Thalemine (https://apps.araport.org/thalemine/) with default parameters [88].

Mass spectrometry analysis

Ground untreated/JA-treated Col-0 and myc2 seedlings tissue was ground and lysed in YeastBuster (71186, Millipore Sigma). Proteins (100 µg per sample) were precipitated using methanol- chloroform. Dried pellets were dissolved in 8 M urea, 100 mM triethylammonium bicarbonate (TEAB), reduced with 5 mM tris (2-carboxyethyl) phosphine hydrochloride (TCEP), and alkylated with 50 mM chloroacetamide. Proteins were then trypsin digested overnight at 37 °C. The digested peptides were labeled with TMT10plex™ Isobaric Label Reagent Set (90309, Thermo Fisher Scientific, lot TE264412) and combined. One hundred micrograms (the pre-enriched sample) was fractionated by basic reverse phase (84868, Thermo Fisher Scientific). Phospho-peptides were enriched from the remaining sample (900 µg) using High-Select Fe-NTA Phospho-peptide Enrichment Kit (A32992, Thermo Fisher Scientific). The TMT labeled samples were analyzed on a Fusion Lumos mass spectrometer (Thermo Fisher Scientific). Samples were injected directly onto a 25 cm, 100 μm ID column packed with BEH 1.7 μm C18 resin (186002350, Waters) and subsequently separated at a flow rate of 300 nL/min on a nLC 1200 (LC140, Thermo Fisher Scientific). Buffer A and B were 0.1% formic acid in water and 90% acetonitrile, respectively. A gradient of 1–20% B over 180 min, an increase to 40% B over 30 min, an increase to 100% B over another 20 min and held at 90% B for a final 10 min of washing was used for 240 min total run time. Column was re-equilibrated with 20 μL of buffer A prior to the injection of sample. Peptides were eluted directly from the tip of the column and Nano sprayed directly into the mass spectrometer by application of 2.8 kV voltage at the back of the column. The Lumos was operated in a data dependent mode. Full MS1 scans were collected in the Orbitrap at 120000 resolution. The cycle time was set to 3 s, and within this 3 s the most abundant ions per scan were selected for CID MS/MS in the ion trap. MS3 analysis with multinotch isolation (SPS3) was utilized for detection of TMT reporter ions at 60000 resolution. Monoisotopic precursor selection was enabled and dynamic exclusion was used with exclusion duration of 10 s.

The raw data were analyzed using MaxQuant version 1.6.3.3 [89]. Spectra were searched, using the Andromeda search engine [90] against the Tair10 proteome file entitled “TAIR10_pep_20101214” that was downloaded from the TAIR website (https://www.arabidopsis.org/download/indexauto.jsp?dir=%2Fdownload_files%2FProteins%2FTAIR10_protein_lists) and was complemented with reverse decoy sequences and common contaminants by MaxQuant. Carbamidomethyl cysteine was set as a fixed modification while methionine oxidation and protein N-terminal acetylation were set as variable modifications. For the phoshoproteome “Phosho STY” was also set as a variable modification. The sample type was set to “Reporter Ion MS3” with “10plex TMT selected for both lysine and N-termini”. Digestion parameters were set to “specific” and “Trypsin/P;LysC”. Up to two missed cleavages were allowed. A false discovery rate, calculated in MaxQuant using a target-decoy strategy [91] less than 0.01 at both the peptide spectral match and protein identification level was required. The ‘second peptide’ option identify co-fragmented peptides was not used. Differentially expressed proteins and phospho-sites were identified using PoissonSeq [92] with a q-value cutoff of 0.1. Sample loading normalization was performed before differential expression analysis.

Transcript quantification and identification of isoform switches

Quantification of transcripts was performed using Salmon v0.8.1 in conjunction with the AtRTD2-QUASI transcript reference [93,94] The quasi mapping-based index was built using an auxiliary k-mer hash over k-mers of length 31 (k=31). For quantification, all parameters of Salmon were kept at default, except that the option to correct for the fragment-level GC biases (“–gcBias”) was turned on. The TSIS R package, which is designed for detecting alternatively spliced isoform switch events in time-series transcriptome data, was used to perform the isoform switch analysis [95] Only transcripts whose average transcript per million (TPM) across all time points was >1 were included in the TSIS analysis. The mean expression approach was used to search interaction points. Significant switch events were identified using the following filtering parameters: (1) probability cutoff >0.5; (2) differences cutoff >1; (3) p-value cutoff < 0.05; (4) min time in interval > 1.

Gene regulatory network (GRN) inference

All GRNs were constructed using the Regression Tree Pipeline for Spatial, Temporal, and Replicate data (RTP-STAR) [96,97]. Prior to GRN inference, genes were clustered based on transcriptome, proteome, or phosphoproteome data using Dynamic Time Warping (DTW) (and the dtwclust package in R [98]. These clusters were then used in the RTP-STAR pipeline. For the transcriptome networks, one network was inferred for genes differentially expressed at each time point (8 networks total), and then the networks were combined in a union. For each network, the biological replicates for that individual time point and the 0 h (control) time point were used to infer the network. The sign (activation/repression) of each edge was inferred using all of the time points in the time course. For the proteome and phosphoproteome networks, one network was inferred for genes differentially expressed in any of the comparisons. The biological replicates for all of the (phospho) proteome samples were used to infer the network. The sign of each edge was not inferred as the (phospho) proteome data only consisted of one time point. After the transcriptome, proteome, and phosphoproteome networks were combined in a union, a Network Motif Score (NMS) [99] was calculated to determine the importance of each gene. Feedback loop, feedforward loop, bi-fan, and diamond motifs were used in this score as they have been previously shown to contain genes important for biological processes [100-102]. All motifs were significantly enriched in the combined network compared to a randomly generated network of the same size. The number of times each gene appeared in each motif was counted, the counts were normalized to a scale of 0 to 1, and the counts were summed to calculate the NMS. The higher the NMS, the more functionally important the gene is. All code for RTP-STAR is available at https://github.com/nmclark2/RTP-STAR. The parameters used for all networks in this paper are provided in Supplementary Table 20.

Source Data Fig. 1 Overview of generated, analyzed and used sequencing data for Figure 1

Source Data Fig. 3 Overview of generated, analyzed and used sequencing data for Figure 3

Source Data Fig. 2 Overview of generated, analyzed and used sequencing data for Figure 2

Source Data Fig. 4 Overview of generated, analyzed and used sequencing and (phospho)proteome data for Figure 4

Source Data Fig. 5 Overview of generated, analyzed and used sequencing data for Figure 5

Source Data Extended Data Fig. 1 Overview of generated, analyzed and used sequencing data for Extended Data Figure 1

Source Data Extended Data Fig. 3 Overview of generated, analyzed and used sequencing data for Extended Data Figure 3

Source Data Extended Data Fig. 2 Overview of generated, analyzed and used sequencing and (phospho)proteome data in Extended Data Figure 2

Source Data Extended Data Fig. 4 Overview of generated, analyzed and used sequencing data for Extended Data Figure 4

Source Data Extended Data Fig. 5 Overview of generated, analyzed and used sequencing data for Extended Data Figure 5

Source Data Extended Data Fig. 6 Overview of generated, analyzed and used sequencing data for Extended Data Figure 6

Source Data Extended Data Fig. 8 Overview of generated, analyzed and used sequencing data for Extended Data Figure 8

Source Data Extended Data Fig. 7 Overview of generated, analyzed and used sequencing data for Extended Data Figure 7

Source Data Extended Data Fig. 9 Overview of generated, analyzed and used sequencing and (phospho)proteome data in Extended Data Figure 9

Source Data Extended Data Fig. 10 Overview of generated, analyzed and used sequencing and (phospho)proteome data in Extended Data Figure 10