Adaptive differentiation of Festuca rubra along a climate gradient revealed by molecular markers and quantitative traits.

Species response to climate change is influenced by predictable (selective) and unpredictable (random) evolutionary processes. To understand how climate change will affect present-day species, it is necessary to assess their adaptive potential and distinguish it from the effects of random processes. This will allow predicting how different genotypes will respond to forecasted environmental change. Space for time substitution experiments are an elegant way to test the response of present day populations to climate variation in real time. Here we assess neutral and putatively adaptive variation in 11 populations of Festuca rubra situated along crossed gradients of temperature and moisture using molecular markers and phenotypic measurements, respectively. By comparing population differentiation in putatively neutral molecular markers and phenotypic traits (QST-FST comparisons), we show the existence of adaptive differentiation in phenotypic traits and their plasticity across the climatic gradient. The observed patterns of differentiation are due to the high genotypic and phenotypic differentiation of the populations from the coldest (and wettest) environment. Finally, we observe statistically significant covariation between markers and phenotypic traits, which is likely caused by isolation by adaptation. These results contribute to a better understanding of the current adaptation and evolutionary potential to face climate change of a widespread species. They can also be extrapolated to understand how the studied populations will adjust to upcoming climate change without going through the lengthy process of phenotyping.

Introduction

Understanding how climate change will affect the evolution of existing populations is one of the greatest challenges of evolutionary ecology at present. Climate change will result in novel environments but also increased environmental variation, both of which will generate new selective pressures [1]. Existing populations will have to adjust to these environmental changes in order to avoid extinction. Predictions of the adaptive response of existing populations to these new environmental conditions have practical application in conservation biology or management of invasive species, but will also help us to understand the fundamental evolutionary and ecological mechanisms of rapid adjustments to environmental variation [2,3].

Three mechanisms are commonly distinguished as possible responses to climate change–migration, genetic adaptation and adaptive plasticity [3-5]. In plants, it has been shown that migration is limited–for many plant species, the migration front does not allow to keep up with the pace of predicted local climate change [6]. Genetic adaptation is the result of selection acting on a population or taxon causing heritable modification of the trait(s) under selection, which will, in turn, increase individual fitness [7]. If a population has sufficient genetic variation upon which selection can act, adaptation can occur on a relatively small time scale [3,8]. Organisms can also respond to climate change through phenotypic plasticity–the ability of a given genotype to modify its phenotype in response to environmental variation. If a plastic genotype has higher fitness across variable environments than a non-plastic one, then plasticity is adaptive [9,10]. Adaptive plasticity can itself be the subject of selection [11,12] and differ between populations [13].

When the existing populations encounter novel climates, the new selective pressures to which they will have to adapt will result in population divergence from their current structure [14]. Understanding and quantifying this divergence is thus a necessary step towards understanding the effects of future climate change on the adaptive evolution of present populations. One simple way to do this is to study the adaptive differentiation of existing populations that are exposed to a climatic gradient, i.e. space for time substitutions [15]. By replacing temporal variation with spatial variation, it is possible to study long-term processes, such as population evolutionary response to climate change, without relying on predictive simulations and selection experiments and when long-term data is not available [16-18]. However, populations can also differentiate in the absence of divergent selection, because of random evolutionary processes such as genetic drift [19]. To study the adaptive response of populations, it is thus necessary to separate population differentiation due to selective and random evolutionary processes. A common method to do this is to compare population differentiation based on quantitative genetic or phenotypic traits (Q) with that based on putatively neutral molecular markers (F) [20-24]. F is the neutral baseline against which population differentiation is measured. If Q = F, then the observed differentiation in phenotypic traits is solely due to random evolutionary processes, whereas if Q > F, then the phenotypic traits have likely differentiated under the effects of divergent selection, which were higher than random evolutionary processes.

In this study, we estimated population differentiation along a temperature and precipitation gradient of eleven populations of F. rubra for phenotypic traits and their plasticity, which are putatively under selection, and for molecular markers, which are supposedly not under selection. Our model species is a perennial grass, Festuca rubra, which has shown significant differentiation in traits and their plasticity on various spatial scales [25-27]. We used plants originating from locations of different temperature and precipitation from a unique natural grassland ‘climate grid’ established in western Norway [22,28]. We examined the relationship between the phenotypic, molecular and environmental data, with the goals to i) detect patterns of adaptive differentiation among the studied populations, ii) test whether these patterns are consistent with adaptation to climate change. If patterns of differentiation are detected, regardless of their causes, we further iii) tested the relationship between neutral molecular markers and phenotypic traits. Based on the space for time substitution approach, these results will provide input about the capacity of F. rubra in western Norway to sustain future climate change.

Materials and methods

Studied species

Festuca rubra L. is a common perennial grass species of temperate grasslands in Europe. F. rubra ssp. rubra used in the experiment is a widespread hexaploid type from the F. rubra complex. It reproduces sexually through outcrossing [29] as well as vegetatively by producing both intravaginal and extravaginal tillers on rhizomes. Clones of F. rubra possess considerable genetic variability and plasticity [25-27].

Study sites and sampling

The study material was sampled in a natural climate grid in Western Norway (the SeedClim project, [28,30]). The 12 natural grassland sites represented a full cross design of four levels of mean annual precipitation (noted 1–4 in increasing order, S1 Fig) and three levels of summer temperature (ALP—alpine, SUB—subalpine and BOR—boreal, in increasing order, S1 Fig). Sites were selected specifically to keep management and history, bedrock, slope, aspect and vegetation types as constant as possible [28,30]. The geographical distance between sites ranged from 0.65 to 175 km (average 15km, S1 Fig). The sampling areas were of variable size, depending on F. rubra density, plot morphology and vegetation homogeneity (S1 Table). Transects were established, and plants were collected until the required number of plants was reached. To avoid sampling clones of the same genotype, the minimum distance between two sampled plants was 1 m. In June 2014, 25 plants were collected at each locality and transported to the experimental garden of the Institute of Botany, Academy of Sciences of the Czech Republic in Průhonice, Czech Republic (49°59'38.972"N, 14°33'57.637"E; mean annual temperature 8.6°C and regular watering during the vegetation season) and planted into pots. The hexaploid status, which is the most widespread cytotype of F. rubra, was confirmed using flow cytometry using the protocol described in [31]. In locality ALP2, no F. rubra hexaploids occurred, so this site was excluded from further analyses. To ensure that each of the collected plants corresponded to a single genotype, plants were reduced to a single ramet at the end of June 2014. The maternal ramets were grown in a greenhouse for nine months to acclimate them to growth in pots and eliminate potential differences due to growth in different environments of origin.

Climatic chamber experiments

In March 2015, four single young ramets originating from the maternal ramets were taken from each plant, and each was placed in a different growth chamber. Each growth chamber simulated the spring and summer conditions in one of the four most extreme natural environments observed in the climate grid (ALP1, ALP4, BOR1, BOR4). The temperature in the growth chamber differed between the cold and warm treatments and changed over the growing season following the course of temperature at the natural localities (Table 1). Soil moisture was monitored continuously during the whole experiment and watering was modified to ensure constant moisture throughout the experiment. In brief, the dry regime plants were watered with about 20 mL of tap water per plant applied to the trays if the soil moisture was lower than 15%. In the wet regime, plants were cultivated under full soil saturation with about 1.5 cm of water in the bottom of the tray. For all the regimes, the same day length and radiation were used, i.e. 16 hours of full light (6 am– 10 pm) and 4 hours of full dark with a gradual change of light availability in the transition between the light and dark period over 2 hours. For an illustration of the experiment see S3 Fig. For a detailed description of the experimental protocol and the growth chamber conditions see [27].

Table 1

Climatic chamber settings throughout the study.

Reproduced with modification from [27].

	Alpine			Boreal
Day	Min (°C)	Max (°C)	Av (°C)	Min (°C)	Max (°C)	Av (°C)
1–4	5	15	9.8	5	16	10.1
5–25	3	12.5	7.5	3	16	9.2
26–46	3	12.5	7.5	3	18.5	10.2
46–67	3	12.5	7.5	3	24.3	12.5
68–88	3	14.5	8.4	3.4	25	12.9
89–176	3	14.7	8.5	5	23.8	14.8

The temperatures in the growth chambers were set to reproduce daily (minimal and maximal) and seasonal temperature variation during the spring in the alpine and boreal localities.

Traits measurement

At the end of August 2015 (21 weeks after placement in the growth chamber) we assessed photosynthetic activity, water potential and stomatal density and size as physiological traits. Photosynthetic activity was assessed through the measurements of the maximum quantum yield of primary PS II and the performance index for energy conservation from photons absorbed by PS II antenna [32]. These traits were measured using FluorPen FP-100 MAX/USB (Photon System Instruments, Brno, Czech Republic). Measurements were made in the dark, at 25°C, at the same time of the day for all plants (9.00–10.30pm). Two hours prior to the measurements plants were kept in the dark by switching off the growth chamber lights. Three measurements were made per plant, and their average was used in further analyses. Stomatal density and size were measured at the end of the experiment, prior to the harvest, using imprints that were observed with an Olympus BX53 microscope with Nikon DS-Fi2 camera for live projection and LIM NIS Elements Ar 4.11 software. Stomatal density was assessed in a 500μm square edge grid in three randomly chosen positions on one leaf. For each of the selected positions, the length of three randomly chosen stomata was measured and then averaged to obtain stomatal size. To measure water potential, all plants were saturated with water one day prior to their harvest. After harvest, about five leaves (or less depending on the volume they occupied) were placed in 1-mL Braun Omnifix-F syringes with a seal from synthetic rubber (B. Braun Melsungen AG, Melsungen, Germany) and frozen to break down the cell walls. After the material was unfrozen, cell cytoplasm was extracted by applying pressure on the plant sample with the syringe. Water potential measurement was made on 10μl of cytoplasm placed on 80 g/m2, wood pulp filter paper, using dew point micro voltmeter Wescor HR33 (Wescor Inc., Logan, USA).

At the end of the experiment, the number of intravaginal and extravaginal ramets and the plant height were assessed. The total number of ramets was calculated as the sum of intravaginal and extravaginal ramets. Finally, the plants were harvested, dried at 60°C to constant weight, and the aboveground (cut at the 3 cm height above the soil surface), and belowground (rhizome and root separately) biomasses were weighed. We calculated the below:aboveground biomass ratio.

We calculated the plasticity of all of the measured traits within each genotype as (|traitMAX−traitmin|)/traitMAX where traitMAX and traitMIN were respectively the maximal and minimal value of the trait measured on the same genotype across the four growth chambers [33]. We chose this index as it is easy to use, robust, widely applied approach that allows comparison among traits [34]. This estimate only characterized the maximal plastic capacity of an individual in variable environments without taking into account the direction of the plastic response or the change in intensity with environmental variation (i.e. reaction norms).

Microsatellite data

All of the studied genotypes were genotyped using standard DNA extraction and amplification methods for the purpose of a different study. For the genotyping we used extra ramets of the clones used in the growth chamber experiment. In short, total genomic DNA was extracted from dry leaves stored in silica gel using a silica-based DNeasy Plant Mini kit (QIAGEN, Dusseldorf, Germany). Quality and quantity of extracted DNA were measured by NanoDrop 2000 (Thermo Scientific, Waltham, USA), and all the samples were normalized to 20 ng/μl for subsequent PCR. Four microsatellite loci were amplified using Festuca and Festuca—Lolium complex specific primers HVM20 and HVM3 [35] and B4-D9 and B3-B8 [36]. Multiplex PCR reaction was designed combining all microsatellite primers to one reaction. Multiplex PCR reaction contained 2.5 μl QIAGEN Multiplex PCR Master Mix (Qiagen, Hilden, Germany), 0.25 μl of HVM 20 primers, 0.05 μl of HVM3, B3-B8, B4-D9 primers (10 μM each in initial volume) 0.7 μl H2O and 20 ng of DNA dissolved in 1 μl TE buffer. An initial denaturation step at 95°C for 15 min was followed by 40 cycles of denaturation (95°C for 30 s), annealing (60°C for 90 s) and extension (72°C for 60 s) steps, and a final extension step at 72°C for 10 min. A 1 μl aliquot of the PCR product was mixed with 11 μl of a 120: 1 solution of formamide and size standard GeneScan 500 LIZ (Life Technologies, Foster City, CA, USA). Fragment lengths were determined by capillary gel electrophoresis with an ABI 3130 Genetic Analyzer using Gene Mapper 4.0 (Life Technologies) and the peaks were scored manually. Overall, the four loci produced 62 alleles (20, 16, 15 and 11 for loci B3-B8, HVM2, B4-D9 and HVM3 respectively, S5 Table).

Data analysis

All analyses were done in R 3.2.5 [37]. The effect of initial ramet size (i.e. at the beginning of the growth chamber experiment) on phenotypic trait values and their plasticity was non-significant when tested in generalized linear models as we attempted to standardize it when setting up the experiment [27]. Therefore we did not include initial ramet size as an explanatory variable in our analyses.

Q-F and P-F comparisons

To estimate the population quantitative genetic differentiation (Q), we first estimated additive genetic variance and between population variance following a similar approach as in [38] using the animal model [39]. Within and between population quantitative genetic variation were estimated using REML approach in ASReml-R[40] as follows: y = Xβ + Zp + Za + ε where X and Z–incidence matrices assigning fixed and random effects to measurements in vector y, β –the vector of the fixed effects of growth chamber climate and any other unspecific chamber effects, p~N(0,Iσp)–the vector of random population effects, with σ –between population variance and I is identity matrix, a~N(0, 2Gσa)–the vector of random individual plant genetic effects, with σ—within population (additive genetic) variance and G–Loisselle’s kinship coefficient [41] matrix estimated in SPAGeDi 1.5 [42] using the four microsatellite markers of this study (S5 Table). Since trait plasticity was calculated as the difference of trait values measured for a single genotype grown in different growth chambers, within and between population quantitative genetic variation of trait plasticity was estimated with the model y = Xβ + Zp + Za + ε Narrow-sense Q was then calculated for trait values and plasticity as with and being the REML estimates of between and within-population genetic variance.

Kinship coefficient estimates based on four loci could be imprecise and downwardly bias the estimates of within-population variance, resulting in an overestimate of Q. We thus estimated the lower bound of Q through the population phenotypic differentiation index (P) as defined in [43]. Within and between population variation were estimated using a general mixed model as follows: y = β + b + ε, where β –the vector of the fixed effects of growth chamber and any other unspecific chamber effects, b~N(0,σ)–the vector of random population effects, with σ the variance of random population effects, and ε~N(0, σ)–the vector of residuals, with σ the residual population variance. We then calculated P for trait values and plasticity as with and being the estimates of within and between population variation, c the proportion of between population variance that is due to additive genetic factors only, and h trait heritability, or the proportion of within-population genetic variance that is due to additive genetic factors. We assume c = h = 1, meaning that all of the observed phenotypic variance is due to additive genetic factors only [44]. The assumption c = 1 plausible in our study, given that all populations are grown in a common environment, which reduces between population differences that can be due to environmental factors. Assuming h = 1 is clearly an overestimate of the additive genetic component of σ, but it is relevant in our case because we want to estimate the lower boundary of P, which is done by maximising the within-population variance.

To test for the effects of selection on population divergence in quantitative traits, we used the method described in [45] as implemented in an R script available from [46]. In brief, a neutral Q-F (P-F) distribution was simulated using the Lewontin-Krakauer distribution [47], and the estimates of F (from molecular markers), and of and for Q (of and for P). The quantile of the observed value compared against the neutral distribution was obtained in order to determine the p-value of the null hypothesis that Q (P) equals F. This method was particularly suitable for our data set, as it gives reliable results when used with relatively few neutral molecular markers, when population differentiation in molecular markers is low, and when the number of populations is relatively high (ten or higher). Given that microsatellite data can have mutation rates that are higher than migration rates, microsatellite based F can be downwardly biased. To avoid this, we also made Q—R comparisons, with R being an F analogue based on allele size, calculated following [48]. R should not be affected by the microsatellite mutation rate, provided that microsatellite size variance is proportional to their genetic distance [49,50].

Coinertia analysis

Coinertia analysis is a multivariate ordination method that measures the concordance between two data sets. The goal of a COA is to find a multidimensional projection of the two data sets which is a compromise between the maximal variance of each data set and the maximal covariance between the two data sets [51,52]. COA can be used to explore the shared structure between genetic diversity as estimated by molecular markers and phenotypic traits [53]. We used COA analyses to simultaneously examine the relationship between differentiation of molecular markers, phenotypic traits (or plasticity). If this relationship was significant, we further tested how the covariation between traits and markers was affected by environmental variation, instead of proceeding by pairwise comparisons between the three data sets. A significant effect of the environment would confirm that the observed joint population differentiation in traits and markers can be at least partly predicted by the environment, and thus the molecular markers could be used as indicators of the phenotypic value of the individuals of the studied populations.

COA was made using the projections of the first two axes of Principal Component Analysis (PCA) for trait values (or plasticity) and molecular markers. For phenotypic trait values PCA we used the Euclidian distance matrices based on the trait values averaged across growth chambers calculated by the R package ade4 [54], and for molecular markers PCA we used Nei-distances calculated by the R package adegenet [55]. The significance of the correlation between the matrix of phenotypic trait values (or plasticities) and that of neutral genetic patterns was assessed by 999 bootstraps. COA with trait means and individual genetic differentiation will be further referred to as COAmean, and COA with trait plasticity and individual genetic differentiation will be referred to as COAplast.

If significant coinertia was detected between trait values (or plasticities) and molecular markers, we then tested how are the individual projections on the main COA axis (COA axis 1) affected by climate variation. For this, we used ANOVA with temperature and precipitation of origin and their interaction as fixed, quantitative, explanatory variables. Significance levels were estimated using Fisher’s F statistic. A significant ANOVA result meant that the molecular marker variation which was associated with variation of trait value (or plasticity) was also influenced by climate variation. If such a relationship existed, we further tested for alleles that were associates with specific trait values, or which were characteristic for some populations (see below). Note that if allelic frequencies were influenced by phenotypic trait values and/or environmental variables, it does not mean that they were under selection, as microsatellite loci are by default considered to be neutral markers. This association could have been caused by pure genetic drift, or non-random mating patterns. Independent of the causes of the association, the alleles associated with phenotypic traits and climatic variables could further be used to identify individuals with specific phenotypes or thriving in specific climates.

Testing for the causes and consequences of population genetic differentiation

We tested for isolation by distance (IBD) and isolation by adaptation (IBA) by estimating the proportion of genetic differentiation between populations explained by geographical and environmental distances, respectively, using Mantel tests (R package vegan). The genetic distance matrix was based on Bruvo distances (R package polysat,[56]), and the geographic and environmental (combining temperature and precipitation data) distances were based on Euclidean distances. Significance testing was made with 999 bootstraps. In addition to testing the effect of geographical and environmental distances on their own, we used partial Mantel test to examine the effect of environmental distance after accounting for geographical structure, and of geographical distances after accounting for environmental structure.

We also tested the relationship between population genetic diversity and population trait means or plasticity using Pearson’s product moment correlations. Since the estimate of genetic diversity in hexaploid organisms is not straightforward because complete information about allelic frequencies cannot be obtained [57,58], several different estimators were used: number of alleles, number of effective alleles, allelic richness and expected heterozygosity and Pons and Petit’s index of population genetic diversity of non-ordered alleles relationship between individual genotypes [59]. All of these estimators were calculated with SPAGeDi 1.5 [42]. Being based on expected rather than observed allele frequencies, these estimators are informative about the population effective size. Thus a significant correlation between genetic diversity and trait means could indicate that populations with a lower effective size are phenotypically highly differentiated, possibly because of the effects of IBD or IBA [60].

Testing the association between molecular markers and traits and their plasticities in the studied populations

The statistical significance of the association between specific molecular marker alleles and phenotypic traits was tested with Multivariate Analysis of Variance (MANOVA). As the total number of alleles was higher than degrees of freedom associated with the phenotypic traits, we first selected the alleles that were most suitable for the analyses based on abundance criteria. We thus excluded rare or private alleles by selecting only alleles that were either present in seven or more populations, or alleles that were present in less than seven populations but they had population frequencies higher than 0.05. In locus HVM3 we excluded two alleles whose frequency was equal or close to 1 in all populations. Then we tested the effect of all alleles within each locus separately by using step-wise selection procedure within MANOVA. All alleles that were retained for each locus after the step-wise selection were combined together as explanatory variables in a new MANOVA. To make sure that the order of introduction of the alleles did not affect the MANOVA results, type 3 testing was carried out by permuting the introduction order of the alleles. To enable comparison between traits with different units, trait means were z-transformed using the scale function in R. Significance levels were assessed using the Pillai statistic. When a significant effect of at least one allele was observed, we performed individual ANOVAs for each trait separately in order to extract the coefficient of the slope associated with each of the tested alleles.

Results

Patterns of divergent selection

The global population genetic differentiation was very low, regardless of the estimator (F: 0.039 ± 0.011 standard error, R: 0.057 ± 0.047 standard error). Q estimates were obtained for ten out of the twelve phenotypic trait values. Q values for below:aboveground biomass and water potential could not be estimated with our model because their variance components were fixed at the boundary. Q estimates for trait values were generally low (0.072 ± 0.024). Six trait values showed Q-F significantly different from the neutral Q-F distribution (Table 2) and could be under divergent selection. Five of these traits were significantly different than the neutral Q-R distribution as well (Table 2). Q estimates for plasticity were obtained for nine phenotypic traits. Q estimates for plasticities were generally much higher than Q for trait values (0.477 ± 0.131), and close to 1 for stomatal density, water potential and number of ramets. However, the precision of the Q estimates for trait plasticity was lower than that for trait values (see standard errors in Table 2). This is likely because the sample size for the Q estimates of plasticity is four times smaller than that of Q for trait values. Both Q-F and Q-R were significantly different from the neutral distribution for plasticity of morphological traits, most of the plasticities of physiological traits, but not for the plasticities of the two resource acquisition traits (Table 2). Thus according to this result plasticity of morphological and physiological traits is likely under divergent selection. The lower bound of Q, estimated as population phenotypic differentiation, P, did not show these patterns of divergent selection (Table 2)–neither trait values nor plasticity showed patterns consistent with the effects of divergent selection except for the value of stomata density.

Table 2

Animal model estimates of population differentiation (Q and P) for phenotypic trait values and their plasticity.

	Trait value			Trait plasticity			Trait value				Trait plasticity
Trait	PST	std. error	PST—FST	PST	std. error	PST—FST	QST	std. error	QST—FST p-value	QST—RST p-value	QST	std. error	QST—FST p-value	QST—RST p-value
			p-value			p-value
Plant Height	0.045	0.021	0.308	0.027	0.02	0.776	0.097	0.056	0.007	0.067	0.592	0.758	<0.0001	<0.0001
Number of ramets	0.033	0.016	0.617	0.011	0.011	0.979	0.059	0.035	0.117	0.41	NA	NA	NA	NA
% extravaginal ramets	0.029	0.014	0.718	0.009	0.011	0.988	0.072	0.043	0.044	0.221	0.019	0.121	0.895	0.953
Aboveground biomass	0.033	0.017	0.598	0.036	0.024	0.534	0.132	0.074	0.0002	0.0101	0.644	0.766	<0.0001	<0.0001
Belowground biomass	0.032	0.016	0.642	0.025	0.016	0.836	0.132	0.082	0.0005	0.101	0.866	1.618	<0.0001	<0.0001
Rhizome biomass	0.020	0.011	0.893	<0.0001	<0.0001	0.997	0.048	0.037	0.2398	0.609	<0.0001	<0.0001	0.997	0.994
Below:aboveground biomass	<0.0001	0.001	0.997	0.018	0.016	0.914	NA	NA	NA	NA	NA	NA	NA	NA
FI.P0	0.009	0.006	0.982	0.007	0.009	0.993	0.014	0.028	0.96	0.974	0.033	0.46	0.67	0.850
Piabs	0.027	0.014	0.756	0.014	0.014	0.960	0.076	0.05	0.031	0.1777	0.335	0.69	<0.0001	<0.0001
Water potential	0.004	0.004	0.993	0.023	0.019	0.824	NA	NA	NA	NA	0.988	1.326	<0.0001	<0.0001
Stomata density	0.110	0.049	0.003	0.055	0.041	0.203	0.261	0.124	<0.0001	<0.0001	NA	NA	NA	NA
Stomata size	0.028	0.018	0.731	0.020	0.024	0.951	0.05	0.063	0.223	0.579	0.337	0.804	<0.0001	<0.0001

p-values were calculated using the Whitlock and Guillaume [45] method for P—F, Q—F and Q—R comparisons. Significant values are in bold. PIABS−photosynthetic index (see main text for explanations).

Effect of environmental variation on quantitative traits and molecular markers

COA results showed significant shared levels of inertia between trait values and molecular markers (COAmean), and between plasticity and molecular markers (COAplast; Table 3). The first axis of COAmean was significantly associated with temperature, precipitation and their interaction. The effect of temperature was due to the differentiation of ALP populations from SUB and BOR populations, whereas the effect of precipitation was not so clear. The interaction effect was mostly due to the differentiation of populations ALP4 and BOR1 whose coinertia patterns were similar to each other and contrasting to the remaining nine populations (Table 3, Figs 1A and S2). The individuals of populations ALP4 and BOR1 were characterised by short leaves, low below- and above-ground biomass; and high below:aboveground biomass ratio, rhizome biomass, proportion of extravaginal ramets, and stomatal density; as well as higher frequencies of the alleles B4-D9 245, B3-B8 281, HVM2 114 and HVM3 159. These qualitative associations between trait values and specific alleles were partly confirmed (and quantified) by the MANOVA results (S2 Table). The first axis of COAplast was affected by temperature and precipitation x temperature, but not precipitation alone (Table 3). Likewise for COAmean, the effect of temperature was due to the differentiation of ALP populations from SUB and BOR populations. The interaction effect was due to the high differentiation of population ALP4 from the remaining ten populations (Table 3, Fig 1B). ALP4 individuals had high plasticity of four growth traits (number of ramets, plant height, aboveground biomass and below ground biomass), and stomata size; low plasticity of resource acquisition and photosynthetic traits; and were associated with the alleles HVM3 117, B4-D9 249, 239 and B3-B8 281 (S2 Fig). These associations were not supported by the MANOVA results. Locality showed a significant effect on the individual projections of COAmean and COAplast (Table 3). The variation explained by locality itself was much higher than by temperature and precipitation combined, suggesting that other factors than climate might contribute to population differentiation.

Table 3

Results of coinertia analysis between phenotypic trait values and molecular markers (COAmean) or phenotypic plasticity and molecular markers (COAplast).

Analysis	R	p-value	Variable	Df	F	p-value	explained variation (%)
COAmean	0.207	0.001	Effect of environment
			Temperature	1	32.8	< 0.001	18.16
			Precipitation	1	5.546	0.02	3.07
			Temperature x Precipitation	1	45.28	< 0.001	25.07
			Effect of locality
			Locality	10	28.35	< 0.001	75.92
COAplast	0.193	0.001	Effect of environment
			Temperature	1	29.49	< 0.001	20.29
			Precipitation	1	1.343	0.249	0.92
			Temperature x Precipitation	1	9.446	0.003	6.50
			Effect of locality
			Locality	10	7.886	< 0.001	44.59

R–correlation coefficient of the coinertia analysis. p-value was obtained by 999 bootstraps. Significant values are shown in bold. The right part of the table shows the results of ANOVA testing the effect of climate variation (temperature and precipitation) or the effect of locality, on the projected values of the principal COA axis. DF–degrees of freedom, F–Fisher’s F. Significant values are in bold.

Fig 1

Boxplot of population differentiation along the first COA axis.

A. COA between phenotypic trait values and molecular markers B. COA with phenotypic plasticity and molecular markers. Localities with the same letters do not differ significantly according to Tukey’s post-hoc test.

Causes and consequences of population genetic differentiation

According to the Mantel tests, there was no significant correlation between genetic and geographic distances between populations (r = 0.010, p-value = 0.255), even in partial tests taking into account environmental distances (r = -0.005, p-value = 0.471). Environmental distances showed significant correlation with genetic distances when tested alone (r = 0.143, p-value ≤ 0.001) or in partial tests with geographic distances (r = 0.143, p-value ≤ 0.001).

All estimators of population genetic diversity (number of alleles, number of effective alleles, allelic richness and expected heterozygosity) showed similar estimates, thus only the correlations with Pons and Petit’s index of population genetic diversity of non ordered alleles are shown (Table 4). The correlation between population genetic diversity and trait means per population was significant for four traits. Proportion of extravaginal ramets and stomatal density were negatively correlated with genetic diversity, and below- and above-ground biomasses were positively correlated with genetic diversity. In line with this, these two pairs of traits show positive correlations within pair and negative correlations between pairs according to the COA (S2 Fig). Regarding mean trait plasticity, population genetic diversity was significantly and negatively correlated with the plasticity of number of ramets and aboveground biomass, and marginally and positively correlated with proportion of extravaginal ramets (Table 4).

Table 4

Results of the Pearson’s product moment correlation of mean population genetic diversity and trait values or phenotypic plasticity mean per population.

	Trait value		Plasticity
Trait	r	p-value	r	p-value
Plant height	0.324	0.330	-0.196	0.563
Number of ramets	-0.263	0.434	-0.628	0.039
% extravaginal ramets	-0.691	0.019	0.593	0.054
Rhizome biomass	-0.263	0.435	-0.348	0.294
Belowground biomass	0.627	0.039	-0.141	0.680
Aboveground biomass	0.828	0.002	-0.613	0.045
Below:aboveground biomass	-0.521	0.100	-0.570	0.067
PI.abs	0.304	0.364	-0.538	0.088
Water potential	0.599	0.052	-0.318	0.341
Chlorophyll fluorescence	0.381	0.248	-0.211	0.533
Stomatal density	-0.682	0.021	0.234	0.489
Stomatal size	0.098	0.774	-0.513	0.106

Discussion

Population differentiation and patterns of selection

Overall, the population genetic differentiation estimated from molecular markers was very low (both F and R). These results need to be taken with precaution, as population genetic differentiation in this study can be somewhat underestimated. First, the of lack of exact allele frequency estimates due to the hexaploid status of F. rubra [57,58] likely overestimates the rare allele frequencies and underestimates the abundant allele frequencies, resulting in an overall underestimate of population differentiation. Second, microsatellite markers are known to have a high mutation rate, which can further underestimate population differentiation [49]. We partly circumvented this second bias using an estimate of population differentiation that takes into account the marker mutation rate, R [48], under the assumption that the microsatellite markers of this study comply to a stepwise mutation model [61]. In spite of these potential biases, testing for isolation by distance in the studied populations using Mantel tests did not show significant structure either, which is consistent with low differentiation of F. rubra populations due to neutral genetic processes.

Population phenotypic differentiation (P) was overall very low, and did not differ significantly from F, except for stomatal density. Stomatal density is a trait regulates carbon uptake and water use efficiency of the plant, and as such is highly responsive to environmental variation [62]. Depending on the species, stomatal density response to environmental variation can be highly adaptive [63], and has often been used as an indicator for climate change based on plant fossil records [64,65]. It is thus not surprising, given the overall adaptive response of this trait across multiple species, that even with the conservative P−F test, stomatal density shows patterns of adaptive differentiation in our study. For all the other traits, based on P−F comparisons, if the selection is acting on these populations, its effects cannot be disentangled from the effects of genetic drift. This result is contradictory to the findings of previous studies using the same approach to estimate the effects of divergent selection [44,66,67] but see [68,69]. P estimates in the previous studies were made on traits measured in natural populations. Thus the observed variation between and within populations could include environmental effects, genotype x environment interactions, and non-additive genetic effects, which can inflate the estimated variances [70]. Since the populations in our study were grown in a common environment, most of the environmental effects on variance estimates were likely cancelled out, resulting in a decreased estimate of between population variance. The effects of the environment were further minimized by growing the maternal ramets in a common garden for nine months prior to the setting of the experiment. This was important as the conditions experienced by the maternal ramets were shown to have strong effects on offspring ramet performance in the same model system [71]. Within-population phenotypic variance, however, comprised non-additive genetic effects, which may increase its value (S3 Table), and result in an overall underestimate of P.

Within-population variance, as calculated for Q, was unsurprisingly lower than within-population variance as calculated for P. This is expected, as for the Q estimates, within-population variance was partitioned into additive variance, explained by the genetic relatedness between individuals, and residual variance, and only the additive variance was used to calculate Q. Thus Q estimates obtained with this approach were overall higher than their respective P. Furthermore, the majority of trait values and plasticities Q were significantly higher than F, consistent with the effects of divergent selection [72,73]. These results are nevertheless to be taken with precaution, as the Q estimates could be inflated. Indeed, the kinship estimates from only four loci could be imprecise, leading to increased “noise” underestimates of the within-population additive variance, σa2.

Based solely on the range of P and Q estimates, we cannot reliably conclude regarding the adaptation of the populations in this study to temperature and precipitation variation. However, the results of coinertia analyses show a significant effect of temperature and precipitation on the shared inertia. Furthermore, a previous study, equivalent to a ‘reciprocal climate common garden experiment’, using the same plant material, showed that the climate of origin has a small, but significant effect on phenotypic trait values [27]. The intensity of the adaptive response is rather low, as suggested by the weak effects of selection observed in these populations (estimated between -0.04 and 0.132, Stojanova et al. pers. obs.), which could explain the absence of clear cut patterns of adaptive population differentiation in our study. Taken altogether, these results cannot quantify the adaptive response of F. rubra populations in western Norway, but they can be considered as qualitative indicators of its existence.

Relationship between molecular markers, quantitative traits and environmental variation

The observed population differentiation in our study this seems to be mainly because populations from the coldest, alpine, climate are highly differentiated from populations from subalpine and boreal climate. Indeed, COA results show that alpine populations, and in particular the one from the locality with the highest precipitation (ALP4), have more effective foraging abilities (high proportion of extravaginal ramets and rhizome biomass) and slower growth (low number of ramets, leaf length, biomass measures). Increased foraging and reduced size is a pattern commonly observed in plants in response to different types of unfavourable environments [74-76]. In line with this, foraging and growth traits both show significant, but opposing correlations of their population means with population genetic diversity. This is probably due to the fact that alpine populations, notably ALP4, have reduced effective population size (indicated by their lower genetic diversity, S4 Table), while at the same time exhibiting adaptive trait combinations to their climate of origin. The causes of this association are discussed in the following section.

Interestingly, our results suggest that trait plasticity could also be under divergent selection. In our study, differentiation in plasticity is solely due to the higher plasticity of alpine populations according to the coinertia analyses, although alpine and boreal populations are facing the same absolute environmental distances (with opposing signs) in the growth chambers. This could be because alpine populations need to efficiently use the short ‘windows of opportunity’ when favourable conditions occur in the harsh alpine environment to rapidly increase their growth, and have thus evolved higher plasticity to cope with local, micro-environmental variations [27]. This hypothesis assumes adaptive phenotypic plasticity, i.e. that plastic individuals would have higher fitness than non-plastic ones in the studied environments [9,10]. If plasticity is non-adaptive for the studied populations, then the high phenotypic variability in alpine populations can also be due to the lack of canalisation [77]. In this case, the inability of the plants originating from extreme climatic environments to produce a stable phenotype when facing environmental variation results in non-adaptive plasticity. Non-adaptive phenotypic plasticity has indeed been observed in alpine plant populations (e.g. [70]). However, to evaluate the (non) adaptive character of plasticity in F. rubra, it is necessary to relate the degree of plasticity of each individual to an estimate of its fitness [9,10,78], which is a complicated task in clonal grasses. Indeed, the high longevity of clonal ramets (up to several hundred years, [79]) makes the lifetime fitness related to clonal reproduction nearly impossible to assess. In addition, our experimental plants flowered very rarely, and early flowering in long-lived perennials may not be an indication of high fitness, but a response to stress [80].

The plasticity Q estimates have rather low precision and non-significant lower bound estimates (P), and thus need to be considered with precaution. Only a few other studies have tested for among-population differentiation in phenotypic plasticity using Q-F comparisons and have found low or non-significant Q for plasticity [46,81-83]. However, the result found in our study is likely due to the reduced sample size available for the estimates, rather than the absence of a population differentiation mechanism. Indeed, it has been shown that plasticity is an evolving trait like any other quantitative trait [84], and can respond to selection imposed by controlled climatic variation [13].

COA results showed significant positive covariation between molecular markers and phenotypic trait values; and between molecular markers and phenotypic plasticity. The association between molecular markers and phenotypes is a general tendency in studies of the divergence of neutral traits and traits under selection [20,21,73]. It can be caused by different evolutionary mechanisms. First, it can be caused by a genetic linkage between molecular markers and loci under selection, resulting in genetic hitch-hiking of the molecular markers [85,86]. In our study, this is unlikely, given that microsatellite markers are putatively neutral, and have in general low genome coverage. Isolation by distance (IBD) and isolation by adaptation (IBA) can also result in simultaneous population differentiation in molecular markers and phenotypic traits. In the case of IBD, the differentiation cause is genetic drift, which acts on neutral markers as well as on loci under selection [20]. In the case of isolation by adaptation, the selective elimination of locally maladapted genotypes in the population causes a barrier to gene flow that will differentiate the whole genome, although the differentiation should be less pronounced in molecular markers than in traits under selection [73]. Under IBA, population genetic differentiation at molecular markers should be positively correlated with genetic differentiation at loci under selection, and both should be positively correlated to selection intensity. The COA and Mantel test results are consistent with this observation, suggesting that IBA might be responsible for the observed patterns of covariation between molecular markers and phenotypic traits. This conclusion should be considered with caution given that our results show only weak divergent selection among populations. Thus the observed patterns can also be partly due to IBD. However, no evidence of IBD was found for these populations based on Mantel tests for correlation between genetic and geographic distances. Furthermore, even low levels of population differentiation can contribute to variation in allelic frequencies. In coniferous tree species, for instance, neutral F below 5% can still generate strong biases in genetic association or environmental association analyses neglecting this structure [87].

Regardless of the causes of the observed associations, the information it provides can be used to identify individuals in the populations of this study bearing traits with specific values without going through the lengthy process of phenotyping. Since microsatellites are neutral molecular markers, a statistical association between adaptive traits and microsatellite markers is not informative about the genetic basis of the adaptive trait. However, establishing a relationship between molecular markers, phenotypic traits and environmental variation would help to easily identify genotypes and phenotypes in the studied populations that can be of interest for the future exploration of the response of the studied populations to climate change, by helping to select for appropriate genotypes (e.g. genotypes with high plasticity, or with specific phenotypic values) to test their response to climate variation. Although the observed marker-trait associations are only relevant for the populations of this study, the same methods can be used to identify analogous associations in other populations, and thus select for new genotypes that have comparable trait associations as those in our study.

Adaptive response of F. rubra to current and future climate change

The alpine populations of this study were overall highly differentiated from boreal and subalpine populations, and among the alpine populations, the population from the wettest habitat was the most genetically and phenotypically distinct. These results are somewhat consistent with an adaptive response to climate, although other environmental factors should not be overlooked. Indeed, we showed that locality explains between 30 and 50% more of the shared inertia than temperature and precipitation combined, so other factors, unaccounted for in this study also contribute and can be more important for shaping the observed patterns.

Finally, according to the climatic predictions for western Norway in the next century, the present populations of F. rubra will face an increasingly hotter and wetter environment [88], which based on the results of other studies on the same system is less selective than colder and drier environments [27]. Alpine populations could, therefore, benefit from the new, more favourable climatic conditions and increase their local abundance. Provided that warmer and wetter climate is favourable even beyond the climate values tested in this study, boreal populations could also benefit from climate change. Although it is also possible that F. rubra has an upper limit tolerance of temperature and moisture given that the species grows in considerably more temperate environments than Norway, it is unlikely that the boreal populations will face climate restrictions in the foreseeable future. Nevertheless, to reliably predict the response of boreal populations, we need data of their performance in temperature and precipitation conditions that are warmer and wetter than those tested in our study.

In conclusion, quantitative genetic differentiation between populations of F. rubra is consistent with the effect of divergent selection that is at least partly due to climate. The variation in adaptive traits and their plasticity is also correlated with variation in neutral molecular markers, and the covariation patterns are themselves influenced by variation in climate. Regardless of the causes of the statistical association between quantitative traits and molecular markers, the information it provides can be used to identify individuals bearing traits with specific values within our studied system without going through the lengthy process of phenotyping.

Map of the populations used in this study.

(PDF)

Click here for additional data file.

Results of coinertia analyses between trait values and molecular markers (A, B, C) and phenotypic plasticity and molecular markers (D, E, F).

A. and D. Population projection on the first two COA axes. The grey arrows represent the PCA projections of trait values (beginning of arrow) and molecular markers (end of arrow). B. and E. COA axes loadings for phenotypic variables. Black lines–growth related trait means, bold lines–resource acquisition trait means, dashed lines–physiological trait means.φP0– maximum quantum yield of primary PS II photochemistry, PIABS–performance index for energy conservation from photons absorbed by PS II antenna. C. and F. COA axes loadings for genotypic structure. For clarity of the representation, only 15 alleles that have the highest loadings for either axis are shown. Full black lines–locus HVM3, dotted black lines– B3-B8, full grey lines– B4-D9, dotted grey lines–HVM2.

(PDF)

Click here for additional data file.

Experimental design.

(PDF)

Click here for additional data file.

Characteristics of the sampled area for the natural populations of the study.

The dimensions are calculated as the distance between the most distant collecting points at two perpendicular directions.

(PDF)

Click here for additional data file.

A. MANOVA results of the multivariate analysis of variance of the effect of selected alleles on trait values and plasticities. B. Regression coefficients of individual ANOVAs testing the effect of allelic presence/absence on trait mean and on plasticity.

(PDF)

Click here for additional data file.

Estimates of within and between population variance for trait values and plasticities between population phenotypic variance (population variance),within population phenotypic variance (residual variance), between population additive genetic variance estimated from the animal model, within population additive genetic variance estimated from the animal model.

(PDF)

Click here for additional data file.

Genetic and genotypic diversity of the studied Festuca rubra populations.

(PDF)

Click here for additional data file.

Original genetic and phenotypic data used in this study.

(XLSX)

Click here for additional data file.