Abstract
Over very short spatial scales, the habitat of a species can differ in multiple abiotic and biotic factors. These factors may impose natural selection on several traits and can cause genetic differentiation within a population. We studied multivariate genetic differentiation in a plant species of a sand dune landscape by linking environmental variation with differences in genotypic trait values and gene expression levels to find traits and candidate genes of microgeographical adaptation.
Maternal seed families of Arabidopsis lyrata were collected in Saugatuck Dunes State Park, Michigan, USA, and environmental parameters were recorded at each collection site. Offspring plants were raised in climate chambers and exposed to one of three temperature treatments: regular occurrence of frost, heat, or constant control conditions. Several traits were assessed: plant growth, time to flowering, and frost and heat resistance.
The strongest trait–environment association was between a fast switch to sexual reproduction and weaker growth under frost, and growing in the open, away from trees. The second strongest association was between the trait combination of small plant size and early flowering under control conditions combined with large size under frost, and the combination of environmental conditions of growing close to trees, at low vegetation cover, on dune bottoms. Gene expression analysis by RNA-seq revealed candidate genes involved in multivariate trait differentiation.
The results support the hypothesis that in natural populations, many environmental factors impose selection, and that they affect multiple traits, with the relative direction of trait change being complex. The results highlight that heterogeneity in the selection environment over small spatial scales is a main driver of the maintenance of adaptive genetic variation within populations.
INTRODUCTION
Environmental variation in both abiotic and biotic factors may result in heterogeneous natural selection over small spatial scales, which in turn may cause genetic differentiation within populations (Hedrick, 2006). Divergent adaptive evolution within populations is most likely in species with limited dispersal of gametes or offspring (Lenormand, 2002). Adaptive differences to local environmental conditions have been suggested to be a crucial aspect underlying community structure and ecosystem functioning (Menge et al., 2002; Fenu et al., 2012). Furthermore, divergent selection over a fine spatial scale may help maintain genetic polymorphism within populations (Hedrick, 2006). Fine-scale adaptive genetic differentiation has been detected in several plant species using common garden studies of populations inhabiting strong environmental gradients, such as shady versus sunny sites (Donohue et al., 2000), different edaphic conditions (Antonovics and Bradshaw, 1970) or between the top and bottom of dunes (Paccard et al., 2013). Studies have mostly focused on adaptive genetic differentiation along single environmental variables, but conditions may actually change in a multivariate manner. Here we assess the relative importance of several abiotic and biotic variables in explaining genotypic trait differentiation in a plant population inhabiting a sand dune landscape, and use transcriptome data to compile a list of genes and pathways probably involved in within-population differentiation.
Well-established theory predicts the conditions under which fine-scale adaptation and adaptive genetic differentiation is possible. Richardson et al. (2014) reviewed four conditions that contribute to microgeographical adaptation. Two refer to the intensity of selection and the amount of gene flow. When microhabitats differ and some genotypes perform better in one but not in the other habitat, genetic differentiation may arise due to divergent selection. However, genetic differentiation may be impeded by gene flow (Lenormand, 2002). If gene flow between two microhabitats is large compared to selection, the evolution of phenotypic plasticity should be favoured (Crispo, 2008). Therefore, adaptive genetic differentiation is more likely to occur when gene flow is limited by ecological or physical barriers. Two other conditions can promote microhabitat adaptation. The first is spatially auto-correlated selection, where microhabitats with similar biotic or abiotic characteristics tend to be spatially aggregated. In this case, one microhabitat mainly receives immigrants from nearby similar microhabitats, which increases differentiation from more distant locations (Urban, 2011). The second is habitat choice, where a particular genotype preferentially migrates into an environment that matches its phenotype (Edelaar et al., 2008). In plants, habitat choice can occur on a very fine scale, by directional growth of foraging organs or clonal growth (Bazzaz, 1991).
In plants, spatial heterogeneity has long been known to cause within-population trait variation (Hutchings et al., 2000). Differences may be due to a plastic response to the local environment or due to past divergent selection and, consequently, genetic differences. Adaptive genetic differences are typically revealed by rearing plants from different origins under standard conditions. For instance, flowering time in common garden experiments has been shown to vary with soil temperature in Agrostis sp. and Mimulus guttatus (Tercek and Whitbeck, 2004; Lekberg et al., 2012), with soil moisture content in Arabidopsis lyrata (Paccard et al., 2013), and with light availability in Impatiens capensis (Schmitt, 1993). Microgeographical adaptation in response to edaphic heterogeneity and soil-water availability has been demonstrated in Plantago lanceolata (Ravenscroft et al., 2014). Because the natural environment of herbaceous plants is often heterogeneous, and gene flow is limited over short spatial scales (Vekemans and Hardy, 2004), phenotypic differences in phenology and stress resistance within populations may be commonly caused by divergent natural selection.
Empirical studies focusing on more complex environmental gradients and their link to multivariate trait variation have so far been rare in the context of microhabitat adaptation. The classic demonstrations of local adaptation on complex gradients in plants have been conducted over relatively large spatial scales and come from forestry research, such as on lodgepole pine, Pinus contorta. When raised in a common garden, pines from a latitudinal gradient differed in a combination of phenotypic traits involving frost resistance, vegetative growth and phenology, and differences were linked to multiple climatic factors of the sites of origin, including temperature and precipitation (Liepe et al., 2016). Transplant experiments on the same species showed that plant height was affected by combinations of climatic factors involving mean annual temperature and soil moisture (Wang et al., 2006, 2010).
While common garden experiments reveal the quantitative genetic contribution to phenotypic differentiation, population genomic studies have focused on determining the genes underlying such differentiation (Stinchcombe and Hoekstra, 2007). Candidate genes underlying functional traits may be associated with changes in nucleotide polymorphism and gene expression level (Mitchell-Olds and Schmitt, 2006). Insight into the genetic basis of phenotypic trait differentiation is typically gained by combining common garden experiments with genome-wide associations or transcriptome analyses (de Villemereuil et al., 2016). For instance, Oleksiak et al. (2002) tested for differential gene expression between individuals within teleost fish populations that differed in some physiological and cold adaptive traits. They found that a relatively small number of genes exhibiting large differences in expression levels were relevant for explaining variation. So far, candidate genes underlying phenotypic trait differentiation at very fine spatial scales remain largely unknown. In principle, studies targeting the genetic basis of fine-scale adaptive genetic differentiation should have a high potential to find candidate genes because when adaptation evolves despite some gene flow, mainly genes of large phenotypic effect should differ, while the genetic background should be similar (Yeaman and Whitlock, 2011). Here, we relate expression differences with variables summarizing the environment and the phenotype, using overlap between gene sets to narrow down the list of candidate genes.
Our study organism was the perennial plant Arabidopsis lyrata ssp. lyrata. We focused on a population inhabiting a partially forested sand dune landscape on the shore of Lake Michigan, USA. First, we performed a common garden experiment to study the link between several environmental parameters and multiple traits depicting aspects of life history and thermal stress performance. Second, we studied gene expression differences among plants that were most extreme with respect to their local environment and their traits linked to environmental factors. Arabidopsis lyrata experiences spatially variable environmental conditions in at least three respects. First, plants in open, unshaded areas are exposed to high irradiation, dryer conditions and greater temperature fluctuation (Maun, 2009). Second, plants on dune crests are exposed to stronger winds, dryer conditions and more temperature variation than plants in dune swales (Leege and Murphy, 2001; Maun, 2009). Finally, plants differ in exposure to competition, depending on the density of conspecifics or grasses and other herbs while others grow with hardly any potential interspecific competitors, which is likely to affect growth and reproductive traits (Miller and Schemske, 1990). Spatial autocorrelation in these micro-habitat variables is not strong. This population is outcrossing and has limited gene flow beyond 10 m (Willi and Määttänen, 2010, 2011; Paccard et al., 2013). These conditions suggest that the association between traits and the environment is the result of past divergent natural selection and adaptive genetic differentiation.
MATERIAL AND METHODS
Sampling design and plant rearing
In June 2010, two healthy fruits of 40 plants of Arabidopsis lyrata ssp. lyrata were collected over an area of 2.5 ha at Saugatuck Dunes State Park, on Lake Michigan, USA (42°42′N, 86°12′W). At the site of each maternal plant, four environmental parameters were recorded (Fig. 1): distance from the ground-projected crown edge (distance from trees, in metres), which depicts the gradient from closed to open habitat; fraction down from the dune top (dune position, value between 0 and 1); coverage by herbs and grasses on a quarter square-metre with the maternal plant in the centre, as an indicator of interspecific competition (vegetation cover; in per cent); and the number of A. lyrata plants on the same surface area as an indicator of intraspecific competition, multiplied by four to get an estimate per square-metre (intraspecific density).
Location of the 40 seed families of Arabidopsis lyrata collected in Saugatuck Dunes State Park, Michigan, USA, with information on the four environmental variables studied: distance from the crown edge of the nearest tree (in metres, m), position on dunes (0 = top of dunes, 1 = bottom of dunes), per cent coverage by herbs and grasses (vegetation cover) and density of A. lyrata plants (per quarter m2).
Seedlings of the 40 seed families were raised under three experimental treatments: frost stress, heat stress and control conditions, with three replicates for each treatment (n = 40 families × 3 treatments × 3 replicates, each in a different block = 360 plants). Two seeds were sown in each of nine pots (7 cm diameter, 5 cm depth; 1: 1 peat/sand ratio) per seed family to ensure that there was at least one seedling per pot. The pots were then split over three blocks within five larger holding trays per block. The three pots of the same family within a block were later exposed to one of the three treatments. Pot position was random over the five trays. Seeds were stratified for 1 week (4 °C in darkness under wet conditions) and then transferred to a growth chamber for germination (Grobank, CLF, Wertingen, Germany); trays of a block were put on a separate shelf in the growth chamber. Initial conditions mimicked short day length: 8 h light/16 h dark; light intensity: 150 µmol m–2 s−1; 18 °C constant; relative humidity (RH): 40–70 %. Germination started after 4 d and was scored every day for 2 weeks. After 10 d when germination rate had reached 80 %, day length and temperature were increased: 12 h light/12 h dark; light intensity: 180 µmol m–2 s−1, 20 °C/18 °C day/night, RH: 40–70 %. At the same time, we haphazardly removed surplus seedlings from each pot and, when possible, transplanted left-over seedlings in pots of the same maternal line with no germination. At the end, we had 353 pots with one seedling each. Plants were watered twice a week for the entire experiment.
Stress treatments
When plants were at the four-leaf stage, 4 weeks after the onset of germination, one of three treatments were applied: regular frost stress, regular heat stress or constant temperature. Stress temperatures were chosen to reflect short frost events experienced by plants in the early morning or short heat events during midday, as they can occur during the phase of vegetative growth in spring and early summer (average temperature in mid-June during the hottest two hours, from 1500 to 1700h: 43.7 °C on a dune top site, 39.6 °C at the edge of a small patch of trees). All pots of a stress treatment were moved to a separate climate chamber a few hours before stress application. Temperature was changed with a ramping rate of 0.5 °C min–1. For frost stress, temperature was first decreased to 0 °C for 1 h then to −3 °C for 1 h, to 0 °C for 1 h, and back to 18 °C; treatment ended 1 h before the lights were turned on. Heat stress began 3 h after the lights were turned on, and temperature was increased to 30 °C for 1 h, then to 47 °C for 1 h, to 30 °C for 1 h and back to 20 °C. At the end of each treatment cycle, plants were placed back in the growth chamber. Stress temperatures were applied on three subsequent days per week, for 3 weeks.
Phenotypic trait assessment
Asymptotic rosette size.
Every week, starting before treatment application and for a duration of 5 weeks, pictures of holding trays were taken. On each picture and plant, the length of the two longest leaves was measured. The time series of mean leaf length for each plant was used to fit seven growth models, of which the two-parameter logistic growth model was best supported by the data (Wos and Willi, 2018). The two parameters of that model are asymptotic rosette size and the rate of exponential growth. We only used asymptotic rosette size as an estimate of growth performance, for the three different treatments (sizecontrol, sizefrost, sizeheat).
Temperature stress resistance.
Thermal stress resistance was estimated by measuring percentage of electrolyte leakage (PEL) at the end of the growth treatments, 7 weeks after germination (Wos and Willi, 2015). Frost resistance was measured on plants raised under regular frost exposure only, and heat resistance was measured on plants raised under regular heat exposure only. For each plant, two discs were excised from the fifth rosette leaf and put in two separate tubes. We exposed the two discs to one of two treatments: (1) control: incubation at 20 °C for 1 h; (2) stress: incubation at −14 °C for 1 h for the frost stress or 47 °C for 1 h for the heat stress. After incubation, we measured conductivity for the first time (Conductivity Meter FiveEasy FE30, Mettler Toledo, Columbus, OH, USA). We then placed the tubes in a boiling bath for 30 min and measured conductivity a second time. We measured PEL as the ratio of the conductivity before boiling to that after the boiling bath (Cornelissen et al., 2003). We calculated resistance as PEL of the control disc minus PEL of the treated disc, so that a low value corresponded to low resistance. For another study on the same plants, thermal stress resistance was also estimated for plants raised under control conditions. Results revealed that thermal stress resistance was significantly increased after stress pre-exposure during growth compared to control plants (Wos and Willi, 2018). Here only thermal resistance of acclimated plants was used for further analyses (frost resistancefrost, heat resistanceheat).
Time to flowering.
We used flowering time of control plants (flowering timecontrol) as another life history trait that was defined as the number of days between germination and the appearance of the first open flower. As flowering is induced by exposure to long day length in A. lyrata, we gradually increased the day length and temperature regime 2 weeks after the end of the growth treatments (16 h/8 h light/dark; 22 °C/18 °C day/night; light intensity: 240 µmol m–2 s−1; RH: 40–70 %). Flowering time was recorded three times per week until the end of the experiment by day 390 after sowing. At that time, 97 out of 120 plants had started flowering. The remaining alive and healthy plants (17 of 23 plants) were assigned a date of flowering 2 months later.
Spatial autocorrelation.
We tested whether maternal line means of the six phenotypic traits and data of the four environmental variables were spatially independent. For this, we performed a Moran’s I test that estimated the similarity between observations to detect spatial autocorrelation (Moran, 1950). We tested for spatial autocorrelation among maternal lines with less than 12 m distance (Paccard et al., 2013). Moran’s I statistic can range between −1 and 1; a value of 0 indicates spatial randomness. Moran’s I for each phenotypic and environmental variable varied between −0.20 and 0.27, indicating no significant relationship among neighbouring maternal lines (Supplementary Data Table S1).
Statistical analysis
Mixed model analysis.
Mixed model analysis tested the relationship between the four environmental variables (x variables): distance from trees, dune position, vegetation cover and intraspecific density; and the six phenotypic variables (y variables): sizecontrol, sizefrost, sizeheat, flowering timecontrol, frost resistancefrost and heat resistanceheat. Three x and two y variables were log-transformed to approach a normal distribution: distance form trees, vegetation cover, intraspecific density, flowering timecontrol and heat resistanceheat. The random effects were plant nested within seed family, and seed family, and the fixed effects were block and the four continuous environmental variables. The mixed model was run with the GLIMMIX procedure of SAS v9.4 (SAS Institute, 2011). We did not correct for multiple testing as the purpose of the analysis was to detect any association between environmental variables and traits, including those with high residual variation.
Canonical correlation analysis.
We used canonical correlation analysis (CCA) to describe the relationship between the four environmental variables depicting spatial environmental heterogeneity and the six phenotypic traits. Canonical correlation analysis is a type of multivariate analysis for studying linear relationships between sets of x and y variables (Hotelling, 1936). Canonical variables are produced based on linear combinations within each set of variables (Hair et al., 2010). The canonical correlation coefficient (Rc) represents the optimal linear combination between canonical variables of the two data sets. CCA assumes independent variables within each set to minimize problems with multicollinearity (Hair et al., 2010). Therefore, we computed Pearson correlation coefficients for pairs of variables of the two data sets (Table S2). For the set of phenotypic variables, two significant correlations were detected, with fairly low coefficients: between flowering timecontrol and sizecontrol (r = 0.44, P < 0.01), and between sizeheat and sizecontrol (r = 0.34, P < 0.05). No significant correlation was detected among the set of environmental variables. Canonical correlation analysis was based on seed family means for phenotypic traits (n = 40) and was therefore based on genotypic trait values. We used the F-approximation of Wilks’ lambda to test for the significance of the canonical correlation coefficients. To determine the relative importance of each original variable to the canonical correlation, we computed the canonical loadings and the canonical cross-loadings. Canonical loadings depict the linear correlation between variables and their respective canonical variables and are considered significant if they are greater than |0.40| (Hair et al., 2010). Canonical cross-loadings, calculated by the product of the canonical loading and the canonical coefficient, are the linear correlations between the variables with the opposite canonical variables and represent a more direct measure of the relationship between the two sets of variables. Cross-loadings above |0.25| are generally considered significant. Analysis were performed with the package CCA in R (González et al., 2008; R Core Team, 2013).
Gene expression analysis
We conducted gene expression analysis to identify candidate genes underlying the associations between the environmental variables and genotypic trait variables revealed by CCA.
RNA extraction, library construction and sequencing.
For plants raised under control conditions, leaf tissue from 16 randomly chosen maternal lines of the first block were collected 7 weeks after germination, prior to PEL measurement. Samples were immediately snap-frozen in liquid nitrogen for RNA extraction (RNeasy Mini Kit protocol, Qiagen, Hilden, Germany). RNA samples were subjected to a DNase treatment (RNase-Free DNase Set, Qiagen). RNA concentration and integrity were checked with a Nanodrop 2000 (Thermo Scientific, Wilmington, DE, USA) and a Bioanalyzer 2100 (Agilent Technologies, Palo Alto, CA, USA). cDNA libraries for each individual were prepared using the Illumina TruSeq Stranded mRNA Sample Preparation Kit (catalogue number: RS-122–2101) and specific TruSeq adapters for individual sequencing. Sequencing was done on an Illumina NextSeq 500 platform (Illumina, San Diego, CA, USA) and produced 75-bp single-end reads. Raw data were filtered to remove low-quality reads. Raw sequence data are available in the EMBL Sequence Read Archive (http://www.ebi.ac.uk/ena) under accession PRJEB13660 (Wos and Willi, 2018).
Alignment of reads and differential gene expression.
Quality of the sequence data was checked with the software FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and TruSeq adapters were trimmed with the software cutadapt (Martin, 2011). We used the read-mapping tool TopHat v2.1.0 (Trapnell et al., 2012) for read alignment on the A. lyrata reference genome (Hu et al., 2011) with the version-2 annotation (Rawat et al., 2015). The number of reads mapped on each gene was counted with HTSeq-count v0.6.1p2 (Anders et al., 2015) and analysed with EdgeR v3.12.0 (Robinson et al., 2010) for gene expression differences. For each significant canonical correlation, we conducted two association analyses, one with the canonical variable of the environment and one with the canonical variable of traits. For each variable, canonical scores of each maternal line – depicting the direction and degree of involvement on the canonical variable – were extracted, and the 16 randomly chosen maternal lines were sorted according to their canonical scores (Table S3). Based on the sorting, two groups were established, one with high and one with low scores. Groups were then compared for expression differences. For each canonical correlation, genes were considered as candidates if they were significant in both association analyses and had the same sign of expression difference. P-values were adjusted for multiple testing with the Benjamini and Hochberg false discovery rate correction (FDR) and the significance level was set at FDR ≤ 0.01. Gene description was obtained from The Arabidopsis Information Resource (Berardini et al., 2015). Gene Ontology (GO) term enrichment analysis for biological processes was performed on the list of candidate genes with BiNGO v3.0.3 software in cytoscape v3.5.1 (Maere et al., 2005). GO terms were considered significantly enriched if the FDR-adjusted P-value was <0.05.
RESULTS
Phenotypic traits and their relationship with environmental variables
Results of univariate mixed model analysis testing the relationship between the environmental variables and the phenotypic traits are shown in Table 1. Several significant relationships were found: Distance from trees was associated with plant size under control conditions, with plants originating away from trees growing to larger size (Fig. 2A). Dune position and intraspecific density were linked to flowering time, with plants originating from dune bottoms and sites with high intraspecific density flowering earlier (Fig. 2B, C). In addition, vegetation cover was associated with frost resistance, with plants from sites with high vegetation cover being more frost resistant (Fig. 2D). Furthermore, four non-significant trends (P < 0.1) were found: two negative trends, between distance from trees and both size after regular frost and frost resistance, and two positive trends, between vegetation cover and both size under control conditions and time to flowering. Size after regular exposure to heat and heat resistance were unrelated to any of the environmental variables.
Results of mixed model analysis testing the relationship between four environmental variables: distance from the crown edge (distance from trees), the fraction down on dunes (dune position), per cent of ground covered by grass and herbs (vegetation cover) and density of Arabidopsis lyrata plants (intraspecific density); and six phenotypic traits: rosette size of plants under control conditions (sizecontrol), under regular frost stress (sizefrost) and under regular heat stress (sizeheat), flowering time of control plants, frost resistance of frost-acclimated plants (frost resistancefrost) and heat resistance of heat-acclimated plants (heat resistanceheat) (n = 117, 116, 115, 113, 111 and 114 respectively)
| Dependent variable | d.f.num | d.f.den | Sizecontrol | Sizefrost | Sizeheat | Flowering timecontrol | Frost resistancefrost | Heat resistanceheat |
|---|---|---|---|---|---|---|---|---|
| F | F | F | F | F | F | |||
| Block | 2 | 69 to 75 | 6.12 (0.003) | 4.23 (0.018) | 3.79 (0.027) | 2.56 (0.085) | 20.07 (<0.001) | 0.73 (0.484) |
| t | t | t | t | t | t | |||
| Distance from trees | 1 | 35 | 2.06 (0.047) | −1.96 (0.058) | −0.32 (0.755) | −1.50 (0.142) | −1.76 (0.087) | 1.67 (0.103) |
| Dune position | 1 | 35 | −0.51 (0.614) | 0.64 (0.524) | −0.18 (0.860) | −2.98 (0.005) | 0.58 (0.566) | 0.00 (1.000) |
| Vegetation cover | 1 | 35 | 1.90 (0.065) | −0.49 (0.626) | 0.43 (0.673) | 1.80 (0.080) | 2.51 (0.017) | −1.21 (0.235) |
| Intraspecific density | 1 | 35 | −0.29 (0.770) | −0.97 (0.341) | 0.96 (0.344) | −2.52 (0.016) | −1.23 (0.228) | −0.69 (0.493) |
| Dependent variable | d.f.num | d.f.den | Sizecontrol | Sizefrost | Sizeheat | Flowering timecontrol | Frost resistancefrost | Heat resistanceheat |
|---|---|---|---|---|---|---|---|---|
| F | F | F | F | F | F | |||
| Block | 2 | 69 to 75 | 6.12 (0.003) | 4.23 (0.018) | 3.79 (0.027) | 2.56 (0.085) | 20.07 (<0.001) | 0.73 (0.484) |
| t | t | t | t | t | t | |||
| Distance from trees | 1 | 35 | 2.06 (0.047) | −1.96 (0.058) | −0.32 (0.755) | −1.50 (0.142) | −1.76 (0.087) | 1.67 (0.103) |
| Dune position | 1 | 35 | −0.51 (0.614) | 0.64 (0.524) | −0.18 (0.860) | −2.98 (0.005) | 0.58 (0.566) | 0.00 (1.000) |
| Vegetation cover | 1 | 35 | 1.90 (0.065) | −0.49 (0.626) | 0.43 (0.673) | 1.80 (0.080) | 2.51 (0.017) | −1.21 (0.235) |
| Intraspecific density | 1 | 35 | −0.29 (0.770) | −0.97 (0.341) | 0.96 (0.344) | −2.52 (0.016) | −1.23 (0.228) | −0.69 (0.493) |
The table shows F-values for the effect of block, t-values for all other effects, and P-values in parentheses. Statistics for the random effects are not shown. Significant P-values (α-level of 0.05) are indicated in bold type.
Results of mixed model analysis testing the relationship between four environmental variables: distance from the crown edge (distance from trees), the fraction down on dunes (dune position), per cent of ground covered by grass and herbs (vegetation cover) and density of Arabidopsis lyrata plants (intraspecific density); and six phenotypic traits: rosette size of plants under control conditions (sizecontrol), under regular frost stress (sizefrost) and under regular heat stress (sizeheat), flowering time of control plants, frost resistance of frost-acclimated plants (frost resistancefrost) and heat resistance of heat-acclimated plants (heat resistanceheat) (n = 117, 116, 115, 113, 111 and 114 respectively)
| Dependent variable | d.f.num | d.f.den | Sizecontrol | Sizefrost | Sizeheat | Flowering timecontrol | Frost resistancefrost | Heat resistanceheat |
|---|---|---|---|---|---|---|---|---|
| F | F | F | F | F | F | |||
| Block | 2 | 69 to 75 | 6.12 (0.003) | 4.23 (0.018) | 3.79 (0.027) | 2.56 (0.085) | 20.07 (<0.001) | 0.73 (0.484) |
| t | t | t | t | t | t | |||
| Distance from trees | 1 | 35 | 2.06 (0.047) | −1.96 (0.058) | −0.32 (0.755) | −1.50 (0.142) | −1.76 (0.087) | 1.67 (0.103) |
| Dune position | 1 | 35 | −0.51 (0.614) | 0.64 (0.524) | −0.18 (0.860) | −2.98 (0.005) | 0.58 (0.566) | 0.00 (1.000) |
| Vegetation cover | 1 | 35 | 1.90 (0.065) | −0.49 (0.626) | 0.43 (0.673) | 1.80 (0.080) | 2.51 (0.017) | −1.21 (0.235) |
| Intraspecific density | 1 | 35 | −0.29 (0.770) | −0.97 (0.341) | 0.96 (0.344) | −2.52 (0.016) | −1.23 (0.228) | −0.69 (0.493) |
| Dependent variable | d.f.num | d.f.den | Sizecontrol | Sizefrost | Sizeheat | Flowering timecontrol | Frost resistancefrost | Heat resistanceheat |
|---|---|---|---|---|---|---|---|---|
| F | F | F | F | F | F | |||
| Block | 2 | 69 to 75 | 6.12 (0.003) | 4.23 (0.018) | 3.79 (0.027) | 2.56 (0.085) | 20.07 (<0.001) | 0.73 (0.484) |
| t | t | t | t | t | t | |||
| Distance from trees | 1 | 35 | 2.06 (0.047) | −1.96 (0.058) | −0.32 (0.755) | −1.50 (0.142) | −1.76 (0.087) | 1.67 (0.103) |
| Dune position | 1 | 35 | −0.51 (0.614) | 0.64 (0.524) | −0.18 (0.860) | −2.98 (0.005) | 0.58 (0.566) | 0.00 (1.000) |
| Vegetation cover | 1 | 35 | 1.90 (0.065) | −0.49 (0.626) | 0.43 (0.673) | 1.80 (0.080) | 2.51 (0.017) | −1.21 (0.235) |
| Intraspecific density | 1 | 35 | −0.29 (0.770) | −0.97 (0.341) | 0.96 (0.344) | −2.52 (0.016) | −1.23 (0.228) | −0.69 (0.493) |
The table shows F-values for the effect of block, t-values for all other effects, and P-values in parentheses. Statistics for the random effects are not shown. Significant P-values (α-level of 0.05) are indicated in bold type.
Relationship between environmental variables measured in the field and family-mean trait values: between distance from trees and rosette size under control conditions (A), dune position and time to flowering under control conditions (0 = top of dunes, 1 = bottom of dunes) (B), density of A. lyrata and time to flowering (C), and per cent coverage by herbs and grass and frost resistance after acclimation (D). Symbols depict family means (n = 40) and lines the linear regression relationship.
CCA on phenotypic traits and environmental variables
The main results of the canonical correlation analysis are presented in Table 2. The first approximate F value of 2.21 corresponded to the test that all four canonical correlations were zero; the hypothesis was rejected at the 0.05 significance level. The first canonical correlation was 0.68 and represents the highest possible correlation between any linear combination of the phenotypic variables and any linear combination of the environmental variables. The second approximate F value of 1.89 corresponded to the test that the second to fourth correlations were zero; this hypothesis was also rejected. All further tests, on the third and fourth canonical correlations, failed to reject the null hypothesis, meaning that only the first and the second canonical correlations were significant.
Statistics for each of the four canonical correlation coefficients
| Rc | R2 | Stat | F-approx. | d.f.num/d.f.den | P | |
|---|---|---|---|---|---|---|
| 1 | 0.682 | 0.465 | 0.243 | 2.206 | 24/105.868 | 0.003 |
| 2 | 0.575 | 0.331 | 0.455 | 1.895 | 15/85.979 | 0.035 |
| 3 | 0.514 | 0.264 | 0.679 | 1.708 | 8/64.000 | 0.114 |
| 4 | 0.278 | 0.077 | 0.923 | 0.918 | 3/33.000 | 0.443 |
| Rc | R2 | Stat | F-approx. | d.f.num/d.f.den | P | |
|---|---|---|---|---|---|---|
| 1 | 0.682 | 0.465 | 0.243 | 2.206 | 24/105.868 | 0.003 |
| 2 | 0.575 | 0.331 | 0.455 | 1.895 | 15/85.979 | 0.035 |
| 3 | 0.514 | 0.264 | 0.679 | 1.708 | 8/64.000 | 0.114 |
| 4 | 0.278 | 0.077 | 0.923 | 0.918 | 3/33.000 | 0.443 |
We used F-approximation of Wilks’ lambda to test for the significance of each coefficient (Menzel, 2012). The table shows the canonical correlation coefficient (Rc), proportion of shared variance (R2), Wilks’ lambda statistic (Stat), F-approximation for Wilks’ lambda (F-approx.), numerator degrees of freedom (d.f.num) and denominator degrees of freedom (d.f.den) and significance probability based on F-approximation. Significant P-values (α-level of 0.05) are indicated in bold type.
Statistics for each of the four canonical correlation coefficients
| Rc | R2 | Stat | F-approx. | d.f.num/d.f.den | P | |
|---|---|---|---|---|---|---|
| 1 | 0.682 | 0.465 | 0.243 | 2.206 | 24/105.868 | 0.003 |
| 2 | 0.575 | 0.331 | 0.455 | 1.895 | 15/85.979 | 0.035 |
| 3 | 0.514 | 0.264 | 0.679 | 1.708 | 8/64.000 | 0.114 |
| 4 | 0.278 | 0.077 | 0.923 | 0.918 | 3/33.000 | 0.443 |
| Rc | R2 | Stat | F-approx. | d.f.num/d.f.den | P | |
|---|---|---|---|---|---|---|
| 1 | 0.682 | 0.465 | 0.243 | 2.206 | 24/105.868 | 0.003 |
| 2 | 0.575 | 0.331 | 0.455 | 1.895 | 15/85.979 | 0.035 |
| 3 | 0.514 | 0.264 | 0.679 | 1.708 | 8/64.000 | 0.114 |
| 4 | 0.278 | 0.077 | 0.923 | 0.918 | 3/33.000 | 0.443 |
We used F-approximation of Wilks’ lambda to test for the significance of each coefficient (Menzel, 2012). The table shows the canonical correlation coefficient (Rc), proportion of shared variance (R2), Wilks’ lambda statistic (Stat), F-approximation for Wilks’ lambda (F-approx.), numerator degrees of freedom (d.f.num) and denominator degrees of freedom (d.f.den) and significance probability based on F-approximation. Significant P-values (α-level of 0.05) are indicated in bold type.
The pair of canonical variables (of traits and environment) for each canonical correlation were investigated for their representation of the different original variables (Table 3). Loading and cross-loading coefficients of the first canonical variable of the trait data set revealed that plant size under frost and flowering under control conditions contributed the most to the canonical correlation. Plants that flowered early had smaller rosettes under frost stress. For the second canonical variable of the trait data set, timing of flowering, plant size under control conditions and size under frost stress contributed most to the canonical correlation, with plants that grew to large size and flowered late under favourable conditions being smaller under frost stress. The coefficients for the environmental variables showed that distance from trees was the most important predictor of the first canonical correlation, and vegetation cover, dune position and distance from tress the most important predictors of the second canonical correlation. In conclusion, the first canonical correlation depicted an association between the trait combination of fast reproductive development combined with reduced size under frost stress and the environmental condition of growing in open habitats, away from trees. The second canonical correlation depicted an association between the trait combination of small size and early flowering under control conditions, but large size under frost stress and proximity to trees where vegetation cover was low, on dune bottoms.
Canonical loadings and cross-loadings for the two first canonical correlations
| 1st canonical correlation | 2nd canonical correlation | |||
|---|---|---|---|---|
| Trait variables | Loadings | Cross-loadings | Loadings | Cross-loadings |
| Sizecontrol | −0.37 | −0.25 | −0.75 | −0.43 |
| Sizefrost | 0.49 | 0.33 | 0.43 | 0.25 |
| Sizeheat | −0.03 | −0.02 | 0.02 | 0.01 |
| Flowering timecontrol | 0.50 | 0.34 | −0.79 | −0.46 |
| Frost resistancefrost | 0.38 | 0.26 | −0.32 | −0.18 |
| Heat resistanceheat | −0.22 | −0.15 | −0.00 | −0.00 |
| Environmental variables | Loadings | Cross-loadings | Loadings | Cross-loadings |
| Distance from trees | −0.79 | −0.54 | −0.48 | −0.27 |
| Dune position | −0.12 | −0.08 | 0.57 | 0.33 |
| Vegetation cover | −0.16 | −0.11 | −0.76 | −0.44 |
| Intraspecific density | −0.35 | −0.24 | 0.41 | 0.23 |
| 1st canonical correlation | 2nd canonical correlation | |||
|---|---|---|---|---|
| Trait variables | Loadings | Cross-loadings | Loadings | Cross-loadings |
| Sizecontrol | −0.37 | −0.25 | −0.75 | −0.43 |
| Sizefrost | 0.49 | 0.33 | 0.43 | 0.25 |
| Sizeheat | −0.03 | −0.02 | 0.02 | 0.01 |
| Flowering timecontrol | 0.50 | 0.34 | −0.79 | −0.46 |
| Frost resistancefrost | 0.38 | 0.26 | −0.32 | −0.18 |
| Heat resistanceheat | −0.22 | −0.15 | −0.00 | −0.00 |
| Environmental variables | Loadings | Cross-loadings | Loadings | Cross-loadings |
| Distance from trees | −0.79 | −0.54 | −0.48 | −0.27 |
| Dune position | −0.12 | −0.08 | 0.57 | 0.33 |
| Vegetation cover | −0.16 | −0.11 | −0.76 | −0.44 |
| Intraspecific density | −0.35 | −0.24 | 0.41 | 0.23 |
The trait variables were: rosette size of plants under control conditions (sizecontrol), under regular frost stress (sizefrost) and under regular heat stress (sizeheat), flowering time of control plants, frost resistance of frost-acclimated plants (frost resistancefrost) and heat resistance of heat-acclimated plants (heat resistanceheat). The environmental variables were: distance from the crown edge (distance from trees), the fraction down on dunes (dune position), per cent of ground covered by grass and herbs (vegetation cover) and density of Arabidopsis lyrata plants (intraspecific density). n = 40. Variables were considered significant (in bold) when both canonical loadings and cross-loadings were greater than |0.40| and |0.25|, respectively.
Canonical loadings and cross-loadings for the two first canonical correlations
| 1st canonical correlation | 2nd canonical correlation | |||
|---|---|---|---|---|
| Trait variables | Loadings | Cross-loadings | Loadings | Cross-loadings |
| Sizecontrol | −0.37 | −0.25 | −0.75 | −0.43 |
| Sizefrost | 0.49 | 0.33 | 0.43 | 0.25 |
| Sizeheat | −0.03 | −0.02 | 0.02 | 0.01 |
| Flowering timecontrol | 0.50 | 0.34 | −0.79 | −0.46 |
| Frost resistancefrost | 0.38 | 0.26 | −0.32 | −0.18 |
| Heat resistanceheat | −0.22 | −0.15 | −0.00 | −0.00 |
| Environmental variables | Loadings | Cross-loadings | Loadings | Cross-loadings |
| Distance from trees | −0.79 | −0.54 | −0.48 | −0.27 |
| Dune position | −0.12 | −0.08 | 0.57 | 0.33 |
| Vegetation cover | −0.16 | −0.11 | −0.76 | −0.44 |
| Intraspecific density | −0.35 | −0.24 | 0.41 | 0.23 |
| 1st canonical correlation | 2nd canonical correlation | |||
|---|---|---|---|---|
| Trait variables | Loadings | Cross-loadings | Loadings | Cross-loadings |
| Sizecontrol | −0.37 | −0.25 | −0.75 | −0.43 |
| Sizefrost | 0.49 | 0.33 | 0.43 | 0.25 |
| Sizeheat | −0.03 | −0.02 | 0.02 | 0.01 |
| Flowering timecontrol | 0.50 | 0.34 | −0.79 | −0.46 |
| Frost resistancefrost | 0.38 | 0.26 | −0.32 | −0.18 |
| Heat resistanceheat | −0.22 | −0.15 | −0.00 | −0.00 |
| Environmental variables | Loadings | Cross-loadings | Loadings | Cross-loadings |
| Distance from trees | −0.79 | −0.54 | −0.48 | −0.27 |
| Dune position | −0.12 | −0.08 | 0.57 | 0.33 |
| Vegetation cover | −0.16 | −0.11 | −0.76 | −0.44 |
| Intraspecific density | −0.35 | −0.24 | 0.41 | 0.23 |
The trait variables were: rosette size of plants under control conditions (sizecontrol), under regular frost stress (sizefrost) and under regular heat stress (sizeheat), flowering time of control plants, frost resistance of frost-acclimated plants (frost resistancefrost) and heat resistance of heat-acclimated plants (heat resistanceheat). The environmental variables were: distance from the crown edge (distance from trees), the fraction down on dunes (dune position), per cent of ground covered by grass and herbs (vegetation cover) and density of Arabidopsis lyrata plants (intraspecific density). n = 40. Variables were considered significant (in bold) when both canonical loadings and cross-loadings were greater than |0.40| and |0.25|, respectively.
Gene expression differences
First canonical correlation.
Gene expression analysis revealed a total of 115 differentially expressed genes (DEGs) associated with the canonical variable of traits. Fifty-eight genes were up-regulated in the plants with early flowering and small size under frost and 57 genes were down-regulated (Table S4). We found 63 DEGs associated with the canonical variable of the environment, with 40 genes being up-regulated in plants away from trees and 23 genes being down-regulated (Table S5). Nineteen genes overlapped between the two subsets of up-regulated genes, and ten genes overlapped between the two subsets of down-regulated genes (Table S6). The GO term enrichment analysis on these 29 genes revealed 15 enriched GO terms (Table S7). The GO term ‘response to stimulus’ included the highest number of genes, seven in total, five of which were also part of ‘response to chemical stimulus’: a down-regulated manganese transporter MTP11, two up-regulated drought and salt responsive genes [LIPID TRANSFER PROTEIN 4, NAC DOMAIN CONTAINING PROTEIN (NAC) 6], an up-regulated pathogen-responsive gene (BASIC PATHOGENESIS-RELATED PROTEIN 1) and an up-regulated (biotic) stress-induced senescence gene [SENESCENCE-ASSOCIATED GENE 12 (SAG12)]. The two genes specific to ‘response to stimulus’ were down-regulated disease resistance genes (AT3G51570 and AT4G19520). Additional GO terms contained the same genes as the GO term ‘response to stimulus’: the additional GO terms were related to immune system, senescence or cell death, response to ethylene and manganese ion homeostasis.
Second canonical correlation.
Gene expression analysis revealed a total of 68 DEGs for the canonical variable of traits. Forty-one genes were up-regulated in the plants of small size and early flowering under control conditions and large size under frost stress; 27 genes were down-regulated (Table S8). We found 104 DEGs associated with the canonical variable of the environment. Forty-one genes were up-regulated in plants close to trees where vegetation cover was low, on dune bottoms, and 63 genes were down-regulated (Table S9). We found that six genes overlapped between the two subsets of up-regulated genes and four genes overlapped between down-regulated genes (Table S10). Among the six up-regulated genes, three had a known function and were involved in RNA processing (RIBOSOMAL PROTEIN S5), DNA or RNA metabolism and biotic stress response (DIR1-LIKE) and response to stimulus [NUDIX HYDROLASE HOMOLOG 9 (NUDIX9]). Among the four down-regulated genes, two had a known function and were linked to signal transduction (RHOGAP FAMILY PROTEIN) and the response to disease (DISEASE RESISTANCE PROTEIN).
DISCUSSION
Our study on microgeographical adaptive differentiation in one population of Arabidopsis lyrata in a partially forested sand dune landscape revealed two main axes of trait–environment association. The first was between fast reproduction under control conditions but weak performance under frost, and growing in the open, away from trees. The second strongest association was between the trait combination of small size and fast reproduction under control conditions and large size under frost, and the environmental variables of proximity to trees, low vegetation cover, on dune bottoms. The genetic basis of these associations was further explored by gene expression analysis, which revealed candidate genes related to signalling pathways.
Relationship between habitat heterogeneity and life history traits and thermal performance
Our study was based on maternal seed families collected over a dune landscape and assessment of phenotypic traits in offspring plants. Working with seed families has the benefit that genotypic trait values can be estimated via family means. However, because seeds had been collected in the field, differences due to maternal effects cannot be excluded. Maternal effects commonly affect traits of early life stages such as seed size but to a lesser extent than those of later life stages (El-Keblawy and Lovett-Doust, 1999). In line with this, a previous study on A. lyrata reported no significant effect of seed size on phenotypic traits such as rosette size and flowering time (Paccard et al., 2013). Also, the seeds were the result of some gene flow as the population studied is outcrossing. This may have weakened trait–environment associations if pollen had come from sites differing in the environment. In this regard, we can say that any association found may actually be stronger in nature due to ongoing selection. Overall, our results established a clear relationship between spatial environmental heterogeneity and timing of flowering, rosette size under benign conditions and coping with short bouts of frost. The results are consistent with local abiotic and biotic conditions selecting for distinct genotypes and contributing to ecotype formation (Turesson, 1922; Hufford and Mazer, 2003).
Spatial heterogeneity in abiotic conditions such as climate, topology or substrate, and in biotic assemblage strongly influence life history traits in most species (Hutchings et al., 2000). The transition from the vegetative to the reproductive phase in response to complex environmental cues represents an especially crucial aspect of the life history of plants. To maximize reproductive output, the timing of flowering needs to be adjusted such that it occurs during favourable environmental conditions, and – in obligate outcrossers – when other plants flower as well (Augspurger, 1981; Putterill et al., 2004). Our results emphasize the genetic control of differentiation in the timing of flowering in a heterogeneous landscape and suggest that the difference in phenology arose from divergent selection experienced in the field. Selection seems to act such that early flowering is favoured at sites some distance from trees, and also at sites offering the combination of conditions of close to trees and low vegetation cover, at dune bottoms. Previous studies on other herbaceous plants suggest that selection regimes represented by these environmental variables may be linked to extreme temperature and drought conditions away from trees, and shading under trees (Schmitt, 1993; Lekberg et al., 2012). Genetically based differences may have evolved as a direct adaptive strategy to escape temperature and drought stress, and to compensate for the environmental effect of shading on flowering in order to maintain synchrony in flowering within the population (countergradient variation; Conover et al., 2009).
We found that frost tolerance was another trait involved in micro-habitat adaptation. Landscape topography, duration of snow cover, or the presence of a canopy creates microhabitats that greatly vary in the severity or duration of extreme freezing events (Sakai and Larcher, 2012). Previous studies reported significant relationships between single environmental variables and freezing resistance. Yellow cedars were studied along a canopy-cover gradient, and trees from colder open-canopy sites were found to be more frost resistant than trees from closed-canopy stands (Schaberg et al., 2005). Furthermore, along a gradient of timing of snow melt, plants of Aciphylla glacialis were more frost resistant when they originated from early snow-melt sites compared to those from late snow-melt sites (Briceño et al., 2014). Topographic depressions may generally have lower temperatures in the winter (Sakai and Larcher, 2012), and therefore plants on dune bottoms with some trees probably experience the coldest conditions during a large fraction of the year. In line with this, plants at forested dune bottoms were more frost tolerant, as their rosettes grew to larger size under regular occurrence of experimental frosts. In contrast, for sites away from trees, we found that frost tolerance was reduced. Conditions away from trees may select for fast reproductive development at the cost of frost tolerance, while conditions at forested dune bottoms seem to select both for fast reproductive development and for more frost tolerance. These complex associations between sets of environmental differences and trait differences emphasize the importance of considering multivariate environmental gradients and the assessment of several aspects of phenotype in the study of microhabitat adaptation.
Finally, we found no significant trait–environment association involving traits in coping with heat stress. Plants inhabiting sand dune habitats frequently experience extremely high soil and air temperatures in summer (Maun, 2009). Therefore, we expected to find adaptive genetic differentiation in coping with heat stress. We worked only with one heat stress regime (and one frost regime), but we think it simulated well an aspect of environmental stress as it occurs in nature. Physiological studies conducted on several sand-dune plants demonstrated that they have high concentrations of sugars and polyols, which helps to maintain enzyme activity and increases cell membrane stability under stress (Smirnoff and Stewart, 1985). Plants specialized to sand dunes may constitutively express high degrees of heat resistance with minor differences between microhabitats. However, only a few studies are available for comparison and more data on micro-climate and local heat performance are needed.
Gene expression differences
Analysis of gene expression differences on the first canonical correlation depicting the link between fast reproduction and weak performance under frost, and growing in open habitat, revealed a total of 29 candidate genes. The GO term enrichment analysis pointed to the most important groups of genes based on functioning: response to stimulus, genes involved in immunity, in senescence and cell death, in response to ethylene and manganese ion homeostasis. In plants away of trees with early flowering and reduced frost tolerance, two drought- and salt-responsive genes were up-regulated (LIPID TRANSFER PROTEIN 4, NAC 6). Furthermore, four genes related to coping with pathogens or biotic stress were either up- or down-regulated (up-regulated: BASIC PATHOGENESIS-RELATED PROTEIN 1, SAG12; down-regulated: AT3G51570, AT4G19520), and one gene related to manganese ion homeostasis was down-regulated (MTP11). The same genes are known to be involved in other functions, such as senescence and cell death. No genes linked to vegetative growth or flowering were detected. The reason may be that we missed some signals due to the timing of leaf collection for expression analysis, which was after the main growth phase and some time before bolting. However, what we detected was a signal of genes involved in copying with drought and a faster life cycle, which matches the ecology of plants living away of trees. They experience drier conditions during late spring and summer (Leege and Murphy, 2001). To escape drought, they seem to accelerate the switch to reproduction (drought escape strategy; Franks, 2011) and as a result may senesce faster. Senescence plays a major role in plant development because of the remobilization of nutrients to other parts of the plants, and one member of the senescence-associated gene family was previously reported as a candidate for local adaptation in A. thaliana and linked to survival (Fournier-Level et al., 2011). Also in A. thaliana, it was found that plants that bolt early also senesce early (Levey and Wingler, 2005), but this result may be linked to the annual life cycle of that species. Beyond this, senescence was shown to be induced by many different factors including some abiotic stressors, such as water, salt or temperature, that may accelerate the process of senescence (Buchanan-Wollaston, 1997; Noodén, 2012).
Transcriptome analysis on the second canonical correlation depicting the link between small size, early flowering under control conditions and large size under frost, and proximity to trees, with low vegetation cover and on dune bottoms, revealed a total of 10 candidate genes, with one down-regulated disease resistance gene (AT4G19520) overlapping with the candidates of the first canonical correlation. The specific candidates are involved in stress response, RNA processing or signal transduction. A second disease resistance gene was up-regulated, DIR1-LIKE, involved in systemic acquired resistance (Champigny et al., 2013). Disease resistance genes often exhibit strong polymorphisms and are an important class of genes involved in the process of local adaptation (Ross-Ibarra et al., 2008). The other candidate genes were: RIBOSOMAL PROTEIN S5 involved in RNA degradation (Tsuzuki et al., 2017); NUDIX9 involved in a wide range of metabolic processes (McLennan, 2006); and the RhoGAP family protein probably involved in cell growth (Tcherkezian and Lamarche-Vane, 2007). The list of genes suggests that trait differences between tree-free dune tops and forested dune bottoms may go back to genes generally important in plant development rather than being specific to a particular pathway or stress.
CONCLUSIONS
Our study described habitat heterogeneity based on several environmental factors and linked it to microgeographical trait differentiation. The focus was on multivariate genotypic trait differences, to determine how selection by the environment targets different life history traits, and thermal stress performance. We found that several environmental variables were linked to genotypic trait differences in Arabidopsis lyrata over a sand dune landscape, suggesting complex divergent selection and adaptive genetic differentiation over a scale of a few metres. We conclude that: (1) in a heterogeneous landscape, a combination of environmental factors imposes selection on traits; (2) the same traits can respond to combinations of environmental variables, here mostly timing of flowering and frost tolerance, but the signs can change between trait variables; and (3) the genetic basis of trait differentiation – even though the same traits are involved – seems to be specific to the environmental gradient.
SUPPLEMENTARY DATA
Supplementary data are available online at https://dbpia.nl.go.kr/aob and consist of the following. Table S1: Spatial autocorrelation analysis on environmental and trait data. Table S2: Matrices of Pearson correlation coefficients for environmental and trait data. Table S3: Canonical scores. Tables S4–S6: Gene expression analysis on the first canonical correlation. Table S7: Gene Ontology (GO) term enrichment analysis on the first canonical correlation. Tables S8–S10: Gene expression analysis on the second canonical correlation.
ACKNOWLEDGEMENTS
We thank Reyhan Sonmez and Olivier Bachmann for their help with measuring plant traits. Collection permit was granted by the Michigan Department of Natural Resources. Library construction and sequencing were performed at the Department of Biosystems Science and Engineering (D-BSSE) of ETH Zürich in Basel and the University of Basel. This work was supported by the Swiss National Science Foundation (PP00P3-123396, PP00P3_146342, 31003A_166322), and the Fondation Pierre Mercier pour la Science, Lausanne.
