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.

Cyanogenesis in T. repens is a polymorphic trait regulated by independently segregating polymorphisms at the Ac and Li 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 Ac based on its historical use in the literature (Corkill 1942; Daday 1954a; Hughes 1991). The second locus (Li) 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., Ac–Li–) 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). 

CALCULATION OF ALLELE FREQUENCIES

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 T. repensbetween 2008 and 2018 (Table S2). In T. repens’ 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 Ac and Li alleles within each population was inferred from cyanogenesis phenotyping with Feigl-Anger assays (Feigl and Anger 1966; Gleadow et al. 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₂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₂O and 80µL of 0.2 EU ml^–1 linamarase (LGC Standards CDX-00012238-100, Teddington, UK), or 50 µL of H₂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 (Ac and Li), the only genotypes we could determine with certainty were the homozygous recessives acac and lili (i.e., from Acli, acLi, and acli cyanotypes). Therefore, we calculated allele frequencies by first determining the frequency of the homozygous recessive genotypes (acac and lili) within populations, and then assumed loci were in Hardy-Weinberg equilibrium to determine the frequencies of Ac and Li. 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 T. repens (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).

HERBIVORY MEASUREMENTS

Herbivory was measured in the native range of T. repens 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, Ac, and Li 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).

CLIMATE DATA

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 raster 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. 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 T. repens 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.

See README.