Intraspecific variation influences performance of moss transplants along microclimate gradients.

Identifying the environmental drivers of population dynamics is crucial to predict changes in species abundances and distributions under climate change. Populations of the same species might differ in their responses as a result of intraspecific variation. Yet the importance of such differences remains largely unexplored. We examined the responses of latitudinally distant populations of the forest moss Hylocomiastrum umbratum along microclimate gradients in Sweden. We transplanted moss mats from southern and northern populations to 30 sites with contrasting microclimates (i.e., replicated field common gardens) within a forest landscape, and recorded growth and survival of individual shoots over 3 yr. To evaluate the importance of intraspecific variation in responses to environmental factors, we assessed effects of the interactions between population origin and microclimate drivers on growth and survival. Effects on overall performance of transplanted populations were estimated using the product of survival and growth. We found differences between southern and northern populations in the response to summer temperature and snowmelt date in one of three yearly transitions. In this year, southern populations performed better in warm, southern-like conditions than in cold, northern-like conditions; and the reverse pattern was true for northern populations. Survival of all populations decreased with evaporation, consistent with the high hydric demands and poikilohydric nature of mosses. Our results are consistent with population adaptation to local climate, and suggest that intraspecific variation among populations can have important effects on the response of species to microclimate drivers. These findings highlight the need to account for differential responses in predictions of species abundance and distribution under climate change.

Introduction

Understanding responses to changes in climate and environmental conditions is key to predict changes in species abundances and distributions. Responses to environmental drivers might vary not only among, but also within, species, because populations are adapted to their local environment (Clausen et al. 1940, Leimu and Fischer 2008, Hereford 2009, Malyshev et al. 2016). For example, populations at higher latitudes have been shown to compensate for shorter growing seasons by faster development, earlier reproduction, or higher survival than populations at lower latitudes (Chapin and Chapin 1981, Laugen et al. 2003, Doak and Morris 2010). As a result of local adaptations, climate change might affect populations across the species latitudinal range differently (Harte et al. 2004, Reich and Oleksyn 2008, Putnam and Reich 2017, Peterson et al. 2018). For instance, warming may affect high‐latitude populations more negatively than low‐latitude populations of vascular plants if the former are locally adapted to low temperatures (De Frenne et al. 2011, Beierkuhnlein et al. 2011, Peterson et al. 2018). Still, the extent to which local adaptation influences species’ response to environmental changes remains poorly understood.

Classical reciprocal transplantation experiments can inform us about the extent of local adaptation among populations by comparing the performance of local vs. foreign genotypes in each genotype’s local environment (Clausen et al. 1940, Kawecki and Ebert 2004, Hereford 2009). However, because reciprocal transplantations typically include few sites representing limited environmental variation, they provide insufficient information about how populations differ in their responses to different environmental drivers. An alternative to reciprocal transplantations is to transplant populations to multiple common gardens along environmental gradients (de Villemereuil et al. 2016). This allows comparison of responses to explicit environmental drivers among populations and enables statistical testing of effects of the interactions between drivers and origin of populations. Another limitation of many transplantation experiments is that they assess population responses only in terms of single fitness components or vital rates (Leimu and Fischer 2008). This is problematic because vital rates might respond in opposing ways to environmental drivers (Doak and Morris 2010, Peterson et al. 2018, Pironon et al. 2018). Assessments of population responses based on integrated fitness measures, such as population growth rate (λ) or the product of survival and growth thus provide more relevant estimates of local adaptation than single fitness components (Joshi et al. 2001, Chen and Schemske 2015, Peterson et al. 2016). However, approaches identifying differences in the effects of explicit environmental drivers on integrated measures of performance among populations are limited.

Environmental drivers affecting population performance need to be measured at relevant spatial scales. Local climate or microclimate, rather than average large‐scale climate, drives population processes, as it represents the conditions that individuals experience in situ (Bramer et al. 2018). Microclimate variation can determine species abundance and distributions by buffering the impact of large‐scale climate change on population extinctions and range shifts (De Frenne et al. 2013, Suggitt et al. 2018). Widening our knowledge about the effects of microclimate on populations’ performance will thus serve to refine predictions of future species abundance and distribution patterns (Bramer et al. 2018, Lembrechts et al. 2019). For this reason, identifying explicit microclimate drivers that affect vital rates and population dynamics is drawing increasing attention (Nicolè et al. 2011, Weegman et al. 2017, Oldfather and Ackerly 2019). Yet, knowledge about differences in responses to microclimate among populations is scarce (but see, e.g., Oldfather and Ackerly 2019).

Bryophytes constitute an excellent group of organisms to examine differences in population responses to microclimate, as they are highly sensitive to changes in their local environment (Busby et al. 1978, Hylander et al. 2002, Proctor 2008). This is because they are poikilohydric; that is, they lack mechanisms to regulate water loss, and their growth and survival depends on their hydration status (Proctor 2008). Understanding bryophyte responses to environmental changes in ecosystems where they play essential ecological roles is particularly relevant, for example, in boreal forests (Fenton et al. 2015). Many forest bryophytes occur in specific environmental conditions. For instance, some species are more common on north‐ than south‐facing slopes, suggesting that they are adapted to certain microclimates (Åström et al. 2007). In fact, microclimate variation strongly influences bryophytes’ abundance and distribution at regional and stand scales (Fenton et al. 2015). Genetic differentiation patterns associated with gradients of disturbance, elevation, or pH have been documented for mosses (Wyatt 1992, Cronberg 2004, Mikulášková et al. 2015), suggesting that populations might also differ in their response to microclimate variation. However, to the best of our knowledge, and despite the ease of translocating and monitoring mosses, transplantation experiments examining differences among populations in responses to microclimate variation are missing (see Jägerbrand et al. 2014, Doherty et al. 2018 for controlled experiments).

We examined intraspecific variation in the response to microclimate variation in the forest moss Hylocomiastrum umbratum. Specifically, we hypothesized that (1) different microclimate drivers affect vital rates and population performance differently, (2) the same microclimate driver may have different effects on different vital rates, and (3) populations across the species range vary in their responses to microclimate. To test these hypotheses, we transplanted moss shoots within mats from six populations from its southern and northern distribution ranges in Sweden to 30 field common gardens with different microclimates. We assessed the effects of microclimate variation on shoot survival and growth, as well as on overall performance estimated as the product of shoot survival and growth, and compared the overall performance of shoots from northern and southern populations along microclimate gradients.

Material and Methods

Study species

Hylocomiastrum umbratum is a long‐lived clonal pleurocarpous moss inhabiting boreal and hemiboreal forests (Hedenäs et al. 2014). Its distribution is incompletely circumpolar, occurring mainly in suboceanic and oceanic parts of Europe, Asia, and North America (Ratcliffe 1968; Koponen 1979, www.gbif.org). In Sweden it is common in the central area, occurring scattered towards the southernmost and northernmost parts Hedenäs et al. 2014; www.artdatabanken.se). It mainly occurs under shaded and moist conditions directly on the forest floor or on logs and boulders (Hedenäs et al. 2014).

Hylocomiastrum umbratum forms mats made of entangled shoots which are relatively easy to monitor, given their size and growth pattern similar to that of Hylocomium splendens (Økland 1995). Each shoot consists of a chain of annually produced segments (clonal modules) developed during the growing season (Hylander et al. 2002; Fig. 1). In Sweden, the growing season extends from late spring to the first snowfall (ca. May–December), but the total length differs by several weeks between northern and southern regions (Appendix S1: Table S1). Conditions enhancing shoot growth, namely combinations of high moisture and light availability, probably concentrate in summer and autumn (K. Hylander, personal observation). One or several new segments emerge in late autumn as small buds (growing points). Growing points grow and become segments during the following growing season (Fig. 1A). Occasionally, the segment instead of forming new buds elongates by apical growth the next season (Fig. 1B). New multiple segments may sometimes also develop from broken old segments (Fig. 1C). Old segments are progressively overgrown by newly produced segments and eventually decompose, although sometimes also old segments produce new segments, increasing the complexity of the branching pattern. Up to four living segment generations can usually be distinguished within a shoot. It is a dioecious species and sporophytes are quite rarely produced (Hedenäs et al. 2014, K. Hylander and C. J. Dahlberg, personal observation). Thus, as in most clonal mosses, local population growth rate is more influenced by vegetative growth (i.e., shoot branching) than by sexual reproduction (During 1979, Økland 1995).

Figure 1

The growth pattern of Hylocomiastrum umbratum illustrated by (A) a typical shoot in late autumn where a new growing point (G, small new segment) emerges from a segment (S) that grew the previous autumn from a previous segment (old S); (B) apical growth of a segment that has continued to grow from the tip (ap) of an old segment; (C) two segments that have grown from a broken old segment (b); (D) a transplant mat at the beginning of the experiment in June 2012; (E) a growing point (G) marked by a PVC ring (black) and measured in the first year; (F) subsequent marking of growing points and measures of segments the next year. The growing point in (E) corresponds to the new segment in (F) after a yearly transition. Thick lines: segments or growing points; thin lines: lateral branches.

Collection and transplantation area

In May 2012, we collected mats of H. umbratum from six populations in mature forests in two distant regions of origin: three populations in the southern distribution range in southern Sweden (boreonemoral zone), and three populations in its northern distribution range in northern Sweden (boreal subzone) (Sjörs 1999; Fig. 2A). Regions of origin were approximately 600 km apart to increase likelihood of finding genetic differentiation between regions. The northern populations experience a shorter growing season (ca. 10–40 d less) and a longer period with snow cover (ca 25–75 d longer) than the southern populations (average for 1961–1990, Swedish Meteorological and Hydrological Institute 2018; Appendix S1: Table S1). Winter, autumn and spring are colder in the northern than in the southern region (2.6°C lower on average; average for 2000–2010, Meineri and Hylander 2017). However, mean summer temperature and annual precipitation are similar between regions (14.9 and 15.4°C for the northern and southern region, respectively [Meineri and Hylander 2017]; and ca. 750 mm in both regions [Swedish Meteorological and Hydrological Institute 2018 ]).

Figure 2

(A) Location of the six source populations of Hylocomiastrum umbratum in two regions of origin (north and south) transplanted in the northern area of the species distribution in Sweden (red square). Environmental data of the source populations are shown in Appendix S1: Table S1. (B) Display of the 30 field common gardens (cross) in the transplantation area. Background maps correspond to elevation (the darker the higher elevation), streams, and water bodies (source: ©Lantmäteriet Gävle 2014 (I2014/00691)).

The transplantation area is situated in central‐north Sweden, within the northern area of the species distribution (southern Ångermanland, latitude 62°47–63°08 N, 1,200 km2; Fig. 2B). Altitude in this hilly area ranges from 23 to 470 m above sea level. The bedrock is composed of gneiss with podzolic soils (the Swedish Mapping, Cadastral and Land Registration Authority; www.lantmateriet.se). The area belongs to the middle boreal subzone (Sjörs 1999) and is dominated by coniferous forests (see Dahlberg et al. 2014 for details). The regional climate of the area is similar to that of the origin of the northern populations (Appendix S1: Table S1).

Transplantation experiment and demographic data collection

We transplanted moss mats in 30 forest sites in a replicated field common garden design (Fig. 2B). To increase variation in climate conditions within the transplantation area we did a stratified selection of sites across topographic gradients. Based on aerial photos, we selected 30 mature Norway spruce (Picea abies [L.] Karst.) forest stands across a range of slope aspects to include as much microclimate variation caused by differences in incoming solar radiation as possible (Appendix S1: Table S2). The forest stands were at least 5 km away from rivers and 20 km from the Bothnian Sea to reduce their climate influence (Vercauteren et al. 2013). Within each forest stand, we selected a 4 × 4 m transplantation site situated at least 50 m from open areas and streams, 25 m from younger forest stands, and 10 m from vertical cliffs higher than 5 m, with mesic soils and homogeneous terrain in a 50‐m radius. The tree layer surrounding the sites was dominated by Norway spruce with sparse Betula spp. trees (Appendix S1: Table S2).

Before transplantation, we stored the moss mats on the ground in shaded moist spruce forests near the transplantation area. We placed six H. umbratum transplants (mats of about 10‐cm diameter; Fig. 1D), one from each population, in each of the 30 sites during 1–17 June 2012. At each site, we chose a 2 × 0.5 m plot on flat ground. To standardize soil conditions among sites we removed the vegetation down to the soil and added a 5‐cm layer of planting peat without additives (pH 3.5–4.5; Hasselfors Solmull/Sphagnum ugødet/Naturtorv) within each plot. We attached the transplants to the soil with two wooden sticks and a U‐shaped rubber‐coated steel wire and placed them 20 cm from each other.

On individual shoots within the transplanted mats, we carried out two types of measurements. First, to estimate shoot survival for two transition intervals we randomly marked seven growing points (i.e., small segments) in different shoots within the central part of each transplant, with PVC rings in June 2012 (Fig. 1E; Økland 1995). After 1 yr, we measured the length of the marked segments and marked the new growing points produced in the same shoots. We subsequently marked new growing points and recorded shoot survival and growth in June 2014 and June 2015 (Fig. 1F). We considered that a shoot was dead when it lacked vitality, that is, was embrowned, or did not produce new segments (did not grow) during two consecutive transition intervals. We assumed that shoots that did not grow in one transition interval but had new segments in the following transition interval were dormant in the first transition interval. We could follow this procedure only for the second transition (2013–2014). For the first transition interval (2012–2013), we marked already growing shoots that had produced new very short segments (the ones that we marked) at the start of this interval in 2012. All marked shoots in the transition 2012–2013 were thus alive according to our definition of survival. For the third transition (2014–2015) we estimated survival in the following way. Because we could not decide whether shoots that did not grow for the first time in 2014–2015 were dead or dormant, we assumed that the proportion of dead and dormant shoots for that transition was the same as in the transition of 2013–2014 (6.2% and 9.1%, respectively). We thus calculated mortality in 2014–2015 as the probability of a shoot of not growing in 2014–2015 multiplied by the proportion of nongrowing shoots that were dead in 2013–2014. Survival was calculated as 1 ‐ mortality. Of the 1,260 shoots initially marked, 59% were impossible to relocate at the end of the experiment because of loss of rings.

Second, we estimated shoot growth as changes in shoot size (i.e., the natural logarithm of the total length of new segments [mm]) between years for three transition intervals based on retrospective measures of segments (Callaghan et al. 1997). We collected five random nonmarked shoots from each transplant at the end of the experiment in June 2015 (n = 15 shoots per region of origin at each site). We determined the length (±2 mm) and age of each segment by counting back from the last transition interval segment (2014–2015) to the time of transplantation (2012–2013) and before (2011–2012), assuming that shoots produced segments every year (Callaghan et al. 1997). We decided to use total shoot new length as a shoot size measure because it was not always possible to assess the length and kinship of individual segments in an unambiguous way. For example, new segments sometimes emerged both from last year’s segments and from older segments. This measure of growth did thus not account for the age of the segments from which new segments grew or the explicit branching pattern, but included branching implicitly by accounting for the length of all new segments. Transplants from each region of origin had similar initial shoot size (Appendix S1: Fig. S1). By using these retrospective measures of growth, we retrieved a large balanced sample with growth data also for the first transition interval. We did not pool retrospective and prospective measurements of growth because we did not have both types of measurements for all transition intervals, and pooling the two data sets, which were collected in different ways, could thus bias comparisons among years.

Microclimate data

We measured the ground surface temperature using one logger (ibuttons, DS1922) per site during the three study years (19 June 2012–10 June 2015). Each logger was placed inside two plastic zip bags underneath a mat of the moss Hylocomium splendens at the center of the plot, and recorded temperature (°C) every 70 min. We calculated the daily, monthly, and seasonal mean maximum temperature for each site. Because of logger failure or loss, we missed values of temperature in some months for some sites (4.4% of values). For these cases, we imputed values using the coefficients from linear regression models with mean maximum temperature as response and site + month × year as predictors. We fitted these maximum temperature models for summer (June, July, and August; F 37, 220 = 16.1, R 2 = 0.73, P < 0.001) and autumn months (September, October, and November; F 37, 215 = 470.2, R 2 = 0.99, P < 0.001) for the 3 yr and for the 30 study sites. We estimated the snowmelt day in spring for each site and year based on the ground temperature measurements (Lundquist and Lott 2008). When the standard deviation of the mean maximum daily temperature of three consecutive days is >0.5℃, it indicates that the logger is uncovered by snow. Using the coefficients of the linear regression model snowmelt day ~ site + year (F 31,52 = 4.75, R 2 = 0.74, P < 0.001), we imputed missing snowmelt day values for the years 2014 and 2015 (6.7% of values in total).

We recorded evaporation at each site by measuring the evaporation of distilled water contained in two narrow (12‐cm diameter), 50‐cm‐long plastic cylinders between mid‐June and mid‐September in 2013. At the first date, we filled the cylinders with water and attached them vertically to a thin wooden stick with their bottoms 10 cm above the ground. Their tops were sealed, and the cylinder bottom was only covered with a filter paper through which the water evaporated. We assessed the evaporation as the amount of lost water in each cylinder (in millimeters) and calculated a mean value from the two cylinders at each site (except in a few cases were only one cylinder worked; Appendix S1: Table S2). The cylinder works in a similar way as a Piche evaporimeter, for which the evaporation correlate to vapor‐pressure deficit in the air (relative humidity × saturated vapor pressure), but also to wind velocity (Stanhill 1962, Papaioannou et al. 1996). We assumed that this evaporation measurement is ecologically relevant for bryophyte performance, because evaporation through a filter paper resembles the passive transpiration from a moss leaf (cf. Proctor 2008).

We measured soil moisture (%) with a soil moisture meter (HH2, Delta‐T devices, UK) at each transplant and retrieved the average value for each site. Measurements were done once during four consecutive days without heavy rain (2–6 June 2015). To estimate canopy openness, we took a canopy digital photograph with a Canon PowerShot G11 camera at the center of each plot from a horizontal position at 50 cm from the ground in June 2012. We digitalized the pictures in ImageJ (1,46f version, Wayne Rasband, National Institutes of Health, Bethesda, Maryland, USA) and calculated the proportion of sky uncovered by tree canopy. For these drivers, measured only once, we assumed that the relative differences in microclimate among sites remained similar between years during our study.

Analyses of vital rates

We did all statistical analyses in R v. 3.4.3 (R Development Core Team 2017). We modeled vital rates using (generalized) linear mixed models (GLMMs) using the package lme4 (Bates et al. 2015). First, we modeled probability of not growing with the data from the monitored marked shoots, using data for all transition intervals, as a function of shoot Size, transition interval, and the microclimate drivers corresponding to each transition interval and their interactions with region of origin. If probability of not growing differed between transition intervals, we modeled it for each interval separately. Second, to detect differences in growth among transition intervals we modeled Size + 1 as a function of Size using the size data from all years from the retrospective shoot size measures, microclimate drivers corresponding to each of the three years at each site, transition interval, and their interactions with region of origin. If growth among transition intervals differed, we modeled growth for each transition interval separately, as described above. Because these models accounted for size in the previous year (size in the year t included as a predictor), differences in size in year t + 1 in the model can be interpreted as differences in growth. To account for the structure of the data we included source population identity and site as random factors in all models, and shoot identity in the growth model with data from all years to account for repeated measures on the same shoots. Binomial error and logit link‐functions were used for probability of not growing models (fitted by maximum likelihood criteria [ML]) and Gaussian error and identity function were used for the growth models (fitted by restricted maximum likelihood criteria [REML]). We also tested for nonlinear relationships between Size + 1 and Size both by including a quadratic term (Size 2) in the linear models, and by performing generalized additive mixed models (GAMMs) that allow nonmonotonic relationships (Dahlgren et al. 2011). Because neither model including quadratic terms or GAMMs were markedly better than models assuming linear relationships (data not shown), we used linear models for simplicity. Variance in growth was independent of size; the predicted vs fitted residuals of the Size + 1 ~ Size model were homoscedastic.

We estimated survival probability as 1‐ (probability of not growing × the proportion of shoots that were dead in the interval 2013–2014; 0.41). We did not consider effects of population density because initial density was similar for all transplants and did not change much during the study period. All continuous predictors were centered and standardized with the “scale” function. To decrease problems associated with multicollinearity of environmental predictors and to keep the number of predictors low, we included only linear effects of uncorrelated microclimate drivers in the models (Appendix S1: Table S2; Pearson correlation tests, r < 0.7 and P > 0.05).

We used model selection based on Akaike Information Criterion corrected for small sample sizes (AICc; Burnham and Anderson 2002) to identify the best combination of the fixed effects of Size, region of origin, microclimate drivers, and their interactions (region of origin × microclimate drivers, and region of origin × Size). For model selection we used the MuMIn package (Barton 2018). For all final models we calculated P values of the fixed effects based on Satterthwaite degrees of freedom approximation using the lmerTest package (Kuznetsova et al. 2017), and the marginal R 2 (R 2m) value to estimate the variance explained by the fixed effects (Nakagawa and Schielzeth 2013). Residuals of the models were checked for normality and homoscedasticity.

Analyses of overall performance of shoots

We used the product of survival and growth as an estimate of the overall performance of transplanted shoots, and used this measure to compare the response among moss transplants along microclimate gradients (cf. Chen and Schemske 2015). This combination of vital rates constitutes a measure of growth weighted by the probability of survival. Although it does not include shoot branching explicitly, it incorporates all effects of shoot branching on total shoot length and survival. Because we were unable to identify dormancy, and thus transitions to and from dormancy in the shoots used for estimating growth, dormant shoots were not included in these estimates. This implies that we implicitly assumed that growth of dormant shoots was equal to that of nondormant shoots. If the size of dormant shoots was smaller after dormancy than before dormancy, then this implies that we are slightly overestimating growth. Yet this bias is unlikely to influence differences in performance estimates across gradients of microclimate conditions, or between regions of origin.

The method we used to calculate overall performance can be described in four points. First, to calculate this integrated estimate of performance, we multiplied the probability of survival and growth of shoots given survival. For estimates of these vital rates, we used the estimates from the linear mixed models in Table 1. To assess effects over a representative sample of shoot sizes, we estimated vital rate products for a cohort of 100 shoots randomly drawn from the original shoot size distribution in the transition interval 2013–2014.

Table 1

Mixed linear models of effects of transition interval (Trans.; e.g., 1314 denotes the transition between 2013 and 2014), shoot size the previous year (Size; ln‐transformed total shoot new length in Year, mm), microclimate drivers (evaporation, summer temperature, snowmelt date) and region of origin of populations (north and south) on the probability of not growing at all (used to estimate survival) and growth (Size + 1, ln‐transformed total shoot new length in Year + 1, mm) of six populations of Hylocomiastrum umbratum transplanted along microclimate gradients in 30 field common gardens in northern Sweden. The data comprise the Wald‐type Z statistic (Z) for logistic models (probability of not growing), t‐test (t) based on Satterthwaite’s method for the Gaussian models (growth), coefficients and standard error (SE) of the predictors included in the models, significance value (P), marginal R 2 (R 2m) and standard deviation (SD) of random effects. Sample size for each model is in brackets.

	Fixed effect	Coefficient (SE)	Z/t	P	R 2m	Random effect
All transitions
Probability of not growing (n = 1,179)	Intercept	−1.88 (0.13)	−14.8	<0.001	0.089	SD intercept = 0.32
	Sizet	−0.53 (0.08)	−6.6	<0.001
	Trans. (1415)	−0.45 (0.20)	−2.2	0.025
Sizet + 1 (n = 2,603)	Intercept	3.22 (0.04)	71.7	<0.001	0.112	SDintercept = 0.30
	Sizet	0.09 (0.02)	5.7	<0.001		SDresiduals = 0.54
	Origin (S)	0.16 (0.05)	3.3	0.007
	Trans. (1314)	0.08 (0.04)	2.2	0.026
	Trans. (1415)	−0.32 (0.04)	‐8.3	<0.001
	Sizet × Origin (S)	0.05 (0.02)	2.0	0.043
	Origin (S) × Trans. (1314)	−0.12 (0.05)	−2.3	0.02
	Origin (S) × Trans. (1415)	−0.07 (0.05)	−1.3	0.193
Transition 2012–2013
Sizet + 1 (n = 870)	Intercept	3.19 (0.05)	67.5	<0.001	0.063	SD intercept = 0.20
	Sizet	0.11 (0.02)	5.9	<0.001		SDresiduals = 0.54
	Origin (S)	0.15 (0.06)	2.8	0.051
Transition 2013–2014
Probability of not growing (n = 748)	Intercept	−2.30 (0.20)	−11.3	<0.001	0.243	SDintercept = 0.29
	Sizet	−0.81 (0.10)	−7.8	<0.001
	Evaporation	0.30 (0.13)	2.3	0.022
	Origin (S)	0.37 (0.24)	1.5	0.13
	Tsummer	0.11 (0.18)	0.6	0.56
	Origin (S) × Tsummer	−0.91 (0.27)	−3.4	<0.001
Sizet + 1 (n = 875)	Intercept	3.32 (0.05)	67.7	<0.001	0.073	SDintercept = 0.20
	Sizet	0.13 (0.02)	6.8	<0.001		SDresiduals = 0.52
	Origin (S)	0.04 (0.06)	0.7	0.26
	Tsummer	0.03 (0.04)	0.8	0.40
	Origin (S) × Tsummer	−0.09 (0.04)	−2.6	0.008
	SnowmeltDate	0.07 (0.04)	1.9	0.07
	Origin (S) × SnowmeltDate	−0.11 (0.04)	−3.0	0.003
Transition 2014–2015
Probability of not growing (n = 471)	Intercept	−2.34 (0.00)	−659	<0.001	0	SDintercept = 0.46
Sizet + 1 (n = 858)	Intercept	2.97 (0.05)	55.3	<0.001	0.077	SDintercept = 0.30
	Sizet	0.18 (0.02)	8.7	<0.001		SDresiduals = 0.57

Origin: region of origin; S: south; Tsummer: maximum summer temperature at each site; SnowmeltDate: day of snow melting in spring at each site.

Second, to compare overall performance between the transition intervals 2013–2014 and 2014–2015, we estimated the overall performance at constant mean values of microclimate drivers using the vital rate effect estimates for each transition interval. If region of origin was significant in the vital rate models for each transition interval, we estimated mean overall performance for each region and averaged these values.

Third, to estimate the effects of a focal microclimate driver and region of origin on overall performance, we included significant effects of the microclimate predictor and its interaction with region of origin from the vital rate models. We did this for the transition interval 2013–2014 only, as no significant effects of microclimate drivers and region of origin on vital rates were detected in the 2014–2015 interval. We used region‐specific estimates from the vital rates models to calculate overall performance for each region of origin over the observed range of one microclimate predictor (100 values distributed evenly), while setting the other predictors constant at their mean values. We also explored differences in performance between the coldest (northern‐like) and warmest (southern‐like) environmental conditions, represented by combinations of observed minimum and maximum values of evaporation (50.5 and 139 mm for coldest and warmest environments, respectively), summer temperature (13.5 and 19.3°C for coldest and warmest environments, respectively) and snowmelt date (99 and 126 d for warmest and coldest environments, respectively) at the transplant sites.

Fourth, to account for uncertainty in estimates of effects of microclimate drivers on overall performance, we randomly sampled estimates of effects of significant microclimate predictors on the respective vital rate using the “rnorm” function, from a distribution of effect estimates based on effect estimates from the statistical models and their associated standard errors. We applied the sampled effect estimate to the cohort of 100 shoots, repeated this sampling procedure 200 times per region of origin, and constructed 95% confidence intervals.

Results

Survival and growth of transplanted H. umbratum shoots differed among transition intervals. Survival was slightly lower in 2013–2014 than in 2014–2015 (93.8% vs. 96.1%; n = 755 and 480 shoots, respectively), and increased with increasing shoot size in 2013–2014 but not in 2014–2015 (Table 1; Appendix S2: Fig. S1A). Shoot growth was significantly lower in 2014–2015 than in the other two transition intervals (Table 1; Appendix S2: Fig. S1B). Mean overall performance of transplants at average microclimate conditions was significantly higher in 2013–2014 than in 2014–2015 (mean [probability of survival × ln (Size + 1)] ± 1 SE = 3.171 ± 0.001, and 2.848 ± 0.002, respectively; analysis of variance [ANOVA], F 1,598 = 444,900, P < 0.001).

Effects of origin and microclimate on vital rates

Models of the data set including all transitions indicated that shoots from southern populations grew more than those from northern populations (significant effect of the interaction Size and region of origin in Table 1). However, the effect of the interaction was not significant in models for single transition intervals.

We found significant effects of microclimate drivers, region of origin, and their interactions on vital rates in the second transition interval (2013–2014), but not in the other two intervals (Table 1). In this interval, survival increased with decreasing evaporation for shoots from both regions (see probability of not growing in Table 1). Survival of shoots from southern populations also increased with increasing summer temperature, whereas survival of northern populations increased with decreasing summer temperature (Table 1, Fig. 3A). The effects of summer temperature and snowmelt date on growth were also in opposite directions in northern and southern populations (Table 1). Growth of southern populations increased with cold summer temperature, whereas growth of northern populations increased with warm summer temperature (Fig. 3B). Growth of southern populations increased with early snowmelt, whereas growth of northern populations increased with late snowmelt (Fig. 3C). Although the effects of the interactions between microclimate and region of origin were large and significant, the R 2m was < 0.25 in all models, implying that much variation in vital rates remained unexplained.

Figure 3

Predicted survival and growth of transplanted northern (purple‐dashed) and southern (orange‐solid) populations of Hylocomiastrum umbratum in the transition interval 2013–2014, for a shoot of average size based on the linear mixed models in Table 1. Survival (A) and growth (ln‐transformed total shoot new length in Year + 1, mm) for different summer temperatures (B), and snowmelt dates (C). The 95% confidence intervals for the significant interactions between fixed effects are shown.

Effects of origin and microclimate on overall performance

The effects of summer temperature and snowmelt date on overall performance differed between populations from different regions of origin (Fig. 4). High summer temperatures had weak positive effects on southern and northern populations (Fig. 4A). Low summer temperatures had more negative effects on transplants from southern populations than on transplants from northern populations (Fig. 4A). Overall performance increased with earlier snowmelt in southern populations but decreased in northern populations (Fig. 4B). Lastly, northern populations performed better in the coldest environment and southern populations performed better in the warmest environment (Fig. 4C).

Figure 4

Predicted overall performance estimated as the product of shoot survival (logit‐transformed) and shoot growth (ln‐transformed total shoot new length in Year + 1, mm). Estimated median with 95% confidence intervals of transplanted northern (purple‐dashed) and southern (orange‐solid) populations of Hylocomiastrum umbratum in the transition interval 2013–2014 for different summer temperatures (A), snowmelt dates (B), and at the observed extremes of the microclimate gradients at the transplantation sites, “cold” and “warm” (C) are shown.

Discussion

Our study demonstrated intraspecific differences in the response of growth, survival, and overall performance to microclimate variation in the forest moss H. umbratum in Sweden. Populations are expected to be adapted to the local climate (Leimu and Fischer 2008, Hereford 2009) and therefore to respond in different ways to environmental factors (Harte et al. 2004, Reich and Oleksyn 2008). In agreement with these expectations, we found that northern (high latitude) populations of H. umbratum performed better under cold, northern‐like conditions, whereas southern (low latitude) populations performed better under warm, southern‐like conditions. Specifically, the vital rate responses of northern and southern populations to differences in summer temperature and snowmelt date were in opposite directions. Such contrasting responses of populations to environmental factors have previously been reported in vascular plants. For instance, southern populations of maple trees were more negatively affected by winter temperature and more positively affected by summer temperature than northern populations (Putnam and Reich 2017). Similarly, survival of southern populations of Arabidopsis thaliana decreased with decreasing winter soil temperature and survival of northern populations was not affected (Ågren and Schemske 2012). Taken together, these results strongly suggest that populations can respond in opposite ways to environmental and climate change.

Survival of populations from both regions decreased with increasing evaporation (i.e., increasing drought), which is consistent with high moisture requirements of bryophytes like H. umbratum (Busby et al. 1978, Hylander et al. 2002, Hedenäs et al. 2014). The differences in the response of survival to summer temperature between populations from the two regions are in agreement with what one would expect based on local adaptation. However, the response of growth differed from what one would expect based on the conditions at the site of origin, as shoots from northern populations grew more at high summer temperature than at low temperatures, whereas southern shoots showed the reversed pattern (Fig. 3). One possible explanation for these patterns is that mosses from northern populations are adapted to maximize growth during the brief interval of the year when high temperatures occur, rather than to grow better at low temperatures. Jägerbrand et al. (2014) also found that mosses from cold areas grew faster than those from warmer areas under similar conditions. These responses are consistent with countergradient variation patterns where populations at northern latitudes evolve higher growth rates at a given temperature than populations at southern latitudes to counteract the effect of colder, shorter growing seasons (Levins 1968, Chapin and Chapin 1981, Conover and Schultz 1995, Laugen et al. 2003). Differences in the effects of environmental drivers among vital rates appear to be relatively common in vascular plants (Doak and Morris 2010, Nicolè et al. 2011, Peterson et al. 2018, Pironon et al. 2018), but have not been documented in bryophytes before. Our results thus emphasize the need to consider effects of environmental drivers on integrated measures of fitness such as overall performance, rather than on single vital rates.

By transplanting populations to multiple common gardens representing an extended set of microclimates and following the transplants over several years, we maximized the probability of detecting differences in responses to environmental drivers among populations. Patterns suggesting differential responses emerged only in the second of the three study years. The expression of genetically based differences among populations may result in performance differences only under environmental conditions appearing only in certain years (Matesanz and Ramírez‐Valiente 2019). Possibly, the unusually wet September in 2013 (data from the two closest weather stations in Kramfors and Sollefteå; not shown), or other aspects of climate conditions differing among the study years led to differences between regions being expressed only during one transition interval. Interestingly, the average performance was higher in the transition interval when these effects were found, suggesting that average conditions were more favorable in that year. On a more general note, we cannot exclude that the observed differences in performance among populations from different regions result from mechanisms other than local population genetic differentiation, such as maternal environment effects on the phenotypic responses (e.g., transgenerational plasticity [Latzel and Klimešová 2010]) or epigenetics (Verhoeven and Preite 2014), or from eco‐physiological acclimation to local environments. However, it is reasonable to assume that mosses quickly acclimatized to the transplantation conditions over the study period, as effects of region of origin were not evident the first year, when ecophysiological responses to transplantation effects should be strongest. In conclusion, our results suggest that effects of local adaptation might only emerge under certain conditions and that transplantation experiments therefore need to be carried out over a sufficiently long period to cover a broad spectrum of environmental conditions.

Using an innovative approach based on transplantation of latitudinally distant populations of a forest moss along multiple microclimate gradients, we showed that intraspecific variation strongly influenced population responses to environmental drivers. This finding has four important broader implications. First, populations across a species range may differ in their vulnerability to altered microclimate conditions because of climate or land‐use change. Second, sites within the distribution range of a species at a large scale are not equally suitable, stressing the need to include microclimate in predictions of species distribution (Lembrechts et al. 2019). Third, bryophytes are suitable study systems to assess effects of local population adaptation in a broad sense. Finally, because populations differ in their responses to microclimate drivers, niche modeling approaches such as species distribution models are likely to overestimate niche breadth if considering species to be homogeneous entities, and predictions of effects of climate change on the distribution and abundance of species may be considerably improved by accounting for intraspecific variation in the response to environmental drivers.

Click here for additional data file.

Click here for additional data file.