DOI: 10.5061/dryad.76hdr7szt
title: Evolution in response to climate in the native and introduced ranges of a globally distributed plant

description: The extent to which species can adapt to spatiotemporal climatic variation
 in their native and introduced ranges remains unresolved. To address this,
 we examined how clines in cyanogenesis (HCN production—an antiherbivore
 defense associated with decreased tolerance to freezing) have shifted in
 response to climatic variation in space and time over a 60-year period in
 both the native and introduced ranges of Trifolium repens. HCN
 production is a polymorphic trait controlled by variation at two Mendelian
 loci (Ac and Li). Using phenotypic assays, we estimated
 within-population frequencies of HCN production and dominant alleles at
 both loci (i.e., Ac and Li) from 10,575 plants
 sampled from 131 populations on 5 continents, and then compared these
 frequencies to those from historical data collected in the
 1950s. There were no clear relationships between changes
 in the frequency of HCN production, Ac,
 or Li and changes in temperature between contemporary
 and historical samples. We did detect evidence of continued evolution to
 temperature gradients in the introduced range, whereby the slope of
 contemporary clines for HCN and Ac in relation to winter
 temperature became steeper than historical clines and more similar to
 native clines. These results suggest that cyanogenesis clines show no
 clear changes through time in response to global warming, but introduced
 populations continue to adapt to their contemporary environments.
description: Cyanogenesis
 in <em>T. repens</em> is a polymorphic trait regulated by
 independently segregating polymorphisms at
 the <em>Ac</em> and <em>Li</em> loci (Mirande
 1912; Armstrong et al. 1913; Corkill 1942; Atwood and Sullivan 1943).
 Expression of this polymorphism is controlled by the epistatic interaction
 between these two loci. The first locus encodes a closely linked
 three-gene cluster for the biosynthesis of the cyanogenic glucosides
 lotaustralin and linamarin (Melville and Doak 1940; Olsen et al. 2008;
 Olsen and Small 2018)—we denote this locus <em>Ac</em> based
 on its historical use in the literature (Corkill 1942; Daday 1954a; Hughes
 1991). The second locus (<em>Li</em>) encodes the hydrolyzing
 enzyme linamarase (Coop 1940; Olsen et al. 2007), which cleaves sugar
 moieties from lotaustralin and linamarin leading to the production of HCN.
 Because cyanogenic glucosides are stored in vacuoles, whereas linamarase
 is found in the apoplast (Kakes 1985; Gruhnert et al. 1994), HCN is only
 produced following tissue lysis when cyanogenic glucosides and linamarase
 have the opportunity to interact (Hughes 1991). HCN interrupts the
 electron transport chain in cellular respiration by inhibiting
 cytochrome-c oxidase (Antonini et al. 1971), making it a potent
 antiherbivore defense (Angeesing 1974; Dirzo and Harper 1982; Thompson and
 Johnson 2016; Santangelo et al. 2019). The presence of a single dominant
 allele at both loci
 (i.e., <em>Ac</em>–<em>Li</em>–) confers a plant
 with the cyanogenic phenotype (cyanotype AcLi: Corkill 1942). Individuals
 homozygous for gene deletions at either locus are acyanogenic (cyanotypes
 Acli, acLi and acli; Olsen et al. 2007, 2008, 2013).    <strong>CALCULATION
 OF ALLELE FREQUENCIES</strong> The classic studies of
 Hunor Daday served as the source for all historical allele frequencies for
 this project (Daday 1954a, 1958). We extracted allele frequency data from
 the appendices of his 1954a and 1958 papers, for all countries shown in
 Fig. 2. In regions where historical sampling was sparse (e.g., New Zealand
 and South America), we chose to keep locations from outside of our
 contemporary sampling to increase statistical power in these sampling
 regions. We excluded historical samples from the UK and Ireland because
 there were many populations (n = 12) not represented in our sampling,
 which could have biased our results. Regardless, our results are
 qualitatively identical even if these samples are included (Table
 S1). Contemporary samples were taken from the native and introduced ranges
 of <em>T. repens</em>between 2008 and 2018 (Table S2).
 In <em>T. repens</em>’ native range, samples were collected in
 summer 2018 from sites spanning continental Europe, from the south of
 Spain to the north of Sweden, and west from the Atlantic coast to eastern
 Europe (Fig. 2). At each site, ~20 ripe infructescences were collected,
 ensuring that plants were at least 3 m apart to avoid sampling the same
 clone. Seeds were returned to the University of Toronto, Canada, scarified
 with sandpaper and germinated in 100 mL pots filled with Pro-Mix LP15
 (Premier Tech, Rivière-du-Loup, Canada) potting soil. Plants were grown in
 Conviron MTPS (Winnipeg, Canada) environmental chambers at 350 µmol of
 light with an 18hr day: 6hr night cycle at 25°C. Each individual received
 3-5 pellets of Nutricote Total 14-13-13 (Arysta, New York, USA) once its
 first true leaf unfurled. Tissue was collected for cyanogenesis assays
 when plants had produced at least 3-4 medium sized leaves. 
             Samples from
 the introduced range were obtained from North America, South America,
 Japan and New Zealand between 2008 and 2018 with the aim of resampling the
 geographic range investigated by Daday (1958) in each region. In North
 America, samples were taken from the Pacific Northwest, and the central
 and eastern half of Canada and the United States, extending south to
 Louisiana and Florida, and north to the coast of Hudson’s Bay in Manitoba
 and Ontario (Fig. 2). In South America, sampling was done from
 south-central Chile to northern Argentina, west to the Pacific coast and
 east to central Argentina (Fig. 2). In Japan, samples were taken from the
 southern tip of Kyushu to northern Hokkaido (Fig. 2). Finally, in New
 Zealand samples were taken from Otago in the south to Pegasus Bay on the
 South Island (Fig. 2). A combination of stolon cuttings and mature
 infructescences were collected from each population, maintaining at least
 3 m between plants. When stolon cuttings were collected, tissue was
 immediately stored in a −80°C freezer and then assayed. When
 infructescences were collected, seeds were removed, germinated and grown
 to the 3-4 leaf stage. From these plants, we collected, stored and assayed
 fresh tissue in an identical manner as described above for stolons
 collected from natural populations. A total of 6,947 stolons were
 collected from 35 populations originating from North America and New
 Zealand (mean ± SE: 212 ± 32.7, range = 31–1032). A total of 3,628
 individuals from 96 populations originating from North America, South
 America, Japan and Europe (mean ± SE: 37 ± 1.18, range = 9–54) were grown
 from seed. In all ranges, sampling of rural sites was prioritized where
 possible since recent research shows that cyanogenesis frequencies evolve
 in response to urbanization (Thompson et al. 2016; Johnson et al. 2018;
 Santangelo et al. 2020).             The frequency
 of
 dominant <em>Ac</em> and <em>Li </em>alleles within each population was inferred from cyanogenesis phenotyping with Feigl-Anger assays (Feigl and Anger 1966; Gleadow et al<em>.</em> 2011), using an optimized method from Thompson et al. (2016). Briefly, leaf tissue from individuals was placed into every second well of 96-well plates, leaving empty wells between each sample to prevent bleeding on the assay paper between plants. We added 80 µL of H<sub>2</sub>O to each well and macerated the tissue within them. A pre-cut filter paper saturated, then dried, with Feigl-Anger reagent was then placed over each 96-well plate and placed in an incubator at 37°C. After 3 hours, individuals were scored for the presence of HCN by visual inspection of a blue colour on the filter paper (cyanotype AcLi). For all acyanogenic individuals, successive tests were conducted by adding 20 µL of H<sub>2</sub>O and 80µL of 0.2 EU ml<sup>–1</sup> linamarase (LGC Standards CDX-00012238-100, Teddington, UK), or 50 µL of H<sub>2</sub>O and 30µL of 10 mM linamarin (Sigma-Aldrich 68264, St-Louis, USA) to identify Acli and acLi cyanotypes, respectively. If an individual failed to produce any positive score after all three assays, its genotype was deemed homozygous recessive at both loci (i.e., acli cyanotype). Previous studies that have compared cyanotype based on Feigl-Anger assays with genotypes using PCR based assays have found near perfect correspondence between the two (Olsen et al. 2007, 2008; Kooyers and Olsen 2012; Thompson et al. 2016). Daday used a comparable but different assay method, utilizing picric acid (Corkill 1940), to assess the presence of dominant alleles. Although both Feigl-Anger and picric acid give comparable results (Hughes 1991), the former is preferred because picric acid is potentially explosive. We compared both tests on a subset of samples (107 individuals) and found 96.3% of plants showed concordant results between methods (Text S1). Since cyanogenesis is controlled by the presence of a single dominant allele at both loci underlying this trait (<em>Ac</em> and <em>Li</em>), the only genotypes we could determine with certainty were the homozygous recessives <em>acac</em> and <em>lili</em> (i.e., from Acli, acLi, and acli cyanotypes). Therefore, we calculated allele frequencies by first determining the frequency of the homozygous recessive genotypes (<em>acac</em> and <em>lili</em>) within populations, and then assumed loci were in Hardy-Weinberg equilibrium to determine the frequencies of <em>Ac</em> and <em>Li</em>. This assumption is reasonable given the obligately outcrossing mating system, the random (non-assortative) movement of pollinators among cyanotypes and the typically large population sizes of <em>T. repens</em> (Burdon 1983; Lachance 2009; Santangelo et al. 2019). In total, 10,575 individuals were assayed as the contemporary sample from 131 populations spanning 5 continents (mean no. individuals per population ± SE: 80.7 ± 10.7, range = 9–1032).   <strong>HERBIVORY MEASUREMENTS</strong> Herbivory was measured in the native range of <em>T. repens</em> during 2018 to assess how the frequency of HCN production within populations is related to the amount of leaf tissue eaten by herbivores. We focused our estimates of herbivory on native populations for logistical reasons, and because we expect that the frequency of HCN production, <em>Ac</em>, and <em>Li</em> is most likely to be in equilibrium due to a long history of selection by herbivores and temperature in the native range. Estimates of consumed leaf tissue (hereafter herbivory) were visually scored to the nearest percentage of leaf damage, which is an accurate method for quantifying herbivory that is comparable to digital methods (Johnson et al. 2016). At each European site, one haphazardly selected leaf from each of ~20 individuals spaced at least 3m apart, was scored for leaf damage to the nearest 1%. Plants used for herbivory estimates were haphazardly selected from the same populations used for the collection of infructescences. Individuals were not scored if leaf damage was found to be the result of mowing (i.e., leaf damage showed a dried and frayed edge). Plants were scored from July 7 to August 28 of 2018, beginning with the southernmost sites in Spain. Sampling progressed northward and eastward before ending in northern Sweden. This approach of sampling in a northward direction followed the methods of Anstett et al. (2014), and allowed us to ensure that plants experienced a roughly comparable amount of the growing season. To control for confounding differences in sampling date, we used growing degree days experienced by each population in 2018 as a covariate in our statistical analyses of herbivory. Growing degree days (GDD) is a metric of thermal accumulation through time and is frequently used to predict plant development and insect population growth (McMaster and Smika 1988). We calculated GDD based on the thermal accumulation above 5ºC with an upper threshold of 30ºC during 2018 until the date of sampling (Mence 1964; McMaster and Wilhelm 1997), extracted from the E-OBS daily gridded temperature dataset (Cornes et al. 2018). In total, 974 individuals from 50 populations were scored for herbivory (mean no. individuals per population ± SE: 19.5 ± 2.13, range = 10–20).   <strong>CLIMATE DATA</strong> Cyanogenesis clines have repeatedly evolved in response to environmental gradients (Daday 1954a,b, 1958; Ganders 1990; Kooyers and Olsen 2012, 2013). We extracted a number of climatic variables to investigate their possible role in generating the observed clines. All climate data were extracted from the CRU TS v4.02 dataset and managed using the <em>raster </em>package in R v. 3.5.2 (Climatic Research Unit, University of East Anglia; Harris et al. 2020; R Core Team 2018; Hijmans 2019). These data are monthly means for a range of climatic variables (e.g., mean temperatures, precipitation, etc.) at 0.5° resolution covering all land areas except for Antarctica. Data were extracted for each population using corresponding longitudes and latitudes. Time calibrated five-year averages were used for each time of collection by taking the average of mean winter temperature (January and July for northern and southern hemispheres, respectively) from 1950–1954 and 2013–2017 for our historical and contemporary data, respectively. Calibrating average mean winter temperature in this way allowed us to account for temporal changes in temperature due to climate change within each population when testing for continued adaptation to spatial gradients in temperature in the introduced range.<em> </em>Mean monthly precipitation and potential evapotranspiration were also extracted and averaged from 1950–1954 and from 2013–2017 for historical and contemporary populations, respectively, and were used to create an aridity index (precipitation/potential evapotranspiration). Extracting temperature data for each population from the time at which it was sampled allowed us to test our research question about whether <em>T. repens</em> has continued to adapt to spatial gradients in temperature in the introduced range while controlling for climate change through time. Data were processed and analyzed in R v3.5.2. All code associated with processing and analysis of these data can be found in this Github repository.
description: See README.