Intraspecific variation in phenotypic plasticity can affect the ability of populations, and thus species, to respond to environmental changes. However, the prevalence and drivers of such variation are not well known. Most proposed explanations for intraspecific variation in phenotypic plasticity involve mechanisms associated with a population’s position within the species’ geographic range or the environmental heterogeneity experienced by the population. To assess the effect of these two drivers, and their potential interaction, we use a combination of germination and greenhouse experiments to measure thermal phenotypic plasticity in traits ranging from germination probability to flower abundance in populations of three Hypericum species sampled across their European ranges. We then relate thermal plasticity to each population's position within the species’ range and to the environmental heterogeneity of the sampling site. Our results revealed that, while average thermal plasticity in several traits was similar among the three tested Hypericum species, it varied among the conspecific populations. Specifically, populations from closer to the range edge tended to be more plastic in germination probability and plant height, while populations from more heterogeneous environments tended to be more plastic in flowering phenology, plant height, and flower abundance. Interestingly, for plasticity in germination phenology, plant height, and flower abundance, we found a substantial interactive effect with accentuated plasticity in heterogeneous sites near the range edge. This suggests that populations in heterogeneous environments at range edges may adjust to environmental change via phenotypic plasticity more effectively than do other conspecific populations. These results support both the tested drivers and reveal important interactive patterns for some of the tested traits. Furthermore, they encourage further research on plasticity that considers both range position and environmental heterogeneity.

Study populations

We acquired seeds of the three species from different parts of their native ranges, totaling 8 H. montanum, 18 H. perforatum, and 12 H. maculatum populations (after filtering as described in Section 2.5), either by collecting them in the field or acquiring seed accessions from European seed banks (Fig. 1, Table S1). The seeds were primarily collected between 2017-2021, with one H. perforatum population collected in 2007 and one H. montanum population collected in 1998. The seed material from seed banks and our field collections was collected according to ENSCONET guidelines (generally originating from at least 50 individuals; ENSCONET, 2009). For the self-collected material, seeds were sampled per mother plant, and an equal number of seeds from each mother were pooled for each temperature treatment.

Germination experiments

The seeds were cleaned from debris with an aspirator (Agriculex DB-1, column seed cleaner CB115009) before sowing. For each population, we sowed a maximum of 200 seeds (depending on availability) 2 mm apart from each other on 8 petri dishes (25 seeds per petri dish, diameter 5.5 cm) filled with 1 % agar (Sigma agar, lots SLBL4283V & SLBX7044). After eight weeks of cold stratification at 4 °C, we placed the petri dishes in growth chambers (Incubator numbers: LMS Cooled COLD 13234/20P4, MEAN 13233/20P4, WARM 13235/20P4, HOT 132636/20P4, United Kingdom), set to four temperature treatments. Daytime temperatures (16 h) were set to 16 °C (cold), 20 °C (medium), 24 °C (warm), and 28 °C (hot). The nighttime temperatures (8 h) were set to 10 °C below the daytime temperature. Photoperiod in all treatments was 16/8 h light/dark. Each treatment was further divided into two replicates that were located on different shelves of the same growth chamber. The choice of temperature treatments was based on data on average summer temperatures at the trailing, core, and leading areas of the study species’ ranges, as well as the predicted thermal conditions at the trailing edge in 2070, using data from WorldClim (Fick & Hijmans, 2017).

Greenhouse experiments

Seeds for the greenhouse experiments were cold-stratified for 4 weeks at 4 °C in dry paper bags. For each population, we sowed a maximum of 200 seeds (depending on availability) in 8 trays (25 seeds per tray) filled with sowing mixture (Kekkilä sowing mixture W HS R8017; KEK31116) covered by a thin layer (1-2 mm) of coarse sand. The trays were placed in greenhouse compartments with four temperature treatments identical to those used in the germination experiment, except that the night-time temperature for the vegetative stage was changed to 8 °C below the daytime temperature. Each treatment was further divided into two replicates placed in distinct greenhouse compartments. After the seeds had germinated, a maximum of 10 seedlings (depending on availability) were randomly chosen from each population and potted into individual 1L pots filled with soil (Kekkilä Professional coarse potting mixture; KEK33933). The pots were placed on a water-retaining rug on growing tables. The plants were watered by an automated watering system by soaking the rug underneath the pots. The watering schedule in each treatment was adjusted to keep the plants equally moist in all treatments. The plants were grown in the greenhouses from December 2021 to May 2022. At the end of March, the watering system broke, leaving the H. perforatum plants dry for some days in replicate A. We accounted for this in the interpretation of the results. The plants were fertilized six weeks after sowing with a 0.075 % solution of Kekkilä Turve Superex (NPK 12–5–27) and subsequently every two weeks with a 0.2 % solution of the same fertilizer. To avoid any effects of differing conditions within the greenhouse compartment, the germination trays and pots were periodically rotated (dates of rotation: Dec 8, 2021, Dec 15, 2021, Dec 22, 2021, Jan 19, 2022, Jan 26, 2022, Feb 2, 2022, Mar 4, 2022, Apr 8, 2022).

Trait measurements

During the germination experiments, we recorded the number of germinated seeds weekly over four weeks, based on seedlings having a radicle longer than 2 mm. Radicle length was defined as the distance from the radicle tip to the last occurrence of visible root hairs. Once a seed had been recorded as germinated, it was removed from the petri dish. In the fifth week, we performed cut-tests to determine the viability of the remaining seeds, i.e,. whether the seed had germinated during the last week of the experiment, or whether it was full, empty, moldy, or infested. At this stage, seedlings with a root radicle shorter than 2 mm were also scored as germinated. Based on this information, we determined the number of viable seeds (excluding empty and infested seeds) following the recommendations of the Millennium Seed Bank (Germination testing: procedures and evaluation. Millennium Seed Bank Partnership, 2022). For estimating germination in further analyses, we included only populations where the proportion of seeds that germinated out of the total number of viable seeds was more than 5%. This was done to exclude accessions where low germination was likely due to incomplete dormancy release. From these measurements, we derived two traits for our analyses: germination proportion, defined as the proportion of viable seeds that germinated, and germination phenology, defined as the number of days from sowing to germination for seeds that emerged during the 4-week monitoring period.

During the greenhouse experiments, we measured three traits: flowering phenology (the number of days from sowing to the onset of flowering), plant height (the length of the longest branch in cm), and flower abundance (the number of flowers produced by the end of the experiment, including withered and open flowers as well as full and empty seed capsules). Plant height and flower abundance were only measured for H. montanum and H. perforatum due to time limitations during data collection.

Data analyses

All data analyses were conducted in RStudio, version 2025.09.2 (Posit team 2025; R Core Team, 2025).

Measuring range position

We measured range position using two alternative metrics: 1) the geographic distance between the population and the edge of the species’ range (DRE), and 2) the climatic distance between the conditions at each source site and the average conditions across the range of the species (DCE). To calculate these metrics, we digitized distribution maps from Hultén & Fries (1986) using QGIS (Transformation type: Thin Plate Spline; Resampling method: Lanczos; Coordinate system: Arctic Polar Stereographic EPSG:3995) (QGIS.org, 2025). For DRE, we calculated the distance (in kilometers) from each study population’s geographic location to the northernmost or southernmost range edge, whichever was located closer to the population (Fig. 2A). To calculate DCE, we extracted and scaled (mean = 0, SD = 1) bioclimatic variables (WorldClim; Fick & Hijmans, 2017) for each ~7 x 7 km (5 minutes) resolution cell that fell within the species’ range and for each sampled population. We then conducted a Principal Component Analysis (PCA; Greenacre et al., 2022) for each species, and formed a convex hull around the first and second principal components to represent the species’ climatic niche. The DCE metric was then calculated as the climatic distance from each study population's position to the nearest edge of the convex hull. Both metrics were finally scaled to zero mean and unit variance (mean = 0, SD = 1).

Measuring environmental heterogeneity

We measured environmental heterogeneity using three alternative metrics, two of which focus on different aspects of land cover heterogeneity, and one describing topographic heterogeneity. Land cover and topography were selected as the landscape features of interest because of their known impact on variation in microclimatic temperatures (Barry & Blanken, 2016; Opedal et al., 2015), which in turn may select for thermal plasticity (Graae et al., 2018). More specifically, we calculated the Shannon diversity of land cover types (SHDI), which represents compositional heterogeneity in land cover, the mean of the perimeter-area ratio of land cover patches (PAR), which represents configurational heterogeneity in land cover, and average roughness in elevation (ARE), which represents topographic heterogeneity. To calculate SHDI and PAR, we used the function calculate_lsm() from the ‘landscapemetrics’ R package (Hesselbarth et al., 2019), and to calculate ARE, we used the sa() function from the ‘geodiv’ R package (Smith et al., 2021). We calculated all three metrics within a 500-m square buffer (i.e., a square extending 500 m from the central point in all directions) around each study population. A square buffer was used to enable including whole grid cells in the buffer. Deciding the size of the square buffer was based on the estimated dispersal distance, according to which the studied Hypericum species would be able to reach areas 500 meters from the seed collection site within 100 years given yearly sexual reproduction (with a maximum dispersal distance of 5 meters, Lososová et al., 2023). For SHDI and PAR, we used S2GLC 2017 land cover data (Malinowski et al., 2020), and for AR, E we used EEA-10 Copernicus DEM from the year 2022 (Copernicus Sentinel data, 2022), both in 10-meter resolution. All three metrics were finally scaled to zero mean and unit variance (mean = 0, SD = 1).

Quantifying the effects of range position and environmental heterogeneity

Estimating thermal plasticity

We fitted generalized linear mixed effects models (GLMMs) implemented with the glmmTMB() function in the ‘glmmTMB’ R package (McGillycuddy et al., 2025) to model each trait as a population-level linear function of temperature. We allowed populations to differ both in their mean deviation and in the slope of their reaction norms by including the interaction between population and treatment as a random effect (i.e., we fitted random-regression models) (Arnold et al., 2019). To account for the nonindependence of individuals grown in the same greenhouse chamber or incubator, we included replicate ID as a random effect. We sea t Gaussian error distribution with an identity link function for all traits except germination, for which we set a binomial error distribution with a logit link function. We treated temperature as a continuous variable, which we scaled to zero mean and unit variance (mean = 0, SD = 1) to facilitate model fitting. We evaluated normality and homoscedasticity of residuals by visually inspecting QQ-plots and by plotting residuals against fitted values and statistically using Shapiro-Wilk and Breusch-Pagan tests, and found minor deviations from normality and homoscedasticity. Flower abundance was square-root-transformed to better meet the assumptions. We extracted the random regression-slope coefficients for each population, converted them to absolute values, and used them as estimates of the magnitude of the population-specific plastic response to temperature in the subsequent analyses. To facilitate interpretation, we derived germination probabilities from the logit-scale model estimates by applying the inverse logit function: p = exp(η)/1+exp(η), where η represents the linear predictor. We then calculated treatment effects as the difference between treatment probability and control probability. We used AIC (Akaike’s Information Criterion) comparisons to assess statistical support for population-specific responses to temperature by comparing models including only a random intercept to models including both a random intercept and a random slope for the population.

Selecting metrics with the highest explanatory power for variation in thermal plasticity

We fitted linear models (LMs) to test the effect of the two range position metrics (DRE and DCE) and the three environmental-heterogeneity metrics (SHDI, PAR, and ARE) on thermal plasticity. For each trait, we fitted twelve models in tota, to test all combinations of one range-position and one environmental-heterogeneity variable, as well as their interaction. Thus, each model included a maximum of two explanatory variables and their interaction. We also tested whether accounting for non-independence between the populations of the same species improved the models, but found no qualitative change in the results. We computed Spearman correlations among explanatory variables using the cor_test() function in the ‘rstatix’ R package (Kassambara, 2025), and all were below 0.7, indicating no problematic multicollinearity (Dormann et al., 2013). We ranked the models based on AIC values, and for each trait selected the model with the lowest AIC value (ΔAIC >2; Burnham & Anderson, 2004; Symonds & Moussalli, 2011). If the highest-ranked model included an interaction, this parameter was also incorporated into the final models. If there were no detectable differences between the models, we selected the one with the lowest AIC value for further investigation and assesse,d case by case, whether any of the competing models provided additional insight. We evaluated normality and homoscedasticity of residuals by visually inspecting QQ-plots and by plotting residuals against fitted values and statistically using Shapiro-Wilk and Breusch-Pagan tests, and found no major deviation from normality and homoscedasticity.

Analyzing the effect of range position and environmental heterogeneity

To analyze the effect of range position and environmental heterogeneity on among-population thermal plasticity, we assessed the parameter estimates of the models selected in the previous step (Section 2.6.2.2). Additionally, we examined the amount of variance explained by each term by calculating their relative importances with the calc.relimp() function from the ‘relaimpo’ R package (Groemping, 2005). Predictions were produced using the ggpredict() function in the ‘ggeffects’ R package (Lüdecke, 2018).

References

• Arnold, P. A., Kruuk, L. E. B., & Nicotra, A. B. (2019). How to analyse plant phenotypic plasticity in response to a changing climate. The New Phytologist, 222(3), 1235–1241. https://doi.org/10.1111/nph.15656
• Barry, R. G., & Blanken, P. D. (2016). Microclimate and Local Climate. Cambridge University
• Burnham, K. P., & Anderson, D. R. (2004). Multimodel Inference: Understanding AIC and BIC in Model Selection. Sociological Methods & Research, 33(2), 261–304. https://doi.org/10.1177/0049124104268644
• Copernicus Sentinel data (2022). https://doi.org/10.5270/ESA-c5d3d65
• Dormann, C. F., Elith, J., Bacher, S., Buchmann, C., Carl, G., Carré, G., Marquéz, J. R. G., Gruber, B., Lafourcade, B., Leitão, P. J., Münkemüller, T., McClean, C., Osborne, P. E., Reineking, B., Schröder, B., Skidmore, A. K., Zurell, D., & Lautenbach, S. (2013). Collinearity: A review of methods to deal with it and a simulation study evaluating their performance. Ecography, 36(1), 27–46. https://doi.org/10.1111/j.1600-0587.2012.07348.x
• ENSCONET (2009). ENSCONET seed collecting manual for wild species. p. 29.
• Fick, S. E., & Hijmans, R. J. (2017). WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology, 37(12), 4302–4315. https://doi.org/10.1002/joc.5086
• Graae, B. J., Vandvik, V., Armbruster, W. S., Eiserhardt, W. L., Svenning, J.-C., Hylander, K.,  Ehrlén, J., Speed, J. D. M., Klanderud, K., Bråthen, K. A., Milbau, A., Opedal, Ø. H.,  Alsos, I. G., Ejrnæs, R., Bruun, H. H., Birks, H. J. B., Westergaard, K. B., Birks, H. H., & Lenoir, J. (2018). Stay or go – how topographic complexity influences alpine plant population and community responses to climate change. Perspectives in Plant                           Ecology, Evolution and Systematics, Special Issue on Alpine and Arctic Plant Communities : A Worldwide Perspective, 30, 41–50. https://doi.org/10.1016/j.ppees.2017.09.008
• Germination testing: procedures and evaluation. Millennium Seed Bank Partnership (2022). https://brahmsonline.kew.org/Content/Projects/msbp/resources/Training/13a-Germination-testing-procedures.pdf
• Greenacre, M., Groenen, P. J. F., Hastie, T., D’Enza, A. I., Markos, A., & Tuzhilina, E. (2022). Principal component analysis. Nature Reviews Methods Primers, 2(1), 100. https://doi.org/10.1038/s43586-022-00184-w
• Groemping, U. (2005). relaimpo: Relative Importance of Regressors in Linear Models (p. 2.2-7) [Dataset]. https://doi.org/10.32614/CRAN.package.relaimpo
• Hesselbarth, M. H. K., Sciaini, M., With, K. A., Wiegand, K., & Nowosad, J. (2019). landscapemetrics: An open‐source R tool to calculate landscape metrics. https://doi.org/10.1111/ecog.04617
• Hultén, E., & Fries, M. (1986). Atlas of North European vascular plants: North of the Tropic of Cancer. Koeltz Scientific Books.
• Kassambara, A. (2025). rstatix: Pipe-Friendly Framework for Basic Statistical Tests (Version 0.7.3) [Computer software]. https://cran.r-project.org/web/packages/rstatix/index.html
• Lososová, Z., Axmanová, I., Chytrý, M., Midolo, G., Abdulhak, S., Karger, D. N., Renaud, J., Van Es, J., Vittoz, P., & Thuiller, W. (2023). Seed dispersal distance classes and dispersal modes for the European flora. Global Ecology and Biogeography, 32(9), 1485–1494. https://doi.org/10.1111/geb.13712
• Lüdecke, D. (2018). ggeffects: Tidy Data Frames of Marginal Effects from Regression Models. Journal of Open Source Software, 3(26), 772. https://doi.org/10.21105/joss.00772
• Malinowski, R., Lewiński, S., Rybicki, M., Gromny, E., Jenerowicz, M., Krupiński, M., Nowakowski, A., Wojtkowski, C., Krupiński, M., Krätzschmar, E., Schauer, P. (2020), Automated Production of a Land Cover/Use Map of Europe Based on Sentinel-2 Imagery https://doi:10.3390/rs12213523
• McGillycuddy, M., Popovic, G., Bolker, B. M., & Warton, D. I. (2025). Parsimoniously Fitting Large Multivariate Random Effects in glmmTMB. Journal of Statistical Software, 112, 1–19. https://doi.org/10.18637/jss.v112.i01
• Opedal, Ø. H., Armbruster, W., & Graae, B. (2015). Linking small-scale topography with microclimate, plant species diversit,y and intra-specific trait variation in an alpine landscape. Plant Ecology & Diversity, 8, 305–315. https://doi.org/10.1080/17550874.2014.987330
• Posit team (2025). RStudio: Integrated Development Environment for R. Posit Software, PBC, Boston, MA. http://www.posit.co/
• QGIS.org (2025). QGIS Geographic Information System. QGIS Association. http://www.qgis.org
• R Core Team (2025). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.
• Smith, A. C., Dahlin, K. M., Record, S., Costanza, J. K., Wilson, A. M., & Zarnetske, P. L. (2021). The geodiv R package: Tools for calculating gradient surface metrics. Methods in Ecology and Evolution, 12(11), 2094–2100. https://doi.org/10.1111/2041-210X.13677
• Symonds, M. R. E., & Moussalli, A. (2011). A brief guide to model selection, multimodel inference, and model averaging in behavioural ecology using Akaike’s information criterion. Behavioral Ecology and Sociobiology, 65(1), 13–21. https://doi.org/10.1007/s00265-010-1037-6