A persistent major mutation in canonical jasmonate signaling is embedded in an herbivory-elicited gene network

We constructed a 26-parent MAGIC RIL population capturing a large fraction of N. attenuata’s total natural variation to identify the genetic bases of various ecological traits. The selection of 19 of the parental lines was based on the accessions’ metabolic diversity in response to M. sexta OS elicitation (SI Appendix, Fig. S3A; see SI Appendix, Materials and Methods) and to dissect the genetic basis of traits involved in plant–herbivore interactions. A Principal component analysis of the Single-nucleotide polymorphisms (SNP) of the parents shows a distinct population structure, partly correlating with their geographic origins (SI Appendix, Fig. S4D), and the two replicate inbred MAGIC RIL populations, L1 and L2 (sixth generation) form a tight cluster, indicating a homogeneous admixture. This genetic homogeneity allowed for the use of high-resolution quantitative GWA trait mapping without population-structure corrections (63) and after correcting for kinship effects, the model delivered robust trait mappings. Given the importance of the JA-signaling cascade, through its most active metabolite, JA-Ile, in orchestrating this plant’s response to herbivory (64), we quantified JA and its conjugates at 1 h after the OS elicitation of a single leaf from two replicates of 650 MAGIC population RILs and parents growing in a field plantation in Arizona in 2019. Interestingly, JA-Ile and JA-Val levels showed greater variation in the RILs compared to the parental lines (SI Appendix, Fig. S3C), consistent with the loss of regulatory activity of the OS-regulated gene networks on these phytohormones, which had been lost in the creation of the RILs. The subsequent GWA analysis of OS-elicited JA-Ile and JA-Valine (JA-Val) bursts imputed a significant QTL at the same region of chromosome 5, consistent with the two amino acids being conjugated to JA by the same biochemical mechanism (Fig. 1 A and B) (65, 66). Results from a replicate glasshouse (GH) screen of the same MAGIC RIL population imputed the same locus but at a marginal level of significance (SI Appendix, Fig. S5), underscoring the importance of screening plants in their native habitats to elicit ecologically relevant signals.

Mining the assembled genome at chromosome 5, we identified candidate genes, and further eQTL analyses of 350 RILs from the same experiment identified jasmonoyl-L-amino acid synthetase (NaJAR4) as a candidate (SI Appendix, Fig. S6), a known enzyme conjugating Ile and Val to JA (66). The two SNPs showing the highest association with the JA-Ile QTL, and JAR4 eQTL lie adjacent to each other and show strong linkage (SI Appendix, Fig. S6E). A phylogenetic analysis of the NaJAR4 region delineating the variation found in the 26 parents of the MAGIC RIL population (SI Appendix, Fig. S4B), revealed that six variant (V) lines form a monophyletic clade, and the variant allele correlated with a strong reduction in NaJAR4 expression (SI Appendix, Fig. S6 B and D). Additionally, phylogenetic analysis of all NaJAR genes in N. attenuata and well-studied taxa in the Solanaceae placed NaJAR4 in a clade distinct from the other NaJARs in the genome (SI Appendix, Fig. S7) and the OS-elicited expression of NaJAR4 was highest among the NaJARs (SI Appendix, Fig. S6F). We then randomly selected six other lines from the other clades with high NaJAR4 expression and used these as non-variant (NV) lines for further comparisons. Sequence analysis and additional PCRs confirmed a deletion in NaJAR4’s 5′ untranslated region (UTR) in the V lines, which renders the gene inactive (Fig. 1C and SI Appendix, Fig. S6 A, B, and D).

We next scrutinized the geographic origins of V and NV lines to evaluate whether the NaJAR4 mutation was a spatially limited phenomenon. Clustering the accessions based on their genetic distance captured only minor associations with their respective geographic origin (SI Appendix, Fig. S4D); the V lines did not form a distinct cluster but were distributed along the tree (Fig. 1D). This indicates a lack of population substructure associated with the NaJAR4 locus and is consistent with considerable gene flow among the different populations. We focused on a small region in SW Utah from where some V lines originate and examined the frequencies of the variant allele (Fig. 1D and SI Appendix, Table S2) to examine how their distributions changed over time in this annual plant. Using DNA extracted directly from individual or bulk seed collections, we quantified allele frequencies from multiple seed collections made over 9 y from this 400 km² area. The allele frequencies were highly variable both spatially and temporally, with no clear patterns of fixation (Fig. 1D and SI Appendix, Table S2). In the Jackson Spring (JS) population, for example, 8 of 11 (72%) individuals sampled in 1990 harbored the variant allele. However, in a 1999 bulk collection from the same location, this ratio dropped to 23%. Populations within the 6 km² Gold Strike Mine (GSM) area also harbored high frequencies of the variant allele, whereas the variant allele was not found in seeds from several years of collections at the DI ranch (the origin of the Utah (UT) genotype) located just 5 km downstream from the GSM. From these results, we infer that balancing selection through spatiotemporal variation (67) is acting on the NaJAR4 loci, which maintains this allele in natural populations (see SI Appendix, Supplemental text for additional discussion of this point).

To evaluate whether this NaJAR4 mutation could be maintained in natural populations through gene duplication/neofunctionalizations of canonical JA-signaling genes, we used a candidate gene approach to evaluate whether other NaJARs or NaCOIs could provide NaJAR4 mutants with compensatory responses (SI Appendix, Supplemental text). As N. attenuata is an annual plant that rapidly transitions from vegetative to reproductive growth as its habitat dries out (68), the priority of defending leaves rapidly switches to one of defending reproductive parts as plants transition into flowering. Prior research had established that constitutive JA-Ile levels in flower buds are typically 10x higher than those of OS-elicited leaves and that floral defense engages a floral-specific NaJAZi repressor that physically interacts with the NINJA-like protein rather than the canonical NINJA (44). Floral JA-Ile levels were dramatically lower (by 20-fold) in a previously characterized NaJAR6-RNAi line (69) compared to those of wild type (WT) plants, indicating additional tissue-specific conjugation activity of NaJAR6 (SI Appendix, Fig. S9). We then asked whether N. attenuata has evolved a tissue-specific JA-Ile receptor through the second COI in its genome, NaCOI2. We tested this hypothesis in floral tissues by silencing NaCOI2 expression by stable RNAi (SI Appendix, Fig. S10A); however, phytohormone and metabolic profiles of irNaCOI2 lines revealed that NaCOI2 did not buffer the likely fitness consequences of the JA-Ile variation (SI Appendix, Supplemental text and Fig. S11).

Taken together, these results establish a causal association between natural genetic variation in the NaJAR4 locus of the N. attenuata MAGIC population with OS-elicited JA-Ile levels in leaf tissues and a genomic deletion in the NaJAR4 5′ UTR that lowers OS-elicited JA-Ile levels in the V lines. Additionally, the effect of silencing NaJAR4 in floral tissues was much less than that of silencing NaJAR6, and thus the loss of NaJAR4 function could potentially be buffered by intact NaJAR6 activity in floral tissues. However, we needed to evaluate both the biochemical and fitness consequences of the NaJAR4 mutation in greater detail.

NaJAR4 adenylates JA for conjugation with Ile and Val (66). To determine whether the V lines have lower Ile conjugation rates compared to the NV lines due to abrogated NaJAR4 activity, we conducted in vivo JA-conjugation kinetic experiments with isotopically labeled ¹³C-Ile. We utilized previously published NaJAR4- and NaJAR6-RNAi lines produced in the UT genetic background (an NV line) as positive controls for the experiments and quantified JA and its conjugates at 45 and 90 min after simulated herbivory, just before and exactly when OS-elicited JA-Ile peaks (Fig. 2 A and B) (35, 65, 69). We reasoned that if the JA-Ile levels are lower in the V lines due to lower conjugation activity and not the lack of Ile substrate, then supplementing puncture wounds with ¹³C-Ile in OS and subsequently quantifying the ¹³C-labeled JA-Ile levels should reflect activity differences between the two groups. Additionally, to rule out the effect of JA variation, we analyzed the JA-Ile levels with respect to the JA levels, estimating conjugation rates from the slopes of JA-Ile vs. JA regressions (SI Appendix, Table S3).

Both ¹³C-JA-Ile and endogenous JA-Ile conjugation rate estimates at 90 min were significantly lower in NaJAR4 and NaJAR4x6-RNAi lines compared to those of EV plants (Fig. 2 G and H). The V lines phenocopied the NaJAR4-RNAi line in both their ¹³C-JA-Ile and endogenous JA-Ile conjugation rate estimates and had lower conjugation rates compared to NV lines (Fig. 2 I and J). The estimated ¹³C-JA-Ile conjugation rates of NaJAR6-RNAi plants were significantly higher than those of EV plants, which suggests that NaJAR4 activity dominates Ile-conjugating activity at 90 min in leaves. Moreover, the conjugation estimates of NaJAR4- and NaJAR4x6-RNAi plants did not differ, suggesting that NaJAR6 activity contributes little to the total Ile-conjugation activity at 90 min. The similarity of the response in the V lines and the NaJAR4-RNAi lines is consistent with the inference that the deletion in the 5′ UTR of their NaJAR4 gene (Fig. 1C) is responsible for their lower conjugation rates.

Interestingly, at 45 min after OS elicitation, ¹³C-JA-Ile conjugation rate estimates of all NaJAR-RNAi lines were significantly lower compared to those of EV; however, no difference was observed in the endogenous JA-Ile conjugation rate estimate at this early time point (Fig. 2 C and D). The opposite pattern was observed in V vs. NV lines (Fig. 2 E and F) where endogenous JA-Ile conjugation estimates differed but not the exogenous ¹³C-JA-Ile conjugation estimates. From these results, we inferred that a larger apoplastic flux of ¹³C-Ile conjugated by NaJAR6, as opposed to NaJAR4, was responsible, as the former is intact in these lines and dominates the conjugation activity at 45 min. This is in contrast to the lower flux of endogenously synthesized Ile (SI Appendix, Fig. S12A), possibly conjugated by basal NaJAR4 activity which showed a significantly lower conjugation rate in V compared to NV lines. We quantified levels of Ile and its precursor, threonine, and found no evidence of substrate limitations at either time point (SI Appendix, Fig. S12B). In summary, silencing NaJAR4 had a stronger effect on OS-elicited JA-Ile accumulations than did silencing NaJAR6, consistent with previous work (69), and the two enzymes contribute to conjugation activity at different times following OS elicitation in leaves.

We next quantified the whole-plant fitness consequences associated with the NaJAR4 natural variant. Since N. attenuata is a self-fertilizing annual plant, we used life-time seed set as an ecologically relevant proxy for its fitness and MeJA-induced trypsin proteinase inhibitor (TPI) activity was used as a proxy for quantifying defense (SI Appendix, Supporting Information Text). Moreover, JA-signaling is preserved in floral tissues, the fitness consequences of NaJAR4 variation can be accurately estimated in this currency. Additionally, after germination in nature, N. attenuata plants compete for ephemeral resources with conspecifics, which can be either full-siblings or unrelated plants, and relative growth rates are another important parameter in determining fitness outcomes, depending on genetic ancestry. Therefore, we used a two-plant/pot competitive growth system that facilitates the quantification of fitness differences under ecologically realistic conditions (70). The two plant/pot experimental setup allowed us to quantify the relative growth/defense trade-off incurred through TPI by the different genotypes associated with NaJAR4 variation (10, 71). We expected that if plants deficient in NaJAR4 activity (Fig. 3 A and B) are trading off their investment in defense toward growth and fitness, they would acquire more resources than slower-growing but well-defended competitors and be clearly distinguishable in the defense-fitness trait space. In the case of plants harboring the intact NaJAR4 copy, their growth/defense trajectory would be the same as that of their well-defended neighbor.

The MeJA-elicited NaJAR4-deficient V lines (X position in Fig. 3) produced significantly greater seed set (by 37%), while the MeJA-induced TPI activity was significantly lower (by 50%) compared with NV lines (Fig. 3A), which is quantitatively consistent with previously reported results (50). In line with these expectations, the competing UT counterparts (inbred NV line used as a standard wild-type in decades of N. attenuata research, Y) showed MeJA-induced trends similar to those of the NV lines in the fitness-defense trait space. However, in unelicited pairs, their realized fitness matched that of their competitors. These data indicate that the realized fitness gain of the V lines did not come at the expense of the competing (Y) plants, as the MeJA-induced fitness of UT competitors did not differ between the V and NV pairs. The fitness benefit of the V lines is rather a direct consequence of the NaJAR4 deficiency, which reduced TPI activity by 50%, indicating a reduced investment in plastic defenses resulting in a 37% gain in seed set. Previous research with this species in a 2 plant/pot competition setup found that when competing with non-kin, plants dramatically increased their growth rates when not induced with MeJA (70), a result consistent with Hamilton’s rule (72) (see SI Appendix, Supporting Information Text for additional discussion of this point). Furthermore, the relative defense/fitness trade-off effect for V lines is likely underestimated due to the inclusion of two NV lines (AZ and P-149) originating from Arizona, which are impaired in their ability to produce TPIs due to an evolutionarily recent Nonsense Mediated Decay pseudogenation of the TPI gene (73).

The NaJAR-RNAi lines did not show significantly lower capsule counts in response to MeJA treatments (Fig. 3B). Both silenced lines accumulated two-fold lower TPI activity (Fig. 3B) compared to NV plants (Fig. 3A). When compared to their isogenic UT competitors, the NaJAR-RNAi lines produced significantly more (by 42%) capsules. Final stalk lengths, providing a quantitative proxy for growth, revealed that V lines phenocopied the NaJAR4-RNAi lines with comparable effect sizes, together with strong growth reductions in MeJA-induced Y competitors in comparison to the JAR-deficient plants in both groups (Fig. 3C). The estimated standardized effect size between uninduced X and Y plants in V and NaJAR4-RNAi lines remained small (0.32 and 0.26) but increased to 0.84 and 1.56 in response to MeJA-induction, respectively. In contrast, the effect size remained negative in the NV group, reducing to −0.45 upon MeJA-induction, indicating strong competition for resources between X and Y plants. In comparison to the NaJAR4-RNAi line, the MeJA-induced growth reduction was stronger in the NaJAR4/6-RNAi group, where under both control and induced conditions, Y plants were significantly shorter than X plants, as revealed by the larger effect sizes that likely result from the additive effect of silencing NaJAR6.

These findings reveal that abrogating NaJAR activity through RNAi in the UT WT genotype disrupts JA-induced regulation of growth/defense responses. The V lines have presumably adapted to their NaJAR deficient JA-signaling and, when MeJA-induced, maximize their seed set by exploiting this deficiency (Fig. 3A). In contrast, MeJA induction of NaJAR-RNAi lines in the NV UT background did not increase seed set despite consistently displaying a NaJAR deficient growth/defense phenotype, suggesting that their metabolic configurations are not adapted to this deficiency (Fig. 3B). Recombination, as evidenced by gene flow patterns (Fig. 1) and sporadic outcrossing in N. attenuata populations (74, 75), both processes that antagonize selective sweeps, makes it unlikely that these adaptations are in linkage disequilibrium with the NaJAR4 loci. It is more likely that the beneficial adaptations are dispersed across the genome. Overall, V lines lacking NaJAR4 activity realize higher fitness when induced in herbivore-free environments compared to NV lines, displaying similar effects on defense/growth-fitness trade-offs in magnitude, but not induction patterns, as the NaJAR4-RNAi line. From these results, we infer that other unknown compensatory responses have evolved in the V lines, enabling their survival in natural populations. To better understand these compensatory responses, we next characterized the OS-elicited coexpression network in the natural accessions.

In higher plants, JA-Ile is the critical bioactive hormone that allows plants to contextualize an impressive array of physiological and developmental processes to environmental contingencies. The canonical model of JA-signaling posits that responses will be proportional to the amounts of JA-Ile elicited by environmental signals above some thresholds, and the diversity of responses elicited emerges from the spatiotemporal diversity of JA-Ile accumulations, their turnover, and specific JAZ/transcription factor complexes that regulate downstream responses (40). The mutation we identified dramatically lowers JA-Ile accumulations in leaves and significantly affects JA-Ile-associated traits. Previous studies showed that coexpression modules can deliver mechanistically insightful genes, regulating specialized metabolic responses (76, 77). Extending our search from a canonical JA-signaling candidate gene approach, we wanted to identify a conserved OS-induced coexpression network and other hub genes within it, that might be buffering the effect of the NaJAR4 variation in natural accessions. This would enable us to gain mechanistic insights at a molecular level into how the mutation is maintained in these populations, which can otherwise be explained by balancing selection at a functional level of analysis (11). We note that the coexpression network structure regulating effector-triggered immunity can be largely conserved even in JA/SA-deficient Arabidopsis plants (78), where only the amplitude of the network was delayed, due to abrogation of the phytohormonal network.

To this end, we analyzed the OS-elicited leaf transcriptomes of 29 N. attenuata natural accessions, including 25 of the 26 the MAGIC parental lines (SI Appendix, Table S4), using gene coexpression analysis and constructed a putative OS-elicited network. The analysis revealed seven major topologically connected modules of coexpressed genes in which both NaJAR4 and NaJAR6 clustered in the same M7 module (Fig. 4A). Gene ontology (GO) term enrichment analysis for each module revealed the M7 module was enriched for terms like “defense response,” whereas another module, M5, was enriched in terms like “jasmonic acid metabolic process,” “modulation of process of another organism,” and “response to biotic stimuli,” indicating that a defense-response-related gene network was captured in these two modules (SI Appendix, Fig. S13), which additionally, based on gene expression similarity, were also functionally closely related (Fig. 4B).

We focused on the M5 and M7 modules (Fig. 4A) to identify potential hub genes that could be targeted for gene silencing to test the effect on JA-Ile-associated traits with respect to the canonical JA-signaling embedded in the gene network. Hub genes were identified by ranking the genes by 11 network connectivity metrics, and further manual curation of the GO term enrichment results led to a selection of four genes, namely, Ethylene Response Factor 1B (NaERF) and F-BoxAt61340 (NaFB61) from the M5 module, and Glutamate receptor (NaGLR), and F-BoxAt67340 (NaFB67) from the M7 module. These selections were based on the observations that ethylene signaling is strongly elicited by OS in N. attenuata and fine-tunes JA-Ile responses (79, 80); glutamate receptors are likely involved in systemic JA signaling through Ca²⁺ signaling (81); the F-Box NaFB67 contains a conserved Tetratrico Peptide Repeats motif which is involved in protein–protein interactions and has been shown to be essential for phytohormone-mediated signaling (82, 83); and the Arabidopsis homologue of NaFB61 is known to regulate multiple JA- and ABA-responsive genes in stress situations (84). This suggested that the two selected F-Box genes could be involved in OS-elicited hormone sensing, albeit by different mechanisms than that of the canonical COI1. Furthermore, none of these genes co-localized with NaJAR4 loci; thus, it is unlikely that their role in the JA-signaling pathway would be confounded by selection pressures acting on the NaJAR4 locus.

To further understand the correlation patterns between gene expression and downstream metabolite accumulations, we conducted an integrative analysis using multiblock sPLSDA of the OS-elicited transcriptomes, metabolomes, and phytohormones of the 29 natural accessions. This dimension reduction (or latent variable approach) to data integration helps to highlight correlations across multiple datasets, parse-out outliers, or batch effects in individual sets, and provides a robust selection (of genes and metabolites) to facilitate further downstream analyses (85). The analysis revealed that major defense metabolites were linearly correlated with the expression of the four selected hub genes and NaJAR6 (Fig. 4C). NaCOI1 in the M2 module was negatively correlated with component 1 and distinct from the two selected F-Box genes that correlated positively with component 1, consistent with the inference that the selected F-Box genes are not functionally similar to the canonical JA-Ile/COI1 receptor complex. Interestingly, NaJAR4 was not strongly correlated with JA or JA-Ile, or with JA-induced metabolites along the first two sPLSDA components, representing linear dimensionality reduction. This suggested that NaJAR4 expression may share a more complex nonlinear relation within the network, not captured by the sPLSDA. Indeed, Uniform Manifold Approximation and Projection (UMAP) clustering of the module gene expression, which is a nonlinear dimensionality reduction technique, revealed that NaJAR4 clusters well within the M7 module (Fig. 4D) in a two-dimensional space. Taken together, these results indicate that although the NaJAR4 variation has major fitness consequences as observed in glasshouse experiments (Fig. 3), the variation itself is embedded in a complex gene network in which other interconnected genes could compensate for the deleterious effects of this major mutation on N. attenuata’s OS-elicited defense responses.

To test the effects of hub genes on OS-elicited defense responses, we silenced the expression of the four hub genes from the M5 and M7 modules by transforming N. attenuata plants with inverted repeat (ir) RNAi constructs. We conducted the transformations in the UT genetic background (an NV genotype) to produce lines that function with an intact NaJAR4 allele and are thus not confounded by compensatory responses from the connected genes in the network (Fig. 3). Two homozygous T₂ plants were generated harboring a single T-DNA insertion effectively silencing the expression (SI Appendix, Figs. S10B, S14, and S15) of NaERF and NaFB61 (M5), NaGLR, and NaFB67 (M7) genes in individual lines, showing consistent elicitation patterns of induced metabolites in both lines. Thus, one replicate per construct was chosen for comparisons with the V, NV, and NaJAR4/6-RNAi lines. Rosette-stage leaves were elicited with M. sexta OS and sampled after 1 h for phytohormone analyses, and untargeted liquid chromatography–mass spectrometry (LC-MS) analyses were used to characterize their defense responses (Fig. 5).

The OS-elicited JA-Ile levels in V lines were significantly lower (two-fold) than those of the NV lines (Fig. 5B). The lower JA conjugation efficiency in V lines was likely responsible for significantly higher JA levels (1.7-fold) compared to NV (Fig. 5A). The same trend was seen in the NaJAR4- and NaJAR4x6-RNAi lines (Fig. 5B) which also had significantly reduced basal (2.6-fold for both lines) and induced JA-Ile levels (4.4-fold and 8.3-fold, respectively) and tended to have higher induced JA levels compared to those of WT plants (Fig. 5A). JA-induced metabolites, CP, DCS (51), and total diterpene glycosides (DTGs), sectors of specialized metabolism which are strongly regulated by JA-Ile (52) (SI Appendix, Fig. S16), were proportionally down-regulated in the NaJAR4x6-RNAi line, but the defense reductions in the V lines were not as dramatic. Moreover, PLS-DA of mass features revealed that OS-elicited V lines clustered differently than those of the NaJAR4 and NaJAR4x6-RNAi lines, despite all of these lines being comparably impaired in their NaJAR4 activity (Fig. 5F).

Both the irNaGLR and irNaFB67 lines, silenced in genes of the topologically proximal M7 module, had significantly reduced basal JA and JA-Ile levels (Fig. 5 A and B), but did not differ from WT in their OS-elicited levels. Consistent with the canonical role of JA-signaling in defense responses, basal levels of CP and DCS were significantly lower in these lines (2.5-fold and 1.5-fold, respectively). Interestingly, OS-elicited levels were also low (by 1.5-fold and 1.7-fold) and hence not proportional to the lines’ induced JA/JA-Ile levels. Moreover, no effect was observed in total DTGs in these lines, with the JA-Ile-elicited DTG, dimalonylated nicotianoside II, also not being affected (SI Appendix, Fig. S16). PLS-DA of the OS-induced mass features and the relative contributions of the targeted compounds confirmed that irNaGLR and irNaFB67 lines had low constitutive and induced metabolite levels that clustered separately from WT along component 1 (Fig. 5G and SI Appendix, Fig. S17), revealing that variations correlated with canonical JA-signaling in uninduced leaves, but not with the expected induced changes.

For the genes of the topologically distal M5 module, both RNAi-silenced lines showed WT or higher JA/JA-Ile levels after OS elicitation (Fig. 5 A and B). OS-elicited JA-Ile in irNaFB61 plants were significantly higher (1.3-fold), and given that the JA-Ile degradation products (OH-JA, OH-JA-Ile, and COOH-JA-Ile) (86) of both lines were dramatically higher than those of WT plants after OS elicitation (SI Appendix, Figs. S18 and S19), it is likely that the flux of oxylipins through JA-signaling was significantly up-regulated in both lines. Surprisingly, basal levels of JA-regulated DCS (by 5.9-fold and 3.3-fold) and DTG (by 2.6-fold and 2.3-fold) were dramatically up-regulated (Fig. 5 D and E) and PLS-DA analysis of the LC-MS mass features revealed that irNaERF and irNaFB61 lines displayed a “constitutively-on” defense phenotype along component 1 (Fig. 5H and SI Appendix, Fig. S16), which was not proportional to constitutive JA/JA-Ile levels, a complete departure from the canonical JA-signaling model. Furthermore, OS elicitation did not increase the levels of these metabolites. To test the putative constitutively defended phenotype of irNaERF and irNaFB61 plants, we conducted a caterpillar bioassay, in which M. sexta larvae were fed razor-excised leaves from these plants (cut at the base of petioles to minimize wound-induced responses) and larval weight was measured as a proxy for herbivore performance over 10 d. The mass of larvae fed irNaERF and irNaFB61 leaves was significantly lower (by 54 and 53%) than those feeding on WT leaves (SI Appendix, Fig. S20). In summary, disrupting hub genes in a topologically distal module resulted in larger disruptions of the canonical model of JA-elicited defenses.

Overall, these results reveal that OS-elicited defense responses in V lines show minor deviations compared to the NV lines, which are likely buffered by the larger OS-elicited gene network. In contrast, the NaJAR4/6-RNAi lines in the NV background, which have not had the opportunity to evolve a compensatory gene network, exhibit significant differences. While both V lines and NaJAR4/6-RNAi lines exhibit significantly reduced OS-elicited JA-Ile levels and V lines phenocopy the NaJAR4-RNAi in specific defense metabolites controlled by NaJAR4, subtle differences between V lines and NaJAR4-RNAi line are evident in the entire OS-induced metabolome (Fig. 5F). Interestingly, these results closely follow the quantitative predictions of the “omnigene” model (18), where the OS-elicited defense responses in V lines are perhaps buffered by the larger gene network, and that the NaJAR4 variation only affects the proximal M7 module, resulting in minor deviations from the NV lines. In contrast, deviations resulting from silencing of the hub genes of the network, particularly those from module M5, topologically most distal (Fig. 4A) from the module harboring NaJAR4, leads to the most aberrant patterns of OS-elicited metabolites. This can be explained by the amplifying effect of the M7 module, where perturbations in hub genes of the distal M5 module cascade through the highly interconnected gene network, causing large cumulative effects directly on the proximal M7 module resulting in a dramatic change in metabolic response (18). Moreover, the flux of oxylipins through JA-signaling became increasingly unregulated with the topological distance of perturbed genes from NaJAR4, as revealed by the OS-elicited OH-JA-Ile and COOH-JA-Ile levels across the 3 groups (SI Appendix, Fig. S19). The M5 lines exhibit the highest JA-Ile degradation levels, followed by irGLR in the M7 and the JAR variant group. Thus, silencing hub genes in the OS-elicited network results in major aberrations in defense responses, proportional to their topological distance in the network from a major mutation in a core gene of the pathway. This result provides insight into how a gene network can buffer OS-induced defense traits through phenotypic plasticity, allowing JA-Ile natural mutants to maintain their chemical defense when faced with herbivory. This interpretation complements the functional explanation at a different level of analysis, where varying degrees of herbivory or other environmental factors facilitate balancing selection regimes maintaining NaJAR4 variants in natural populations.

This work provides an example of how a single gene mutation in a major signaling pathway can be buffered by a robust GRN and provides the outlines of a mechanistic explanation for how the mutants are able to persist in nature. In such cases, the GRN likely evolve toward a stable state where a single gene mutation becomes less likely to affect individual fitness (87). This robustness however does not come at a cost of polygenic selection, as it can facilitate accumulation of other loss-of-function mutations, thereby increasing genetic diversity and providing additional mechanisms for adaptation in changing environments, elegantly described as “evolutionary capacitance” (88). Thus, a mechanistic understanding of a trait in the context of its larger GRN can provide evolutionary insights into organismal fitness. For example, when grown under laboratory conditions, Arabidopsis JAR1-mutants display an early flowering phenotype but are also hypersensitive to drought, which suggests the hypothesis that the mutation could maximize fitness under drought conditions (89). Whether N. attenuata V lines are similarly adapted to dryer sites, and whether the GRN minimizes the V line’s likely drought susceptibility, are testable hypotheses that could inform the analysis of potential balancing selection regimes.

The melding of evolutionary and molecular biology approaches is leading to a new era of consilience between these disciplines, where the focus is shifting from gene-first to interaction-first models (14). Molecular biology has primarily focused on understanding the genetic bases of a trait, often using model systems to characterize signaling pathways. We have learned that these systems are often complex, with multiple levels of fine-tuning and cross-talk among signaling pathways (90). This complexity can be amortized through an omnigene model, which emphasizes the interconnectedness of genes and their role in shaping complex phenotypes (18). The integration of molecular and evolutionary approaches can furthermore provide a more nuanced understanding of the relationship between genes and phenotypes and the genetic bases underlying complex traits in the context of evolutionary outcomes. This framework accommodates both mechanistic and functional explanations for trait evolution in an organismic context (14). However, it is still critical to examine alternative hypotheses in the same level of analysis to holistically understand biological systems (11). Our inference of both spatiotemporal balancing selection and the OS-elicited gene network facilitating the maintenance of the NaJAR4 mutation are not alternative hypotheses as they are posed at different levels of analysis. Thus, the analysis of the four hub genes from the OS-elicited gene network presented here should be viewed in the context of this broader perspective, as part of a larger network of interacting players that contribute to plant defense and fitness.

The founder lines for 26-parent MAGIC RIL population were selected by screening 422 plants from 74 different natural accessions for various ecologically relevant traits, focusing on metabolic response to herbivory (99). The 24 plants with the most extreme phenotypes (SI Appendix, Fig. S3A) plus the well-characterized Utah (UT, 30× inbred line) and Arizona (AZ, 22× inbred line) accessions were chosen as MAGIC founder lines. This is a common selection criterion for founder lines for multiparent plant population, when the primary focus is to dissect the genetic basis of agronomic traits and improve them (100) and not specifically to sample spatially stratified populations to conduct population genetic analyses. Five rounds of systematic intercrossing ensured the resulting population of 325 plants had genetic contribution from all 26 parents and were subsequently inbred for six generation to achieve approximately 99% homozygosity. Detailed methods and technical considerations are given in SI Appendix, Materials and Methods. A graphical overview of the creation of this MAGIC population is given in ref. 64. All the plants for all crosses were grown in GH of Max Planck Institute for Chemical Ecology (MPI-CE), Jena, Germany, under the conditions described earlier.

Two replicates of the 650 MAGIC RILs population and the 26 parental lines were planted at the Walnut Creek Center for Education and Research (WCCER) in Prescott, AZ. Smoke- and GA₃-treated N. attenuata seeds were germinated in Jiffy-7® peat pellets hydrated with a native soil extract to provide plants with an appropriate microbiome during germination (101) and then planted in the field plot. A standardized M. sexta-specific response in a kinetically defined manner was elicited in the second and third stem leaves of every plant by wounding and applying 20 µL of 1:5 diluted M. sexta OS. The first (lowest) stem leaf was left untreated and harvested as an unelicited control leaf immediately before the elicitation. Leaf samples were harvested at 0 h, 1 h, and 72 h, respectively, and immediately frozen on dry ice, transported on dry ice in dry-ice shipping containers to the MPI-CE, where they were stored in −80 °C freezers until further processing. Ground frozen leaf tissue (1 h OS-elicited) aliquots from 350 plants were shipped to Novogene (HK) Company Ltd. for RNA-Seq and sequenced on Illumina Novaseq.

The details of the field watering and planting scheme as well as experimental and harvesting protocols along with library preparations for sequencing are given in SI Appendix.

R.R., R.H., M.C.S., and I.T.B. designed research; R.R., R.H., K.G., S.M.L., M.C.S., and I.T.B. performed research; K.G., N.R., and I.T.B. contributed new reagents/analytic tools; R.R., R.H., S.M.L., M.C.S., and I.T.B. analyzed data; and R.R. and I.T.B. wrote the paper, with contributions from R.H., K.G., and M.C.S.

The authors declare no competing interest.