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

Rishav Ray, Rayko Halitschke, Klaus Gase, Sabrina Leddy, Meredith

Schuman, Nathalie Rodde, Ian Baldwin

► To cite this version:

Rishav Ray, Rayko Halitschke, Klaus Gase, Sabrina Leddy, Meredith Schuman, et al.. A persistent major mutation in canonical jasmonate signaling is embedded in an herbivory-elicited gene network. Proceedings of the National Academy of Sciences of the United States of America, 2023, 120 (35), 10.1073/pnas.2308500120. hal-04609746

HAL Id: hal-04609746

https://hal.inrae.fr/hal-04609746

Submitted on 12 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# A persistent major mutation in canonical jasmonate signaling is embedded in an herbivory-elicited gene network

Rishav Ray, Rayko Halitschke, Klaus Gase, Sabrina Leddy, Meredith

Schuman, Nathalie Rodde, Ian Baldwin

► To cite this version:

Rishav Ray, Rayko Halitschke, Klaus Gase, Sabrina Leddy, Meredith Schuman, et al.. A persistent major mutation in canonical jasmonate signaling is embedded in an herbivory-elicited gene network. Proceedings of the National Academy of Sciences of the United States of America, 2023, 120 (35), 10.1073/pnas.2308500120. hal-04609746

HAL Id: hal-04609746

https://hal.inrae.fr/hal-04609746v1

Submitted on 12 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

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

Rishav Rayª ID, Rayko Halitschkeª ID, Klaus Gaseᵇ ID, Sabrina M. Leddyᶜ, Meredith C. Schumanᵈ,ᵉ ID, Nathalie Roddef, and Ian T. Baldwinª,¹ ID

Edited by Anne Osbourn, John Innes Centre, Norwich, United Kingdom; received May 21, 2023; accepted July 19, 2023

When insect herbivores attack plants, elicitors from oral secretions and regurgitants (OS) enter wounds during feeding, eliciting defense responses. These generally require plant jasmonate (JA) signaling, specifically, a jasmonoyl-L-isoleucine (JA-Ile) burst, for their activation and are well studied in the native tobacco *Nicotiana attenuata*. We used intraspecific diversity captured in a 26-parent MAGIC population planted in nature and an updated genome assembly to impute natural variation in the OS-elicited JA-Ile burst linked to a mutation in the JA-Ile biosynthetic gene *NaJAR4*. Experiments revealed that *NaJAR4* variants were associated with higher fitness in the absence of herbivores but compromised foliar defenses, with two *NaJAR* homologues (4 and 6) complementing each other spatially and temporally. From decade-long seed collections of natural populations, we uncovered enzymatically inactive variants occurring at variable frequencies, consistent with a balancing selection regime maintaining variants. Integrative analyses of OS-induced transcriptomes and metabolomes of natural accessions revealed that *NaJAR4* is embedded in a nonlinear complex gene coexpression network orchestrating responses to OS, which we tested by silencing four hub genes in two connected coexpressed networks and examining their OS-elicited metabolic responses. Lines silenced in two hub genes (*NaGLR* and *NaFB67*) co-occurring in the *NaJAR4/6* module showed responses proportional to JA-Ile accumulations; two from an adjacent module (*NaERF* and *NaFB61*) had constitutively expressed defenses with high resistance. We infer that mutations with large fitness consequences can persist in natural populations due to compensatory responses from gene networks, which allow for diversification in conserved signaling pathways and are generally consistent with predictions of an omnigene model.

jasmonate signaling | plant defense | MAGIC populations | *NaJAR4* | natural mutations with large fitness effects

Fisher in his seminal 1930 book proposed the geometric model of adaptation to reconcile Darwinian and Mendelian views on genetics, predicting that most mutations are pleiotropic and that evolutionary adaptations, which drive populations to realize fitness optima, result from changes in many genes of small fitness effects (1–3). Yet recent empirical research has identified several examples of single-gene mutations with large fitness effects in natural populations, from genes that influence coat color and burrowing behavior of mice and in flower color that mediates pollination syndromes (4–6). In natural accessions of the model plant, *Arabidopsis thaliana*, major mutations in several genes in the aliphatic glucosinolate defense pathway (7) and in the salicylate (SA)-mediated pathogen resistance locus, *ACD6* (8), are known to influence plant fitness in context- and season-specific manners. Balancing selection regimes, mediating temporally and spatially varying and conflicting selection pressures from different biotic and abiotic selective agents, are thought to maintain these major mutations in natural populations. For example, attack from seasonally variable herbivores and pathogens engage different and sometimes incompatible defense responses (9), and the growth costs associated with some of these defense responses can impair competitive abilities to produce frequency-dependent selection regimes (10). While direct empirical evidence of these balancing selection regimes, particularly from natural environments, is challenging to obtain, patterns of allelic diversity at mutated loci provide compelling population genetic evidence for the existence of these balancing selection regimes, as elegantly shown for the *ACD6* locus (8).

While balancing selection provides an explanation for the maintenance of major mutations at an evolutionary level of analysis (11), other processes operating at a mechanistic level of analysis can buffer the fitness effects of these major mutations at the level of the individual. This phenotypic buffering is commonly explored by examining the potential compensatory effects of paralogs of the mutated loci, or when the mutation occurs in well-characterized pathways, other elements of the pathway are explored for potential fitness-buffering effects. Genome wide association (GWA) studies have consistently found

## Significance

Single-gene mutations with large fitness effects are found in natural populations, contrary to theoretical expectations: understanding how these fitness effects are buffered would complement the functional explanations of balancing selection for their maintenance. Jasmonate (JA) signaling mediates many environmental responses in higher plants; the formation of its central ligand, jasmonoyl-L-isoleucine (JA-Ile), was thought to be invariable. In natural populations of native tobacco plants, we uncover a major mutation in forming JA-Ile; demonstrate its molecular basis, biochemical function, fitness consequences, and find it maintained in several populations over a decade of sampling. When highly connected hub genes were silenced from a coexpression network, JA-mediated defense responses varied with distance from the mutation, consistent with its embedding in a gene regulatory network.

The authors declare no competing interest.

This article is a PNAS Direct Submission.

Copyright © 2023 the Author(s). Published by PNAS. This open access article is distributed under Creative Commons Attribution License 4.0 (CC BY).

Although PNAS asks authors to adhere to United Nations naming conventions for maps (https://www.un.org/geospatial/mapsgeo), our policy is to publish maps as provided by the authors.

¹To whom correspondence may be addressed. Email: baldwin@ice.mpg.de.

This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2308500120/-/DCSupplemental.

Published August 22, 2023.
<figure><img src="image_2.png" /></figure>

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

Rishav Rayª ID, Rayko Halitschkeª ID, Klaus Gaseb ID, Sabrina M. Leddyc, Meredith C. Schuman d,e ID, Nathalie Roddef, and Ian T. Baldwina,¹ ID

Edited by Anne Osbourn, John Innes Centre, Norwich, United Kingdom; received May 21, 2023; accepted July 19, 2023

When insect herbivores attack plants, elicitors from oral secretions and regurgitants (OS) enter wounds during feeding, eliciting defense responses. These generally require plant jasmonate (JA) signaling, specifically, a jasmonoyl-L-isoleucine (JA-Ile) burst, for their activation and are well studied in the native tobacco *Nicotiana attenuata*. We used intraspecific diversity captured in a 26-parent MAGIC population planted in nature and an updated genome assembly to impute natural variation in the OS-elicited JA-Ile burst linked to a mutation in the JA-Ile biosynthetic gene *NaJAR4*. Experiments revealed that *NaJAR4* variants were associated with higher fitness in the absence of herbivores but compromised foliar defenses, with two *NaJAR* homologues (4 and 6) complementing each other spatially and temporally. From decade-long seed collections of natural populations, we uncovered enzymatically inactive variants occurring at variable frequencies, consistent with a balancing selection regime maintaining variants. Integrative analyses of OS-induced transcriptomes and metabolomes of natural accessions revealed that *NaJAR4* is embedded in a nonlinear complex gene coexpression network orchestrating responses to OS, which we tested by silencing four hub genes in two connected coexpressed networks and examining their OS-elicited metabolic responses. Lines silenced in two hub genes (*NaGLR* and *NaFB67*) co-occurring in the *NaJAR4/6* module showed responses proportional to JA-Ile accumulations; two from an adjacent module (*NaERF* and *NaFB61*) had constitutively expressed defenses with high resistance. We infer that mutations with large fitness consequences can persist in natural populations due to compensatory responses from gene networks, which allow for diversification in conserved signaling pathways and are generally consistent with predictions of an omnigene model.

jasmonate signaling | plant defense | MAGIC populations | *NaJAR4* | natural mutations with large fitness effects

Fisher in his seminal 1930 book proposed the geometric model of adaptation to reconcile Darwinian and Mendelian views on genetics, predicting that most mutations are pleiotropic and that evolutionary adaptations, which drive populations to realize fitness optima, result from changes in many genes of small fitness effects (1–3). Yet recent empirical research has identified several examples of single-gene mutations with large fitness effects in natural populations, from genes that influence coat color and burrowing behavior of mice and in flower color that mediates pollination syndromes (4–6). In natural accessions of the model plant, *Arabidopsis thaliana*, major mutations in several genes in the aliphatic glucosinolate defense pathway (7) and in the salicylate (SA)-mediated pathogen resistance locus, *ACD6* (8), are known to influence plant fitness in context- and season-specific manners. Balancing selection regimes, mediating temporally and spatially varying and conflicting selection pressures from different biotic and abiotic selective agents, are thought to maintain these major mutations in natural populations. For example, attack from seasonally variable herbivores and pathogens engage different and sometimes incompatible defense responses (9), and the growth costs associated with some of these defense responses can impair competitive abilities to produce frequency-dependent selection regimes (10). While direct empirical evidence of these balancing selection regimes, particularly from natural environments, is challenging to obtain, patterns of allelic diversity at mutated loci provide compelling population genetic evidence for the existence of these balancing selection regimes, as elegantly shown for the *ACD6* locus (8).

While balancing selection provides an explanation for the maintenance of major mutations at an evolutionary level of analysis (11), other processes operating at a mechanistic level of analysis can buffer the fitness effects of these major mutations at the level of the individual. This phenotypic buffering is commonly explored by examining the potential compensatory effects of paralogs of the mutated loci, or when the mutation occurs in well-characterized pathways, other elements of the pathway are explored for potential fitness-buffering effects. Genome wide association (GWA) studies have consistently found

## Significance

Single-gene mutations with large fitness effects are found in natural populations, contrary to theoretical expectations: understanding how these fitness effects are buffered would complement the functional explanations of balancing selection for their maintenance. Jasmonate (JA) signaling mediates many environmental responses in higher plants; the formation of its central ligand, jasmonoyl-L-isoleucine (JA-Ile), was thought to be invariable. In natural populations of native tobacco plants, we uncover a major mutation in forming JA-Ile; demonstrate its molecular basis, biochemical function, fitness consequences, and find it maintained in several populations over a decade of sampling. When highly connected hub genes were silenced from a coexpression network, JA-mediated defense responses varied with distance from the mutation, consistent with its embedding in a gene regulatory network.

The authors declare no competing interest.

This article is a PNAS Direct Submission.

Copyright © 2023 the Author(s). Published by PNAS. This open access article is distributed under Creative Commons Attribution License 4.0 (CC BY).

Although PNAS asks authors to adhere to United Nations naming conventions for maps (https://www.un.org/geospatial/mapsgeo), our policy is to publish maps as provided by the authors.

¹To whom correspondence may be addressed. Email: baldwin@ice.mpg.de.

This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2308500120/-/DCSupplemental.

Published August 22, 2023.
low heritabilities for most traits, even those associated with large fitness effects, such as complex diseases (12), likely due to pleiotropic effects of the mutated loci. These pleiotropic effects are frequently explored within the context of presumed balancing selection regimes, such as for growth-related traits in defense pathway mutants (7, 8). However, when examined more globally, the pleiotropy is often found to extend more widely, often to traits not clearly related to the proposed balancing selection regime, presumably through gene regulatory networks (GRNs) (7, 13, 14).

Recently, the “omnigene” model (15), framed at the mechanistic level of analysis, was proposed to reconcile the low heritabilities observed in most GWA studies, by proposing that large-effect mutations in a core pathway have additional nonzero effects from “peripheral” genes in the larger GRNs in which traits are embedded. These highly interconnected networks have the “small world” property in which most nodes are connected by few steps, even the “hub” genes in these nodes (16, 17). These tightly coregulated gene networks have modular topologies that distribute phenotypic effects to peripheral nodes and could account for the low heritabilities (18). The omnigene model has been applied to explain genetic and epigenetic bases of several human diseases (19–21) and tested in populations of *Populus nigra* trees (22) using transcriptomic analyses, but manipulative tests of the model are rare. The fine-grained modularity of GRN architecture was demonstrated with perturbation experiments of the different subunits of the conserved Mediator transcriptional coregulator of *Saccharomyces cerevisiae* (23). We note that knocking down the expression of *ACD6* in different *Arabidopsis* accessions results in phenotypes that extend beyond the protein’s canonical function as a transmembrane, ankyrin-repeat domain protein regulating the pathogenesis phytohormone, SA, to another phytohormone, jasmonic acid (JA) (8, 24). Hence, in studies that provide some of the strongest support for balancing selection-maintaining mutations with large fitness effects, the known pleiotropic interactions are generally consistent with the predictions of the omnigene model. Here, we identify a major mutation in the JA signaling pathway in a non-model plant and take inspiration from the omnigene model to explore how this major mutation’s fitness effects could be buffered.

In natural environments, plants are continuously attacked by heterotrophs and compete with each other for fitness-limiting resources; these primary environmental challenges have selected for growth-defense tradeoffs largely mediated by the JA phytohormonal signaling systems that plants use to recognize herbivore attack (25–30). When lepidopteran insect herbivores attack the native tobacco plant, *Nicotiana attenuata*, fatty acid-amino acid conjugates in the herbivores’ oral secretions and regurgitants (OS) are introduced into wounds during feeding to elicit the JA-signaling cascade (31–34), through its most active metabolite, JA-Ile (35, 36) acting as a ligand for coronatine-insensitive 1 (COI1), an FBox protein (37–39) which forms a receptor complex with hypervariable jasmonate ZIM-domain (JAZ) proteins and inositol pentakisphosphate (40–45). Most higher plants deploy this signaling system in response to many environmental challenges (26, 45–48). JA-induced responses activate the production of a plant’s arsenal of diverse specialized metabolites that function as direct and indirect defenses and are well studied in the *N. attenuata-Manduca sexta* natural model plant-herbivore system (48–52) (*SI Appendix, Supplemental text*).

At the heart of JA signaling, sits the conserved JA-Ile COI1 receptor system, the evolutionary history of which has illuminated the GRN architecture of JA signaling. The nonvascular liverwort, *Marchantia polymorpha*, neither produces nor perceives JA-Ile but uses a JA-Ile precursor, dinor-cis/iso-12-oxo-10,15(Z)-phytodienoic acid, as a ligand for its MpCOI1 receptor (53, 54). A single-residue substitution in MpCOI1 switches ligand specificity to JA-Ile, a

more polar metabolite that is more readily transported through the vasculature (53), a trait lacking in liverworts. The JA-Ile receptor complex is thought to be highly conserved among vascular land plants, consistent with the observation that many biotrophic pathogens, such as the model plant pathogen, *Pseudomonas syringae*, produce coronatine, a structural mimic of JA-Ile that activates the COI1-JAZ receptor system and functions a key virulence factor by dampening SA signaling responses (55, 56).

Two decades ago, when we released *NaLOX3*-silenced, JA-deficient *N. attenuata* plants into natural habitats, they attracted an entirely new herbivore species, an *Empoasca* leaf hopper (57), which subsequent research revealed to be a natural insect “bloodhound” for JA mutants (58). This insect guild is prevented from using tobacco plants with intact JA signaling as hosts by a combination of direct and indirect defensive chemistry that is elicited in response to leaf-hopper probing and provides a rarely documented mechanism for non-host resistance against insect herbivores (59). Furthermore, when *N. attenuata* plants, silenced in *NaMYC2* expression, a major JA-responsive hub gene, were released during two field seasons that differed dramatically in their herbivory pressures, the *NaMYC2*-silenced plants realized a significant growth advantage in the year when herbivore pressure was low but were severely compromised during the season of high herbivory pressure and did not survive to reproductive maturity (60). These natural history interactions suggested that natural populations of *N. attenuata* would harbor natural JA-signaling mutants that balance the JA-mediated growth/defense trade-off. Yet our understanding of how the fitness effects of such putative JA-signaling mutants could have been buffered by OS-elicited GRNs was hampered because previous attempts to understand OS-mediated signaling in the *Nicotiana* genus using closely related species did not resolve JA signaling within the OS-elicited GRN (61, 62). We hypothesized that the failure of this previous attempt using interspecific variation was due to the lack of shared OS-mediated JA-signaling GRNs among these closely related species; in other words, we had underestimated the speed with which OS-elicited GRNs evolve (62).

Here, we revisit the OS-elicited gene network with an intraspecific approach harnessing natural variation within *N. attenuata*. Together with a chromosome-length contiguous genome assembly (N50 82 Mb), and a 26-parent Multiparent Advanced Generation Inter Cross (MAGIC) Recombinant Inbred Line (RIL) population of *N. attenuata* (*SI Appendix, Materials and Methods*, Figs. S1 and S2, and Table S1), we imputed a locus associated with variation in OS-elicited JA-Ile bursts. We screened for candidate genes using RNAi-silenced lines of genes in the pathway and, with stable-isotope labeling experiments, confirmed *NaJAR4* as the causal gene. We estimated the frequency of this variant locus in natural populations and characterized the whole-plant fitness consequences of harboring these variants in environments lacking herbivores. After failing to find compensatory buffering effects in jasmonate resistant (JAR) and COI paralogs, we undertook a broader multiomics integrative approach and identified a putative GRN from a gene coexpression network in which the mutation resides, and by manipulating gene expression of four hub genes in two topologically connected modules, we characterized the OS-induced metabolic responses to examine predictions inferred from the omnigene framework (15, 18).

## Results and Discussion

**MAGIC Population and NaJAR4 Natural Variation.** 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
low heritabilities for most traits, even those associated with large fitness effects, such as complex diseases (12), likely due to pleiotropic effects of the mutated loci. These pleiotropic effects are frequently explored within the context of presumed balancing selection regimes, such as for growth-related traits in defense pathway mutants (7, 8). However, when examined more globally, the pleiotropy is often found to extend more widely, often to traits not clearly related to the proposed balancing selection regime, presumably through gene regulatory networks (GRNs) (7, 13, 14).

Recently, the “omnigene” model (15), framed at the mechanistic level of analysis, was proposed to reconcile the low heritabilities observed in most GWA studies, by proposing that large-effect mutations in a core pathway have additional nonzero effects from “peripheral” genes in the larger GRNs in which traits are embedded. These highly interconnected networks have the “small world” property in which most nodes are connected by few steps, even the “hub” genes in these nodes (16, 17). These tightly coregulated gene networks have modular topologies that distribute phenotypic effects to peripheral nodes and could account for the low heritabilities (18). The omnigene model has been applied to explain genetic and epigenetic bases of several human diseases (19–21) and tested in populations of *Populus nigra* trees (22) using transcriptomic analyses, but manipulative tests of the model are rare. The fine-grained modularity of GRN architecture was demonstrated with perturbation experiments of the different subunits of the conserved Mediator transcriptional coregulator of *Saccharomyces cerevisiae* (23). We note that knocking down the expression of *ACD6* in different *Arabidopsis* accessions results in phenotypes that extend beyond the protein's canonical function as a transmembrane, ankyrin-repeat domain protein regulating the pathogenesis phytohormone, SA, to another phytohormone, jasmonic acid (JA) (8, 24). Hence, in studies that provide some of the strongest support for balancing selection-maintaining mutations with large fitness effects, the known pleiotropic interactions are generally consistent with the predictions of the omnigene model. Here, we identify a major mutation in the JA signaling pathway in a non-model plant and take inspiration from the omnigene model to explore how this major mutation's fitness effects could be buffered.

In natural environments, plants are continuously attacked by heterotrophs and compete with each other for fitness-limiting resources; these primary environmental challenges have selected for growth-defense tradeoffs largely mediated by the JA phytohormonal signaling systems that plants use to recognize herbivore attack (25–30). When lepidopteran insect herbivores attack the native tobacco plant, *Nicotiana attenuata*, fatty acid-amino acid conjugates in the herbivores' oral secretions and regurgitants (OS) are introduced into wounds during feeding to elicit the JA-signaling cascade (31–34), through its most active metabolite, JA-Ile (35, 36) acting as a ligand for coronatine-insensitive 1 (COI1), an FBox protein (37–39) which forms a receptor complex with hypervariable jasmonate ZIM-domain (JAZ) proteins and inositol pentakisphosphate (40–45). Most higher plants deploy this signaling system in response to many environmental challenges (26, 45–48). JA-induced responses activate the production of a plant's arsenal of diverse specialized metabolites that function as direct and indirect defenses and are well studied in the *N. attenuata-Manduca sexta* natural model plant-herbivore system (48–52) (*SI Appendix, Supplemental text*).

At the heart of JA signaling, sits the conserved JA-Ile COI1 receptor system, the evolutionary history of which has illuminated the GRN architecture of JA signaling. The nonvascular liverwort, *Marchantia polymorpha*, neither produces nor perceives JA-Ile but uses a JA-Ile precursor, dinor-cis/iso-12-oxo-10,15(Z)-phytodienoic acid, as a ligand for its MpCOI1 receptor (53, 54). A single-residue substitution in MpCOI1 switches ligand specificity to JA-Ile, a

more polar metabolite that is more readily transported through the vasculature (53), a trait lacking in liverworts. The JA-Ile receptor complex is thought to be highly conserved among vascular land plants, consistent with the observation that many biotrophic pathogens, such as the model plant pathogen, *Pseudomonas syringae*, produce coronatine, a structural mimic of JA-Ile that activates the COI1-JAZ receptor system and functions a key virulence factor by dampening SA signaling responses (55, 56).

Two decades ago, when we released *NaLOX3*-silenced, JA-deficient *N. attenuata* plants into natural habitats, they attracted an entirely new herbivore species, an *Empoasca* leaf hopper (57), which subsequent research revealed to be a natural insect “bloodhound” for JA mutants (58). This insect guild is prevented from using tobacco plants with intact JA signaling as hosts by a combination of direct and indirect defensive chemistry that is elicited in response to leaf-hopper probing and provides a rarely documented mechanism for non-host resistance against insect herbivores (59). Furthermore, when *N. attenuata* plants, silenced in *NaMYC2* expression, a major JA-responsive hub gene, were released during two field seasons that differed dramatically in their herbivory pressures, the *NaMYC2*-silenced plants realized a significant growth advantage in the year when herbivore pressure was low but were severely compromised during the season of high herbivory pressure and did not survive to reproductive maturity (60). These natural history interactions suggested that natural populations of *N. attenuata* would harbor natural JA-signaling mutants that balance the JA-mediated growth/defense trade-off. Yet our understanding of how the fitness effects of such putative JA-signaling mutants could have been buffered by OS-elicited GRNs was hampered because previous attempts to understand OS-mediated signaling in the *Nicotiana* genus using closely related species did not resolve JA signaling within the OS-elicited GRN (61, 62). We hypothesized that the failure of this previous attempt using interspecific variation was due to the lack of shared OS-mediated JA-signaling GRNs among these closely related species; in other words, we had underestimated the speed with which OS-elicited GRNs evolve (62).

Here, we revisit the OS-elicited gene network with an intraspecific approach harnessing natural variation within *N. attenuata*. Together with a chromosome-length contiguous genome assembly (N50 82 Mb), and a 26-parent Multiparent Advanced Generation Inter Cross (MAGIC) Recombinant Inbred Line (RIL) population of *N. attenuata* (*SI Appendix, Materials and Methods*, Figs. S1 and S2, and Table S1), we imputed a locus associated with variation in OS-elicited JA-Ile bursts. We screened for candidate genes using RNAi-silenced lines of genes in the pathway and, with stable-isotope labeling experiments, confirmed *NaJAR4* as the causal gene. We estimated the frequency of this variant locus in natural populations and characterized the whole-plant fitness consequences of harboring these variants in environments lacking herbivores. After failing to find compensatory buffering effects in jasmonate resistant (JAR) and COI paralogs, we undertook a broader multiomics integrative approach and identified a putative GRN from a gene coexpression network in which the mutation resides, and by manipulating gene expression of four hub genes in two topologically connected modules, we characterized the OS-induced metabolic responses to examine predictions inferred from the omnigene framework (15, 18).

## Results and Discussion

**MAGIC Population and NaJAR4 Natural Variation.** 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

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** GWA trait mapping of OS-elicited JA-Ile and JA-Val in MAGIC RIL populations and locations of the *NaJAR4* mutation in *N. attenuata* natural populations (A and B) GWA analysis of OS-elicited JA-Ile and JA-Val bursts. (C) The structures of the JAR4 gene showing Whole Genome Sequence (WGS) alignments of V and NV lines and the deletion in the 5' UTR region of *NaJAR4* in the V lines. (D) Genetic distance of natural accessions of *N. attenuata* originating in the SW USA. Blue accessions indicate V lines and orange accessions indicate NV lines. Locations of seed collections in the SW Utah region over 9y from 32 natural populations that were genotyped for the *NaJAR4* variant (blue): Pie diagrams are genotype ratios of either bulked or averages of individual genotypes in each location represented by the black circles, while points reflect each genotyped collection. Details of each collection are listed in SI Appendix, Table S2.</figcaption></figure>
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

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** GWA trait mapping of OS-elicited JA-Ile and JA-Val in MAGIC RIL populations and locations of the *NaJAR4* mutation in *N. attenuata* natural populations (A and B) GWA analysis of OS-elicited JA-Ile and JA-Val bursts. (C) The structures of the JAR4 gene showing Whole Genome Sequence (WGS) alignments of V and NV lines and the deletion in the 5' UTR region of *NaJAR4* in the V lines. (D) Genetic distance of natural accessions of *N. attenuata* originating in the SW USA. Blue accessions indicate V lines and orange accessions indicate NV lines. Locations of seed collections in the SW Utah region over 9y from 32 natural populations that were genotyped for the *NaJAR4* variant (blue): Pie diagrams are genotype ratios of either bulked or averages of individual genotypes in each location represented by the black circles, while points reflect each genotyped collection. Details of each collection are listed in *SI Appendix*, Table S2.</figcaption></figure>
*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.

**OS-Elicited JA-Ile Conjugation Kinetics.** *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. 2A 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. 2G 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. 2I 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. 2C and D). The opposite pattern was observed in V vs. NV lines (Fig. 2E 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.

**Fitness Consequences of the NaJAR4 Variation.** 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
*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.

**OS-Elicited JA-Ile Conjugation Kinetics.** *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. 2A 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. 2G 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. 2I 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. 2C and D). The opposite pattern was observed in V vs. NV lines (Fig. 2E 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.

**Fitness Consequences of the NaJAR4 Variation.** 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
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** JA-Ile synthesis and accumulation kinetics in OS-elicited leaves. (A) *NaJAR4* and/or *NaJAR6* catalyzes the conjugation step of JA-Ile synthesis, conjugating either free endogenously produced isoleucine (Ile) or exogenously applied ¹³C-labeled Ile to JA in response to OS elicitation. (B) Schematic representation of typical JA and JA-Ile bursts in *N. attenuata* leaves immediately after simulated herbivory (OS elicitation). (C and D) Endogenous JA-Ile and ¹³C-JA-Ile vs. endogenous JA accumulations at 45 min after OS elicitation in three RNAi knockdown lines and empty vector (EV) as controls (E and F) Endogenous JA-Ile and ¹³C-JA-Ile vs. JA accumulation in 6 V and 6 NV MAGIC parental accessions at 45 min after OS elicitation. (G and H) Endogenous JA-Ile and ¹³C-JA-Ile vs. JA accumulations at 90 min after OS elicitation in the different RNAi lines. (I and J) Endogenous JA-Ile and ¹³C-JA-Ile vs. JA accumulations in 6 V and 6 NV MAGIC parental accessions at 90 min after OS elicitation. The values are the differences in the slopes of regression lines calculated with the emmeans package in R and are listed in SI Appendix, Table S3; an asterisk denotes significant differences.</figcaption></figure>

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
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** JA-Ile synthesis and accumulation kinetics in OS-elicited leaves. (A) *NaJAR4* and/or *NaJAR6* catalyzes the conjugation step of JA-Ile synthesis, conjugating either free endogenously produced isoleucine (Ile) or exogenously applied ¹³C-labeled Ile to JA in response to OS elicitation. (B) Schematic representation of typical JA and JA-Ile bursts in *N. attenuata* leaves immediately after simulated herbivory (OS elicitation). (C and D) Endogenous JA-Ile and ¹³C-JA-Ile vs. endogenous JA accumulations at 45 min after OS elicitation in three RNAi knockdown lines and empty vector (EV) as controls (E and F) Endogenous JA-Ile and ¹³C-JA-Ile vs. JA accumulation in 6 V and 6 NV MAGIC parental accessions at 45 min after OS elicitation. (G and H) Endogenous JA-Ile and ¹³C-JA-Ile vs. JA accumulations at 90 min after OS elicitation in the different RNAi lines. (I and J) Endogenous JA-Ile and ¹³C-JA-Ile vs. JA accumulations in 6 V and 6 NV MAGIC parental accessions at 90 min after OS elicitation. The values are the differences in the slopes of regression lines calculated with the emmeans package in R and are listed in SI Appendix, Table S3; an asterisk denotes significant differences.</figcaption></figure>

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. 3A 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
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Growth and fitness consequences of the *NaJAR4* natural variation. (A and B) An ecologically realistic competition experiment was used to estimate fitness-defense consequences of the *NaJAR4* variant. The middle schematics describes two size-matched 10-d-old seedlings that were planted in a single pot (n = 10), of which one (Y-plant) was always a UT-WT plant, with the other competing neighbor (X-plant) being, either (A) a non-isogenic V or NV WT plant; or (B) an isogenic RNAi-line silenced in either *NaJAR4* or both *NaJAR4* and *NaJAR6* (*NaJAR4/6*). Bivariate plots describe the fitness/defense consequences of methyl jasmonate (MeJA)-elicitation during reproductive maturation in fitness (lifetime capsule numbers: Y axis) and leaf defense (fitness-costly trypsin proteinase inhibitor (TPI) defenses: X axis). Traits of the X competitors are shown in the left-side panels while those of the Y competitors are shown in the right-side panels. (C) Flowering stalk lengths of X (solid lines) and Y plants (dashed lines) measured at the final harvest, 37 d after treatments. Error bars indicate ± 95% CIs. The numbers on top of error bars indicate standardized effect sizes estimated by Cohen's D method.</figcaption></figure>

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.

**Multi-omic Integration of OS-Elicited *N. attenuata* 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
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Growth and fitness consequences of the *NaJAR4* natural variation. (A and B) An ecologically realistic competition experiment was used to estimate fitness-defense consequences of the *NaJAR4* variant. The middle schematics describes two size-matched 10-d-old seedlings that were planted in a single pot (n = 10), of which one (Y-plant) was always a UT-WT plant, with the other competing neighbor (X-plant) being, either (A) a non-isogenic V or NV WT plant; or (B) an isogenic RNAi-line silenced in either *NaJAR4* or both *NaJAR4* and *NaJAR6* (*NaJAR4/6*). Bivariate plots describe the fitness/defense consequences of methyl jasmonate (MeJA)-elicitation during reproductive maturation in fitness (lifetime capsule numbers: Y axis) and leaf defense (fitness-costly trypsin proteinase inhibitor (TPI) defenses: X axis). Traits of the X competitors are shown in the left-side panels while those of the Y competitors are shown in the right-side panels. (C) Flowering stalk lengths of X (solid lines) and Y plants (dashed lines) measured at the final harvest, 37 d after treatments. Error bars indicate ± 95% CIs. The numbers on top of error bars indicate standardized effect sizes estimated by Cohen's D method.</figcaption></figure>

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.

**Multi-omic Integration of OS-Elicited *N. attenuata* 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 founded 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

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Coexpression analysis of OS-elicited transcriptomes, metabolomes, and phytohormones of *N. attenuata* natural accessions (A) Coexpression analysis of OS-elicited transcriptomes of *N. attenuata* natural accessions revealed seven gene modules. Each node in the dendrogram represents a single gene, and the number of genes in each module is M1 = 499, M2 = 93, M3 = 392, M4 = 214, M5 = 187, M6 = 89, and M7 = 259. (B) Module similarity was computed using the (1-similarity) matrix obtained in the previous step. (C) Correlation plot of integrative analysis of the transcriptome, metabolome, and phytohormones performed using multiblock sparse partial least squares discriminant analysis (sPLSDA), showing the correlation of each feature against two latent variables, namely components 1 and 2. The two vectors C and W+OS are diametrically opposite, indicating that these two components adequately captured the variation due to the OS elicitation. (D) UMAP clustering was performed on the OS-induced transcriptomes of *N. attenuata* natural accessions, and colored based on the modules returned by the coexpression analysis in panel A.</figcaption></figure>
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 founded 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

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Coexpression analysis of OS-elicited transcriptomes, metabolomes, and phytohormones of *N. attenuata* natural accessions (A) Coexpression analysis of OS-elicited transcriptomes of *N. attenuata* natural accessions revealed seven gene modules. Each node in the dendrogram represents a single gene, and the number of genes in each module is M1 = 499, M2 = 93, M3 = 392, M4 = 214, M5 = 187, M6 = 89, and M7 = 259. (B) Module similarity was computed using the (1-similarity) matrix obtained in the previous step. (C) Correlation plot of integrative analysis of the transcriptome, metabolome, and phytohormones performed using multiblock sparse partial least squares discriminant analysis (sPLSDA), showing the correlation of each feature against two latent variables, namely components 1 and 2. The two vectors C and W+OS are diametrically opposite, indicating that these two components adequately captured the variation due to the OS elicitation. (D) UMAP clustering was performed on the OS-induced transcriptomes of *N. attenuata* natural accessions, and colored based on the modules returned by the coexpression analysis in panel A.</figcaption></figure>
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.

**Metabolomic Analysis of Silenced Hub Genes.** 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. 5A 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. 5A 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. 5D 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
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.

**Metabolomic Analysis of Silenced Hub Genes.** 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. 5A 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, dimalkonylated 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. 5A 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. 5D 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
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Metabolic characterization of four RNAi lines silenced in the expression of hub genes, two *NaJAR*-RNAi lines, V and NV plants (A–E) JA, JA-Ile, caffeoylputrescine (CP), dicaffeoylspermidine (DCS), and total diterpene glycoside (DTG) levels in unelicited (0 h) control plants and 1 h after OS elicitation in *irNaERF*, *irNaFB61*, *irNaGLR*, and *irNaFB67*, *NaJAR4/6*-RNAi and the NV and V lines. Metabolite levels were normalized against the common UT genotype that was included in all experiments and used as a genotypic control for all RNAi lines. The module dendrogram on the top is a subset of the coexpression clustering from Fig. 4A (with reversed branch ordering for better visualization), indicating the relationships of the genes within their modules and their respective silenced lines. (F) Partial least square discriminant analysis (PLS-DA) projection in two components of control and OS-elicited metabolomes of the *NaJAR4* variant group. The WT lines are circled within the larger NV group. (G) PLS-DA biplot of control and OS-elicited module 7 lines. (H) PLS-DA biplot of control and OS-elicited metabolomes of module 5 lines. See SI Appendix, Fig. S16 for quantifications of the JA-regulated sector of the potent HG-DTG defense metabolites. Error bars are SEs. Significance was calculated by ANOVAs followed by Tukey's Honestly significant difference tests comparing against WT in the M5 and M7 groups and the *NaJAR*-RNAi group, and NV in natural variant, across two treatment groups (*P*-values < 0.05; ** < 0.01; *** < 0.001; **** < 0.0001).</figcaption></figure>

genes of the distal M5 module cascade through the highly inter-
connected 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
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Metabolic characterization of four RNAi lines silenced in the expression of hub genes, two *NaJAR*-RNAi lines, V and NV plants (A–E) JA, JA-lle, caffeoylputrescine (CP), dicaffeoylspermidine (DCS), and total diterpene glycoside (DTG) levels in unelicited (0 h) control plants and 1 h after OS elicitation in *irNaERF*, *irNaFB61*, *irNaGLR*, and *irNaFB67*, *NaJAR4/6*-RNAi and the NV and V lines. Metabolite levels were normalized against the common UT genotype that was included in all experiments and used as a genotypic control for all RNAi lines. The module dendrogram on the top is a subset of the coexpression clustering from Fig. 4A (with reversed branch ordering for better visualization), indicating the relationships of the genes within their modules and their respective silenced lines. (F) Partial least square discriminant analysis (PLS-DA) projection in two components of control and OS-elicited metabolomes of the *NaJAR4* variant group. The WT lines are circled within the larger NV group. (G) PLS-DA biplot of control and OS-elicited module 7 lines. (H) PLS-DA biplot of control and OS-elicited metabolomes of module 5 lines. See SI Appendix, Fig. S16 for quantifications of the JA-regulated sector of the potent HG-DTG defense metabolites. Error bars are SEs. Significance was calculated by ANOVAs followed by Tukey's Honestly significant difference tests comparing against WT in the M5 and M7 groups and the *NaJAR*-RNAi group, and NV in natural variant, across two treatment groups (*P*-values < 0.05; ** < 0.01; *** < 0.001; **** < 0.0001).</figcaption></figure>

genes of the distal M5 module cascade through the highly inter-
connected 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-elicitèd 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-elicitèd 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.

**Proximate and Ultimate Levels of Explanation of Major Mutations.**

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-elicitèd 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-elicitèd 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.

**Conclusions**

The identification of natural variation in JA-signaling in *N. attenuata* brings full circle two-decades of research with this system in which these natural mutants were anticipated based on the herbivore community's response to field releases of lines from a reverse-genetics toolbox (57). The dramatic response of leafhoppers to field releases of JA-deficient RNAi lines, combined with sporadic observations of leafhopper damage in wild populations (seemingly corresponding to JA deficiency in damaged plants, (58)) suggested that natural populations of *N. attenuata* must harbor JA-signaling mutants. This is what we report here by developing a forward-genetics toolbox that captured most of the species genetic variability in ecological traits. An updated contiguous reference genome allowed natural variation in OS-elicitèd JA-Ile bursts in field plantings of the MAGIC RIL population to be imputed to the *NaJAR4* locus. We mined our extensive seed collections of natural *N. attenuata* populations to describe the spatial and temporal persistence of the mutation in natural populations and its fitness effects in controlled environments. We further show how this large-effect single-gene mutation could be maintained in natural populations as a consequence of being embedded in stabilizing layers of an OS-induced coexpression network, which, together with overlapping homolog functions, allow mutants to maximize their fitness in varying environments.

It is likely that natural populations of other plant species harbor such major mutations in other canonical signaling cascades, despite large fitness effects of these mutations. Indeed, two commonly used models in plant biology are consistent with this inference: The laboratory strain of *N. benthamiana*, a workhorse for plant biotechnology and synthetic biology, is a natural *RdR1* mutant, highly deficient in the RNAi mechanisms required for viral resistance (91–93); the reference strain of *A. thaliana*, Col-0, is a natural hydroperoxide lyase mutant (94, 95) deficient in green leaf volatile production. The coarse-grained perturbation analysis presented here of an herbivory-induced gene network stabilizing a natural mutation in canonical JA-signaling reveals how much more remains to be discovered in the complex signaling systems that plants use to optimize their fitness in heterogeneous natural environments.

The results of this work also harmonize with recent work on *Caenorhabditis elegans* natural populations (96). Interindividual metabolic variation in amino acid (again, Ile and Val) conjugation to 3-hydroxypropionate in the species, resulting from mutations in propionyl-CoA carboxylase activity, which in humans causes a rare, but severe, metabolic disorder—propionic acidemia—is maintained on a small spatial scale in Hawaiian island populations of worms. It is likely that the compensatory transcriptomic and metabolic responses result from the pleiotropic effect of “peripheral” genes, which maintain these large-effect mutations in the Ile/Val conjugation steps in both natural populations of plants and nematodes. In nematodes, the biochemical details of these compensatory responses, elegantly characterized as “shunts within shunts” are known; in plants, these mechanisms remain unresolved in the complicated regulatory circuits of the 446-gene OS-elicitèd gene network that buffers canonical JA-signaling from perturbations. Clearly, more fine-grained perturbations of network topologies, as pioneered in *S. cerevisiae* (23), are needed to evaluate whether and how these complex gene networks and metabolic networks maintain large-effect mutations in natural populations.

The commonalities of metabolic variation uncovered in natural populations of nematodes and plants suggest that the concepts of precision medicine (97) apply equally well to animals and plants. In plants, such concepts have been in practice through “precision” pest
OS-elicted 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 ir*GLR* in the M7 and the JAR variant group. Thus, silencing hub genes in the OS-elicted 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.

**Proximate and Ultimate Levels of Explanation of Major Mutations.**

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-elicted 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-elicted 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.

**Conclusions**

The identification of natural variation in JA-signaling in *N. attenuata* brings full circle two-decades of research with this system in which these natural mutants were anticipated based on the herbivore community's response to field releases of lines from a reverse-genetics toolbox (57). The dramatic response of leafhoppers to field releases of JA-deficient RNAi lines, combined with sporadic observations of leafhopper damage in wild populations (seemingly corresponding to JA deficiency in damaged plants, (58)) suggested that natural populations of *N. attenuata* must harbor JA-signaling mutants. This is what we report here by developing a forward-genetics toolbox that captured most of the species genetic variability in ecological traits. An updated contiguous reference genome allowed natural variation in OS-elicted JA-Ile bursts in field plantings of the MAGIC RIL population to be imputed to the *NaJAR4* locus. We mined our extensive seed collections of natural *N. attenuata* populations to describe the spatial and temporal persistence of the mutation in natural populations and its fitness effects in controlled environments. We further show how this large-effect single-gene mutation could be maintained in natural populations as a consequence of being embedded in stabilizing layers of an OS-induced coexpression network, which, together with overlapping homolog functions, allow mutants to maximize their fitness in varying environments.

It is likely that natural populations of other plant species harbor such major mutations in other canonical signaling cascades, despite large fitness effects of these mutations. Indeed, two commonly used models in plant biology are consistent with this inference: The laboratory strain of *N. benthamiana*, a workhorse for plant biotechnology and synthetic biology, is a natural *RdR1* mutant, highly deficient in the RNAi mechanisms required for viral resistance (91–93); the reference strain of *A. thaliana*, Col-0, is a natural hydroperoxide lyase mutant (94, 95) deficient in green leaf volatile production. The coarse-grained perturbation analysis presented here of an herbivory-induced gene network stabilizing a natural mutation in canonical JA-signaling reveals how much more remains to be discovered in the complex signaling systems that plants use to optimize their fitness in heterogeneous natural environments.

The results of this work also harmonize with recent work on *Caenorhabditis elegans* natural populations (96). Interindividual metabolic variation in amino acid (again, Ile and Val) conjugation to 3-hydroxypropionate in the species, resulting from mutations in propionyl-CoA carboxylase activity, which in humans causes a rare, but severe, metabolic disorder—propionic acidemia—is maintained on a small spatial scale in Hawaiian island populations of worms. It is likely that the compensatory transcriptomic and metabolic responses result from the pleiotropic effect of “peripheral” genes, which maintain these large-effect mutations in the Ile/Val conjugation steps in both natural populations of plants and nematodes. In nematodes, the biochemical details of these compensatory responses, elegantly characterized as “shunts within shunts” are known; in plants, these mechanisms remain unresolved in the complicated regulatory circuits of the 446-gene OS-elicted gene network that buffers canonical JA-signaling from perturbations. Clearly, more fine-grained perturbations of network topologies, as pioneered in *S. cerevisiae* (23), are needed to evaluate whether and how these complex gene networks and metabolic networks maintain large-effect mutations in natural populations.

The commonalities of metabolic variation uncovered in natural populations of nematodes and plants suggest that the concepts of precision medicine (97) apply equally well to animals and plants. In plants, such concepts have been in practice through “precision” pest
control and disease resistance for eons. These similarities also provide valuable guidance to the conservation stewards of a planet rapidly losing its unexplored biodiversity (98): Identifying and preserving hotspots of genetic variance should be of the highest priority.

## Materials and Methods

The description of plant materials, cultivation, plant treatment and sampling in the GH together with subsequent transcriptome sequencing and analysis, metabolite analyses, genetic, genomic analyses, and larval growth assays are all based on published procedures and detailed in *SI Appendix*. As the MAGIC population and the field screening and its analysis motivated our focus on *NaJAR4* variation and its characterization, those methods are briefly described here, with additional details in *SI Appendix*.

**Creation of MAGIC RIL Population.** 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.

**Field Experiments and RNA-Seq from Field Samples.** 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 unelicitd 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*.

**Data, Materials, and Software Availability.** All raw sequence and optical man data were deposited in the NCBI SRA database. The Whole Genome project has been deposited at DDBJ/ENA/GenBank under the accession JAPQWU000000000 with BioProject id PRJNA900029 (102). The version described in this paper is version JAPQWU010000000. The MAGIC parents and RIL WGS are available under BioProject id PRJNA907539 (103). The two transcriptome datasets are available under BioProject PRJNA898122(104) and PRJNA911202 (105).

**ACKNOWLEDGMENTS.** We thank all the members of the Department of Molecular Ecology who contributed to the decade-long MAGIC project. Detailed acknowledgments with individual contributions are listed in *SI Appendix*. This work was supported by the Max Planck Society's core funding of I.T.B.'s department, Advanced grant no. 293926 of the European Research Council to I.T.B., and the Deutsche Forschungsgemeinschaft-SFB 1127/2 ChemBioSys-239748522 grant to I.T.B.

Author affiliations: *Department of Molecular Ecology, Max Planck Institute for Chemical Ecology, 07745 Jena, Germany; *Department of Natural Product Biosynthesis, Max Planck Institute for Chemical Ecology, 07745 Jena, Germany; *Department of Molecular Biology and Genetics, Cornell University, Ithaca, NY 14850; *Department of Geography, University of Zurich, 8006 Zurich, Switzerland; *Department of Chemistry, University of Zurich, 8006 Zurich, Switzerland; and *Institut national de recherche pour l'agriculture, l'alimentation et l'environnement, Centre National de Resources Génomiques Végétales, French Plant Genomic Resource Center, Castanet Tolasan F-31326, France

Author contributions: 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.

1. R. A. Fisher, *The Genetical Theory of Natural Selection* (Oxford University Press, 1930).

2. C. Darwin, *On the Origin of Species* (Routledge, 1859).

3. N. H. Barton, A. M. Etheridge, A. Vêber, The infinitesimal model: Definition, derivation, and implications. *Theor. Popul. Biol.* **118**, 50–73 (2017).

4. R. D. Barrett et al., Linking a mutation to survival in wild mice. *Science* **363**, 499–504 (2019).

5. C. R. Linnen et al., Adaptive evolution of multiple traits through multiple mutations at a single gene. *Science* **339**, 1312–1316 (2013).

6. A. E. Berardi et al., Complex evolution of novel red floral color in *Petunia*. *Plant Cell* **33**, 2273–2295 (2012).

7. R. Kerwin et al., Natural genetic variation in *Arabidopsis thaliana* defense metabolism genes modulates field fitness. *eLife* **4**, e05604 (2015).

8. M. Todesco et al., Natural allelic variation underlying a major fitness trade-off in *Arabidopsis thaliana*. *Nature* **465**, 632–636 (2010).

9. C. Ponzio, B. T. Weldegergis, M. Dicke, R. Gols, Compatible and incompatible pathogen-plant interactions differentially affect plant volatile emissions and the attraction of parasitoid wasps. *Funct. Ecol.* **30**, 1779–1789 (2016).

10. T. Züst, A. A. Agrawal, Trade-offs between plant growth and defense against insect herbivory: An emerging mechanistic synthesis. *Annu. Rev. Plant. Biol.* **68**, 513–534 (2017).

11. P. W. Sherman, The levels of analysis. *Anim. Behav.* **36**, 616–619 (1988).

12. T. A. Manolio et al., Finding the missing heritability of complex diseases. *Nature* **461**, 747–753 (2009).

13. I. Barrio-Hernandez et al., Network expansion of genetic associations defines a pleiotropy map of human cell biology. *Nat. Genet.* **55**, 1–10 (2023).

14. M. Fagny, F. Austerlitz, Polygenic adaptation: Integrating population genetics and gene regulatory networks. *Trends Genet.* **37**, 631–638 (2021).

15. E. A. Boyle, Y. I. Li, J. K. Pritchard, An expanded view of complex traits: From polygenic to omnigenic. *Cell* **169**, 1177–1186 (2017).

16. D. J. Watts, S. H. Strogatz, Collective dynamics of "small-world" networks. *Nature* **393**, 440–442 (1998).

17. S. H. Strogatz, Exploring complex networks. *Nature* **410**, 268–276 (2001).

18. X. Liu, Y. I. Li, J. K. Pritchard, Trans effects on gene expression can drive omnigenic inheritance. *Cell* **177**, 1022–1034.e6 (2019).

19. E. O. Paull et al., A modular master regulator landscape controls cancer transcriptional identity. *Cell* **184**, 334–351 (2021).

20. J. Peedicyayil, D. R. Grayson, An epigenetic basis for an omnigenic model of psychiatric disorders. *Theor. Biol.* **443**, 52 (2018).

21. P. Van Der Harst, N. Verweij, Identification of 64 novel genetic loci provides an expanded view on the genetic architecture of coronary artery disease. *Circ. Res.* **122**, 433–443 (2018).

22. A. Chateigner et al., Gene expression predictions and networks in natural populations supports the omnigenic theory. *BMC Genomics* **21**, 1–16 (2020).

23. P. Kemmeren et al., Large-scale genetic perturbations reveal regulatory networks and an abundance of gene-specific repressors. *Cell* **157**, 740–752 (2014).

24. D. N. Rate, J. V. Cuenca, G. R. Bowman, D. S. Guttman, J. T. Greenberg, The gain-of-function *Arabidopsis acd6* mutant reveals novel regulation and function of the salicylic acid signaling pathway in controlling cell death, defenses, and cell growth. *Plant Cell* **11**, 1695–1708 (1999).

25. M. Erb, D. J. Kliefenstein, Plant secondary metabolites as defenses, regulators, and primary metabolites: The blurred functional trichotomy. *Plant Physiol.* **184**, 39–52 (2020).

26. A. Kessler, I. T. Baldwin, Plant responses to insect herbivory: The emerging molecular analysis. *Annu. Rev. Plant. Biol.* **53**, 299–328 (2002).

27. D. Cipollini, D. Walters, C. Voelckel, "Costs of resistance in plants: from theory to evidence" in *Annual Plant Reviews, Insect-Plant Interactions*, C. Voelckel, G. Jander, Eds. (John Wiley & Sons, Ltd, 2014), pp. 263–307.

28. P. D. Coley, J. P. Bryant, F. S. Chapin III, Resource availability and plant antiherbivore defense. *Science* **230**, 895–899 (1985).

29. D. A. Herms, W. J. Mattson, The dilemma of plants: To grow or defend. *Q. Rev. Biol.* **67**, 283–335 (1992).

30. M. Dicke, I. T. Baldwin, The evolutionary context for herbivore-induced plant volatiles: Beyond the "cry for help". *Trends Plant Sci.* **15**, 167–175 (2010).

31. N. Yoshinaga et al., Fatty acid-amino acid conjugates diversification in lepidopteran caterpillars. *J. Chem. Ecol.* **36**, 319–325 (2010).

32. R. Halitschke, U. Schittko, G. Pohnert, W. Boland, I. T. Baldwin, Molecular interactions between the specialist herbivore *Manduca sexta* (Lepidoptera, Sphingidae) and its natural host *Nicotiana attenuata*. *J. Insect. Physiol.* **55**, 171–177 (2009).
control and disease resistance for eons. These similarities also provide valuable guidance to the conservation stewards of a planet rapidly losing its unexplored biodiversity (98): Identifying and preserving hotspots of genetic variance should be of the highest priority.

## Materials and Methods

The description of plant materials, cultivation, plant treatment and sampling in the GH together with subsequent transcriptome sequencing and analysis, metabolite analyses, genetic, genomic analyses, and larval growth assays are all based on published procedures and detailed in *SI Appendix*. As the MAGIC population and the field screening and its analysis motivated our focus on *NaJAR4* variation and its characterization, those methods are briefly described here, with additional details in *SI Appendix*.

**Creation of MAGIC RIL Population.** 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.

**Field Experiments and RNA-Seq from Field Samples.** 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

1. R. A. Fisher, *The Genetical Theory of Natural Selection* (Oxford University Press, 1930).

2. C. Darwin, *On the Origin of Species* (Routledge, 1859).

3. N. H. Barton, A. M. Etheridge, A. Vêber, The infinitesimal model: Definition, derivation, and implications. *Theor. Popul. Biol.* **118**, 50–73 (2017).

4. R. D. Barrett et al., Linking a mutation to survival in wild mice. *Science* **363**, 499–504 (2019).

5. C. R. Linnen et al., Adaptive evolution of multiple traits through multiple mutations at a single gene. *Science* **339**, 1312–1316 (2013).

6. A. E. Berardi et al., Complex evolution of novel red floral color in *Petunia*. *Plant Cell* **33**, 2273–2295 (2012).

7. R. Kerwin et al., Natural genetic variation in *Arabidopsis thaliana* defense metabolism genes modulates field fitness. *eLife* **4**, e05604 (2015).

8. M. Todesco et al., Natural allelic variation underlying a major fitness trade-off in *Arabidopsis thaliana*. *Nature* **465**, 632–636 (2010).

9. C. Ponzio, B. T. Weldegergis, M. Dicke, R. Gols, Compatible and incompatible pathogen-plant interactions differentially affect plant volatile emissions and the attraction of parasitoid wasps. *Funct. Ecol.* **30**, 1779–1789 (2016).

10. T. Züst, A. A. Agrawal, Trade-offs between plant growth and defense against insect herbivory: An emerging mechanistic synthesis. *Annu. Rev. Plant. Biol.* **68**, 513–534 (2017).

11. P. W. Sherman, The levels of analysis. *Anim. Behav.* **36**, 616–619 (1988).

12. T. A. Manolio et al., Finding the missing heritability of complex diseases. *Nature* **461**, 747–753 (2009).

13. I. Barrio-Hernandez et al., Network expansion of genetic associations defines a pleiotropy map of human cell biology. *Nat. Genet.* **55**, 1–10 (2023).

14. M. Fagny, F. Austerlitz, Polygenic adaptation: Integrating population genetics and gene regulatory networks. *Trends Genet.* **37**, 631–638 (2021).

15. E. A. Boyle, Y. I. Li, J. K. Pritchard, An expanded view of complex traits: From polygenic to omnigenic. *Cell* **169**, 1177–1186 (2017).

16. D. J. Watts, S. H. Strogatz, Collective dynamics of “small-world” networks. *Nature* **393**, 440–442 (1998).

17. S. H. Strogatz, Exploring complex networks. *Nature* **410**, 268–276 (2001).

18. X. Liu, Y. I. Li, J. K. Pritchard, Trans effects on gene expression can drive omnigenic inheritance. *Cell* **177**, 1022–1034.e6 (2019).

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 unelicitied 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*.

**Data, Materials, and Software Availability.** All raw sequence and optical man data were deposited in the NCBI SRA database. The Whole Genome project has been deposited at DDBJ/ENA/GenBank under the accession JAPQWU000000000 with BioProject id PRJNA900029 (102). The version described in this paper is version JAPQWU010000000. The MAGIC parents and RIL WGS are available under BioProject id PRJNA907539 (103). The two transcriptome datasets are available under BioProject PRJNA898122(104) and PRJNA911202 (105).

**ACKNOWLEDGMENTS.** We thank all the members of the Department of Molecular Ecology who contributed to the decade-long MAGIC project. Detailed acknowledgments with individual contributions are listed in *SI Appendix*. This work was supported by the Max Planck Society's core funding of I.T.B.'s department, Advanced grant no. 293926 of the European Research Council to I.T.B., and the Deutsche Forschungsgemeinschaft-SFB 1127/2 ChemBioSys-239748522 grant to I.T.B.

Author affiliations: *Department of Molecular Ecology, Max Planck Institute for Chemical Ecology, 07745 Jena, Germany; *Department of Natural Product Biosynthesis, Max Planck Institute for Chemical Ecology, 07745 Jena, Germany; *Department of Molecular Biology and Genetics, Cornell University, Ithaca, NY 14850; *Department of Geography, University of Zurich, 8006 Zurich, Switzerland; *Department of Chemistry, University of Zurich, 8006 Zurich, Switzerland; and *Institut national de recherche pour l'agriculture, l'alimentation et l'environnement, Centre National de Resources Génomiques Végétales, French Plant Genomic Resource Center, Castanet Tolasan F-31326, France

Author contributions: 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.

19. E. O. Paull et al., A modular master regulator landscape controls cancer transcriptional identity. *Cell* **184**, 334–351 (2021).

20. J. Peedicyayil, D. R. Grayson, An epigenetic basis for an omnigenic model of psychiatric disorders. *Theor. Biol.* **443**, 52 (2018).

21. P. Van Der Harst, N. Verweij, Identification of 64 novel genetic loci provides an expanded view on the genetic architecture of coronary artery disease. *Circ. Res.* **122**, 433–443 (2018).

22. A. Chateigner et al., Gene expression predictions and networks in natural populations supports the omnigenic theory. *BMC Genomics* **21**, 1–16 (2020).

23. P. Kemmeren et al., Large-scale genetic perturbations reveal regulatory networks and an abundance of gene-specific repressors. *Cell* **157**, 740–752 (2014).

24. D. N. Rate, J. V. Cuenca, G. R. Bowman, D. S. Guttman, J. T. Greenberg, The gain-of-function *Arabidopsis* *acd6* mutant reveals novel regulation and function of the salicylic acid signaling pathway in controlling cell death, defenses, and cell growth. *Plant Cell* **11**, 1695–1708 (1999).

25. M. Erb, D. J. Kliefenstein, Plant secondary metabolites as defenses, regulators, and primary metabolites: The blurred functional trichotomy. *Plant Physiol.* **184**, 39–52 (2020).

26. A. Kessler, I. T. Baldwin, Plant responses to insect herbivory: The emerging molecular analysis. *Annu. Rev. Plant. Biol.* **53**, 299–328 (2002).

27. D. Cipollini, D. Walters, C. Voelckel, "Costs of resistance in plants: from theory to evidence" in *Annual Plant Reviews, Insect-Plant Interactions*, C. Voelckel, G. Jander, Eds. (John Wiley & Sons, Ltd, 2014), pp. 263–307.

28. P. D. Coley, J. P. Bryant, F. S. Chapin III, Resource availability and plant antiherbivore defense. *Science* **230**, 895–899 (1985).

29. D. A. Herms, W. J. Mattson, The dilemma of plants: To grow or defend. *Q. Rev. Biol.* **67**, 283–335 (1992).

30. M. Dicke, I. T. Baldwin, The evolutionary context for herbivore-induced plant volatiles: Beyond the "cry for help". *Trends Plant Sci.* **15**, 167–175 (2010).

31. N. Yoshinaga et al., Fatty acid-amino acid conjugates diversification in lepidopteran caterpillars. *J. Chem. Ecol.* **36**, 319–325 (2010).

32. R. Halitschke, U. Schittko, G. Pohnert, W. Boland, I. T. Baldwin, Molecular interactions between the specialist herbivore *Manduca sexta* (Lepidoptera, Sphingidae) and its natural host *Nicotiana attenuata*. *III.* Fatty acid-amino acid conjugates in herbivore oral secretions are necessary and sufficient for herbivore-specific plant responses. *Plant Physiol.* **125**, 711–717 (2001).