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.