1. Ecology
  2. Plant Biology
Download icon

Genetic variation, environment and demography intersect to shape Arabidopsis defense metabolite variation across Europe

  1. Ella Katz
  2. Jia-Jie Li
  3. Benjamin Jaegle
  4. Haim Ashkenazy
  5. Shawn R Abrahams
  6. Clement Bagaza
  7. Samuel Holden
  8. Chris J Pires
  9. Ruthie Angelovici
  10. Daniel J Kliebenstein  Is a corresponding author
  1. Department of Plant Sciences, University of California, Davis, United States
  2. Gregor Mendel Institute, Austrian Academy of Sciences, Vienna Biocenter (VBC), Austria
  3. Department of Molecular Biology, Max Planck Institute for Developmental Biology, Germany
  4. Division of Biological Sciences, Bond Life Sciences Center, University of Missouri, United States
  5. Division of Biological Sciences, Interdisciplinary Plant Group, Christopher S. Bond Life Sciences Center, University of Missouri, United States
  6. DynaMo Center of Excellence, University of Copenhagen, Denmark
Research Article
  • Cited 2
  • Views 1,673
  • Annotations
Cite this article as: eLife 2021;10:e67784 doi: 10.7554/eLife.67784

Abstract

Plants produce diverse metabolites to cope with the challenges presented by complex and ever-changing environments. These challenges drive the diversification of specialized metabolites within and between plant species. However, we are just beginning to understand how frequently new alleles arise controlling specialized metabolite diversity and how the geographic distribution of these alleles may be structured by ecological and demographic pressures. Here, we measure the variation in specialized metabolites across a population of 797 natural Arabidopsis thaliana accessions. We show that a combination of geography, environmental parameters, demography and different genetic processes all combine to influence the specific chemotypes and their distribution. This showed that causal loci in specialized metabolism contain frequent independently generated alleles with patterns suggesting potential within-species convergence. This provides a new perspective about the complexity of the selective forces and mechanisms that shape the generation and distribution of allelic variation that may influence local adaptation.

eLife digest

Since plants cannot move, they have evolved chemical defenses to help them respond to changes in their surroundings. For example, where animals run from predators, plants may produce toxins to put predators off. This approach is why plants are such a rich source of drugs, poisons, dyes and other useful substances. The chemicals plants produce are known as specialized metabolites, and they can change a lot between, and even within, plant species. The variety of specialized metabolites is a result of genetic changes and evolution over millions of years.

Evolution is a slow process, yet plants are able to rapidly develop new specialized metabolites to protect them from new threats. Even different populations of the same species produce many distinct metabolites that help them survive in their surroundings. However, the factors that lead plants to produce new metabolites are not well understood, and it is not known how this affects genetic variation.

To gain a better understanding of this process, Katz et al. studied 797 European variants of a common weed species called Arabidopsis thaliana, which is widely studied. The investigation found that many factors affect the range of specialized metabolites in each variant. These included local geography and environment, as well as genetics and population history (demography). Katz et al. revealed a pattern of relationships between the variants that could mirror their evolutionary history as the species spread and adapted to new locations.

These results highlight the complex network of factors that affect plant evolution. Rapid diversification is key to plant survival in new and changing environments and has resulted in a wide range of specialized metabolites. As such they are of interest both for studying plant evolution and for understanding their ecology. Expanding similar work to more populations and other species will broaden the scope of our ability to understand how plants adapt to their surroundings.

Introduction

Continuous and dynamic change in a plant’s habitat/environment creates a complex system to which a plant must adapt. Central to this adaptation are the production and accumulation of different metabolites ranging from signaling hormones, primary metabolites to a wide array of multi-functional specialized metabolites (Erb and Kliebenstein, 2020; Hanower and Brzozowska, 1975; Hayat et al., 2012; Kim et al., 2012; Kliebenstein, 2004; Malcolm, 1994; Thakur and Rai, 1982; Wolters and Jürgens, 2009; Yang et al., 2000). The complete blend, chemotype, of these metabolites helps to determine the plants’ survival and development, but the creation of any blend is complicated by the fact that individual specialized metabolites can have contrasting effects on the plant. For example, individual specialized metabolites can defend the plant against some stressors while simultaneously making the plant more sensitive to other biotic or abiotic stresses (Agrawal, 2000; Bialy et al., 1990; Erb and Kliebenstein, 2020; Futuyma and Agrawal, 2009; Hu et al., 2018; Lankau, 2007; Opitz and Müller, 2009; Uremıs et al., 2009; Züst and Agrawal, 2017). These opposing effects create offsetting ecological benefits and costs for individual metabolites. Integrating these offsetting effects across dynamic environments involves multiple selective pressures that might contribute to shaping the genetic and metabolic variation within a species (Fan et al., 2019; Kerwin et al., 2015; Malcolm, 1994; Sønderby et al., 2010; Szakiel et al., 2011; Wentzell and Kliebenstein, 2008; Züst et al., 2012).

Significant advances have been made in recent decades to identify genetic sources contributing to metabolic variation. A common finding of these studies is that the metabolic variation within and between species is the result of structural variation at the enzymes responsible for the chemical structures, or variation at the expression levels of these enzymes, which contributes to the quantitative variation in specialized metabolism (Chan et al., 2011; Chan et al., 2010; Fan et al., 2019; Kroymann et al., 2003; Moore et al., 2019; Schilmiller et al., 2012). These structural and regulatory variants and the resulting chemical variation strongly influence plant fitness in response to a broad range of biotic interactions including herbivores and other plant species (Bednarek and Osbourn, 2009; Brachi et al., 2015; Kerwin et al., 2017; Kerwin et al., 2015; Lankau and Kliebenstein, 2009; Lankau and Strauss, 2007; Lankau, 2007). The potential for these genetic variants influencing plant chemical variation is derived from the enhanced proportion of gene duplication in enzyme encoding genes for specialized metabolism, both at the local and whole genome level (Kliebenstein et al., 2001c; Moghe and Last, 2015). Many mechanistic studies of natural variation in specialized metabolism have focused on biallelic phenotypic variation linked to loss-of-function variants. However, it is not clear if biallelic phenotypic variation is created by biallelic genetic causation when investigating a large collection of individuals from wide-ranging populations within a species.

If selective pressures are sufficiently non-linear, it is possible to have repeated and independent generation of structural variants creating the same metabolic variation in processes that are akin to parallel and convergent evolution are used to describe interspecific variation. Specifically, parallel and convergent evolution describe independent evolution of the same trait that differs depending on the beginning state of the organisms. In parallel evolution, the lineages begin from the same state and in parallel evolve to the same new state, while in convergent evolution the lineages start at different states and independently converge on the same new state (Figure 1). In this context, we are focusing the analogy on the fact that the two processes differ where they begin, same or different state. This raises the possibility for chemical variation within a species to exhibit parallel evolution, wherein independent new haplotypes with identical metabolic consequences arise multiple times from single-core haplotype. Equally it may be possible to find within-species convergent evolution, where genotypes with the same metabolic profile actually contain completely different haplotypes that themselves arose from distinct haplotypic lineages. These genetic processes and interplay between genetics and selection overlap with neutral demographic processes like gene flow. Thus, it is necessary to understand how the intersection of environmental pressure, demography and genomic complexity gives rise to the pattern of metabolic variation across a plant species.

Parallel and convergent evolution.

The schema describes our use of parallel (A) and convergent (B) evolution for within-species chemotypic variation. The letters in the blue box represent the state of the source/ancestral haplotypes. The letters within the yellow box represent the newly derived haplotypes that arose by genetic mutation in the source haplotype. Finally, X, Y and Z show the chemotypes that arise from each haplotype. Blue and red arrows represent parallel or convergent genetic changes (respectively), while mustard arrows represent the enzymatic result.

To better understand how genomic variation, demography and environmental pressure shape the variation of specialized metabolism within a species, we used the Arabidopsis glucosinolate (GSL) pathway as a model. GSLs are a diverse class of specialized metabolites that display extensive variation across the order Brassicales, which includes the model plant Arabidopsis (Arabidopsis thaliana) (Bakker et al., 2008; Benderoth et al., 2006; Brachi et al., 2015; Chan et al., 2010; Daxenbichler et al., 1991; Halkier and Gershenzon, 2006; Kerwin et al., 2015; Kliebenstein et al., 2001a; Kliebenstein et al., 2001b; Kliebenstein et al., 2001c; Rodman et al., 1981; Rodman, 1980; Sønderby et al., 2010; Wright et al., 2002). GSLs consist of a common core structure with a diverse side chain that determines biological activity in defense, growth, development and abiotic stress resistance (Beekwilder et al., 2008; Hansen et al., 2008; Hasegawa et al., 2000; Katz et al., 2020; Katz et al., 2015; Malinovsky et al., 2017; Salehin et al., 2019; Yamada et al., 2003). The Arabidopsis-GSL system is an optimal model to study the species-wide processes driving specialized metabolite variation because the identity of the whole biosynthetic pathway is known, including the major causal loci for natural variation (Benderoth et al., 2006; Brachi et al., 2015; Chan et al., 2011; Chan et al., 2010; Hansen et al., 2007; Kliebenstein et al., 2001a; Kliebenstein et al., 2002b; Kliebenstein et al., 2002a; Kroymann and Mitchell-Olds, 2005; Pfalz et al., 2007; Sønderby et al., 2010; Wentzell et al., 2007). These major loci have been proven to influence Arabidopsis fitness and can be linked to herbivore pressure (Brachi et al., 2015; Hansen et al., 2008; Jander et al., 2001; Kerwin et al., 2017; Kerwin et al., 2015; Züst et al., 2012). Beyond the major causal loci, there is also evidence from genome-wide association (GWA) studies for highly polygenic variation in the genetic background that contributes to modulating GSL variation (Chan et al., 2011). The public availability of over 1000 widely distributed accessions with genomic sequences facilitates phenotyping GSL variation across a large spatial scale and analyses of causal haplotypes at the major GSL causal loci.

In Arabidopsis and other Brassicas, the main GSLs are methionine-derived, aliphatic, GSLs. Variation in the structure of aliphatic GSL is controlled by natural genetic variation at three loci: GS-Elong, GS-AOP and GS-OH. The specific alleles at these three loci combine to determine a predominant chemical structure and define chemically distinct aliphatic GSL chemotypes. In addition to these large-effect loci, there is a large suite of loci that can quantitatively alter the total accumulation and relative concentrations of GSLs within each chemotype (Brachi et al., 2015; Chan et al., 2011; Chan et al., 2010). GS-Elong differentially elongates the methionine side chain by the methylthioalkylmalate synthase enzymes (MAM). The elongation of the side chain by one methylene group is the result of one cycle that includes three steps: deamination of the methionine to create a ω-methylthio-2-oxoalkanoic-acid, condensation of the ω-methylthio-2-oxoalkanoic-acid with acetyl-CoA, and then isomerization and oxidative decarboxylation. The one carbon longer outcome can then undergo additional cycles of elongation (Benderoth et al., 2006; Graser et al., 2000; Kroymann et al., 2001; Textor et al., 2007). In Arabidopsis, MAM2 catalyzes the addition of one carbon to the side chain, creating GSLs with three carbon side chains. MAM1 catalyzes the addition of two carbons to make GSLs with four carbon side chains (Figure 2). MAM3 (also known as MAM-L) catalyzes the additions up to six carbons (Kliebenstein et al., 2001c; Kroymann et al., 2003; Mithen et al., 1995; Textor et al., 2007). The core pathway leads to the creation of methylthio GSL (MT). Then, the MT is converted to a methylsulfinyl (MSO) with a matching number of carbons (Giamoustaris and Mithen, 1996; Hansen et al., 2007). Structural variation at the GS-AOP locus leads to differential modification of the MSO by differential expression of a family of 2-oxoacid-dependent dioxygenases (2ODD). The AOP2 enzyme removes the MSO moiety leaving an alkenyl sidechain, while AOP3 leaves a hydroxyl moiety. Previous work has suggested three alleles of GS-AOP: the OHP allele that expresses only AOP3 and accumulates terminal OH containing GLS in the leaves and seeds; an alkenyl allele expressing AOP2 in the leaf and AOP2 and AOP3 in the seed leading to solely alkenyl GLS in the leaf and both alkenyl and OH aliphatic GLS in the seed; and a final allele containing a null mutation in the AOP2 gene that accumulates MSO GLS in the leaf and enhanced MSO and OH GLS in the seed. (Figure 2; Chan et al., 2010; Kliebenstein et al., 2001b; Kliebenstein et al., 2001c; Mithen et al., 1995). The C4 alkenyl side chain can be further modified by adding a hydroxyl group at the 2C via the GS-OH 2-ODD (Figure 2; Hansen et al., 2008). In spite of the evolutionary distance, independent variation at the same three loci influences the structural diversity in aliphatic-GSLs within Brassica, Streptanthus and Arabidopsis (Kliebenstein and Cacho, 2016; Lankau and Kliebenstein, 2009). For example, the MAMs responsible for C3 GSLs in Arabidopsis and Brassica represent two independent lineages, same as the MAMs responsible for C4 GSLs; in fact, the MAM locus contains at least three independent lineages that recreate the same length variation (Abrahams et al., 2020). This indicates repeated evolution across species, but it is not clear how frequently these loci are changing within a single species or how ecological or demographic processes may shape within-species variation at these loci.

Aliphatic glucosinolate (GSL) biosynthesis pathway.

Short names and structures of the GSLs are in black. Genes encoding the causal enzyme for each reaction (arrow) are in gray. GS-OX is a gene family of five or more genes. OH-But: 2-OH-3-Butenyl.

In this work, we described GSL variation in seeds of a collection of 797 A. thaliana natural accessions collected from different locations mainly in and around Europe. The amounts of GSLs can vary across different tissues and life stages, but there is a strong correlation in the type of aliphatic GSL produced across tissues (Brown et al., 2003; Kliebenstein et al., 2001a; Kliebenstein et al., 2001b; Petersen et al., 2002). Thus, in most cases the chemotype of the seeds is the same as the leaves. The seeds have the highest level of GSLs in Arabidopsis and are stable at room temperature until germination, which makes the seeds a perfect tissue to survey variation. Further, GSLs are known to be important for seed defenses against herbivores and pathogens (Raybould and Moyes, 2001). By combining GSL seed measurements with prior whole-genome sequencing in a European collection of accessions, we show that all three major causal loci controlling GSL metabolic diversity contain multiple independently derived alleles that recreate the same phenotypes using a combination of single nucleotide polimorphism (SNPs) and structural variation. Using these causal genotypes and chemotypes in combination with their geographic distribution provided evidence that the distribution of GSL metabolic diversity across Europe is influenced by a combination of demography and ecological factors. The ecological relationships to chemotype suggested a potential for variation in selective processes across the geographic regions studied. Future work will be needed to identify the specific biotic and/or abiotic factors shaping this distribution.

Results

GSL variation across Europe

To investigate the genetic, environmental and demographic parameters influencing the distribution of Arabidopsis GSL chemotypes, we measured GSLs from seeds of a collection of 797 A. thaliana natural accessions ( The 1001 Genomes Consortium, 2016). These Arabidopsis accessions were collected from different geographical locations, mainly in and around Europe (Figure 3A). 23 different GSLs were detected and quantified, identifying a wide diversity in composition and amount among the natural accessions with a median heritability of 83%, ranging from 34% to 93% (Supplementary file 1). To summarize the GSL variation among the accessions, we performed principal component analyses (PCA) on the accumulation of all the individual GSLs across the accessions as an unbiased first step. The first two PCs only captured 33% of the total variation with PC1 describing GSLs with four and seven carbons and PC2 mainly capturing GSLs with eight carbons in their side chain (Figure 3—figure supplement 1). Previous work using a collection of predominantly central European accessions had suggested a simple continental gradient chain-elongation variation from the south-west (SW) (that is enriched with alkenyl and hydroxyalkenyl GSLs) to the north-east (NE) (Brachi et al., 2015; Züst et al., 2012). To assess if this was still apparent in this larger collection, we plotted the accessions based on their geographical locations and colored them based on their PC1 and PC2 scores that are linked to chain elongation variation (Figure 3A and Figure 3—figure supplement 2A, respectively). This larger collection shows that there is not a single gradient shaping GSL diversity across Europe (Figure 3A). Instead, the extended sampling of accessions around the Mediterranean Basin in this collection shows that the SW to NE pattern reiterates within the Iberian Peninsula. In each of these areas (Iberian Peninsula and Central Europe), the SW is enriched with C4 GSLs, and the NE with C3 GSLs (Figure 3—figure supplement 1).

Figure 3 with 3 supplements see all
Glucosinolate variation across Europe is dominated by two loci.

(A) The accessions are plotted on the map based on their collection site and colored based on their principal component (PC)1 score. (B) Manhattan plot of genome-wideassociation analyses using PC1. Horizontal lines represent 5% significance thresholds using Bonferroni (red) and Benjamini–Hochberg (blue).

To test which of the major causal loci are detectable in this collection and to identify new genomic regions that are associated with the observed GSL variation, we performed GWA (with EMMAX algorithms) analyses using the PC1 and PC2 values. This collection of natural accessions presents a dense variant map that is 3× larger than previous GSL GWA mapping populations and includes 6,973,565 SNPs. In spite of the large population size, both PC1- and PC2-based analyses identified the same two major peaks covering two of the known causal gene clusters controlling GSL diversity (Figure 3B for PC1 GWA analyses, Figure 3—figure supplement 2B for PC2 GWA analyses) (Brachi et al., 2015; Chan et al., 2011; Chan et al., 2010). The largest peak in both cases is the GS-Elong locus on chromosome 5, containing the MAM1 (AT5G23010), MAM2 (that is not present in Col-0 plants) and MAM3 (AT5G23020) genes.

The peak on chromosome 4 is the GS-AOP locus containing the AOP2 and AOP3 genes (AT4G03060 and AT4G03050, respectively). Applying a more permissive cutoff did not result in the detection of any other related genes (Supplementary file 2). Previous QTL mapping and molecular experiments have shown that the genes within GS-AOP and GS-Elong loci are the causal genes for GSL variation within these regions (Benderoth et al., 2006; Brachi et al., 2015; Chan et al., 2011; Chan et al., 2010; Kliebenstein et al., 2001a; Kliebenstein et al., 2002a; Kliebenstein et al., 2002a; Kroymann and Mitchell-Olds, 2005; Pfalz et al., 2007; Wentzell et al., 2007). Surprisingly, none of the other known natural variants within the GSL biosynthetic pathway (listed in Supplementary file 2) were identified by GWA including three that were found with 96 accessions and three that were found with 595 accessions using PC1 and 2 (Brachi et al., 2015; Chan et al., 2011; Chan et al., 2010; Kliebenstein, 2009). Performing GWA studies using the accumulation of each of the 23 individual GSL detected in this collection resulted in an identical result, no additional known GSL-related genes were detected, while a few additional unknown genes were found (Figure 3—figure supplement 3 and Supplementary file 2). One explanation for that is that the dense sampling in this collection is available for mainly the Iberian Peninsula, the southern coast of Sweden and the south-western coast of Italy, and is still insufficient for Central Europe. Another possibility is that allelic heterogeneity for the other loci, and more complex patterns of interaction, may hamper their detection and influenced this high false-negative error rate where ~80% of prior validated natural variants found using multiple RIL populations were missed.

Complex GSL chemotypic variation

One potential complicating factor is that GSL chemotypic variation is best described as a discrete multimodal distribution involving the epistatic interaction of multiple genes which PCA’s linear decomposition cannot accurately capture (Figure 2). To test if PCA was inaccurately describing GSL chemotypic variation, we directly called the specific GSL chemotypes in each accession. Using Arabidopsis QTL mapping populations and GWA, we have shown that the GS-AOP, Elong and OH loci determine seven discrete chemotypes, 3MSO, 4MSO, 3OHP, 4OHB, Allyl, 3-Butenyl, 2-OH-3-Butenyl (Figure 2), that can be readily assigned from GSLs’ phenotypic data (Brachi et al., 2015; Chan et al., 2011; Chan et al., 2010; Kliebenstein et al., 2001a). The presence and amounts of these seven chemotypes provide a reliable indication about the existence and activity of each of the major GSL loci. Using accessions with previously known chemotypes and genotypes, we developed a phenotypic classification scheme to assign the chemotype for each accession (Figure 4; for details, see Methods and Figure 4—figure supplements 13; for structures, see Figure 2 and Supplementary file 1). Since the aliphatic GSLs’ composition in the seeds reliably indicates the GSL structural composition in the other plant’s life stages and tissues, assigning a chemotype for each accession based on the seeds’ composition is expected to be highly stable across tissues of the same accession (Brown et al., 2003; Chan et al., 2011; Chan et al., 2010; Kliebenstein et al., 2001a; Kliebenstein et al., 2001b). Most accessions were classified as 2-OH-3-Butenyl (27%) or Allyl (47%) with lower frequencies for the other chemotypes. Mapping the chemotypes onto Europe showed that the PCA decomposition was missing substantial information on GSL chemotype variation (Figure 4). Instead of a continuous distribution across Europe, the chemotype classifications revealed specific geographic patterns. Central and parts of Northern Europe (like north Germany and Poland) were characterized by a high variability involving the co-occurrence of individuals from all chemotypes. In contrast, southern Europe, which presents a dense sampling, including the Iberian Peninsula, Italy and the Balkan, has two predominant chemotypes, Allyl or 2-OH-3-Butenyl, that are separated from each other by a clear and sharp geographic partitioning (Figure 4 and Figure 4—figure supplement 4). Uniquely, Swedish accessions displayed a striking presence of almost solely Allyl chemotypes. Deeper sampling is required to test if this is or is not mirrored on the eastern side of the Baltic Sea as the few accessions from that region are almost solely 3OHP chemotypes (Finnish, Lithuanian, Latvian or Estonian accessions). Directly assigning GSL variation by discrete chemotypes provided a more detailed image not revealed by PCA decomposition. Further, the different chemotypic to geographic patterns suggest that there may be different pressures shaping GSL variation particularly when comparing Central and Southern Europe.

Figure 4 with 4 supplements see all
Phenotypic classification based on glucosinolate (GSL) content.

(A) Using the GSL accumulation, each accession was classified to one of seven aliphatic short-chained GSL chemotypes based on the enzyme functions as follows: MAM2, AOP null: classified as 3MSO dominant, colored in yellow. MAM1, AOP null: classified as 4MSO dominant, colored in pink. MAM2, AOP3: classified as 3OHP dominant, colored in green. MAM1, AOP3: classified as 4OHB dominant, colored in light blue. MAM2, AOP2: classified as Allyl dominant, colored in blue. MAM1, AOP2, GS-OH non-functional: classified as 3-Butenyl dominant, colored in black. MAM1, AOP2, GS-OH functional: classified as 2-OH-3-Butenyl dominant, colored in red. The accessions were plotted on a map based on their collection sites and colored based on their dominant chemotype. (B) The coloring scheme with functional GSL enzymes in the aliphatic GSL pathway is shown with the percentage of accessions in each chemotypes (out of the total 797 accessions) shown in each box.

Figure 4—source data 1

Environmental conditions differentially associate with MAM status in the north versus the south.

Linear model for MAM status (carbon side chain length) was conducted with the indicated environmental parameters, for the northern and southern collection, separately (for more details, see Methods). The tables show p values for each term from the linear model. For the interaction with geography, the linear model was run using the total dataset, and the geography parameter (north or south) was added to the model.

https://cdn.elifesciences.org/articles/67784/elife-67784-fig4-data1-v2.pptx

Geography and environmental parameters affect GSL variation

Because GSL chemotypes may be more reflective of local environment, we proceed to test if they are associated with weather parameters and landscape conditions. Further, given the difference in chemotype occurrence in Central and Southern Europe, we hypothesized that these environmental connections may change between Central and Southern Europe. For these tests, we chose environmental parameters that capture a majority of the environmental variance and by that may describe the type of ecosystem (Ferrero-Serrano and Assmann, 2019). We assigned each accession the environmental value based on its location. These environmental parameters include geographic proximity (distance to the coast), precipitation descriptors (precipitation of wettest and driest month) and temperature descriptors (maximal temperature of warmest month and minimal temperature of coldest month) and capture major abiotic pressures as well as provide information about the type of ecosystem in which each accession exists. Because demography and environment can be confounded, we included demography in our models using the previously assigned genomic groupings as components of the model (The 1001 Genomes Consortium, 2016). Further, we included specific geographic information by assigning the accessions to a northern or a southern collection, based on their location in relation to the following chain of mountains: the Pyrenees, the Alps and the Carpathians (Figure 4—figure supplement 4). We then ran a linear model for each geographic area separately (north and central vs. south) to check if the environmental parameters and the genomic population group associate with specific chemotypes. To directly test for an interaction of environment and geography, we ran the model with all accessions and incorporated the geography parameters and genomic population group. As the most frequent chemotypes in the collection are Allyl and 2-OH-3-Butenyl (Figure 4—figure supplement 4B), we focused the models on these chemotypes. The models showed that the environmental conditions have different relationships to the chemotypes that shift by geographic areas. Moreover, two of the parameters (min temp of coldest month and precipitation of wettest month) have a significant interaction with geography, suggesting that the relationship of these environmental parameters to specific GSL chemotypes is different between Northern and Southern Europe (Table 1; for details on the models, see Methods). This suggests that the relationship of GSL chemotype to environmental parameters varies across geographic regions of Europe rather than fitting a simple linear model.

Table 1
Environmental conditions differentially associate with major chemotypes across geographic location.

Linear model for the two major chemotypes, Allyl and 2-OH-3-Butenyl, was conducted with the indicated environmental parameters, for the northern and southern collection, separately (for more details, see Methods). The table shows p values for each term from the linear model. For the interaction with geography, the linear model was run using the total dataset, and the geography parameter (north or south) was added to the model.

Environmental parameterEffect on chemotype – northEffect on chemotype – southInteraction with geography
Genomic group<0.0001<0.0001<0.0001
Max temperature of warmest month0.0382<0.00010.3574
Min temperature of coldest month0.0007<0.00010.0049
Precipitation of wettest month0.16450.00030.0094
Precipitation of driest month0.06650.20260.47425
Distance to the coast0.27810.026800.1279

As the two main chemotypes in the collection differ by the length of the carbon chain (C3 for Allyl, C4 for 2-OH-3-Butenyl), we created a linear model to further check the interaction between each environmental condition to geography in respect to the carbon chain length. As was shown by the chemotypes models, most of the environmental parameters (min temp of coldest month, precipitation of wettest month and distance to the coast) significantly interacted with geography, showing again that the relationship of environment to GSL alleles changes across Europe (Figure 4—source data 1; for details on the models, see Methods). Conducting this analysis for each of the geographic areas separately highlighted this by showing that these parameters have different effects on the carbon chain length in each of the areas (Figure 4—source data 1).

The genetic architecture of GSL variation

The presence of different GSL chemotype to environmental relationships across Europe raises the question of how these chemotypes are generated. Are these chemotypes from locally derived alleles or obtained by the intermixing of widely distributed causal alleles? Further, if there are multiple alleles, do they display within-species convergent or parallel signatures? We focus on the GS-AOP, GS-Elong and GS-OH loci, the causal genes creating Arabidopsis GSL chemotypes, and use the available genomic sequences in all these accessions to investigate the allelic variation in these genes to map the allelic distribution and test the potential for convergent and/or parallel evolution within each locus.

GS-Elong: Because the variation in the GS-Elong locus is caused by complex structural variation in MAM1 and MAM2 that is not resolvable using the available data from short-read genomic sequence, we used the MAM3 sequence within this locus to ascertain the genomic relationship of accessions at the causal GS-Elong locus (Kroymann et al., 2003). We aligned the MAM3 sequence from each of the accessions, rooted the tree with the Arabidopsis lyrata orthologue (MAMb) and colored the tree tips based on the accessions-dominant chemotype.

The accessions were distributed across eight distinctive clades with each clade clustering accessions having either a C3 or C4 phenotype (Figure 5A and Figure 5—figure supplement 1 for bootstrap support). The clades C3/C4 status altered across the tree with three of the clades expressing C3 (MAM2) and five clades expressing the C4 (MAM1). The use of MAM3 clades as proxy for C3/C4 GLS chemotypes is supported by prior genomic sequencing of the GS-Elong region from 15 accessions (Figure 5B; Kroymann et al., 2003). To test for potential within locus recombination that may influence the overarching patterns, we compared the MAM3 tree to a tree obtained using MYB37, which is on the opposite end of the MAM locus from MAM3 (Figure 5—figure supplement 1F). We found that while the order of the clades in the MYB37 tree is different than their order in the MAM3 tree, the accessions’ classification to clades was similar among the two trees. This suggests that while there are potentially individual instances of within-locus recombination they are not influencing the overall genotype to chemotype linkage from MAM3, and MAM3 can be used as a reliable reflection of the structural variation in this locus.

Figure 5 with 5 supplements see all
MAM3 phylogeny.

(A) MAM3 phylogeny of Arabidopsis thaliana accessions, rooted by Arabidopsis lyrata MAMb, which is not shown because of distance. Tree tips are colored based on the accession chemotype. (B) The genomic structure of the GS-Elong regions in the previously sequenced accessions is shown based on Kroymann et al., 2003. The structures in the box are based on sequences obtained in this work. The numbers left to the structures indicate the number of sequenced accessions in this work (left) or by Kroymann et al., 2003 (right). The numbers are colored based on their clades. Bright gray arrows represent MAM1 sequences, and dashed arrows represent MAM2 sequences. Dark gray arrows represent MAM3 sequences. The number to the right of the genomic cartoon represents the number of carbons in the side chain. (C) Collection sites of the accessions colored by their clade classification (from section A) and shaped based on the side chain length of the aliphatic short-chained glucosinolates (circles for C3, triangles for C4).

Six of the clades in the MAM3 tree include accession/s with a previously sequenced MAM locus (Figure 5B; Kroymann et al., 2003), while two clades (clades 6 and 7) did not include any accession with a previously determined structure. We obtained long-read-based sequencing of 11 additional accessions from the 1001 Genome project for the MAM locus that included accessions in all clades including clades 6 and 7 (Figure 5—figure supplement 2—source data 1 for sequences). This showed that clade 6 are accessions that have a haplotype that contains a previously described chimeric MAM gene that combines the 5′ of MAM2 with the 3′ of MAM1 (Figure 5B and Figure 5—figure supplement 2; Benderoth et al., 2006; Kroymann et al., 2003). In these accessions, the chimeric gene leads to predominantly C3 GSLs. Clade 7 has a haplotype that is highly similar to clade 2 with a single copy of MAM1 leading to C4 GSLs. Comparing transposable elements in the two clades shows that they are different configurations.

The new sequenced accessions present in the clades with existing genomic haplotypes predominantly agreed with these previously published haplotypes. There were only three accessions with differences, two with a local duplication of a truncated MAM1 pseudogene in clades 1 and 2 (PHW-34 and TAL 07, respectively), and a second with a local duplication of a MAM1 pseudogene in clade 2 (Qar-8a, Figure 5—figure supplement 2).

The bootstrap support and smaller trees raised the possibility that clade 2 could be considered as two distinct clades (Figure 5—figure supplement 1, clades 2a and 2b). The chemotypes and haplotype in the accessions do not provide a clear mechanistic basis for separating this clade into two (Kroymann et al., 2003; Figure 5). Comparing the accessions across the main split in this clade suggested that one group of accessions (clade 2b) has lower total GSLs and a higher fraction of short-chain GSLs in comparison to the longer chain structures. Future work involving populations solely focused on this question would be needed to resolve the mechanistic basis of this difference and if this represents two distinct MAM loci.

One complication in interpreting the potential for parallel vs. convergent evolution in this locus is that the relationship between the major chemotype/haplotype groups is not resolvable with very low bootstraps (Figure 5—figure supplement 1). Functional parsimony would suggest that clade 4, by having both MAM1 and MAM2, may represent a single haplotype that can give rise to the other functional haplotypes via independent mutations akin to parallel evolution. Supporting this potential is the observation that MAM1 and MAM2 are likely derived via a tandem duplication with ensuing divergence since the separation from A. lyrata (Figure 5—figure supplement 3; Benderoth et al., 2009; Benderoth et al., 2006). Fully resolving this would require collecting more accessions to identify additional alleles that may contain the information necessary to better resolve the relationships amongst the haplotypes.

Using this phylogeny, we investigated the presence of the different GS-Elong haplotypes across Europe to ask if each region has a specific allele/clade or if the alleles are distributed across the continent. Specifically, we were interested if the strong C3/C4 partitioning in Southern Europe was driven by the creation of local alleles or if this partitioning might contain a wide range of alleles. If the latter is true, this can argue for a selective pressure shaping this C3/C4 divide. To understand the patterning of the C3/C4 haplotypes and chemotypes in Iberia, we plotted the accessions on the map and colored them based on their GS-Elong clade (Figure 5C and Figure 5—figure supplement 4). As expected given that genetic variation in Iberia results from a series of range expansions from Central Europe and Africa (Lee et al., 2017; Durvasula et al., 2017), there is extensive mixing of nearly all major European GS-Elong haplotypes in Iberia, except of clade 3 that is not present. In contrast, there is a sharp partitioning between the C3/C4 chemotypes created by these haplotypes. The strong geographic separation between the two chemotypes involving nearly all causal haplotypes (Figure 5—figure supplement 4) raises the possibility that the strong geographic partitioning of the C3/C4 chemotypes in Iberia may be driven by selective pressure enhancing the partitioning of the chemotypes, and not solely neutral demographic processes. The presence of a few accessions in Iberia that disagree with the sharp C3/C4 partition (Figure 5—figure supplement 4A) suggests that a new configuration of this loci arose in this area and is reflected in a few accessions. However, this requires further assessment.

Shifting focus to all of Europe showed that while most clades were widely distributed across Europe there were a couple of over-arching patterns (Figure 5C and Figure 5—figure supplement 5). GS-Elong clades 1 and 6 provide an example of potential gene flow between Iberia and Central Europe. In contrast, the absence of clade 3 in Iberia is more parsimonious with this haplotype having a glacial refugium in the Balkans followed by a northward flow wherein it mixed with the other clades. Other clades do not present evidence of a gene flow to the north as they are exclusive to the south as shown by clades 5 and 7. While these are both C4 clades, other C4 clades like clades 2 and 8 present a case of a gene flow to the north (Figure 5—figure supplement 5). This suggests that there are either differences in their GSL chemotype influencing their distribution or there are neighboring genes known to be under selection in Arabidopsis like FLC (AT5G10140) that may have influenced their distribution. In combination, this suggests that a complex demography is involved in shaping the chemotype’s identity with some regions, Iberia, showing evidence of local selection while other regions, Central Europe, possibly showing a blend requiring further work to delineate (Figure 5—figure supplement 5).

GS-AOP: Side chain modification of the core MSO GSL is determined by the GS-AOP locus. Most of the accessions contain a copy of AOP2 and a copy of AOP3, but only one of them will be functionally expressed (Chan et al., 2010), while in some cases both will be non-functional. To better understand the demography and evolution of the GS-AOP locus, we separately aligned the AOP2 and AOP3 sequences, rooted each tree with the A. lyrata orthologue and colored the trees tips based on the accessions-dominant chemotype (Figure 6—figure supplement 1).

The phylogenetic trees shared a very similar topology, yielding a clear separation between alkenyl (AOP2 expressed) and hydroxyalkyl (AOP3 expressed) accessions. Alkenyl expressing accessions like Cvi-0 with an expressed copy of the AOP2 enzyme formed a single continuous cluster (Figure 6A and Figure 6—figure supplement 1). In contrast, hydroxyalkyl (AOP3 expressed) accessions clustered into two separate groups with one group of 3OHP-dominant accessions partitioning from the rest of the accessions at the most basal split in the tree (Figure 6—figure supplement 1A, AOP2 tree). This haplotype is marked by having an inversion swapping the AOP2 and AOP3 promoters as shown in bacterial artificial chromosome sequencing of the Ler-0 accession (Figure 6D; Chan et al., 2010). The AOP3 tree also identified a second group of 3OHP-dominant accessions located among the alkenyl accessions. Analyzing the sequences of these accessions reveals that this small group of 3OHP accessions has a complete deletion of AOP2 and contains only AOP3 (Figure 6E). Thus, there are at least two independent transitions from alkenyl to hydroxyalkyl GSLs within Arabidopsis, neither of which are related to the alkenyl to hydroxyalkyl conversion within A. lyrata. This indicates that there are multiple alkenyl to hydroxyalkyl GSL conversions both within and between Arabidopsis species.

Figure 6 with 1 supplement see all
AOP genomic structure.

The genomic structure and causality of the major AOP2/AOP3 haplotypes are illustrated. Pink arrows show the AOP2 gene while yellow arrows represent AOP3. The black arrows represent the direction of transcription from the AOP2 promoter as defined in the Col-0 reference genome. Its position does not change in any of the regions. A-F represent the different structures. The black lines in C and F represent theoretical positions of independent variants creating premature stop codons. The GSL chemotype for each haplotype is listed to the right with the number of the accessions in brackets. The maps show the geographic distribution of the accessions from each structure.

The null accessions (MSO-dominant chemotypes) were identifiable in all the major clades on the tree (Figure 6—figure supplement 1, middle column of heatmap), suggesting that there are independent LOF mutations that abolish either AOP2 or AOP3. Deeper examination of the sequences of these accessions identified three convergent LOF alleles leading to the MSO chemotype. Most of the null accessions harbor a 5 bps deletion in their AOP2 sequence, which causes a frameshift mutation. This mutation arose within the alkenyl haplotype and was first reported in the Col-0 reference genome (Figure 6B; Kliebenstein et al., 2001c). In addition, there are additional independent LOF events arising in both the alkenyl haplotype (e.g., Sp-0, Figure 6C) and within the Ler-0 inversion haplotype (e.g., Fr-2, Figure 6F). Thus, GS-AOP has repeated LOF alleles arising within several of the major AOP haplotypes, suggesting convergent evolution of the MSO chemotype out of both the alkenyl and hydroxyalkyl chemotypes.

Using the combined chemotype/genotype assignments at GS-AOP, we investigated the distribution of the alleles across Europe. The alkenyl haplotype is spread across the entire continent. In contrast, the hydroxyalkyl haplotypes are geographically more restricted. The Ler-like 3OHP haplotype is present in only Central and North Europe (Figure 6D), while the other 3OHP haplotype, possessing only AOP3, is limited to Azerbaijan, along the Caspian Sea (Figure 6E). In contrast to the distinct hydroxyalkyl locations, the distribution of the independent LOF null haplotypes overlaps with all of them being located within Central and North Europe (Figure 6B, C and F). The fact that these independently derived LOF alleles are all contiguous suggests that they may be beneficial or neutral in Central Europe.

GS-OH: The final major determinant of natural variation in Arabidopsis GSL chemotype is the GS-OH enzyme that adds a hydroxyl group to the carbon 2 on 3-butentyl GSL to create 2-OH-3-Butenyl GSL. Previous work had suggested two GS-OH alleles measurable in the seed, a functional allele in almost all accessions and a non-functional allele caused by active site mutations represented by the Cvi-0 accession (Hansen et al., 2008). Because of functional epistasis, we can only obtain functional phenotypic information from accessions that accumulate the GS-OH substrate, 3-Butenyl GLS. This identified 11 accessions with a non-functional GS-OH. Surveying these 11 accessions in the polymorph database (The 1001 Genomes Consortium, 2016) identified multiple independent LOF events. One of these 11 accessions has the Cvi active site mutations, two accessions have a shared nonsense SNP that introduces premature stop codons and two accessions have a complete loss of this gene (Table 2). The other six accessions with a loss of enzyme activity had an unidentified lesion due to sequence quality for this locus.

Table 2
GSOH structure.

The structures of GS-OH in the 3-Butenyl accessions are illustrated. Gray boxes represent exons, and blue lines represent introns. Black line represents a mutation, and gray lines represent unknown lesions in hypothetical locations. The above mutations create non-functional GS-OH alleles. The fractions of the mutated accessions were calculated out of the total number of 3-Butenyl and OH-3-Butenyl accessions. The observed frequencies were calculated as the ratio between number of accessions with the specific mutation in non-C4 Alkenyl accessions and all the non-C4 Alkenyl accessions, as mentioned in the parentheses.

AccessionType of mutationAllele structureFraction (out of C4 Alkenyl accessions)Observed frequency (out of non-C4 Alkenyl accessions)
Sorbo, PienPolymorphism at SNP108313020.009 (2/226)0.067 (38/564)
Cvi-0Active site mutation0.004 (1/226)0.025 (14/564)
IP-Mot-0, IP-Tri-0Gene deletion0.009 (2/226)0.055 (31/564)
Multiple accessions (T670, FlyA-3, Ting-1, T880, T710, T850)Unidentified mutations0.026 (6/226)Unknown

All these independent GS-OH LOF alleles are found in accessions that do not accumulate 3-Butenyl GSL, for example, three carbon or non-alkenyl accessions, suggesting that functional epistasis may be influencing the maintenance of these alleles in nature. Thus, we searched for the accessions that do not accumulate 3-Butenyl GLS and carry GS-OH LOF events (Table 2). In all cases, the LOF allele is more frequent in the non four carbon-alkenyl accessions than expected by random chance. This suggests that there is selection against 3-Butenyl GSL synthesis since LOF alleles are more frequent when the GS-OH gene is cryptic by functional epistasis. This agrees with the fact that the 3-Butenyl chemotype is the most sensitive to generalist lepidopteran herbivory (Hansen et al., 2008). Thus, these mutations may represent ongoing pseudogenization of the GS-OH gene when it is functionally hidden by epistasis at the GS-AOP and GS-Elong loci. These LOF events would then only be displayed upon rare admixture with 2-OH-3-Butenyl accessions.

Discussion

Understanding the genetic, demographic and environmental factors that shape variation within a trait in a population is key to understanding trait evolution. In this work, we used a family of specialized metabolites, aliphatic GSLs, and measured their amounts in seeds of A. thaliana to query how genetics, geography, environment and demography intersect to shape chemotypic variation across Europe. We found that environmental conditions, together with geography, affect the presence and distribution of chemotypes within the accessions. This was demonstrated by specific traits that were associated with specific environmental conditions, and this association was shifted across the continent. Comparing the associations of traits to specific environmental conditions in Central Europe versus the south revealed different behaviors. This demonstrated that chemotypic variation across Europe is created by a blend of all these processes that differ at the individual loci. This implies that a simultaneous analysis of both genotype and phenotype is required to fully interpret these processes and relationships. The above analysis is extensively using abiotic factors because of their availability while aliphatic GSLs have mainly been linked to biotic interactions. GSLs have been mainly linked to influencing biotic interaction with herbivores considered the primary drivers shaping GSL genetic diversity. However, because climate drives the distribution of biotic factors like herbivores, it is likely that climatic factors might appear as indirectly associated with GSLs. Interestingly, an aliphatic GSL was recently mechanistically linked to drought resistance in Arabidopsis, suggesting a potential role for abiotic factors to influence GSL diversity (Salehin et al., 2019). More work is needed to dissect all the potential components of an environment that influence selection on GSL chemotypes.

All three major aliphatic GSL loci display extensive allelic heterogeneity that is shaped by a blend of evolutionary events at each locus reminiscent of either parallel or convergent evolution. For this analogy, we are defining parallel evolution to be when a new chemotype arises two or more independent times by independent mutations from shared ancestral haplotype. Conversely, we are considering convergent evolution of a chemotype to occur when independent mutations in independent ancestral haplotypes derive the same chemotype. GS-OH provided clear evidence of parallel evolution where a single functional haplotype gave rise to at least four independent LOF alleles all with similar phenotypic consequence. GS-AOP suggested the potential for convergent evolution-style events leading to MSO chemotype arising from LOF events at the GS-AOP locus. The first characterized LOF event was in the Col-0 accession that has a 5 bp frameshift indel in the AOP2 gene arising within the alkenyl AOP2-dominant GS-AOP haplotype. In this study, we identified additional parallel AOP2 LOF events in this haplotype. More critical, we could identify multiple independent LOF events arising in the AOP3 gene within the AOP3-dominant inversion GS-AOP haplotype. Thus, the same GSL chemotype, MSO, arises from independent LOF alleles in two different genes, AOP2 and AOP3, that represent two different ancestral haplotypes (Figure 6). Thus, it appears that parallel-style events occur at all the loci and at least at the GS-AOP locus there is potential for a convergent-style evolutionary event leading to a single chemotype.

In addition to independent evolutionary events, three-way epistasis is shaping the allelic heterogeneity at these loci and their evolutionary potential. For example, the multiple independent GS-OH LOF variants all appear to have arisen in lineages where the GS-AOP and GS-Elong loci had haplotypes that epistatically combined to block the formation of the but-3-enyl GSL precursor for GS-OH (Figure 2). Thus, the parallel evolution of GS-OH LOF alleles is epistatically conditioned by GS-AOP and GS-Elong. A similar epistatic contingency also exists between the two independent AOP3 alleles of GS-AOP (Figure 6D and E) and GS-Elong. Both of AOP3 alleles of GS-AOP are coordinated with C3 haplotypes in GS-Elong. The AZE allele of GS-AOP is associated with a novel geographically limited GS-Elong allele in clade 8 (Figure 5—figure supplement 5) while the Central European AOP3 allele is limited to accessions containing the clade 1, 3 or 6 C3 haplotypes of GS-Elong. There is also within-locus epistasis in the GS-Elong locus wherein a functional MAM1 leads to the creation of C4 GSLs regardless of the functional state of the MAM2 gene and C3 haplotypes are marked by the loss of MAM1. All of these within- and between-loci interactions create a directional arrow for most loci where the haplotypes do not have equal evolutionary potential. For example, the clade 4 haplotype of GS-Elong can equally mutate to create a C3 or C4 chemotype because it has both MAM1 and MAM2. However, the remaining clades like clade 3 have lost one or the other gene limiting their ability to create alternative chemotypes. In this case, the loss of MAM1 likely prevents the ability for this C3 lineage to recreate the C4 chemotype. Thus, the potential evolutionary trajectory of a haplotype/allele at one or even within one of these GSL loci may be epistatically conditioned by the allelic state all the loci within a specific lineage.

This level of allelic diversity at these loci raises a question of how do pathways with this level of diversity and structural variation pass through speciation boundaries (Figure 5—figure supplement 4). Confounding this further is the observation that Brassica ssp. have genetic variation at GS-Elong and GS-AOP creating the exact same chemotypic variation found in Arabidopsis (Heidel et al., 2006; Ramos-Onsins et al., 2004; Windsor et al., 2005). There is similar within-species (GS-AOP) and between-species (GS-AOP and GS-Elong) variation between the closely related Arabidopsis sister species, A. lyrata, A. petraea and A. halleri. However, the underlying genetic basis is independent events at the AOP and MAM loci showing that the variation did not go through the speciation boundary. Instead, this suggests that this variation has been recreated repeatedly in these species This raises the possibility that there may be a class of loci that are being repeatedly sampled by pangenomic variation across species within a family. To test this possibility would need a deeper phylogenetic sampling within and between species, particularly for understanding the intersection of ecology and evolution (Göktay et al., 2021; Durvasula et al., 2017).

Previous work on other biotic interactions genes like pathogen resistance gene-for-gene loci had indicated a predominant model of having two moderate-frequency ancient alleles creating the phenotypic variation within the species (Atwell et al., 2010; Corrion and Day, 2001; MacQueen et al., 2016). In contrast to R-gene loci characterized by old/stable biallelic variation, the GSL loci are characterized by a blend of structural and SNP-based variation with numerous alleles that appear young. In other cases, alleles of genes involved in biotic defense can present more complex patterns, for example, natural variation in the immune gene ACCELERATED CELL DEATH 6 (ACD6) is caused by a rare allele causing an extreme lesion phenotype. It is not yet clear what selective pressures influence ACD6 genetic variation (Todesco et al., 2010; Zhu et al., 2018). Thus, loci controlling resistance to diverse biotic traits under natural conditions have diverse genetic architectures and further work is needed to assess the range of allelic heterogeneity in these adaptive loci.

The allelic diversity at the GSL loci illustrates the benefit of simultaneously tracking the phenotype and genotype when working to understand the distribution of trait variation. For example, the Iberian Peninsula and the Mediterranean Basin had low variability in aliphatic GSL chemotypes, which show strong geographic structure. By contrast, Central/North Europe had high aliphatic GSL diversity with chemotypes showing overlapping geographic distributions. At first glance, this contrasts with previous work showing that the Iberian Peninsula and the Mediterranean Basin are genetically diverse. However, this discrepancy was caused by one of the causal loci. Specifically, the GS-AOP locus was largely fixed as the Alkenyl allele in Iberia/Mediterranean Basin with the alternative GS-AOP alleles enriched in Central Europe. In contrast to GS-AOP, Iberia and the Mediterranean Basin were highly genetically diverse for the GS-Elong locus and appear to contain almost all the variation in GS-Elong found throughout Europe (The 1001 Genomes Consortium, 2016). Thus, the chemotypic divergence from genomic variation expectations was driven by just the GS-AOP locus. This indicates that the high level of chemotypic variation in Central Europe is a blend of alleles that emerged in the south (GS-Elong) and alleles that possibly arose locally (GS-AOP, both nulls and AOP3). Further, the chemotypes found in any one region appear to be created by a combination of alleles as a result of a gene flow across the continent, local generation of new polymorphisms and local selective pressures.

Another challenge potentially caused by allelic heterogeneity and differential selective pressures, as displayed within this system, is detecting the known and validated causal natural variants within a population. Specifically, the GWA with this collection of 797 accessions was unable to find 80% of the known causal loci including one of the three major effect loci, GS-OH. Maximizing the number of genotypes and the SNP marker density was unable to overcome the complications imposed by the complex pressures shaping the distribution of these traits, potentially due to unequal dense sampling from the different areas. In this system, the optimal path to identifying the causal polymorphisms has instead been a small number of Recombinant Inbred Line populations derived from randomly chosen parents. In complex adaptive systems, the optimal solution to identifying causal variants is likely a blend of structured mapping populations and then translating the causal genes from this system to the GWA results and tracking the causal loci directly.

In this work, we combined different approaches to uncover some of the parameters shaping the aliphatic GSL content across Europe. Widening the size of the population will enable us to deepen our understanding on the evolutionary mechanisms shaping a phenotype in a population.

Materials and methods

Plant material

Request a detailed protocol

Seeds for 1135 Arabidopsis (A. thaliana) genotypes were obtained from the 1001 genomes catalog of A. thaliana genetic variation (https://1001genomes.org/). All Arabidopsis genotypes were grown at 22°C/24°C (day/night) under long-day conditions (16 hr of light/8 hr of dark). Two independent replicates were performed, each of them included the full set of genotypes. The replicates obtained from independent maternal plants were grown in randomized fashion. In the analyses, only accessions from Europe and around Europe were included (Figure 3A), resulting in an analysis of 797 accessions. A list of the accessions can be found in Supplementary file 1.

GSL extractions and analyses

Request a detailed protocol

GSLs were measured as previously described (Kliebenstein et al., 2001a; Kliebenstein et al., 2001b; Kliebenstein et al., 2001c). Briefly, ~3 mg of seeds were harvested in 200 μL of 90% methanol. Samples were homogenized for 3 min in a paint shaker, centrifuged, and the supernatants were transferred to a 96-well filter plate with DEAE sephadex. The filter plate with DEAE sephadex was washed with water, 90% methanol and water again. The sephadex-bound GSLs were eluted after an overnight incubation with 110 μL of sulfatase. Individual desulfo-GSLs within each sample were separated and detected by HPLC-DAD, identified, quantified by comparison to standard curves from purified compounds and further normalized to the weight. A list of GSLs and their structure is given in Supplementary file 1A. Raw GSLs data are given in Supplementary file 1B.

Statistics, heritability and data visualization

Request a detailed protocol

Statistical analyses were conducted using R software (https://www.R-project.org/) with the RStudio interface (http://www.rstudio.com/). For each independent GLS, a linear model followed by ANOVA was utilized to analyze the effect of accession, replicate and location in the experiment plate upon the measured GLS amount. Broad-sense heritability (Supplementary file 1C) for the different metabolites was estimated from this model by taking the variance due to accession and dividing it by the total variance. Estimated marginal means (emmeans) for each accession were calculated for each metabolite from the same model using the package emmeans (CRAN, 2021a; Supplementary file 1D). PCAs were done with FactoMineR and factoextra packages (Abdi and Williams, 2010). Data analyses and visualization were done using R software with tidyverse (Wickham et al., 2019) and ggplot2 (Kahle and Wickham, 2013) packages.

Maps were generated using ggmap package (Kahle and Wickham, 2013).

Phenotypic classification based on GSL content

Request a detailed protocol

For each accession, the expressed enzyme in each of the following families was determined based on the content (presence and amounts) of short-chained aliphatic GSLs.

MAM enzymes: The total amount of three carbon GSLs and four carbon GSLs was calculated for each accession. Three carbon GSLs include 3MT, 3MSO, 3OHP and Allyl GSL. Four carbon GSLs include 4MT, 4MSO, 4OHB, 3-Butenyl and 2-OH-3-Butenyl GSL (for structures and details, see Supplementary file 1). Accessions that the majority of aliphatic short-chained GSL contained three carbons in their side chains were classified as MAM2 expressed (Figure 4—figure supplement 1). Accessions that the majority of aliphatic short-chained GSL contained four carbons in their side chains were classified as MAM1 expressed (Figure 4—figure supplement 1). The accessions were plotted on a map based on their original collection sites (Figure 4—figure supplement 1).

AOP enzymes: The relative amount of alkenyl GSL, alkyl GSL and MSO GSL was calculated in respect to the total short-chained aliphatic GSL as follows:

AlkenylGSL(AOP2expressed)=Allyl +2OH3butenyl + 3butenylTotal short chained GSL
AlkylGSL(AOP3expressed)=3OHP + 4OHBTotal short chained GSL
MSOGSL(AOPnull)=3MSO + 4MSOTotal short chained GSL

The expressed AOP enzyme was determined based on those ratios: Accessions with majority alkenyl GSL were classified as AOP2 expressed. Accessions with majority of alkyl GSL were classified as AOP3 expressed. Accessions with majority of MSO GSL were classified as AOP null. The accessions were plotted on a map based on their original collection sites (Figure 4—figure supplement 2).

GS-OH enzyme: The ratio between 2-OH-3-Butenyl GSL to 3-Butenyl GSL was calculated only for MAM1-expressed accessions (accessions that the majority of GSLs contain four carbons in their side chain). Accessions with high amounts of 2-OH-3-Butenyl GSL were classified as GS-OH functional. Accessions with high amounts of 3-Butenyl GSL were classified as GS-OH non-functional. The accessions were plotted on a map based on their original collection sites (Figure 4—figure supplement 3).

Each accession was classified to one of seven aliphatic short-chained GSLs based on the combination of the dominancy of the enzymes as follows: MAM2, AOP null: classified as 3MSO dominant. MAM1, AOP null: classified as 4MSO dominant. MAM2, AOP3: classified as 3OHP dominant. MAM1, AOP3: classified as 4OHB dominant. MAM2, AOP2: classified as Allyl dominant. MAM1, AOP2, GS-OH non-functional: classified as 3-Butenyl dominant. MAM1, AOP2, GS-OH functional: classified as 2-OH-3-Butenyl dominant. The accessions were plotted on a map based on their original collection sites and colored based on their dominant chemotype (Figure 4).

Environmental and demographic data

Request a detailed protocol

Environmental and demographic data (referred to as ‘genomic group’) were obtained from the 1001 genomes website (https://1001genomes.org/, for geographical and demographic data) and from the Arabidopsis CLIMtools (http://www.personal.psu.edu/sma3/CLIMtools.html, Ferrero-Serrano and Assmann, 2019) for environmental data. We chose the five variables that captured a majority of the variance in this dataset based on PCA using different combinations of variables. The chosen variables are maximal temperature of warmest month (WC2_BIO5), minimal temperature of coldest month (WC2_BIO6), precipitation of wettest month (WC2_BIO13), precipitation of driest month (WC2_BIO14) and distance to the coast (in km). Each one of the above variables (including genomic group) was assigned to each one of the accessions.

Environmental models

Request a detailed protocol

Linear models to test the effect of geographical and environmental parameters (Figure 3—figure supplement 1 and Figure 4—source data 1) were conducted using dplyr package (CRAN, 2021b) and included the following parameters:

Figure 3—figure supplement 1 linear models for collection sites: PC score ~ Latitude + Longitude + Latitude * Longitude.

Table 1 and Figure 4—source data 1 for all the data: C length (C3 and C4) or the chemotypes (Allyl and 2-OH-3Butenyl) ~ Genomic group + Geography (north versus south) + Max temperature of warmest month + Min temperature of coldest month + Precipitation of wettest month + Precipitation of driest month + Distance to the coast + Geography * Genomic group + Geography * Max temperature of warmest month + Geography * Min temperature of coldest month + Geography * Precipitation of driest month + Geography * Precipitation of wettest month + Geography * Distance to the coast.

For the north and the south: C length (C3 and C4) or the chemotypes (Allyl and 2-OH-3Butenyl) ~ Genomic group + Geography (north versus south)+ Max temperature of warmest month + Min temperature of coldest month + Precipitation of wettest month + Precipitation of driest month + Distance to the coast.

Genome-wide association studies

Request a detailed protocol

The phenotypes for GWA studies were each accession value for PC1 and 2. GWA was implemented with the easyGWAS tool (Grimm et al., 2017) using the EMMAX algorithms (Kang et al., 2010) and a minor allele frequency (MAF) cutoff of 5%. The results were visualized as Manhattan plots using the qqman package in R (Turner, 2014).

Phylogeny

Request a detailed protocol

Genomic sequences from the accessions for MAM3 – AT5G23020, AOP2 – Chr4, 1351568 until 1354216, AOP3 – AT4G03050.2, GS-OH – AT2G25450 and MYB37 – AT5G23000 were obtained using the Pseudogenomes tool (https://tools.1001genomes.org/pseudogenomes/#select_strains).

Multiple sequence alignment was done with the msa package (default settings) in R using the ClustalW, ClustalOmega and Muscle algorithms (Bodenhofer et al., 2015). Phylogenetic trees were generated with the ‘ape’ package (neighbor-joining tree) (Paradis and Schliep, 2019) and were visualized with ggtree package in R (Yu, 2020). Each tree was rooted by the genes matching A. lyrata’s functional orthologue or closest homologue.

Bootstrap analyses (Bootstrap = 100) was done with ‘ape’ package in R (Paradis and Schliep, 2019), with the same tree inference method as described before. For MAM3 bootstrap analysis, the accessions with low-quality sequencing were excluded.

Amino acid phylogenies: Sequences were taken from Abrahams et al., 2020, which uses A. thaliana Col-0 genome and the MAM2 amino acid sequence 1006452109 from the Arabidopsis Information Resource (TAIR) database. Alignments were run using MAFFT (Katoh et al., 2017; Kuraku et al., 2013) and cleaned using Phyutility at a 50% occupancy threshold (Smith and Dunn, 2008). RAxML was used for phylogenetic inference (Stamatakis, 2014) with the PROTCATWAG model (Bootstrap = 1000).

Sequencing

Request a detailed protocol

PacBio long read-based de novo genome assemblies of the relevant accession were generated as part of the 1001 Genomes Plus project. The genomes were assembled with Canu (v1.71) (Koren et al., 2017) and polished using the long reads followed by a second polishing step with PCR-free short reads.

Data availability

All data generated or analyzed during this study are included in the manuscript and supporting files.

References

    1. Abdi H
    2. Williams LJ
    (2010) Principal component analysis
    Wiley Interdisciplinary Reviews: Computational Statistics 2:433–459.
    https://doi.org/10.1002/wics.101
    1. Benderoth M
    2. Pfalz M
    3. Kroymann J
    (2009) Methylthioalkylmalate synthases: Genetics, ecology and evolution
    Phytochemistry Reviews: Proceedings of the Phytochemical Society of Europe 8:255–268.
    https://doi.org/10.1007/s11101-008-9097-1
  1. Book
    1. Corrion A
    2. Day B
    (2001) Pathogen resistance signalling in plants
    In: Hetherington AM, editors. ELS. Chichester, UK: John Wiley & Sons, Ltd. pp. 1–14.
    https://doi.org/10.1002/9780470015902.a0020119.pub2
    1. Giamoustaris A
    2. Mithen R
    (1996)
    Genetics of aliphatic glucosinolates. IV. Side-chain modification in Brassica oleracea
    TAG. Theoretical and Applied Genetics. Theoretische Und Angewandte Genetik 93:1006–1010.
    1. Kliebenstein DJ
    2. Figuth A
    3. Mitchell-Olds T
    (2002a)
    Genetic architecture of plastic methyl jasmonate responses in Arabidopsis thaliana
    Genetics 161:1685–1696.
    1. Kliebenstein D
    2. Pedersen D
    3. Barker B
    4. Mitchell-Olds T
    (2002b)
    Comparative analysis of quantitative trait loci controlling glucosinolates, myrosinase and insect resistance in Arabidopsis thaliana
    Genetics 161:325–332.
    1. Uremıs I
    2. Arslan M
    3. Sangun MK
    4. Uygur V
    5. Isler N
    (2009)
    Allelopathic potential of rapeseed cultivars on germination and seedling growth of weeds
    Asian Journal of Chemistry 21:2170–2184.

Decision letter

  1. Meredith C Schuman
    Senior and Reviewing Editor; University of Zurich, Switzerland
  2. Arthur Korte
    Reviewer; Vienna Biocenter, Austria

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

The manuscript by Katz and colleagues makes key advances in understanding how glucosinolate chemotype variation across Arabidopsis ecotypes is linked to the underlying genetic variation. It demonstrates the complexity of the interplay between genetic and environmental factors and gives insights into how complex traits could evolve in natural populations. This article reveals the complexities involved in conducting genome-wide association studies on geographically structured populations which cannot be eliminated even in a very well studied system. It thus bears important general lessons for researchers in the genomic era.

Decision letter after peer review:

[Editors’ note: the authors submitted for reconsideration following the decision after peer review. What follows is the decision letter after the first round of review.]

Thank you for submitting your work entitled "Genetic variation, environment and demography intersect to shape Arabidopsis defense metabolite variation across Europe" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and a Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Arthur Korte (Reviewer #1); Juergen Kroymann (Reviewer #3).

Our decision has been reached after consultation between the reviewers. Based on these discussions and the individual reviews below, we regret to inform you that your work will not be considered further for publication in eLife. However, the reviewers and editors agree that both the topic and the general approach are of interest. eLife would consider a substantially revised version of the manuscript which addresses the points raised by the individual reviewers and summarized below. Because we expect that these revisions represent substantial work, we would consider a revised manuscript as a new submission.

Specifically, the reviewers and editors agree that the study could provide an important new insight into patterns and processes of evolution. In the words of Reviewer 1, it is valuable "to present the complexity of the interplay between genetic and environmental factors and its contribution to trait evolution in natural populations." It is perhaps especially valuable to do so in such a well-studied system as Arabidopsis thaliana, as this indicates fundamental complexities which cannot be remedied by large volumes of data alone, but require well-informed approaches to analyze the data and in some cases, perhaps, cannot be resolved. Thus although we are not convinced that any one aspect of this study is very novel, we do think that the study as a whole presents an unusually holistic perspective and thus could be a very valuable contribution. Nevertheless, we also agree that there are substantial gaps and pitfalls in the analyses as they are currently presented, and these must be remedied in order for the conclusions of the study to be well supported:

1. We are concerned that MAM3 is not a good proxy for the GS-Elong locus, and there is evidence of recombination within MAM3; and that resequencing of this locus in the accessions studied is required to support the conclusions that the authors wish to draw in this study. Please see detailed comments from Reviewer 3.

2. Several sections of the study require deep revision in terms of writing and analysis. See major comments from Reviewer 2, as well as the other reviewers.

Reviewer #1:

The authors present a detailed and well-written analyses of the variation in the glucosinolate chemotype of different Arabidopsis ecotypes. These differences were linked to the underlying genetic variation. They could show many complex patterns that are related to the distinct chemotypes. The analyses highlights the complexity of the interplay between genetic and environmental factors and gives insights into how traits evolve in natural populations.

I only have one point that I was wondering about.

23 different GSL have been detected and quantified and GWAS has been performed on PCs that summarize the GSL variation. The GWAs on PC1 and PC2 do only identify the two known loci AOP and NAM and no additional loci are significant. As this first two PCs explain only 33 % of the total phenotypic variation, I wonder if GWAS on the individual GSL concentration, especially as all of these show a high heritability, could identify additional associations? It would be nice to have an overview of these potential associations.

Reviewer #2:

This manuscript describes the contribution of environment and demography to shape the genetic variation underlying the glucosinolate (GLS) composition of Arabidopsis thaliana across Eurasia. To this end, the authors carried out a very extensive characterization of different GLSs in a large collection of natural Arabidopsis accessions. Genetic analyses of these data show that the epistasis underlying causal loci hampers the identification of new loci. More importantly, authors compare these data with environmental factors to show significant environmental associations, which support a differential adaptive role of this variation across Europe. In addition, thorough genetic analyses of the major underlying genes identified multiple independent mutations in several duplicated genes, which demonstrate a repeated evolution pattern with differential geographic histories. Overall, this study shows a complex adaptive history of Arabidopsis plant defenses mediated by GLSs at a continental scale and brings forward our current understanding of adaptation of plant defense mechanisms.

Although most sections of the manuscript are written in a comprehensive manner, several sections of results and discussion are rather confusing due to lack of information on particular analyses. The following sections should be revised:

1. The section entitled "Geography and environmental parameters affect GLS variation" (pages 13-14) is very confusing due to an unclear distinction between demography and environment. Given the strong geographic structure of Arabidopsis genetic diversity across Europe (The 1001 Genomes Consortium, 2016), most genetic diversity can be expected to correlate with environment due to Arabidopsis demographic history but not reflecting true adaptive association. For this reason, to distinguish true environmental associations from correlations related with demographic history, analyses must take into account the genetic structure (genetic or genomic relationships among accessions). Usually, the effect of demographic factors on the distribution of genetic diversity is inferred by the genetic structure. According to Supp. Figure 7 and Materials and Methods, authors seem to have included such genetic structure in their MANOVA analyses, which is referred to as genomic group. However, authors should describe what are those genomic groups, and if they were included also in Figure 4B. As shown by Supp. Figures 7A and 8, the genetic structure explains most of the diversity in MANOVA and Random forest models, indicating that most of the geographic patterns can be explained by the demographic history. However, both Figures also show that there is still a significant association with environmental factors, which supports that the environment contributes to maintain the genetic diversity of chemotypes. The text describing these results should be carefully rewritten to explain how the genetic structure (demography) has been taken into account and separating clearly demographic and environmental components. The genetic structure is only mentioned in line 259 where it is confusingly named as "ancestral state".

2. The section describing the genetic diversity of MAM3 across Europe (lines 302-326) shows certain overinterpretation of the action of natural selection, which should be revised and rewritten in a softer tone. For instance, the sentence in line 306 stating "If the latter is true, this would argue for a selective pressure shaping this C3/C4 divide". The two scenarios described in paragraph of lines 302-306 could be explained by natural selection or demography, and one should be cautious interpreting geographic patterns. In particular, the geographic pattern observed in Iberia is just showing that this region contains a large diversity for C3/C4 chemotypes and MAM3, as shown previously for most genes in this region (The 1001 Genomes Consortium). Therefore, the sentence in lines 310-313 ("This suggests that the strong geographic partitioning of the C3/C4 chemotypes in Iberia may be driven by selective pressure causing the partitioning of the chemotypes rather than neutral demographic processes") or lines 324-326 ("…,Iberia, showing evidence of local selection while other regions, central Europe, possibly showing a blend requiring further work to delineate") are not fully supported by data. In addition, this section on MAM3 could include the recent study on the genetic diversity of Europe, which suggests that Iberia could have also been colonized from central Europe after the last glaciation (Lee et al., 2017; Nature Communications 8:14458).

3. The section describing GS-OH gene and Table 1 (lines 377-390) is rather confuse. In line 377 is stated: "We could not identify the causal LOF in the other six accessions due to sequence quality". However, Table 1 shows the structure of GS-OH gene with apparent location of six mutations. It is unclear which are those six accessions of Table1, what is represented in the structure of the genes, and the meaning of symbols on those genes (A and the vertical bars). Are those accessions of Table 1 the same described in the text? It should also made clear the type of mutations of those six accessions instead of the confusing term "Independent mutation". Furthermore, this section refers to Supplementary Table 2 (line 383), but I could not find such table.

Reviewer #3:

The manuscript is based on an impressive data set of (seed) glucosinolate profiles across Europe and beyond, with a focus on methionine-derived glucosinolates. However, the manuscript reads in many aspects like a 'best of' of previous research, with limited net gain in insight. In addition, there are substantial weaknesses (in particular relating to GS-Elong) that prevent me from recommending this manuscript for publication in eLife its current state.

1. The representation of the MAM reaction is misleading and incorrect (lines 101ff). MAMs catalyze condensation of 2-oxoacids with acetylCoA. Two further reaction steps, catalyzed by an isomerase and a dehydrogenase, are necessary to complete chain extension. This extension consists formally in the addition of a methylene group per reaction cycle. With MAM2, the net result is an extension by one (not two!) methylene group, with MAM1 it is mostly two; note that methionine has already two methylene groups in the side chain! I suggest replacing the reference Abrahams et al., 2020, with Graser et al., 2000 (Archives of Biochemistry and Biophysics, 378, 411-419). Furthermore, Kroymann et al., 2001 (Plant Phys. 127, 1077-1088), Benderoth et al., 2006 and Textor et al., 2007 (Plant Phys. 144, 60-71) should be added as original references for the MAM1, 2 and 3 reactions.

2. The use of MAM3 as a proxy for inferring the configuration of the GS-Elong locus is highly debatable. Even though recombination appears to be suppressed at and around GS-Elong, an additional gene (or additional genes) upstream of GS-Elong should be used to construct (a) gene tree(s) and compare with the MAM3 gene tree for congruency. In any case, without sequencing of the MAM genes at the GS-Elong locus, it is not possible to conclude that 'local' configurations of MAM1/MAM2 have not arisen in Iberia (or elsewhere).

3. The A. lyrata ortholog of Arabidopsis MAM3 is MAMb (not MAMc). The MAM3 gene tree does not provide any statistical support for the branching order. Therefore, it remains unclear whether the presented tree reflects the true branching order and whether the 'most basal' clade in the tree is indeed the most basal. Note also that the Cal-0 gene at the MAM2 position is actually mostly MAM2, except for a converted stretch towards the 3' end of the gene; see: Figure 2 in Kroymann et al., 2003. Hence, the most parsimonious interpretation of available data is that the ancestral state of GS-Elong in Arabidopsis is indeed: MAM2, MAM1 and MAM3. Note also that the Arabidopsis genes in Abrahams et al., 2020, are from TAIR 10, i.e. from Col-0 (see Table 1 in Abrahams et al., 2020). Because Col-0 does not have MAM2, this gene is not included in the gene trees presented by Abrahams et al., 2020.

4. Despite the use of an impressive number of accessions, sufficiently dense sampling appears to be only available for the Iberian peninsula, the southern coast of Sweden and the south-western coast of Italy. In this regard, the authors overstate their findings (e.g., lines 219ff).

5. While the link between glucosinolates and biotic stress is well documented, evidence for a major role of abiotic factors in shaping glucosinolate profiles is less well substantiated. Now, the authors seem to find that certain chemotypes are favored by high precipitation in 'the North' but that the same chemotypes are favored by low precipitation in 'the South'. This is interesting, but what could explain this unexpected pattern? Any idea? Furthermore, with additional sampling on the eastern coast of Italy and the eastern coast of the Adriatic, would you expect to obtain similar patterns?

6. Please, make sure that the citation style is uniform. Please also make sure that the tense is uniform throughout the manuscript. Please, simplify the introduction in a way that it is accessible by non-specialists. Please improve Figure 5D – this is a mess.

[Editors’ note: further revisions were suggested prior to acceptance, as described below.]

Thank you for submitting your article "Genetic variation, environment and demography intersect to shape Arabidopsis defense metabolite variation across Europe" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Meredith Schuman as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Arthur Korte (Reviewer #2).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential Revisions:

While we find the overall story to be insightful and advancing the field, we do think that some of the analyses still need to be redone, while the eco-evolutionary framework needs to be better defined.

1. We have serious concerns regarding the environmental association analyses. Since there is a general correlation between genetic/genomic/phenotypic diversity and environmental factors due to Arabidopsis demographic history, the overall genetic structure needs to be included in the statistical models testing the associations between traits and environmental factors. You should produce tables similar to that of Supplementary Figure 8A for the chemotypes of Figure 4 and derive statistical significances from these tables, in order to provide convincing evidence of the differential climatic patterns in both European regions.

2. The conceptual framework of repeated evolution that differentiates between parallel and convergent evolution was traditionally developed for interspecific evolution, not for intraspecific evolution where it refers to parallel evolution by multiple independent mutations in the same gene. The differentiation between parallel and convergent evolution for intraspecific evolution is not so straightforward. Therefore we want you to provide a more precise definition of intraspecific convergent and parallel chemotype evolution and outline the role of different or the same genes from a tandem repeat cluster, as well as the role of epistasis between duplicated genes, therein. To this end, the LD analysis between MAM genes and flowering genes should also be removed from the argumentation, because this analysis clearly shows that LD between these genes is not larger than average LD between random pairs of genes. Therefore, the conclusion drawn from this LD analysis is not supported.

We are optimistic that regardless of the outcome of the corrected environmental analysis (point 1), this study represents a valuable contribution, especially if the revised framework for intraspecific repeated evolution – helping to conceptually unite studies of population and functional genetics in an apt system – can be articulated and illustrated more clearly in the revised analysis (point 2).

Reviewer #1 (Recommendations for the authors):

In this new version of the Manuscript by Katz et al., authors address numerous important questions raised in the previous version, and rewrite many of the paragraphs that were confusing or incorrect. In this version authors provide better support for some of their conclusions in a more comprehensive manuscript. However, this Reviewer still has serious concerns regarding the environmental association analyses carried out in this study to support their conclusions. In addition, for clarity and comprehension I also suggest addressing an important conceptual issue, as well as the numerous editorial points described below.

1A. As explained in previous Reviewer 2 report to Authors, in A. thaliana there is a general correlation between genetic/genomic/phenotypic diversity and environmental factors due to Arabidopsis demographic history. For this reason, the overall genetic structure needs to be included in the statistical models testing the associations between traits and environmental factors. Authors now provide more detailed description of the statistical tests applied for such associations. However, as stated in Materials and methods, the statistical models used to build Figure 4 and Supplemental Figure 8B do not include the genomic groups. Therefore, significant associations described in these Figures and throughout the manuscript are mainly reflecting the association of genomic groups (demographic history) with climate. These two figures are the main support for the differential environmental associations in Northern and Southern Europe, but they just show that different genomic groups are present in each geographic area. It seems that such differential association between climate and chemotypes might be true, as shown by the table of Supplementary Figure 8A, where authors include the genomic groups to test the differential association specifically for C3 and C4 status. However, in this table the variable "Precipitation for the driest moth" appears as non-significantly associated. In spite of this, authors choose this variable in Supplementary Figure 8B (as well as in Figure 4A) to illustrate that this climatic variable is associated when applying a simple t-test to compare both groups. This test does not support such environmental association because it is not controlling for the genomic group (demographic history). By contrast, the table of Supplementary Figure 8 shows a stronger differential association in Northern and Southern Europe for minimum temperature, which is associated only in the south but not in the north. Tables similar to that of Supplementary Figure 8A should be provided for the chemotypes of Figure 4, and statistical significances should be derived from them, in order to provide convincing evidence of their differential climatic patterns in both European regions.

1B. Related with the environmental associations, authors include random forest analyses. However, in the way these analyses are presented, they are hard to interpret because variables with the lowest mean decrease Gini (the probability of the variable of being wrong according to Authors) mostly have the lowest mean decrease accuracy (suitability of the variable as predictor).

2. A major goal of this manuscript is the distinction between convergent and parallel evolutions in GSL chemotypes, which are described as two different patterns of repeated evolution. However, these concepts have been classically developed to characterize interspecific evolution, and their translation into the intraspecific variation is not straight forward, thus requiring a precise definition of both alternatives. Authors define parallel and convergent evolution as traits derived from either the same or different "genetic backgrounds". In this definition, the concept of genetic background is very vague and it is unclear what is meant. Authors apply these two concepts to the particular case of repeated evolution of GSL chemotypes evolved specifically from two tandem duplicated genes (MAM and AOP loci). Given the epistasis between these pairs of duplicated genes, the same trait can be derived from the same or different gene within. Then, it seems to this Reviewer that Authors mean by parallel and convergent evolution, the evolution occurring from the same or different gene of the cluster, respectively. In my opinion, this dissection of repeated evolution for the specific case of duplicated genes with epistasis is a major contribution of this study, but it is hard to understand. The manuscript and the audience would benefit by making these concepts clearer. I suggest that Authors add this in discussion, (e.g. paragraph of lines 550-554), explaining the role of different or the same gene from tandem repeat cluster, as well as the role of epistasis between duplicated genes for the distinction of their concepts of parallel and convergent evolution.

Reviewer #2 (Recommendations for the authors):

The authors present a substantial revised version of their manuscript that improved in clarity. All my concerns and suggestions are addressed in this version and I don't have any further comments.

https://doi.org/10.7554/eLife.67784.sa1

Author response

[Editors’ note: the authors resubmitted a revised version of the paper for consideration. What follows is the authors’ response to the first round of review.]

1. We are concerned that MAM3 is not a good proxy for the GS-Elong locus, and there is evidence of recombination within MAM3; and that resequencing of this locus in the accessions studied is required to support the conclusions that the authors wish to draw in this study. Please see detailed comments from Reviewer 3.

2. Several sections of the study require deep revision in terms of writing and analysis. See comment from Reviewer 2, as well as the other reviewers.

We have worked extensively to address both of these concerns in the revised manuscript. There are more details below about the specifics. In summary, we have worked on the entire manuscript to work on the wording and clarification in all noted sections and in the entire manuscript. Specifically, we have tried to make the holistic aspect of the study clearer. For the MAM locus, we have redone the analysis with the flanking gene to test for within locus recombination and how this may influence the analysis. Further, we have obtained 11 new sequences for this locus to further test the use of MAM3 as a proxy. Both approaches support the use of the MAM3 as a proxy to derive general conclusions. Finally, we have reworked all discussion on this locus, the enzymes chemistry and its potential evolution to reflect what is and is not supported by the phylogenetics.

Reviewer #1:

The authors present a detailed and well-written analyses of the variation in the glucosinolate chemotype of different Arabidopsis ecotypes. These differences were linked to the underlying genetic variation. They could show many complex patterns that are related to the distinct chemotypes. The analyses highlights the complexity of the interplay between genetic and environmental factors and gives insights into how traits evolve in natural populations.

I only have one point that I was wondering about.

23 different GSL have been detected and quantified and GWAS has been performed on PCs that summarize the GSL variation. The GWAs on PC1 and PC2 do only identify the two known loci AOP and NAM and no additional loci are significant. As this first two PCs explain only 33 % of the total phenotypic variation, I wonder if GWAS on the individual GSL concentration, especially as all of these show a high heritability, could identify additional associations? It would be nice to have an overview of these potential associations.

To test the reviewer’s hypothesis, we have performed GWAS using the accumulation of each of the 23 individual GSL to test if this identifies more associations both with known GSL genes and unknown genes. This analysis did find a few additional SNPs in unknown genes, it failed to detect additional GSL related genes similar to the use of principal components. We added a supplemental figure of the 23 Manhattan plots showing it, and a supplementary table that summarize the number of SNPs in each of the expected GSL genes (Figure 3 - figure supplement 3 and Supplementary file 2).

Reviewer #2:

This manuscript describes the contribution of environment and demography to shape the genetic variation underlying the glucosinolate (GLS) composition of Arabidopsis thaliana across Eurasia. To this end, the authors carried out a very extensive characterization of different GLSs in a large collection of natural Arabidopsis accessions. Genetic analyses of these data show that the epistasis underlying causal loci hampers the identification of new loci. More importantly, authors compare these data with environmental factors to show significant environmental associations, which support a differential adaptive role of this variation across Europe. In addition, thorough genetic analyses of the major underlying genes identified multiple independent mutations in several duplicated genes, which demonstrate a repeated evolution pattern with differential geographic histories. Overall, this study shows a complex adaptive history of Arabidopsis plant defenses mediated by GLSs at a continental scale and brings forward our current understanding of adaptation of plant defense mechanisms.

Although most sections of the manuscript are written in a comprehensive manner, several sections of results and discussion are rather confusing due to lack of information on particular analyses. The following sections should be revised:

1. The section entitled "Geography and environmental parameters affect GLS variation" (pages 13-14) is very confusing due to an unclear distinction between demography and environment. Given the strong geographic structure of Arabidopsis genetic diversity across Europe (The 1001 Genomes Consortium, 2016), most genetic diversity can be expected to correlate with environment due to Arabidopsis demographic history but not reflecting true adaptive association. For this reason, to distinguish true environmental associations from correlations related with demographic history, analyses must take into account the genetic structure (genetic or genomic relationships among accessions). Usually, the effect of demographic factors on the distribution of genetic diversity is inferred by the genetic structure. According to Supp. Figure 7 and Materials and methods, authors seem to have included such genetic structure in their MANOVA analyses, which is referred to as genomic group. However, authors should describe what are those genomic groups, and if they were included also in Figure 4B. As shown by Supp. Figures 7A and 8, the genetic structure explains most of the diversity in MANOVA and Random forest models, indicating that most of the geographic patterns can be explained by the demographic history. However, both Figures also show that there is still a significant association with environmental factors, which supports that the environment contributes to maintain the genetic diversity of chemotypes. The text describing these results should be carefully rewritten to explain how the genetic structure (demography) has been taken into account and separating clearly demographic and environmental components. The genetic structure is only mentioned in line 259 where it is confusingly named as "ancestral state".

We apologize for not making it clearer that we had included population structure in these analyses and we have worked to clarify this and to ensure that it is clear for each analysis how population structure/demography was included. For the MANOVA and Random Forest models, population structure was included using the genomic groupings as previously defined in The 1001 Genomes Consortium, 2016.

2. The section describing the genetic diversity of MAM3 across Europe (lines 302-326) shows certain overinterpretation of the action of natural selection, which should be revised and rewritten in a softer tone. For instance, the sentence in line 306 stating "If the latter is true, this would argue for a selective pressure shaping this C3/C4 divide". The two scenarios described in paragraph of lines 302-306 could be explained by natural selection or demography, and one should be cautious interpreting geographic patterns. In particular, the geographic pattern observed in Iberia is just showing that this region contains a large diversity for C3/C4 chemotypes and MAM3, as shown previously for most genes in this region (The 1001 Genomes Consortium). Therefore, the sentence in lines 310-313 ("This suggests that the strong geographic partitioning of the C3/C4 chemotypes in Iberia may be driven by selective pressure causing the partitioning of the chemotypes rather than neutral demographic processes") or lines 324-326 ("…,Iberia, showing evidence of local selection while other regions, central Europe, possibly showing a blend requiring further work to delineate") are not fully supported by data. In addition, this section on MAM3 could include the recent study on the genetic diversity of Europe, which suggests that Iberia could have also been colonized from central Europe after the last glaciation (Lee et al., 2017; Nature Communications 8:14458).

We rewrote this section, to both soften and better explain our hypothesis. We emphasized the fact that though there is high diversity in the haplotypes across Iberia, there is a clean delineation of the two opposing chemotypes controlled by these haplotypes. Further, the haplotypes within a chemotype are randomly distributed across the portions of Iberia containing that specific chemotype. This suggests that the haplotypes within a chemotype are more random across the landscape than the chemotype providing evidence supportive of selective delineation but not definitive. We also added the necessary citation.

3. The section describing GS-OH gene and Table 1 (lines 377-390) is rather confuse. In line 377 is stated: "We could not identify the causal LOF in the other six accessions due to sequence quality". However, Table 1 shows the structure of GS-OH gene with apparent location of six mutations. It is unclear which are those six accessions of Table1, what is represented in the structure of the genes, and the meaning of symbols on those genes (A and the vertical bars). Are those accessions of Table 1 the same described in the text? It should also made clear the type of mutations of those six accessions instead of the confusing term "Independent mutation". Furthermore, this section refers to Supplementary Table 2 (line 383), but I could not find such table.

We prepared a new version for table 1. In this table we presented a clearer description of what the mutations are (and which ones were not identifiable from the sequence), we added the genes structure, and emphasized that six accessions with a loss of enzyme activity had an unidentified lesion due to sequence quality for this locus. We added an additional sup table presenting the frequency calculations of the mutated accessions.

Reviewer #3:

The manuscript is based on an impressive data set of (seed) glucosinolate profiles across Europe and beyond, with a focus on methionine-derived glucosinolates. However, the manuscript reads in many aspects like a 'best of' of previous research, with limited net gain in insight. In addition, there are substantial weaknesses (in particular relating to GS-Elong) that prevent me from recommending this manuscript for publication in eLife its current state.

1. The representation of the MAM reaction is misleading and incorrect (lines 101ff). MAMs catalyze condensation of 2-oxoacids with acetylCoA. Two further reaction steps, catalyzed by an isomerase and a dehydrogenase, are necessary to complete chain extension. This extension consists formally in the addition of a methylene group per reaction cycle. With MAM2, the net result is an extension by one (not two!) methylene group, with MAM1 it is mostly two; note that methionine has already two methylene groups in the side chain! I suggest replacing the reference Abrahams et al., 2020, with Graser et al., 2000 (Archives of Biochemistry and Biophysics, 378, 411-419). Furthermore, Kroymann et al., 2001 (Plant Phys. 127, 1077-1088), Benderoth et al., 2006 and Textor et al., 2007 (Plant Phys. 144, 60-71) should be added as original references for the MAM1, 2 and 3 reactions.

We apologize for the error and have rewritten this section, corrected the mistakes, and updated the citations.

2. The use of MAM3 as a proxy for inferring the configuration of the GS-Elong locus is highly debatable. Even though recombination appears to be suppressed at and around GS-Elong, an additional gene (or additional genes) upstream of GS-Elong should be used to construct (a) gene tree(s) and compare with the MAM3 gene tree for congruency. In any case, without sequencing of the MAM genes at the GS-Elong locus, it is not possible to conclude that 'local' configurations of MAM1/MAM2 have not arisen in Iberia (or elsewhere).

To address this concern, we have generated an additional phylogenetic tree using the sequences of MYB37, the gene on the opposite end of the GS-Elong locus (At5g23000), from all the 797 accessions. To test for within locus recombination that may be confusing us, we compared this tree to the MAM3. The vast majority of accessions showed the same major clade memberships in both the MAM3 and MYB37 trees arguing against a high level of within locus recombination. We also reanalyzed our analysis using just the 637 accessions with perfect agreement in the two trees and these showed the same message. Thus, our overarching patterns are not confounded by recombination shuffling this locus.

We realized that we had been inaccurate when discussing a local Iberian haplotype as we’d meant to mean a single haplotype explaining all C3/C4 diversity Iberia. The reviewers comment led us to search for possible true local haplotypes by using the few accessions that have an incongruency in the MYB37/MAM3/C3/C4 data. This identified a collection of 6 accessions in a very narrow region of Iberia that appear to have developed a new local haplotype that requires further sequencing. We added a supp. Figure that describes the detection of these accessions and added this comment to the manuscript. Additionally, we have obtained new genomic sequences of the MAM locus for 11 accessions that resample the existing haplotypes and sample new haplotypes from the clades previously not sampled. These sequences support the classification of the accessions to the clades based on the MAM3 sequences, and the use of the MAM3 gene as a proxy to the MAM1/MAM2 haplotypes. In combination this suggests that while there is at least one local haplotype in Iberia, the vast majority of accessions in Iberia don’t have unique events, and the presence of these few accessions will not change the overall message.

3. The A. lyrata ortholog of Arabidopsis MAM3 is MAMb (not MAMc). The MAM3 gene tree does not provide any statistical support for the branching order. Therefore, it remains unclear whether the presented tree reflects the true branching order and whether the 'most basal' clade in the tree is indeed the most basal. Note also that the Cal-0 gene at the MAM2 position is actually mostly MAM2, except for a converted stretch towards the 3' end of the gene; see: Figure 2 in Kroymann et al., 2003. Hence, the most parsimonious interpretation of available data is that the ancestral state of GS-Elong in Arabidopsis is indeed: MAM2, MAM1 and MAM3. Note also that the Arabidopsis genes in Abrahams et al., 2020, are from TAIR 10, i.e. from Col-0 (see Table 1 in Abrahams et al., 2020). Because Col-0 does not have MAM2, this gene is not included in the gene trees presented by Abrahams et al., 2020.

The MAM3 tree was indeed rooted by A. lyrata MAMb, and not MAMc as was written in the manuscript. This was a typo that was corrected. As for the gene tree presented in Abrahams et al., – we teamed up with them to recreate the tree, this time including MAM2 (supp. Figure x).

For the MAM3 tree, we have worked in several ways to determine the proper support for this tree and the claims. The tree with all accessions has some sequences that are not maximal quality which influences the bootstrapping but not the clade assignment. To improve our precision, we created a tree using 637 accessions that have the highest quality accessions for MAM3 (supp. Figure 10). The accessions in the new tree classified to the same clades as in the original MAM3 tree. The only difference between the trees (original tree and the cleaner tree) was the order of the clades. We then created multiple smaller trees, each contain a different combination of 200300 accessions. In each of these trees the accessions classified to the same clades as the original MAM3 tree, with bootstrap support. This agrees with the reviewer’s assessment that the evolutionary relationships between the clades are not resolvable but that the clades are decently defined. As we don’t have a strong support to the orders of the clades, we soften our tone in the manuscript in the part talking about the clades order and the ancestral state of the locus. We have discussed the most functionally parsimonious model of the ancestral state being MAM2/1/3 but also noted that this model requires a MAM1/2 duplication and diversification following the split from A. lyrata and ensuing structural variation to create the haplotypes.

4. Despite the use of an impressive number of accessions, sufficiently dense sampling appears to be only available for the Iberian peninsula, the southern coast of Sweden and the south-western coast of Italy. In this regard, the authors overstate their findings (e.g., lines 219ff).

We added this clarification to several locations along the text (results and discussion) and soften our claims.

5. While the link between glucosinolates and biotic stress is well documented, evidence for a major role of abiotic factors in shaping glucosinolate profiles is less well substantiated. Now, the authors seem to find that certain chemotypes are favored by high precipitation in 'the North' but that the same chemotypes are favored by low precipitation in 'the South'. This is interesting, but what could explain this unexpected pattern? Any idea? Furthermore, with additional sampling on the eastern coast of Italy and the eastern coast of the Adriatic, would you expect to obtain similar patterns?

We think what is happening is that the abiotic factors are capturing the biotic environment but there are no universal databases describing the biotic environment. We have added the following to try and describe this “The above analysis is extensively using Abiotic factors because of their availability while Aliphatic GSLs have mainly been linked biotic interactions. Our best hypothesis is that the abiotic factors used in this model are capturing variation in the biotic environments. However, an Aliphatic GSL was recently mechanistically linked to drought resistance in Arabidopsis (Salehin et al., 2019). This suggests a need to capture the biotic as well as the abiotic components of any given environment to distinguish what may be driving these patterns.”

6. Please, make sure that the citation style is uniform. Please also make sure that the tense is uniform throughout the manuscript. Please, simplify the introduction in a way that it is accessible by non-specialists. Please improve Figure 5D – this is a mess.

We have worked to correct tense and citation style. We re-organized the whole figure to make it cleaner and clearer. We moved part D to the supplementary so part C (the map) will have more space. Hopefully it helped to make things less messy. We have worked on the introduction and have run it past five individuals who are not glucosinolate specialists and they have helped greatly.

[Editors’ note: what follows is the authors’ response to the second round of review.]

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential Revisions:

While we find the overall story to be insightful and advancing the field, we do think that some of the analyses still need to be redone, while the eco-evolutionary framework needs to be better defined.

1. We have serious concerns regarding the environmental association analyses. Since there is a general correlation between genetic/genomic/phenotypic diversity and environmental factors due to Arabidopsis demographic history, the overall genetic structure needs to be included in the statistical models testing the associations between traits and environmental factors. You should produce tables similar to that of Figure 4—source data 1 for the chemotypes of Figure 4 and derive statistical significances from these tables, in order to provide convincing evidence of the differential climatic patterns in both European regions.

We apologize as we had run all the models with and without genomic group in them and did not include this table. We have rerun the models to include the genomic group to account for population structure. As suggested, for Figure 4 (now Table 1) we created a table in the same manner as the one in supp figure 8 that now accounts for structure. All the p-values in the analyses are derived from models that include the genomic group. We rewrote this section of the results, discussion and the methods accordingly.

2. The conceptual framework of repeated evolution that differentiates between parallel and convergent evolution was traditionally developed for interspecific evolution, not for intraspecific evolution where it refers to parallel evolution by multiple independent mutations in the same gene. The differentiation between parallel and convergent evolution for intraspecific evolution is not so straightforward. Therefore we want you to provide a more precise definition of intraspecific convergent and parallel chemotype evolution and outline the role of different or the same genes from a tandem repeat cluster, as well as the role of epistasis between duplicated genes, therein. To this end, the LD analysis between MAM genes and flowering genes should also be removed from the argumentation, because this analysis clearly shows that LD between these genes is not larger than average LD between random pairs of genes. Therefore, the conclusion drawn from this LD analysis is not supported.

We rewrote the sections concerning parallel/convergent evolution throughout the manuscript (more details are below) and added an introductory figure that will hopefully facilitate a better conveyance of how we are attempting to incorporate these terms in this work. The figure and the section concerning the trans-LD analysis were removed from the manuscript according to the reviewers’ suggestion.

We are optimistic that regardless of the outcome of the corrected environmental analysis (point 1), this study represents a valuable contribution, especially if the revised framework for intraspecific repeated evolution – helping to conceptually unite studies of population and functional genetics in an apt system – can be articulated and illustrated more clearly in the revised analysis (point 2).

Reviewer #1 (Recommendations for the authors):

In this new version of the Manuscript by Katz et al., authors address numerous important questions raised in the previous version, and rewrite many of the paragraphs that were confusing or incorrect. In this version authors provide better support for some of their conclusions in a more comprehensive manuscript. However, this Reviewer still has serious concerns regarding the environmental association analyses carried out in this study to support their conclusions. In addition, for clarity and comprehension I also suggest addressing an important conceptual issue, as well as the numerous editorial points described below.

1A. As explained in previous Reviewer 2 report to Authors, in A. thaliana there is a general correlation between genetic/genomic/phenotypic diversity and environmental factors due to Arabidopsis demographic history. For this reason, the overall genetic structure needs to be included in the statistical models testing the associations between traits and environmental factors. Authors now provide more detailed description of the statistical tests applied for such associations. However, as stated in Materials and methods, the statistical models used to build Figure 4 and Supplemental Figure 8B do not include the genomic groups. Therefore, significant associations described in these Figures and throughout the manuscript are mainly reflecting the association of genomic groups (demographic history) with climate. These two figures are the main support for the differential environmental associations in Northern and Southern Europe, but they just show that different genomic groups are present in each geographic area. It seems that such differential association between climate and chemotypes might be true, as shown by the table of Supplementary Figure 8A, where authors include the genomic groups to test the differential association specifically for C3 and C4 status. However, in this table the variable "Precipitation for the driest moth" appears as non-significantly associated. In spite of this, authors choose this variable in Supplementary Figure 8B (as well as in Figure 4A) to illustrate that this climatic variable is associated when applying a simple t-test to compare both groups. This test does not support such environmental association because it is not controlling for the genomic group (demographic history). By contrast, the table of Supplementary Figure 8 shows a stronger differential association in Northern and Southern Europe for minimum temperature, which is associated only in the south but not in the north. Tables similar to that of Supplementary Figure 8A should be provided for the chemotypes of Figure 4, and statistical significances should be derived from them, in order to provide convincing evidence of their differential climatic patterns in both European regions.

We rerun the models, and now all the models include the genomic group. As suggested, for Figure 4 (now Table 1) we created a table in the same manner as the one in supp figure 8. All the p-values in the analyses are derived from models that include the genomic group. We rewrote the corresponding sections accordingly.

1B. Related with the environmental associations, authors include random forest analyses. However, in the way these analyses are presented, they are hard to interpret because variables with the lowest mean decrease Gini (the probability of the variable of being wrong according to Authors) mostly have the lowest mean decrease accuracy (suitability of the variable as predictor).

We had intended the random forest to convey the differential ranks of the different genetic and environmental parameters across geography. We agree that this is not well shown in the current representation and we removed this analysis from the manuscript.

2. A major goal of this manuscript is the distinction between convergent and parallel evolutions in GSL chemotypes, which are described as two different patterns of repeated evolution. However, these concepts have been classically developed to characterize interspecific evolution, and their translation into the intraspecific variation is not straight forward, thus requiring a precise definition of both alternatives. Authors define parallel and convergent evolution as traits derived from either the same or different "genetic backgrounds". In this definition, the concept of genetic background is very vague and it is unclear what is meant. Authors apply these two concepts to the particular case of repeated evolution of GSL chemotypes evolved specifically from two tandem duplicated genes (MAM and AOP loci). Given the epistasis between these pairs of duplicated genes, the same trait can be derived from the same or different gene within. Then, it seems to this Reviewer that Authors mean by parallel and convergent evolution, the evolution occurring from the same or different gene of the cluster, respectively. In my opinion, this dissection of repeated evolution for the specific case of duplicated genes with epistasis is a major contribution of this study, but it is hard to understand. The manuscript and the audience would benefit by making these concepts clearer. I suggest that Authors add this in discussion, (e.g. paragraph of lines 550-554), explaining the role of different or the same gene from tandem repeat cluster, as well as the role of epistasis between duplicated genes for the distinction of their concepts of parallel and convergent evolution.

We rewrote this section to emphasize the above points. For the discussion we split this into two paragraphs, one on epistasis and one on convergent/parallel to better convey these concepts. We have also reworked the initial paragraph within the introduction to try and translate the interspecific concept to the intraspecific analogy we are contemplating. Specifically that genetic background is a specific haplotype and that if parallel evolution would be a new chemotype arising by independent mutations in a single haplotype to give the new chemotype multiple times. Convergent evolution would be where a new chemotype arose by independent mutations in different haplotypes both giving rise to a single new chemotype. To make sure that our idea is clear we added an introductory figure (Figure 1) that presents how we are trying to incorporate these terms in this work.

https://doi.org/10.7554/eLife.67784.sa2

Article and author information

Author details

  1. Ella Katz

    Department of Plant Sciences, University of California, Davis, Davis, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Visualization, Writing – original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1619-5597
  2. Jia-Jie Li

    Department of Plant Sciences, University of California, Davis, Davis, United States
    Contribution
    Data curation
    Competing interests
    No competing interests declared
  3. Benjamin Jaegle

    Gregor Mendel Institute, Austrian Academy of Sciences, Vienna Biocenter (VBC), Vienna, Austria
    Contribution
    Data curation
    Competing interests
    No competing interests declared
  4. Haim Ashkenazy

    Department of Molecular Biology, Max Planck Institute for Developmental Biology, Tübingen, Germany
    Contribution
    Data curation, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-5079-4684
  5. Shawn R Abrahams

    Division of Biological Sciences, Bond Life Sciences Center, University of Missouri, Columbia, United States
    Contribution
    Data curation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1749-2040
  6. Clement Bagaza

    Division of Biological Sciences, Interdisciplinary Plant Group, Christopher S. Bond Life Sciences Center, University of Missouri, Columbia, United States
    Contribution
    Data curation
    Competing interests
    No competing interests declared
  7. Samuel Holden

    Division of Biological Sciences, Interdisciplinary Plant Group, Christopher S. Bond Life Sciences Center, University of Missouri, Columbia, United States
    Contribution
    Data curation
    Competing interests
    No competing interests declared
  8. Chris J Pires

    Division of Biological Sciences, Bond Life Sciences Center, University of Missouri, Columbia, United States
    Contribution
    Writing – review and editing
    Competing interests
    No competing interests declared
  9. Ruthie Angelovici

    Division of Biological Sciences, Interdisciplinary Plant Group, Christopher S. Bond Life Sciences Center, University of Missouri, Columbia, United States
    Contribution
    Writing – review and editing
    Competing interests
    No competing interests declared
  10. Daniel J Kliebenstein

    1. Department of Plant Sciences, University of California, Davis, Davis, United States
    2. DynaMo Center of Excellence, University of Copenhagen, Frederiksberg, Denmark
    Contribution
    Conceptualization, Data curation, Formal analysis, Funding acquisition, Visualization, Writing – original draft
    For correspondence
    kliebenstein@ucdavis.edu
    Competing interests
    Reviewing editor, eLife
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5759-3175

Funding

National Science Foundation (MCB 1906486)

  • Daniel J Kliebenstein

United States - Israel Binational Agricultural Research and Development Fund (FI-560-2017)

  • Ella Katz
  • Daniel J Kliebenstein

National Science Foundation (IOS 1655810)

  • Daniel J Kliebenstein

National Science Foundation (IOS 1754201)

  • Ruthie Angelovici

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work was supported by the National Science Foundation, Directorate for Biological Sciences, Division of Molecular and Cellular Biosciences (grant no. MCB 1906486 to DJK) and Division of Integrative Organismal Systems (grant no. IOS 1655810 to DJK, and IOS 1754201 to RA), and by the United States-Israel Binational Agricultural Research and Development Fund (to DJK and EK, grant no. FI-560-2017). We thank Dr. Allison Gaudinier (Department of Plant and Microbial Biology, University of California, Berkeley), Dr. Tobias Züst (Institute of Plant Sciences, University of Bern), Dr. Christopher W Wheat (Department of Zoology, Stockholm University) and Dr. Daniel Runcie (Department of Plant Sciences, University of California Davis) for critical reading of the manuscript. We thank the 1001 Genomes Plus project for access for newly assembled A. thaliana accession genomes. The 1001 Genomes Plus project is funded by ERA-CAPS through BBSRC, DFG and FWF to Paul Kersey, Detlef Weigel and Magnus Nordborg, respectively.

Senior and Reviewing Editor

  1. Meredith C Schuman, University of Zurich, Switzerland

Reviewer

  1. Arthur Korte, Vienna Biocenter, Austria

Publication history

  1. Received: February 23, 2021
  2. Accepted: May 2, 2021
  3. Accepted Manuscript published: May 5, 2021 (version 1)
  4. Version of Record published: June 15, 2021 (version 2)

Copyright

© 2021, Katz et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,673
    Page views
  • 277
    Downloads
  • 2
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Ecology
    Lei Guo et al.
    Research Article Updated

    The Varroa destructor mite is a devastating parasite of Apis mellifera honeybees. They can cause colonies to collapse by spreading viruses and feeding on the fat reserves of adults and larvae. Amitraz is used to control mites due to its low toxicity to bees; however, the mechanism of bee resistance to amitraz remains unknown. In this study, we found that amitraz and its major metabolite potently activated all four mite octopamine receptors. Behavioral assays using Drosophila null mutants of octopamine receptors identified one receptor subtype Octβ2R as the sole target of amitraz in vivo. We found that thermogenetic activation of octβ2R-expressing neurons mimics amitraz poisoning symptoms in target pests. We next confirmed that the mite Octβ2R was more sensitive to amitraz and its metabolite than the bee Octβ2R in pharmacological assays and transgenic flies. Furthermore, replacement of three bee-specific residues with the counterparts in the mite receptor increased amitraz sensitivity of the bee Octβ2R, indicating that the relative insensitivity of their receptor is the major mechanism for honeybees to resist amitraz. The present findings have important implications for resistance management and the design of safer insecticides that selectively target pests while maintaining low toxicity to non-target pollinators.

    1. Ecology
    2. Evolutionary Biology
    Hanna M Bensch et al.
    Research Article

    Living with relatives can be highly beneficial, enhancing reproduction and survival. High relatedness can, however, increase susceptibility to pathogens. Here, we examine whether the benefits of living with relatives offset the harm caused by pathogens, and if this depends on whether species typically live with kin. Using comparative meta-analysis of plants, animals, and a bacterium (nspecies = 56), we show that high within-group relatedness increases mortality when pathogens are present. In contrast, mortality decreased with relatedness when pathogens were rare, particularly in species that live with kin. Furthermore, across groups variation in mortality was lower when relatedness was high, but abundances of pathogens were more variable. The effects of within-group relatedness were only evident when pathogens were experimentally manipulated, suggesting that the harm caused by pathogens is masked by the benefits of living with relatives in nature. These results highlight the importance of kin selection for understanding disease spread in natural populations.