Materials and Methods
Study site
We utilized the large climatic gradients of the Cascade Range in
Washington, USA to transplant seeds along each of 2 large transects
spanning ~1200 m on the West and East side of the
Cascade Crest on the traditional lands of the Nlaka’pamux, Nooksack,
Okanagan, and Methow peoples. Both transects are characterized by
topographically complex terrain and differ substantially in their
macroclimatic characteristics. The West transect (Mount Baker National
Forest) is warmer and wetter than the East transect (Okanagan National
Forest) (Supporting Information). Both transects are characterized by
montane to subalpine species common in the Pacific Northwest, with an
abundance of cedars, firs, heathers, and understory forbs and grasses.
We selected 15 sites along each of our 2 transects (n = 30 sites) using
satellite images of accessible areas, and identifying areas (i.e.
blocks) that had both low and high tree canopy openness. In the field,
we established two blocks per site (n = 60 blocks) to encompass
different levels of canopy openness, with one block in the relatively
most open area and the other block in the relatively most closed canopy.
Species data
We sowed 25 species encompassing tree (Abies grandis, A.
lasiocarpa, Picea sitchensis, P. engelmannii, Pinus ponderosa, P.
contorta ), shrub (Mahonia nervosa, M. aquifolium, Rubus ursinus,
R. spectabilis, Sambucus cerulea, S. racemosa, Sorbus sitchensis,
Vaccinium parvifolium, V. deliciosum ), forb (Eriophyllum lanatum,
Anemone occidentalis, Erigeron perigrinus, Lupinus latifolius,
Maianthemum dilatum, M. racemosum, Tolmiea menziesii, Tellima
grandiflora ), and graminoid (Carex stipata, Carex spectabilis )
growth forms. We chose congeneric and confamilial pairs of species that
have different regional distributional characteristics (e.g. lower vs.
higher elevation, or wetter vs. drier side of the Cascade Crest), thus
aiming to capture a range of species potential responses. To facilitate
field identification of seedlings, we sowed these pairs onto separate
quadrats within each block for three plots (n = 180 plots) of three 0.25
x 0.25 m quadrat replicates (n = 540 replicates), leaving the third
quadrat unmanipulated to control for background recruitment. We also
included unpaired species with large seeds (L. latifolius ) or
with high regional prevalence (S. sitchensis, A. occidentalis ).
We opportunistically sourced seeds from nearby areas in 2016 and
purchased native seeds for those species for which we had no, or not
enough, locally sourced seeds.
We homogenized all seed sources for a given species and sowed seeds in a
sand mixture in September-October 2017. We recorded recruitment (i.e.
survived to end of first growing season) and yearly survival of
seedlings during the growing season (May-September) in 2018, 2019, and
2020, and recorded only surviving seedlings in 2022. We surveyed sites
three times during the growing season in 2018 and 2019, and visited
sites once at the end of the growing season in 2020 and 2022. We were
not able to access 10 of our sites in 2018 (wildfire closures) or any of
our field sites in 2021 (closed USA border).
Microhabitat data
We measured a suite of abiotic and biotic parameters to quantify
microhabitat (Table S1) and categorize these parameters as being
directly influenced by climate change or not. We measured some of these
variables annually throughout the first three years of the study period
and others in just the last year of the study (due to logistical and
financial constraints). Here, we assume that the rank order of
microhabitat differences among sites remained constant across years.
We measured aspects of microhabitat directly related to climate by
recording soil and air temperature, duration of snow cover, and soil
moisture at each block with two different data loggers (TMS-4 data
logger, TOMST, Prague Czech Republic; Wild et al., 2019; HOBO 64K
Pendant Temperature/Alarm data loggers, Onset, Bourne, Massachusetts
USA). We calculated seasonal variables from these data loggers, using
some functions from package ‘myClim’ (Man et al., 2023; Man, Kalčík,
Macek, Wild, et al., 2023) to generate the following biologically
meaningful microclimate explanatory variables: summer maximum soil
temperature, winter minimum soil temperature, spring days of snow
coverage, summer minimum soil moisture, and summer minimum plant-height
temperature (Table S1, Supporting Information).
We also measured aspects of microhabitat not directly related to climate
by quantifying abiotic and biotic soil microhabitat, namely soil fungal,
bacterial, organic carbon and nitrogen content, and water holding
capacity at each replicate to capture variability within each site. We
also measured canopy openness at each block. Upon individual site snow
melt-out in 2022, we took soil cores to conduct in situmeasurements of fungus:bacteria ratio (F:B) with a microBIOMETER® Test
Kit (Prolific Earth Sciences, Montgomery, New York, USA). We then
quantified soil organic carbon to nitrogen ratio (C:N) and water holding
capacity in the lab (Supporting Information).
Statistical analyses: microhabitat
variation
To answer how variation within sites compares to among sites of each of
our microhabitat parameters, we calculated the variance partition
coefficient from intercept-only linear mixed models (LMMs; package
‘lmerTest’; Kuznetsova et al., 2017) that included block nested within
site as a random effect (microhabitat variable ~ 1 +
(1|site/block)). In cases where LMMs did not converge, we
changed the random effect to include only site. To test if microhabitat
parameters covary with elevation, we fit LMMs with quadratic elevation,
transect, and their interaction and included site as a random effect
(microhabitat variable ~ poly(elevation, 2) * transect +
(1|site)). We also used a principal component analysis plot to
identify any clustering in the microhabitat data. We conducted all data
processing and analyses in R version 4.2.3 (R Core Team, 2023).
Statistical analyses: species
establishment
To interpret our results in terms of likelihood of recruitment, recruit
numbers, and seedling survival, we first transformed our species-level
yearly recruitment and seedling survival data to binary or proportional
responses by calculating: binary recruitment, relative recruit counts in
the first three years (2018, 2019, 2020), 1-year (2018-2019 or
2019-2020), 2-year (2018-2020), and multi-year (survived to 2022)
survival of seedlings. We controlled for background recruitment by
subtracting the recruit counts in control plots from paired seed
addition plots. From this, we then calculated relative recruit counts:
recruit countsspecies i, site-rep j / maximum recruit
countsspecies i. We calculated the proportion of
surviving seedlings in all replicates, including control replicates,
using all non-zero recruitment data and not allowing survival
probability > 1.
While we have general a priori hypotheses of how certain
microhabitat conditions affect species establishment (Table S1),
previous studies have shown that recruitment and seedling survival vary
greatly by species and we have limited prior knowledge of which
microhabitat parameters are important for each species. We thus chose a
data exploratory framework and used an information-theoretic approach
with model averaging (‘MuMIn’ package, Bartoń, 2022) to compare
ecologically meaningful microhabitat variables to identify the most
suitable ones for each species (Tredennick et al., 2021). We fit
separate binomial generalized linear models (GLMs) for each species’
life stage to determine which of the microhabitat variables described
above (uncorrelated with Pearson’s correlation coefficient <
0.7) are important for establishment (i.e. recruitment and seedling
survival). We only analyzed species for which we observed a response in
at least 8 plots for any given life stage.
Each of our global GLMs included all 8 microhabitat variables described
above, plus any quadratic effects identified as important by AICc in a
reduced model (‘AICcmodavg’ package; Mazerolle, 2020; Supporting
Information). To account for soil temperature effects on recruitment and
plant-height temperature on seedling survival, we fit recruitment models
with soil temperature and seedling survival models with plant-height
temperature. To account for microhabitat variation not captured by any
of these parameters, we further included canopy openness and transect as
fixed effects (life stage success ~ microhabitat
variables + canopy openness + transect, family = binomial(link =
‘logit’)). In 1-year survival models, we also included year as a fixed
effect to account for differences among years in overall seedling
survival and different frequencies of site visits (life stage success
~ microhabitat variables + canopy openness + transect +
year, family = binomial(link = ‘logit’)).
We removed variables with a variance inflation factor > 5
(‘car’ package; Fox et al., 2023) to reduce parameter collinearity, and
restricted the models used in model selection to have N/10 maximum
parameters to avoid fitting overly complex models. Where there were
multiple GLMs within 2 AICc points of the best model, we calculated the
full model averaged coefficients (‘MuMIn’ package; Bartoń, 2022), and
otherwise we selected the model with the lowest AICc as the best model.
Because of the link function in our GLMs, our results can be interpreted
as increasing or decreasing the log odds of recruitment or seedling
survival, corresponding to lower or higher likelihoods, respectively.
Statistical analyses: microhabitat
suitability
To answer how microhabitat suitability changes at and beyond species’
range edges, we used the model results from the species establishment
analyses above to predict likelihood of recruitment, relative recruit
counts, and 1-year as well as 2-year seedling survival at each plot
(predict(model, type = ‘response’)) (‘MuMIn’ package; Bartoń, 2022). We
considered these predicted values as proxies for microhabitat
suitability, with increasing suitability at and beyond range edge
indicating greater likelihood of range expansion and decreasing
suitability indicating a lower likelihood for range expansion. For these
analyses, we used only the species where sites extend beyond their
thermal range limit (Fig. S2) and with an observed response in at least
8 plots.