Linking ecological niche models and common garden experiments to predict phenotypic differentiation in stressful environments: Assessing the adaptive value of marginal populations in an alpine plant.

Environmental variation within a species' range can create contrasting selective pressures, leading to divergent selection and novel adaptations. The conservation value of populations inhabiting environmentally marginal areas remains in debate and is closely related to the adaptive potential in changing environments. Strong selection caused by stressful conditions may generate novel adaptations, conferring these populations distinct evolutionary potential and high conservation value under climate change. On the other hand, environmentally marginal populations may be genetically depauperate, with little potential for new adaptations to emerge. Here, we explored the use of ecological niche models (ENMs) linked with common garden experiments to predict and test for genetically determined phenotypic differentiation related to contrasting environmental conditions. To do so, we built an ENM for the alpine plant Silene ciliata in central Spain and conducted common garden experiments, assessing flowering phenology changes and differences in leaf cell resistance to extreme temperatures. The suitability patterns and response curves of the ENM led to the predictions that: (1) the environmentally marginal populations experiencing less snowpack and higher minimum temperatures would have delayed flowering to avoid risks of late-spring frosts and (2) those with higher minimum temperatures and greater potential evapotranspiration would show enhanced cell resistance to high temperatures to deal with physiological stress related to desiccation and heat. The common garden experiments revealed the expected genetically based phenotypic differentiation in flowering phenology. In contrast, they did not show the expected differentiation for cell resistance, but these latter experiments had high variance and hence lower statistical power. The results highlight ENMs as useful tools to identify contrasting putative selective pressures across species ranges. Linking ENMs with common garden experiments provides a theoretically justified and practical way to study adaptive processes, including insights regarding the conservation value of populations inhabiting environmentally marginal areas under ongoing climate change.

INTRODUCTION

Understanding key factors determining species distributions has long been of major interest in biology and now holds greater urgency given ongoing global environmental changes (Brown, 1984; Chen et al., 2011; Grinnell, 1917; Hutchinson, 1957; Parmesan, 2006; Root et al., 2003; Walther et al., 2002). One key concept, the center‐periphery hypothesis, posits that geographically disparate populations across a species’ range may have important differences, including variation in demographic, genetic, and phenotypic characteristics (Abeli et al., 2014; Gaston, 2009; Pironon et al., 2017; Sexton et al., 2009; Soule, 1973; Zakharov & Hellmann, 2008). Understanding the effects of such variation holds clear relevance under rapidly changing climate, land cover, and other environmental conditions because it may be crucial to facilitate adaptation and overcome the resulting threats (Abeli et al., 2014; Abeli & Orsenigo, 2018; Eckert et al., 2008; Gaston, 2009; Holt & Keitt, 2005; Kawecki, 2008; Pennington et al., 2021; Sexton et al., 2009; Swab et al., 2015). In this context, better predictions of future biodiversity patterns at and below the species level depend on two overlapping issues: (1) the relationship between geographic peripherality and ecological marginality, which directly influences patterns of genetic and phenotypic variation and (2) the potential adaptive value (or cost) of ecologically marginal populations under environmental change.

First, it is often assumed that environmental conditions are optimal in the center of a species’ range and suboptimal at the periphery (Brown, 1984; Hoffmann & Blows, 1994; Holt & Keitt, 2005), with geographically peripheral populations also being ecologically marginal. Accordingly, it has been proposed that peripheral populations exhibit lower individual fitness, have lower abundance, and experience declining demographic trends (e.g., fertility, mortality) (Angert, 2006; Brown, 1984; Hengeveld & Haeck, 1982; Herlihy & Eckert, 2005; Villellas et al., 2013). However, geographical and ecological gradients do not necessarily coincide (Hardie & Hutchings, 2010; Pironon et al., 2015; Soley‐Guardia et al., 2014; Soule, 1973), and the corresponding patterns regarding fitness, abundance, and demographic trends are not always consistent (Abeli et al., 2014; Hellmann et al., 2008; Holt, 2020; Osorio‐Olvera et al., 2019; Pironon et al., 2017). Hence, environmentally marginal areas can exist not only in the geographical periphery of a species’ range but also within it (Pironon et al., 2017). Defining environmental marginality is thus key for studying the ecological niches, adaptive processes, and conservation relevance of populations in such areas (Kawecki, 2008).

Second, opposing viewpoints have been presented regarding the genetic patterns and properties of geographically peripheral and/or environmentally marginal populations (Eckert et al., 2008; Kawecki, 2008). On one hand, lower survival and reproductive rates of such populations would lead to reduced effective population size and lower gene flow—and consequently, reductions in neutral and adaptive genetic diversity through genetic drift and inbreeding (Eckert et al., 2008; Kawecki, 2008). As a result, ecologically marginal populations may experience a significant genetic load due to maladaptation via low heterosis and an accumulation of deleterious alleles. On the other hand, strong selective pressures caused by environmentally stressful conditions may yield novel adaptations, conferring the population distinct evolutionary potential (Bontrager et al., 2021; Kawecki, 2008; Rolland et al., 2015). Indeed, the conservation value (or cost) of populations inhabiting environmentally marginal areas constitutes a long‐standing debate now of high importance given rapid global environmental changes (Abeli et al., 2014; Abeli & Orsenigo, 2018; Channell, 2004; Hunter & Hutchinson, 1994; Lesica & Allendorf, 1995; Millar & Libby, 1991) because they may represent a source of valuable adaptations. Therefore, identifying genetically based phenotypic variation related to contrasting environmental conditions across species ranges is a key step in assessing local adaptive processes and their potential to influence range shifts under changing climates (O’Neil et al., 2014; Razgour et al., 2019).

Correlative ecological niche models (ENMs) (or species distribution models; hereafter ENMs) have the potential to guide experiments aimed at testing for genetically based phenotypic differentiation. Correlative ENMs estimate habitat suitability based on associations between occurrence records of the species and predictor variables that characterize aspects of its environment (Elith & Leathwick, 2009; Guevara et al., 2018; Phillips & Dudík, 2008). Ecological niche models model environmental suitability via observed occurrence patterns, which may show mismatches with related population‐level phenomena such as abundance (Dallas & Hastings, 2018), productivity (Serra‐Diaz et al., 2013), and functional traits (Bohner & Diez, 2020; Pagel et al., 2020; Pironon et al., 2018). Hence, they do not necessarily represent an unbiased characterization of the species’ underlying niche, which defines the environmental requirements of its full life cycle (Bohner & Diez, 2020; Holt, 2009; Pagel & Schurr, 2012; Pagel et al., 2020; Pironon et al., 2018; Soberón, 2007).

Nevertheless, under clear assumptions ENMs can be used to study many important issues in ecology, evolution, biogeography, and global change (Anderson, 2013). These include the climatic drivers of species distributions (Guisan & Thuiller, 2005; Peterson et al., 2011), including adaptive processes and responses to climate change (Banta et al., 2012; Bush et al., 2016; Chardon et al., 2020; Csergő et al., 2017; Ikeda et al., 2017; Lee‐Yaw et al., 2016; Lira‐Noriega & Manthey, 2014; Marcer et al., 2016; Rolland et al., 2015; Swab et al., 2015). Such studies typically assume no functional genetic variation across the species’ range for factors related to its distributional ecology (Anderson, 2013; Hällfors et al., 2016; Valladares et al., 2014; but see Capblancq et al., 2020; Fitzpatrick & Keller, 2015). However, an adaptive response can occur if contrasting selective pressures exist across the species’ range and relevant genes hold enough variation (Blanquart et al., 2013; Kawecki & Ebert, 2004). In addition to indicating low‐suitability areas, ENMs can also identify the particular environmental conditions associated with them, or putative limiting factors. By doing so, they can aid in testing specific hypotheses regarding local selection pressures, including predictions of the genetically based phenotypic differentiation that would be expected (Dixon & Busch, 2017).

Proposed approach and study system

Here, we propose an integrated framework aimed at identifying putative environmental stressors using ENMs and testing for expected patterns of phenotypic differentiation via common garden experiments. We implement this framework using a temperate‐latitude alpine plant. We first use ENMs to characterize environmentally optimal and marginal areas and identify particular putative stressors for the latter. We then undertake common garden experiments to test for the genetically based phenotypic differentiation predicted under each of the stressors for populations inhabiting marginal areas. We consider as alpine those montane areas above treeline (Körner, 2003). As alpine ecosystems are locally highly variable (i.e., small differences in elevation can cause large changes in temperature, humidity, and exposure) (Körner, 2003), selective pressures can differ greatly over short distances (Herrera & Bazaga, 2008). Accordingly, in applying this framework, we take advantage of a well‐studied Mediterranean cushion plant, Silene ciliata, which is especially suited for evaluating evolutionary responses to ongoing climatic warming (e.g., García‐Fernández et al., 2012, 2013; Giménez‐Benavides et al., 2011b, 2018). We do so in an alpine ecosystem with a set of particularly challenging environmental conditions that can prompt a wide variety of phenotypic adaptations (Körner, 2003).

Substantial natural history and experimental information exists for S. ciliata, informing general putative selection pressures. Adult plants of S. ciliata face critical issues regarding annual phenological timing for flowering and seed production. Plants must avoid flowering too early after spring snow‐thaw due to the risk of losing their flowers to late‐spring frosts (Wheeler et al., 2014), yet they must flower quickly enough to produce fruits and set seeds before summer drought and early‐autumn frosts. Additionally, the species is dependent on water availability to complete the various stages of its life cycle, including flowering and fruiting as well as subsequent germination and establishment (Giménez‐Benavides et al., 2007a, 2018). Silene ciliata plants also confront particular environmental stressors during their first stages of development. They must survive the dry conditions typical of summer and early autumn as well as the occasional early autumn frosts. Furthermore, those plants go on to experience an extended cold period during the winter and late spring. The species also must deal with the extremely hot conditions and in situ heat‐related damage to plants in open alpine environments in the summer (Gauslaa, 1984; Körner, 2003; Körner & Larcher, 1988; Neuner et al., 2000). Thus, the leaf tissues of young plants must deal with physiological stress caused by drought and extreme temperatures (García‐Fernández et al., 2013; Giménez‐Benavides et al., 2007a, 2018).

Based on existing knowledge of environmental pressures acting in Mediterranean alpine areas, we identified an initial set of environmental variables. These best represent the extreme environmental conditions experienced by the species during its different life stages, show associations with its distributional patterns, and characterize intra‐range environmental variation (minimum temperature of coldest month, precipitation of driest and wettest months, mean snowpack in thaw months, and annual potential evapotranspiration, see Appendix 1 for further detail). Taking into account the species’ natural history (García‐Fernández et al., 2014; Giménez‐Benavides et al., 2007a, 2007b, 2018; Lara‐Romero et al., 2014), we then hypothesized that: (1) the populations experiencing contrasting values of these five environmental variables would exhibit different genetically determined phenological patterns and (2) populations experiencing contrasting values of these environmental variables would have leaf tissues in juvenile plants that exhibit genetically determined differentiation in cell resistance to extreme temperatures.

Guided by these general hypotheses for S. ciliata, we couple ENMs and common garden experiments to identify contrasting putative selection pressures across the species’ distribution and test for them. We posit that the ENM developed for this species will be a valuable tool to define environmental marginality across geography and identify specific stressors that may cause phenotypic differentiation of populations. Additionally, in interpreting the results in the context of other available information for the species, we make inferences regarding the adaptive potential of populations inhabiting environmentally marginal areas. Altogether, this approach could help future research aimed at disentangling species’ phenotypic differentiation and adaptive processes, as well as improve the accuracy of climate change forecasts of biodiversity.

MATERIALS AND METHODS

Research design

This research design links ENMs and common garden experiments following three major phases: pre‐modeling, modeling, and post‐modeling (Figure 1). We first formulated general hypotheses that identify potential environmental stressors likely to cause genetically based phenotypic differentiation between populations of the species (Step 1, Introduction). Subsequently, we selected the environmental variables to be used in ecological niche modeling according to these hypotheses (Step 2, Introduction). Next, we built an optimized ENM to predict habitat suitability across the landscape and delimitated environmentally marginal areas based on these values (Steps 3 and 4, Section 2.3). We then identified the environmental variables that held the strongest responses in the ENM (and thus would be most likely to exert local selective pressures) and assessed their response curves (Step 5, Section 2.3). Based on the hypotheses formulated in Step 1 and the modeled responses of the environmental variables identified via the ENM in Step 5 (in particular for the conditions present in the ecologically marginal areas), we made more specific predictions of expected genetically based phenotypic differentiation between populations regarding flowering phenology and cell resistance to extreme temperatures (Step 6, Section 3.2). We assessed these specific predictions via common garden experiments using populations from environmentally marginal and optimal areas to detect any genetically based phenotypic differentiation associated with relevant functional traits (Step 7, Section 2.4). Finally, we conducted statistical tests to determine if the results match the patterns expected (Step 8, Section 2.4). Taken together, these analyses tested whether populations occupying environmentally marginal and optimal environmental conditions showed the expected phenotypes linked to potentially adaptive genetic differences that could be of conservation value under climate change.

FIGURE 1

Workflow for linking an ecological niche model (ENM) and common garden experiments, illustrated for analyses of the alpine plant Silene ciliata in the Sistema Central mountain range of the Iberian Peninsula. The steps fall in three major phases: pre‐modeling (blue), modeling (yellow), and post‐modeling (light red). Step 1: Generate hypotheses of potential environmental stressors that may act as selective pressures for the species (such as snowpack, precipitation regimes, or temperatures), based on previous knowledge. Step 2: Select fine‐scale environmental variables that characterize the hypothesized stressors, to use as predictors in an ENM. Step 3: Build an optimized ENM aimed at characterizing the environmental gradients within the species’ distribution in this region, as well as modeled responses to the variables. Step 4: Define environmental marginality (optimal and marginal areas) based on the habitat suitability values estimated by the ENM. The map shows the habitat suitability values from the ENM in the Gredos subrange and indicates the location of the particular populations (in optimal and marginal environments) used in later steps of the study. Step 5: Identify the environmental variables with the strongest effects on the ENM, assessing their response curves to determine the values associated with marginal areas. The plot shows the response curve for snowpack of thaw months, indicating the position of optimal and marginal populations. Step 6: Formulate specific predictions of genetically based phenotypic differentiation between optimal and marginal populations that would be expected under selection pressures associated with the identified variables. The image shows the predicted phenotypic differentiation in flowering phenology related with differences in snow cover. Step 7: Undertake common garden experiments to assess specific predictions of genetically based phenotypic differentiation. The photo shows the common garden experiment assessing flowering phenology. Step 8: Undertake statistical tests to determine if the results match the expected patterns. The plot shows differences in flowering rate across time between optimal and marginal populations. Photos in Steps 1 and 7 by Javier Morente‐López

Study species

Silene ciliata Pourret (Caryophyllaceae) is a dwarf cushion, long‐lived perennial plant (Giménez‐Benavides et al., 2011) that is one of the dominant species of Mediterranean alpine habitats (i.e., those above tree line; Figure 1, Step 1). The study was undertaken in the Sistema Central of the Iberian Peninsula, a 500 km long mountain range that comprises the southernmost limit of the species’ distribution (see Kyrkou et al., 2015; Tutin et al., 1964) (Appendix 2A). More specifically, we studied populations in each of its three main mountain subranges: Béjar (BJR), Gredos (GRD), and Guadarrama (GDM) (Figure A2.1). Silene ciliata experiences specific environmental stressors of the Mediterranean alpine ecosystem over its life stages and responds to stark environmental gradients over small geographic distances. In this region, populations that grow at lower elevations (1850 m) start and end the flowering period (late July until late August) earlier than those occurring at higher elevations (2400 m) (Morente‐López et al., 2019). Seeds require exposure to cold conditions to break dormancy and germinate (i.e., cold stratification) (Baskin & Baskin, 1998; Giménez‐Benavides et al., 2005). Effective pollen and seed dispersal distances are low and relatively invariant across populations (Lara‐Romero et al., 2014, 2016).

Ecological niche modeling

Input data

We used species occurrence data in the Sistema Central obtained from the Global Biodiversity Information Facility (www.GBIF.org accessed October 2013), complemented with 10 new records from our fieldwork (with specimens deposited in the Real Jardín Botánico de Madrid herbarium; Appendix 2B). Silene ciliata is easily distinguishable from the other sympatric species of the same genus in the alpine areas of the Sistema Central (Silene boryi Boiss., which inhabits crag environments). Occurrence records with >10 m precision in geographic coordinates were removed. After cleaning the data, 120 total occurrence records were used for model building (Appendix 2C). To account for sampling biases (Boria et al., 2014; Merow et al., 2014) and associated artifactual spatial autocorrelation and environmental biases (Radosavljevic & Anderson, 2014), we used the R package spThin 0.1.0.1 to spatially thin our occurrence records by 0.75 km, yielding 87 records (Aiello‐Lammens et al., 2015; Appendix 3). A model of the species’ full range and inhabited conditions might produce different variable importance values and modeled response curves. However, the populations occurring in central Spain may differ in niche from those in other parts of the species’ range. Because of this, and to focus on identifying and characterizing fine‐grain environmental gradients within this area (i.e., those most likely to drive local adaptation), we selected one mountain range (Sistema Central) and characterized it as accurately as possible.

Our aim for modeling was explanatory rather than predictive, as we sought to identify putative driving variables and environmentally differentiated portions of the range rather than just characterize suitability (Araújo et al., 2019; Araújo & Guisan, 2006; Elith & Leathwick, 2009; Merow et al., 2014). Thus, we selected as model predictors the five environmental variables that we hypothesized could best explain this species’ distributional patterns and characterize intra‐range variation and extremes for environmental conditions (see Section 1 and Appendix 1). All variables had a 200 m spatial resolution and were obtained or derived from the climatic digital atlas of the Iberian Peninsula (Ninyerola et al., 2005, accessed October 2016). Mean snowpack in thaw months was calculated following the methodology proposed by López‐Moreno et al. (2007), reflecting the presence and potential magnitude of snow accumulation in a variable that ranges from 0 to 1. Environmental variables were constructed and managed using the R packages raster 2.8–19 (Hijmans & van Etten, 2012) and envirem 1.4 (Title & Bemmels, 2018). The pairwise correlation values of all variables were below 0.7, indicating relatively low collinearity (see Appendix 1 for detailed information on variable creation and selection).

Ecological niche model parameterization, evaluation, and selection

We used the algorithm Maxent 3.3.3k (Phillips et al., 2006), a presence‐background machine learning technique, to estimate the ecological niche of S. ciliata in the study region. We explored a range of settings affecting model complexity (combinations of feature classes and regularization multipliers) with the R package ENMeval 0.3.0 (Muscarella et al., 2014). To optimize model complexity, we first selected the candidate models with the lowest value of the Akaike Information Criterion corrected for small sample sizes (AICC), specifically those within 2 units of the best model (Warren & Seifert, 2011). We then evaluated these models via cross‐validation performance metrics using a spatial checkerboard partition method (“checkerboard2” from the ENMeval R package 0.3.0; aggregation factor 2) (Muscarella et al., 2014). Threshold‐dependent (omission rate calculated on withheld validation data after applying the 10 percent training omission threshold, or OR10) and threshold‐independent metrics (AUC calculated on withheld validation data, or AUCTEST in ENMeval 0.3.0) were used to evaluate the models in terms of overfitting and discrimination, respectively (see Anderson & Gonzalez, 2011). Details of model development and evaluation including background selection, model settings used (such as combinations of feature classes and regularization multipliers), data partitioning, and performance metrics can be found in Appendix 4 (Anderson & Gonzalez, 2011; Anderson & Raza, 2010; Araújo & Guisan, 2006; Galante et al., 2017; Guevara et al., 2018; Merow et al., 2014; Phillips & Dudík, 2008; Saupe et al., 2012).

Environmental variable importance: Defining environmental marginality, identifying putative selective pressures, and developing specific predictions

We visualized the optimal model in the full study area to indicate habitat suitability (Maxent's logistic transformation, bounded by 0 and 1; Phillips et al., 2006; Figure 1 and Figure A5.2), which we used to distinguish the degree of environmental marginality of different areas. We classified portions of the study area as “optimal” and “marginal” based on the predicted suitability values for the spatially thinned localities. The optimal category was assigned to areas with habitat suitability values at or above the highest 33rd percentile of the suitability value distribution of ranked pixels, whereas the marginal category was assigned to those at or below the lowest 33rd percentile. As we focused on contrasting environmental conditions, we did not consider intermediate suitability values further. We then plotted the occurrence records in environmental space defined by the variables with non‐zero coefficients in the model (see Appendix 5) to visualize environmental differences among the populations selected for addressing phenotypic differentiation (see below). We also identified which mountain subrange each of these occurrence records were associated with to determine if there was any environmental differentiation among subranges within the Sistema Central.

We then used the ENMs to develop specific predictions of the phenotypic patterns that would be expected under genetic differentiation of environmentally marginal populations. Based on the variable importance (Maxent's “permutation importance”) and parameter weights (i.e., lambda values) from the ENM, we identified and focused on the environmental variables that contributed most to explaining the species’ distribution (Galante et al., 2017; Phillips, 2005; Searcy & Shaffer, 2016); we did not consider further the variables that showed little importance for the model. We then plotted marginal response curves (model predictions for one variable when all other variables are set to their means) to determine the direction of the species’ modeled responses to the environmental variables, which potentially describe local selective pressures. This allowed us to pinpoint the particular environmental conditions predicted to be stressful for populations inhabiting environmentally marginal localities (in contrast to optimal ones). Harnessing this information and existing knowledge of the crucial life stages of the species, we developed specific predictions regarding genetically based phenotypic differentiation and outlined expected patterns for phenological and physiological traits, as described in the next section.

Genetically based phenotypic differentiation analysis

To assess associations of phenotype with habitat suitability values and environmental variables identified by the ENM, we selected functional traits with the potential to indicate phenological and physiological adaptations. First, we examined flowering phenological traits, which have been shown to be strongly affected by adaptation to changing environmental conditions (Bradshaw & Holzapfel, 2008; Chuine, 2010; Hoffmann & Sgró, 2011). Second, we characterized physiological traits associated with leaf cell resistance to cold and hot temperatures, representing the species’ capacity to withstand extreme environmental pressures (Nobel et al., 2002).

Common garden experiments and functional trait measurements

We then intersected the ENM results with a set of S. ciliata populations considered in an ongoing project addressing adaptive patterns associated with genetically based phenotypic differentiation (Morente‐López et al., 2020). Genetically based phenotypic differentiation can be detected through the variation in phenotypes found between different groups of individuals grown under the same controlled environmental conditions (since the differences that emerge are related to genetic composition of the groups). Specifically, we characterized the examined populations of S. ciliata as either “optimal” or “marginal” according to the habitat suitability values predicted by the ENM. In each of the three mountain subranges, we considered two populations: one located in an area predicted as “optimal” (high habitat suitability) and one predicted as “marginal” (low habitat suitability but still sufficient to support reproductive populations, see Section 2.3 for thresholding implementation) (Table 1; Figure 1 – Steps 3 and 4, Figure A5.2).

TABLE 1

Six Silene ciliata populations were chosen for phenotypic differentiation analyses under common garden conditions (Population and ID), corresponding to one environmentally marginal and one optimal population for each mountain subrange (Environment; Figure 1). Environmental classification in optimal and marginal areas was made based on habitat suitability (HS) values from the optimal Maxent ecological niche model (logistic transformation, see Methods). Geographic coordinates and elevation are provided for each population. Phenology (flowering onset, peak, end and duration) and leaf cell resistance to extreme temperatures (lethal temperature, LT50) were measured in each of the populations for given numbers of individuals and genets under common garden conditions (see Methods and Appendix 6)

Population	ID	Mountain subrange	Environment	HS values	Latitude	Longitude	Elevation (m)	Phenology trait analysis		LT50 traits analysis
								№ Individuals	№ Genets	№ Individuals	№ Genets
Canchal Negro	NEG	Béjar	Optimal	0.80	40°20'19.97''N	5°41'22.27''W	2360	54	31	6	6
Pico El Aguila	AGI	Béjar	Marginal	0.51	40°21'12.36''N	5°41'46.52''W	1950	79	46	6	6
Altos del Morezón	ZON	Gredos	Optimal	0.85	40°14'57.5''N	5°16'8.3''W	2380	57	29	6	6
El Sestil	SES	Gredos	Marginal	0.43	40°16'24.45''N	5°14'54.93''W	1900	103	33	6	6
Pico de Peñalara	PEÑ	Guadarrama	Optimal	0.87	40°51'2.11''N	3°57'24.02''W	2400	65	35	6	6
Najarra baja	NAJ	Guadarrama	Marginal	0.29	40°49'23.46''N	3°49'52.53''W	1850	74	43	6	6

Flowering phenology

In the late summer of 2013, we collected a minimum of 30 randomly selected adult plants from each of the six populations. From these, we obtained rooted cuttings to generate a sufficient number of genetically identical individuals to overcome potential losses during rooting and growth. The entire collection was allowed to grow in the Universidad Rey Juan Carlos CULTIVE laboratory (690 m, https://urjc‐cultive.webnode.es) under common garden conditions for 7 months to minimize carry‐over effects from the original environment (Appendix 6A). The diameter of all selected individuals of each population (Table 1) was measured on the date when the first plant started flowering (15 April 2014). From then until 25 September 2014, when the last new flower appeared, we counted the total number of new flowers of each plant weekly (n = 23). From these data, we calculated four individual‐based flowering traits as representations of the temporal patterns of reproductive investment: flowering onset, peak, end, and duration (details in Appendix 6A).

Leaf cell resistance to extreme temperatures

The purpose of this study was not to obtain accurate measures of absolute temperature tolerances of individual plants in their natural habitats, but rather to assess whether relative differences existed in assays of temperature tolerance performed under shared artificial conditions, between plants originating from populations in optimal versus marginal environments. With this goal, and recognizing that laboratory conditions will not match the complex environmental sequences experienced in the field, we designed a simple experiment to evaluate possible patterns of genetically based differentiation with all plants exposed to one common state. To do so, we collected seeds from 30 randomly selected plants from each of the six populations. Seeds were germinated and seedlings grown for 8 months (juvenile plants for this species) under common garden conditions in the CULTIVE laboratory. Although leaves may change physiologically over the seasons, the simplified experiment conducted here used two treatments that were subjected to all plants at a single growth stage. Tolerance of juvenile plant leaf cells to extreme temperatures was determined in May–June 2014 by measuring electrolyte leakage, an indicator of cell membrane integrity (Didden‐Zopfy & Nobel, 1982; Drennan, 2009; Soto‐Correa et al., 2013). The lethal temperature, defined as the temperature at which half of the maximum electrolyte leakage was reached, was estimated for each plant separately for both low (hereafter LT50LOW) and high (hereafter LT50HIGH) incubation temperature sequences (see Nobel et al., 2002 detailed description). LT50 values have been found to be good predictors of extreme temperature tolerance in the field (Drennan, 2009; Nobel et al., 2002). Each low (5 to −20°C) and high (30 to 55°C) sequence was composed of six temperatures (5°C intervals), for which the initial temperature acts as a baseline of the electrolytic conductivity for undamaged tissue. The temperatures selected can be reached in this ecosystem under natural conditions (see Section 1 and Appendix 6). Six plants per population were randomly selected for the analysis, and six 2‐mm leaf discs were taken from each plant for testing at multiple temperatures (see Appendix 6 for detailed information).

Data analysis

Categorical approach

We used two different statistical approaches to analyze whether the optimal and marginal environmental categories identified through ENMs were associated with differences among populations in the functional traits. For flowering phenology, we used mixed effects Cox proportional hazards regression models implemented in the coxme function from the R package coxme 2.2–7 (Therneau, 2018). We modeled rate of flowering over time with plant size as a covariate, environmental class as a fixed factor with two levels (optimal/marginal), and population of origin as a random factor. In the case of leaf cell resistance to extreme temperatures, we tested the association between physiological traits (LT50HIGH and LT50LOW) and environmental classes using linear mixed models implemented in the lmer function from the R package lme4 1.1–17 (Bates et al., 2014) with the same model design explained above for phenological traits (Appendix 7).

Continuous approach

We used linear models (R Core Team, 2013) to assess whether habitat suitability values predict genetically based phenotypic differentiation. We assessed the association of: (1) flowering phenology (onset, peak, duration, and end of flowering) and (2) cell resistance to extreme temperatures (LT50HIGH and LT50LOW) with continuous habitat suitability values and plant size (Appendix 7).

RESULTS

Ecological niche model selection, evaluation, and spatial projection

Four ENMs were identified as co‐optimal. They all had delta AICc values <2, similar values for other metrics of performance (averages of test omission rates and AUC, each calculated on data withheld in cross‐validations), and the same feature class combination (Table A8.1). For further analysis, we selected the model with the lowest AICc value, which was also one of the simplest overall (see model 1 in Table A8.1). When we projected this ENM to the full background spatial extent, the highest habitat suitability values were concentrated along mountain crests around 2400 m, and lower elevations corresponded to decreasing habitat suitability values that were <0.1 below 1800 m (Figure 1 Step 3; Figure A5.2).

Environmental variables and specific predictions

The modeled responses for the predictor variables differed in weights and complexity. Snowpack constituted the most important variable, followed by minimum temperature of the coldest month, potential evapotranspiration, and precipitation of the driest month (Table 2). Marginal response curves showed the shapes of the relationships between each of the environmental variables and habitat suitability (logistic transformation, Figure 2; raw outputs in Appendix 8B, Figure A8.1). Minimum temperature of the coldest month and snowpack showed the strongest contrast between optimal and marginal populations, making these variables the most likely to distinguish populations experiencing contrasting environmental conditions (Figure 2). As minimum temperature of the coldest month increased, habitat suitability decreased, dropping rapidly around −8°C and approaching zero at around 3°C (Figure 2a). Habitat suitability was close to zero until snowpack values were positive (i.e., presence of snow cover), where habitat suitability spiked suddenly and reached a maximum at snowpack values of 1 (Figure 2b). Potential evapotranspiration and precipitation of the driest month showed responses similar to each other (Figure 2c and d), with gradually decreasing habitat suitability as each variable increased. However, whereas potential evapotranspiration showed a consistent trend regarding the habitat suitability of optimal vs. marginal populations, precipitation of the driest month showed no clear contrast between them (Figure 2d). Tellingly, populations selected as “optimal” and “marginal” fell in different subsets of the environmental space defined by the two most important variables, snowpack and minimum temperature of coldest month (Appendix 5, Figure A5.1a). No major environmental differentiation was found among mountain subranges (Appendix 5, Figure A5.1b).

TABLE 2

The five input environmental variables used to develop the final Maxent ecological niche model (ENM) for Silene ciliata in the Sistema Central mountain range of the Iberian Peninsula. These variables reflect factors considered likely to affect habitat suitability based on prior ecological knowledge of the species. These include snowpack (mean snow accumulation ratio in the thaw months); minimum temperature (annual minimum temperature of the coldest month); annual potential evapotranspiration (PET); precipitation driest (precipitation of the driest month); and precipitation wettest (precipitation of the wettest month). These environmental variables were then ranked based on their permutation importance obtained from the ENM. Feature classes used and lambda values (coefficient weights) from the ENM are provided for each variable (see Section 2 and Appendix 4)

Variable	Permutation importance	Feature class	Lambda value coefficient
Snowpack	76.89	Linear	0.88
		Hinge	−2.48
Minimum temperature	10.32	Linear	−4.52
PET	7.19	Quadratic	−3.27
Precipitation driest	5.60	Quadratic	−2.84
Precipitation wettest	0	Linear	0

FIGURE 2

Response curves for four environmental variables from the Maxent ecological niche model (ENM) for Silene ciliata in the Sistema Central mountain range of the Iberian Peninsula: (a) minimum temperature of coldest month, (b) snowpack in thaw months, (c) annual potential evapotranspiration, and (d) precipitation of driest month were the four variables identified based on the permutation importance of the final ENM (Table 2). These marginal response curves show the shape of the relationship between the environmental variable and habitat suitability values ranging from 0.0 (lowest) to 1.0 (highest; Maxent's logistic transformation, see Section 2 and Appendix 8B) when all other variables are set to their means. This allows identification of the environmental conditions predicted to be stressful for populations inhabiting environmentally marginal localities. Based on this and existing knowledge regarding crucial life stages of the species, we generated specific predictions of expected genetically based phenotypic differentiation of phenological and physiological traits (see Section 2 and Figure 3). Vertical bars represent values for each of the six populations selected for the common garden experiments (Figure 1 and Table 1; for environmental and habitat suitability values associated with each population, see Appendix 11, Table A11.1). Note that: As minimum temperature of the coldest month increases (a), habitat suitability values decrease. Habitat suitability values are close to zero until snowpack in thaw months is positive (b; presence of snow cover). Precipitation of the driest month (c) and potential evapotranspiration (d) showed similar responses, with gradually decreasing habitat suitability values as each variable increased. Minimum temperature of the coldest month (a) and mean snowpack in the thaw months (b) showed the greatest contrast in suitability between optimal and marginal populations. Thus, these two variables were crucial to generate specific expectations regarding phenotypic differentiation

FIGURE 3

Differences in flowering phenology in common garden experiments between populations of Silene ciliata from the Sistema Central of the Iberian Peninsula from environmentally marginal and optimal conditions (categorical approach). Environmental classification into optimal and marginal areas was based on the habitat suitability values of the ecological niche model (Table 1; Figure A5.2). Samples from populations located in marginal areas displayed significantly later flowering onset (a), peak (b), and end (c) than those from populations inhabiting optimal areas; furthermore, flowering season duration (d) was significantly longer for samples from populations in marginal areas (*p < .05, ***p < .001). These results matched predictions made for all four phenological traits measured, based on the modeled response curves and existing knowledge of the species (Figure 2)

The response curves led to specific predictions for functional traits related to flowering phenology and leaf cell resistance to extreme temperatures, refinements of the original general hypotheses (see Section 2.2). Considering the environmental variables with high importance for the model, response curves indicated that lack of snowpack, higher minimum temperature of the coldest month, and higher potential evapotranspiration led to lower suitability. Translation of these modeled patterns into specific predictions for genetically based phenotypic differentiation produced the following expected differences between populations occurring in contrasting environments:

Flowering phenology

Environmentally marginal areas defined by the model have higher temperatures of the coldest month (which prompt a more active metabolism) and lack of snowpack (which would have provided effective frost protection). This combination can result in higher mortality risk for flowers when late‐spring frosts occur (Hennon et al., 2012; Kreyling, 2010; Schaberg et al., 2008). Therefore, we predicted that populations in such areas experience selective pressure to delay flowering. If populations have evolved according to this pressure, we would expect them to have a later flowering onset, peak, and end, as well as a longer flowering duration.

Leaf cell resistance to extreme temperatures

Environmentally marginal and optimal areas identified by the model differed substantially in temperature and potential evapotranspiration (which is highly correlated with maximum temperature of the warmest month, see Appendix 1). Marginal areas correspond to warmer conditions and higher potential evapotranspiration values (which prompt physiological cell stress to high temperatures and desiccation) (Zandalinas et al., 2018). In contrast, optimal ones are consistently colder and with lower values of PET (with modeled suitability that remains high even for the coldest temperatures within the study region). High temperatures and higher values of PET, especially in late summer and early autumn when water availability is very low, can lead to high physiological stress related to desiccation and greater mortality risk for plants (García‐Fernández et al., 2013). Therefore, we predicted that populations in environmentally marginal areas would experience selective pressure for leaf cells that are resistant to higher temperatures—but that populations in optimal areas would not experience analogous pressure for greater resistance to extreme cold. If populations have evolved according to this pressure, we would expect marginal populations to have greater values of LT50HIGH but that the optimal populations would not have lower values of LT50LOW (see Methods 2.5).

Functional traits

Results of both the categorical approach (environmental classification into optimal and marginal areas) and the continuous one (directly using habitat suitability values) matched the predictions for flowering phenology (Section 3.2). The effect of environmental classification based on the ENM was significant in the expected direction for all phenological variables (Figure 3; Appendix 9, Table A9.1). Samples from populations located in environmentally marginal areas displayed a later flowering onset, peak, and end than those inhabiting optimal areas (Figure 3a–c; Table A9.1). Additionally, flowering season duration was longer in samples from populations in marginal areas (Figure 3d; Table A9.1). Similarly, a significant negative relationship existed between habitat suitability values and each of the four flowering phenology variables (Table A9.1): as suitability values increased, the onset, peak, and end of flowering were all earlier, and the duration of flowering decreased (Figure 4). Parameter estimates of all developed models are in Appendix 10 (Tables A10.1 and A10.2).

FIGURE 4

Phenological traits related to flowering based on common garden experiments for Silene ciliata in the Sistema Central of the Iberian Peninsula, in relation to habitat suitability values from the optimal ecological niche model (continuous approach). Habitat suitability had a significant negative relationship with each of the four phenology variables (***p < .001). As habitat suitability decreased, the onset (a), peak (b), and end (c) of flowering were all later, and the duration (d) of flowering increased. These results matched predictions made for all phenological traits measured, based on the modeled response curves and existing knowledge of the species (Figure 2)

Results for both the environmental classification into optimal and marginal areas (categorical approach) and the habitat suitability values (continuous approach) matched predictions for LT50LOW but not for LT50HIGH (Section 3.2). As expected, no significant associations were found between environmental classification or habitat suitability values and cell resistance to extreme low temperatures (LT50LOW); however, contrary to expectations for cell resistance to extreme high temperatures, no significant associations existed for either characterization of suitability (LT50HIGH, Figure 5; Appendix 9, Table A9.2). Parameter estimates of all developed models are in Appendix 10 (Tables A10.3 and A10.4).

FIGURE 5

Traits related to leaf cell resistance to extreme high and low temperatures (LT50HIGH and LT50LOW, see Section 2.5) based on common garden experiments for Silene ciliata in the Sistema Central of the Iberian Peninsula. Results are shown in relation to habitat suitability values from the final ecological niche model (a and c; continuous approach), and from environmentally marginal and optimal conditions (b and d; categorical approach). Environmental classification into optimal and marginal areas was based on the habitat suitability values of the ecological niche model (Table 1, Figure A5.2). Neither habitat suitability values nor environmental classification had any significant relationship with LT50HIGH or LT50LOW (p > .005)

DISCUSSION

Overall findings and utility of this framework

The ENM characterized specific putative selective pressures in environmentally marginal areas, and common garden experiments confirmed most of the expected differences in phenotypic traits. Specifically, the weights and responses of the model allowed refining the original hypotheses and developing specific predictions of expected phenotypic differentiation between marginal and optimal populations in flowering phenology and leaf cell resistance to extreme temperatures. The common garden experiments confirmed the expectation that marginal populations delay and elongate the flowering period, as predicted in response to selective pressure to avoid the negative effects of late‐spring frost. Experiments regarding leaf cell tolerance to extreme cold and hot temperatures matched the expectation for cell resistance for optimal populations but not for the marginal ones, based on predicted selective pressures to avoid the negative effects of high temperatures and desiccation. In relevant sections below, we discuss caveats to these interpretations (including variation and statistical power), many of which point toward fruitful directions for future research.

Taken together, the results support the proposal that at least in some cases ENMs can be useful for identifying areas within a species’ range that are environmentally marginal, as well as for developing refined predictions regarding expected functional genetic differentiation in such areas. Hence, ENMs can guide experiments with the potential for identifying adaptive population divergence driven by divergent selective pressures in marginal areas (Hällfors et al., 2016). Extrapolating the results found in this study to other areas or species should be undertaken with caution, due to the species‐ and context‐dependency of functional differentiation (Bohner & Diez, 2020; Pagel et al., 2020; Pironon et al., 2018; Thuiller et al., 2010; Wittmann et al., 2016). Indeed, this framework aims to guide the identification of mismatches between suitability and underlying genetics and/or functional performance. Hence, it adds to a growing body of research addressing interesting—and important—mismatches between suitability and factors such as survival, reproductive success, population abundance and growth rates, and both adaptive and neutral genetic patterns (Bohner & Diez, 2020; Dallas & Hastings, 2018; Pagel et al., 2020; Pironon et al., 2018; Serra‐Diaz et al., 2013). The methodological framework proposed here could be applicable across species and study systems that vary widely in the amount of previous knowledge available. When little is known, ENMs can be useful for identifying novel hypotheses for exploration. In contrast, if previous knowledge exists and general hypotheses can be created a priori as in this study, the models can be useful for selecting which hypothesized stressors are most likely and for refining their corresponding predictions.

Interpretation of responses to putative selective pressures

In this study, the key factor leading to low modeled habitat suitability (i.e., the strongest putative driver of selective pressure) was low snowpack in thaw months. Snowpack in thaw months is driven by spring temperature and precipitation (López‐Moreno et al., 2007) and is related to the timing and success of crucial life stages such as plant flowering in S. ciliata (Giménez‐Benavides et al., 2007a, 2011) and other temperate alpine plant species (Dunne et al., 2003; Giménez‐Benavides et al., 2018; Inouye, 2008; Inouye et al., 2002; Lara‐Romero et al., 2014). Moreover, strong evidence has been found for the role of snowpack in local adaptative processes, particularly in alpine plants (Anderson & Wadgymar, 2020; Inouye, 2020). Snowpack thaw directly affects vulnerability to spring freezing of various plant tissues and the growth patterns of alpine plant species (Jonas et al., 2008; Wheeler et al., 2014). In this context, the late‐flowering strategy (seed‐risk sensu Molau, 1993) followed in general by S. ciliata would be especially prevalent in populations inhabiting areas with shorter snow cover duration and, consequently, higher risk of late‐spring frosts (Wheeler et al., 2014). This result highlights the ability of ENMs to identify one or more particular environmental factors strongly associated with great intra‐range differences in habitat suitability, even for species inhabiting landscapes with high local environmental heterogeneity. Notably, snowpack is greatly affected by microtopographic variation, especially in alpine areas where wind scouring occurs on exposed ridges (DeBeer & Pomeroy, 2009); this could be a source of uncertainty in estimating habitat suitability values, especially at fine grains. In any case, marked ecological niche differences related with fine‐grain environmental variables also have been found in other alpine species between populations inhabiting central and marginal areas (Giménez‐Benavides et al., 2007b; Papuga et al., 2018).

Indeed, genetic differentiation was found in S. ciliata for phenological traits related to flowering, as expected for adaptation to avoid late‐spring frosts in areas with low snowpack and in concordance with earlier studies. When grown in common garden settings, samples from populations in marginal areas exhibited later flowering start and peak and flowered for a longer time than those from populations inhabiting optimal ones. It is important to consider the possibility that the phenological variation found could be due at least in part to effects of the environments of origin because plants used in the common garden experiment derive from cuttings obtained from individuals that grew in natural conditions. However, the patterns found are congruent with other analyses of selection, based on the effects of phenological trait differences in fitness for the same species at the same sites, which showed that flowering onset and duration are under selection in populations at different elevations (Giménez‐Benavides et al., 2011b, see also Morente‐López et al., 2019). The patterns of phenological differentiation found are also congruent with results obtained in subsequent experiments implemented with plants grown from seeds derived from the cuttings used in the present study (unpublished results). Moreover, a parallel transcriptomic study with the same populations identified SNPs (single nucleotide polymorphisms) with unusually high allele‐frequency differentiation between marginal and optimal environments, for candidate genes in the flower development pathway (Sacristán‐Bajo et al., 2019). In addition, the phenotypic differentiation found in this study adds to recent evidence of the adaptive potential of marginal populations at early life stages (i.e., seedlings). In a reciprocal field seed‐sowing experiment involving populations used in this study (Morente‐López et al., 2020), seeds from marginal populations showed higher seedling survival rate when sown in both optimal and marginal areas. Additionally, in that study gene flow between marginal populations (by artificial pollination) improved seedling survival, but gene flow from optimal to marginal populations did not.

Although common garden experiments such as this one only can demonstrate the existence of genetically based differentiation (not adaptation per se), the information from other studies mentioned above provides additional evidence toward the conclusion that the phenological differentiation found here likely is the result of divergent selection (Giménez‐Benavides et al., 2011; Morente‐López et al., 2019; Sacristán‐Bajo et al., 2019). Alternatively, such differentiation could result from neutral processes unrelated to selection, such as genetic drift (Garcia‐Ramos & Kirkpatrick, 1997). Therefore, future research to corroborate the adaptive nature of the differentiation found here could involve the use of genomic approaches to quantify the magnitude of genetic differentiation and the nature of the genes involved, as well as reciprocal in situ common garden experiments with fitness measurements.

Similar patterns also have been observed in co‐occurring species in the same ecosystems (Lara‐Romero et al., 2014), suggesting that phenological reproductive traits are likely to be under strong selection. In fact, they have been shown to be critical for plant fitness (Fox, 2003; Levin, 2006), capable of evolving rapidly under climatic variation (Franks et al., 2007), and influencing species distributions and potential responses to climate change (Banta et al., 2012; Chuine, 2010). This also agrees with studies showing that these traits generally show high evolvability (Kawakami et al., 2011; Méndez‐Vigo et al., 2013; Volis, 2011) and with a meta‐analysis indicating consistent selection patterns (Munguía‐Rosas et al., 2011).

Overall, the environmentally marginal populations of S. ciliata that show genetically determined phenotypes expected to perform better with low snowpack may be advantageous under climate change. In the Iberian Peninsula, temperature is expected to increase, and precipitation decrease, causing a reduction in snow cover in the coming decades. Nevertheless, as discussed above, extrapolating the patterns found in this study to other parts of the distribution of S. ciliata or to other species should be undertaken only with caution, because functional differentiation of plant populations related with abiotic factors is commonly species‐ and context‐dependent (Leuschner et al., 2009; Papuga et al., 2018; Wagner et al., 2011).

In contrast to the differentiation found for flowering phenology, the lack of expected association between low habitat suitability values and functional differentiation in the cell resistance of leaves may reflect the absence of strong selective pressures concerning high temperatures per se, but other non‐exclusive possibilities also exist. Notably, the importance of temperature and PET to the model (which led to the expectations regarding cell stress to extreme temperatures and desiccation) was much lower than that of snowpack (which more directly led to the expectations regarding flowering phenology). Previous research in the same ecosystem demonstrated no differences in leaf functional traits linked to water‐use efficiency (e.g., carbon and nitrogen isotope ratios of organic material) or economic spectrum (e.g., specific leaf area and leaf dry matter content), along an elevational gradient for S. ciliata and co‐occurring species (Pescador et al., 2015). Other studies in both plants and animals have shown that phenological traits are more sensitive to divergent selection than physiological ones (Bradshaw & Holzapfel, 2008; Hoffmann & Sgró, 2011). For instance, evolved differences in thermal resistance in organisms are usually weak or absent among populations occupying even a broad range of climatic conditions, which may be related to a lack of genetic variation for these traits at the population level (Bradshaw & Holzapfel, 2008; Hoffmann & Sgró, 2011, but see Sinclair et al., 2012). A previous study performing a drought experiment using the same species, although limited to one population per habitat, observed that plants derived from marginal populations had greater survival than those from optimal ones (García‐Fernández et al., 2013). This fact highlights that decoupling the effects of drought and high temperature on plants is very difficult, particularly for seedlings and juvenile plants. Additionally, the statistical power of the analyses developed in the present study regarding temperature tolerances was hampered by the limited sample size (number of plants studied per population) and the unexpectedly high variance of the results. To reach more conclusive results, additional experiments would be required with a larger sample size. In addition, it would be relevant to test resistance to extreme temperatures at different developmental stages of the plant (e.g., seedlings, mature plants). Indeed, for traits more related to drought resistance physiological differences may be more likely to arise in earlier crucial life stages (such as germination or seedling establishment), (García‐Fernández et al., 2013, 2014; Giménez‐Benavides et al., 2007a, 2008), whereas other traits may differ at later developmental stages.

CONCLUSIONS AND FUTURE DIRECTIONS

Using ecological niche modeling in conjunction with laboratory and field experiments represents a fruitful way to study adaptive processes, including the evolutionary capabilities of species confronted with climate change (Dixon & Busch, 2017; Swab et al., 2015). Our findings demonstrated experimentally that ENMs have the potential not only to define environmental marginality (Lira‐Noriega & Manthey, 2014; Madeira de Madeiros et al., 2018), but also to predict genetically determined variation in a functional trait among populations (Manel et al., 2012; Thuiller et al., 2010; Wittmann et al., 2016). Moreover, this study highlights the value of conserving marginal populations during a time of rapidly changing global climate (Channell, 2004) based on their adaptive potential, which can derive from ecologically relevant genetic uniqueness (Papuga et al., 2018). Although the present results correspond to a single species, the implications and applicability of the proposed methodological framework for evolutionary biology—especially in the context of climate change—should be generalizable for other species and spatial scales.

In addition to correlative models, mechanistic ones may improve predictions of species’ responses to climate change. Such models can consider factors such as survival and reproductive rates, dispersal capability, competition–facilitation relationships, and measured physiological responses of species. However, they require much more data and the establishment of larger‐scale and more complex common garden experiments implemented across a variety of environmental conditions (Benito Garzón et al., 2011; Benito‐Garzón et al., 2019; Kearney & Porter, 2009; Valladares et al., 2014). Although the approach we propose does not include the data necessary to parameterize a mechanistic model, it can provide fundamental information and insights that can lead to further experiments that better facilitate mechanistic models. Either approach—or a combination—likely will prove feasible for taxa conducive to common garden treatments (e.g., plants or arthropods, Pelini et al., 2009; Williams et al., 2012). Paired with landscape genetics and genomic approaches, studies such as this one that combine modeling and experiments should lead to a greater unified understanding of adaptive processes under changing climates and other anthropogenic alterations of the environment (Exposito‐Alonso et al., 2018, 2019; Fournier‐Level et al., 2011).

CONFLICT OF INTEREST

The authors have no conflict of interest to declare.

AUTHOR CONTRIBUTIONS

JML, JMI, RPA, and JMK designed the framework integrating niche models and phenotypic assays. JML, JMI and CLR designed the experimental study, collected the samples, and conducted the field work. JSC undertook the leaf cell resistance experiment laboratory work. JML developed the niche models in collaboration with JMK, RPA, and JMSD. JML analyzed phenotypic data with the help of CLR. JML wrote the paper with the help of RPA, JMI, JMK, CLR, and JMSD. All authors reviewed the paper and approved the final manuscript.

Supplementary Material

Click here for additional data file.