Variable species establishment in response to microhabitat indicates diﬀerent likelihoods of climate-driven range shifts

Climate change is causing geographic range shifts globally, and understanding the factors that inﬂuence species’ range expansions is crucial for predicting future changes in biodiversity. A common, yet untested, assumption in forecasting approaches is that species will shift beyond current range edges into new habitats as they become macroclimatically suitable, even though microhabitat variability could have overriding eﬀects on local population dynamics. We aim to better understand the role of microhabitat in range shifts through its impacts on establishment by i) examining microhabitat variability along large macroclimatic gradients, ii) testing which of these microhabitat variables explain plant recruitment and seedling survival, and iii) predicting microhabitat suitability beyond species range limits. We transplanted seeds of 25 common tree, shrub, forb, and graminoid species across and beyond their current elevational ranges in the Washington Cascade Range, USA, along a large elevational gradient spanning a broad range of macroclimates. Over ﬁve years, we recorded recruitment, survival, and microhabitat characteristics rarely measured in biogeographic studies. We asked whether microhabitat variables correlate with elevation, which variables drive species establishment, and whether microhabitat variables important for establishment are already suitable beyond leading range limits. We found that only 30% of microhabitat parameters covaried in the expected way with elevation. We further observed extremely low recruitment and moderate seedling survival in our study system, and these were generally only weakly explained by microhabitat. Moreover, species and life stages responded in contrasting ways to soil biota, soil moisture, temperature, and snow duration. Microhabitat suitability predictions suggest that distribution shifts are likely to be species-speciﬁc, as diﬀerent species have diﬀerent suitabilities, and availabilities, of microhabitat beyond their present ranges, thus calling into question large-scale macroclimatic projections that will miss such complexities. We encourage further research on species responses to microhabitat and the inclusion of microhabitat in range shift forecasts.


Introduction
The most pressing environmental issues of our times are understanding the ecological effects of ongoing climate change and predicting the ensuing implications for maintaining biodiversity (Lovejoy & Hannah, 2019).Foundational to these issues are an understanding of species range limits, which are the geographic limits to species' spatial distributions.Projection inaccuracy of how species ranges will shift with climate change is alarming (Urban, 2019), and it is unclear whether our inability to project species range shifts comes from our poor knowledge of which climatic variables are important for which species or from the many other factors that can influence range shifts (e.g.microhabitat variation, dispersal limitation, species interactions, and many others).Even explaining the mechanisms underlying observed shifts is difficult, and it is clear that responses are highly species-specific (Freeman et al., 2018;Rumpf et al., 2018).
To address this issue, ecologists are urgently seeking to understand species niche limits in order to project species distribution shifts with climate change.These are commonly predicted with correlative Species Distribution Models (SDMs; Elith et al., 2010;Pacifici et al., 2015), which correlate species occurrences to current macro-climate in order to predict the relative probability of occurrence in space and time (Wiens et al., 2009).Such models have become an integral part of biogeographical research and progress has been made to address some of their shortcomings (Franklin, 2023).While these models implicitly assume that they include sufficient niche variables to accurately describe a species distribution, few models actually test what niche components are necessary to quantify population growth rate and shape species distributions (Pulliam, 2000).
A main driver of population growth rate is the fine-scale microhabitat that mediates the success of individual plants through abiotic (e.g., canopy cover, above-ground temperature, soil moisture) and biotic (e.g., soil biota, soil nutrients) conditions (Kephart & Paladino, 1997;Zurbriggen et al., 2013;Oldfather & Ackerly, 2019;Tanner et al., 2021;Sanczuk et al., 2023;Allsup et al., 2023;Kemppinen et al., 2023).However, studies that quantify population response to microhabitat variation along species ranges are rare, impeding an understanding of which aspects of fine-scale microhabitat are important in shaping species broad-scale distributions (but see Tourville et al., 2022).Despite this, recent work has still shown that incorporating microhabitat into SDMs can provide substantially improved predictions (Maclean & Early, 2023).
How microhabitat and microclimate themselves vary along large elevational, and therefore macroclimatic, gradients that characterize species ranges is poorly understood, and research is increasingly finding that microclimate is often decoupled from macroclimate (Scherrer & Körner, 2010;Ford et al., 2013;Lembrechts, 2023).It is, however, now widely accepted that microclimate differs from macroclimate due to factors such as forest canopy cover (De Frenne et al., 2019;Haesen et al., 2021) and topography (Lawson et al., 2014).This can modify species response to habitat use (Lawson et al., 2014) and create climate change refugia (Dobrowski, 2011;Pradhan et al., 2023) that mediate species response to climate change (Morelli et al., 2020).Very few studies have been able to document how variable plant-relevant environmental drivers of performance are across macroscale gradients, although substantial progress is underway in including microclimate as part of biogeographical research (Kemppinen et al., 2023).
Testing how sensitive species' life stages are to variations in microhabitat can address many of the sources of bias that distribution models face and accounting for demographic processes can greatly improve our understanding of range dynamics (Fig. 1; Normand et al., 2014;Copenhaver-Parry et al., 2020).For plant species, how early life stages (i.e.recruitment and seedling survival; henceforth establishment) respond to microhabitat may be the key to determining species' ability to expand their ranges (Kroiss & HilleRisLambers, 2015), as establishment is essential for a range expansion to occur.Understanding how early life stages respond to variations in microhabitat is thus critical to understanding microhabitat suitability for range shifts and population persistence, especially since life stages can respond differently to environmental conditions (Goodwin & Brown, 2023).Yet, few studies test how microhabitat influences establishment, and important microhabitat variables, such as soil moisture and temperature, are often left out of such work (Schurr et al., 2012).
To comprehensively understand how microhabitat, and other factors, influence species distributions, transplant experiments beyond species' ranges have been suggested as a promising approach (Lee-Yaw et al., 2016;Morris & Ehrlén, 2015).Transplant experiments show either congruence with SDM predictions (Sanczuk et al., 2022) or markedly different responses of populations than predicted by SDMs (Greiser et al., 2020), highlighting the importance of these experiments.Such experiments are particularly valuable at the leading edge of a species range (i.e.edge expanding with climate change), where novel species interactions are likely to be found (Thuiller et al., 2008).Transplant experiments often find that species are establishment, not dispersal, limited, and this can be due to unfavorable microsite conditions (Clark et al., 2007;Davis & Gedalof, 2018).How microhabitat suitability is distributed beyond species current leading edge therefore likely determines species' ability to shift their ranges (Tourville et al., 2022), yet few studies quantify microhabitat along and above species distributions.We utilized the large elevational macroclimatic gradients of the West and East sides of the Washington Cascade Range, USA for a seed transplant experiment encompassing common grasses, forbs, shrubs, and trees to ask: 1. How does microhabitat variation within sites compare to among-site macroclimate variation, and does microhabitat covary expectedly along elevational gradients?2. Does microhabitat explain establishment of common species in our system, and if so, are responses consistent among species?3.If microhabitat variables explain establishment, does the distribution of microhabitats along macroclimatic gradients suggest range shifts will be facilitated, constrained, or unaffected by the availability of suitable microhabitat?

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 counts species i, site-rep j / maximum recruit counts species 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')) ('Mu-MIn' 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.

Microhabitat variability
Microhabitat variation within sites was higher than among site variation for half of the microhabitat parameters for which we could fit nested LMMs.Out of these 9 microhabitat parameters that we measured, only 3 varied predictably with elevation (canopy openness, winter minimum soil temperature, and spring days of snow).The parameters that followed an elevational pattern are both expected to be directly (winter minimum temperature, spring days of snow) as well as indirectly (canopy openness) affected by climate (Fig. 2).
Our principal component analysis plot shows that the microhabitat variables we measured cluster broadly by soil characteristics (carbon:nitrogen, fungus:bacteria, water holding capacity), soil moisture (summer minimum soil moisture, spring days of snow), light availability (canopy openness), soil temperature (summer maximum soil temperature), and above-ground temperature (summer minimum plant-height temperature) (Fig. S3).An additional parameter describing soil temperature (winter minimum soil temperature) clusters in the middle of the microhabitat space.

Species establishment
Out of the 25 species that we sowed, 84% (21/25) recruited in at least 1 plot and 56% (14/25) recruited in 8 plots or more (Table S2).Of the species that recruited, seedlings (i.e. at least one seedling) survived for one year for 57% (8/14) of species and survived for two years for 43% (6/14) of species.Almost all species that had sites beyond their range edge recruited and survived beyond their range edge, a pattern seen for species with either high-or low-elevation range edges (Table S2).Almost all of our models met the GLM assumption of independence of residuals.However, many of our models did not meet homogeneity of variance, variance did not equal the mean, or the link function was sometimes inappropriate (> 0.5 difference from 1 for slope of link function).Together with the low proportion of variance (Table S3) from each set of candidate models, we are therefore cautious in interpreting our results.
Overall, we found microhabitat parameters both related and not related to climate identified as the most important parameters in model selection of the effects of microhabitat on establishment, but models had low explanatory power (Tables 1, S3).At the community level, the main patterns we found were that certain microhabitat variables had largely negative or positive effects, but the direction of these effects changed with different life stages (Fig. 3).At the species level, we found that some species only had few microhabitat parameters chosen in model selection, whereas others had many (Fig. S4).Within the same species, parameters usually had the same directionality of effect on likelihood of recruitment, recruit counts, and seedling survival except for T. grandiflora.For two species, A. lasiocarpa andE.lanatum , the directionality of effects was consistent across life stages.

Microhabitat suitability
We generated predictions of microhabitat suitability to better understand if the distribution of microhabitats along climatic gradients suggests range shifts will likely be facilitated (i.e.where suitability increases with elevation), constrained (i.e.decreases), or unaffected (i.e. is constant).We note that due to poor model fit in the original models, our predictions cannot be used to make definitive forecasts of which species are likely to shift their ranges with climate change, but rather use these predictions to assess how microhabitat suitability changes across a large elevational gradient.Our predictions indicate that only range shifts for L. latifolius are likely to be facilitated with increasing microhabitat suitability with elevation (Fig. 4), although our experimental sites extended only just beyond the range limit for the species.For all other species, microhabitat suitability either declines or has no pattern with elevation and thus range shifts are likely to either be constrained or unaffected, respectively (Figs.S5, S6, S7).We also found that suitability patterns were modified by transect and life stage.
Together with these suitability predictions, the microhabitat parameters that significantly vary with elevation can give further insight on what aspects of microhabitat can facilitate or constrain range shifts.For example, the increasing microhabitat suitability for L. latifolius (Fig. 4a) together with the positive effects of canopy openness on the species (Table 1) and canopy openness increasing with elevation (Fig. 2a) point to a likely facilitated range expansion for the species.Spring snow cover shows the same pattern, and thus is a further microhabitat variable that will likely facilitate range expansion for this species.However, winter minimum soil temperature decreases with elevation in the Okanagan transect but positively affectsL.latifolius recruit counts and seedling survival, so this microhabitat variable could constrain range expansion.

Discussion
While many species ranges are on average shifting in the direction predicted by SDMs, there is large variation among species that is difficult to explain.Microhabitat variability is increasingly recognized as influencing population dynamics across species ranges (Oldfather & Ackerly, 2019) and may be a key factor in understanding and predicting range shifts (Lembrechts, Nijs, et al., 2019;Maclean & Early, 2023;Stickley & Fraterrigo, 2023).However, how microhabitat suitability, and variability therein, is distributed along species ranges is vastly understudied, impeding an understanding of how microhabitat can facilitate or constrain range shifts.We found that microhabitat is highly variable across species ranges and is largely decoupled from the macro-scale at which SDMs are typically constructed.We further found variable, albeit weak, effects of microhabitat for different species and life stages, suggesting that the drivers of establishment are complex and difficult to detect.Our microhabitat suitability predictions show microhabitat will either constrain range expansions, or have no effect, for most species, with the range shifts of just one species (L.latifolius ) likely facilitated by microhabitat.

Complex microhabitat patterns across elevation gradients
Surprisingly, most of our microhabitat variables, even ones related to climate (summer soil moisture, summer soil temperature, summer plant height temperature) did not follow an elevational pattern representative of the macroclimate (Fig. S8) commonly used in SDMs.This is consistent with many studies finding that microhabitat is often decoupled from elevation and regional climate (Ford et al., 2013;Lembrechts, Lenoir, et al., 2019).We further found that none of the soil composition parameters we measured (fungus:bacteria, carbon:nitrogen, water holding capacity) followed an elevational pattern, even though other alpine studies have found soil microbial community (Hiiesalu et al., 2023) and soil carbon:nitrogen (Weintraub et al., 2016) to covary with elevation.
We further found that only canopy openness, winter minimum soil temperature, and spring snow days covary with elevation.Despite higher elevation areas characterized by more open forest canopies, the high variability at each site (Fig. 2) is in line with the large body of work showing that forest canopies can provide climate refugia (De Frenne et al., 2019;Haesen et al., 2021).As species light requirements can be variable, light availability and canopy cover can shape species ranges (Muñoz Mazon et al., 2023) and therefore can be an important component of range shifts (Tourville et al., 2022).

Low and variable effects of microhabitat on establishment
Our results agree with other studies, which find overall low recruitment in seed addition sites (Clark et al., 2007) and decreasing fitness of transplants beyond the range (Lee-Yaw et al., 2016;Stanton-Geddes et al., 2012).Plant establishment beyond current range edges is necessary for plant species to shift their ranges upward, however this is complicated by plants usually being recruitment limited.This means that they are limited by some aspect of their environment (e.g., microhabitat unsuitability, seed predation) and not by dispersal (Clark et al., 2007).While we did not control for seed predation, herbivory, or fungal infections, we captured mortality with our methods.Even though seed predation rates decline in higher elevation areas (Hargreaves et al., 2019), we still removed fleshy fruits around seeds to deter predation.We also did not collect data to test for effects of community composition on establishment, but with low recruitment rates, plants in our system would still be considered recruitment limited.
We did not find that establishment responses to any particular microhabitat were restricted to certain growth forms, with species-specific relationships to microhabitat.This is not surprising, as species have unique environmental requirements (Table S1) and species characteristics even impact predictive power in SDMs (Guisan et al., 2007).We also did not find differences in effects of the microhabitat aspects that are expected to be directly versus indirectly affected by climate change.Interestingly, we found that some parameters, such as winter soil temperature, increased likelihood of recruitment but these effects switched to decreasing likelihood of seedling survival.One exception was spring days with snow, which generally had positive effects, and this is similar to findings by Davis & Gedalof (2018), who found positive effects of winter snow cover on recruitment.The importance of microbial community in plant dynamics is also increasingly acknowledged (Castro et al., 2022), although soil fungus:bacteria ratio was not selected more than other variables in our system.Microbial communities could even be an important factor in determining species range shifts, such as favorable soil microbial composition mediating climate tolerance to promote tree seedling survival (Allsup et al., 2023).

Predicting microhabitat suitability
We found that microhabitat suitability beyond the leading range edge was species-and life stage-specific, with most species showing either decreased suitability or no pattern with elevation.We posit that the role that microhabitat may play in facilitating or constraining range shifts for any given species is closely tied to how microhabitat itself varies with elevation.For example, if a favorable microhabitat parameter for establishment decreases with elevation, establishment will overall not be favored and a range shift could be constrained.We found this in our system, with increasing likelihood of L. latifoliusrecruit counts and seedling survival with warmer winter minimum soil temperatures, but soil temperature decreases with elevation.If a favorable microhabitat parameter increases with elevation, however, this could lead to a facilitated range shift.This is the case for canopy openness, which increases with elevation and also increases likelihood of L. latifolius establishment.Since microhabitat suitability increases beyond the range for this species, these cases highlight the complex ways in which microhabitat parameters may interact to mediate range shifts.The alternative scenarios, where an unfavorable microhabitat parameter either declines with elevation to inversely favor establishment or increases with elevation to inhibit establishment, can also facilitate or inhibit range shifts, respectively.
Macroclimatic parameters vary predictably with elevation (Fig. S8) in our study system, however only a third of our microhabitat variables show an elevational pattern.Since species distributions can be driven by microhabitat, understanding how microhabitat parameters vary across species' ranges can give insights into how microhabitat may facilitate or inhibit range expansion (Lembrechts et al., 2017;Tourville et al., 2022;Kemppinen et al., 2023).Incorporating microhabitat suitability into ecological studies is not new, including quantifying habitat suitability in the field (Jabis & Ayers, 2014), with remote sensing (Falco et al., 2019), and at the leading range edge (Mamet & Kershaw, 2013).Transplant studies often find decreasing macroclimate suitability beyond the range (Lee-Yaw et al., 2016), but to our knowledge no studies exist that assess microhabitat suitability beyond species range edges.Despite increasing work showing the necessity of including microhabitats in SDMs (Lembrechts, Nijs, et al., 2019), the vast majority of SDMs still focus on macroclimatic variables and this may lead to the mismatches between predicted and observed range shifts, or lack thereof.
In our system, we found that suitable microhabitat is unchanged or reduced beyond the leading edge for many species, yet we also find that macroclimatic conditions enable an upwards shift in species' recruitment optima at the community level (unpublished data ).However, this pattern is only evident at the community level, with high species-level variability and no effect of canopy cover.This matches the large variability, and low explanatory power, seen in our microhabitat suitability predictions, and suggests that species can find pockets of suitable microhabitat to recruit beyond their macroclimatic cold range edge.Microhabitat suitability may be just high enough to allow for successful establishment, but is overall very low and the complicated ways in which microhabitat acts to affect species establishment makes it difficult to detect patterns.We also found continued species recruitment at the warm edge of the range (unpublished data ) and this buffering of contractions at the warm edge could be due to microhabitat refugia (De Lombaerde et al., 2022) created by the range of microhabitat found at the lower elevation sites (Fig. 2)

Limitations
Our experimental design is not without limitations.Seed provenance might cause different recruitment responses across microhabitat gradients, however variation in germination rates was not biologically meaningful (Fig. S1).Since we sowed all seeds in the same year, we also cannot test if we sowed in favorable versus unfavorable years and aimed to capture a large germination window by recording recruitment over three years.Finally, we collected some of our data on species responses and microhabitat measurements asynchronously and therefore emphasize that only differences between sites, not absolute values, should be interpreted for their effects.Finally, we strongly urge future studies to incorporate increased replication at each microhabitat value for higher statistical power in explaining results.

Conclusion
To our knowledge, our work is one of the first studies to measure key demographic life stages of a large group of species along a microhabitat gradient both within and beyond their current range limit.As such, this work yields a more comprehensive understanding of the mechanisms that set species range limits and indicates ways in which species shift their distributions with climate change.Our results suggest that complex ways in which microhabitat parameters influence early life stages are complex and difficult to detect, which complicates range shift predictions.A greater understanding of the role of macrohabitat in shaping species' distribution limits will ultimately improve predictions of how species distributions will shift with climate change.We emphasize that predictions accounting for complex microhabitat drivers are necessary to create tailored conservation and management decisions in order to mitigate ongoing biodiversity loss.We show within (σ 2 within ) and among (σ 2 among ) site percent variation explained.These variance partitioning results are from Linear Mixed Models (LMMs) of the microhabitat parameter on the y-axis and block nested within site (microhabitat parameter ˜1 + (1|site/block).In cases of model fitting issues, we ran the LMM with only (1|site) as a random effect and thus only report among-site percent variation explained.We also fit LMMs to test the effects of quadratic elevation, transect, and their interaction on the microhabitat parameter indicated on the y-axis, with a random effect of site (microhabitat parameter ˜poly(elevation, 2) * transect + (1|site)).We show one fitted line for models for significant (alpha = 0.05) effects of elevation (line is shown for MB transect only), and two fitted lines for significant elevation*transect interaction even if elevation itself was not significant.Note that fitted lines do not include random effects.Legend for all plots is as in (a) with different colors indicating the warmer, wetter West (MB) transect or the cooler, drier East (RP) transect.Microhabitat suitability predicted for each plot across the elevation gradient of our study system for three representative species with sites beyond their current thermal range limit (below range limit = solid points; beyond = hollow points; see Fig. S2) varies by species and life stage.Predictions are generated from model averaging results testing the effects of microhabitat on likelihood of recruitment (a-c), relative recruit counts (d-f), and 1-year seedling survival (g-i).Note that these predictions are not based on elevation, and that y-axis scales differ.Points are colored by the warmer, wetter West (green) or cooler, drier East (black) transects and size of points corresponds to observed recruitment (a-c), relative recruit counts (d-f), or proportion of seedlings that survived to one year (g-i) in a plot.Lighter colors for both transects indicate plots without recruitment and thus this size does not reflect an observed value.Full predictions are shown in Fig. S7.
Tables Table 1.Model averaging results (or best model, when no other models were within Δ AIC < 2, indicated with '*') for assessing importance of different microhabitat variables in explaining likelihood of recruitment (i.e. did seeds germinate or not), relative recruit counts (i.e.how many seedlings recruited relative to maximum recruit counts for the species), and 1-or 2-year seedling survival.'+' indicates a positive parameter estimate and '" indicates a negative parameter estimate.Columns are colored by the expectation of those microhabitat variables being directly (blue) or indirectly (green) affected by climate change.Blank cells indicate that the parameter was not chosen in model selection.Because many of our soil loggers were compromised by animal disturbance, we only used those variables for species that showed no response at

Figure 2 .
Figure2.Both microhabitat parameters that are not expected to be directly affected by climate change (a-d) or expected to directly respond to climate change (e-i) differed in their within-site variability and show variable relationships with elevation.We show within (σ 2 within ) and among (σ 2 among ) site percent variation explained.These variance partitioning results are from Linear Mixed Models (LMMs) of the microhabitat parameter on the y-axis and block nested within site (microhabitat parameter ˜1 + (1|site/block).In cases of model fitting issues, we ran the LMM with only (1|site) as a random effect and thus only report among-site percent variation explained.We also fit LMMs to test the effects of quadratic elevation, transect, and their interaction on the microhabitat parameter indicated on the y-axis, with a random effect of site (microhabitat parameter ˜poly(elevation, 2) * transect + (1|site)).We show one fitted line for models for significant (alpha = 0.05) effects of elevation (line is shown for MB transect only), and two fitted lines for significant elevation*transect interaction even if elevation itself was not significant.Note that fitted lines do not include random effects.Legend for all plots is as in (a) with different colors indicating the warmer, wetter West (MB) transect or the cooler, drier East (RP) transect.

Figure 3 .
Figure 3. Microhabitat parameters identified as most important vary for likelihood of recruitment (a), relative recruit counts (b), 1-year seedling survival (c), and 2-year seedling survival (d).Percentages calculated from the total times a parameter was used in model selection for all species, with pink indicating a negative effect and green indicating a positive effect.Legend in all panels is as in (b).Note that each panel corresponds to a different number of models (see Tables 1, S3).Summer temperature and minimum soil moisture were only used in model selection for two species because of missing data at many sites.Parameter abbreviations are C:N -soil carbon:nitrogen; Canopy -canopy openness; F:B -soil fungus:bacteria; Region -transect; S Soil Moist -summer minimum soil moisture; S Temp -summer maximum soil temperature for (a), (b) and summer minimum plant-height temperature for (c), (d); Snow -spring days of snow cover; W Soil Temp -winter minimum soil temperature; WHC -water holding capacity, with 'ˆ2' indicating a quadratic effect.

Figure 4 .
Figure4.Microhabitat suitability predicted for each plot across the elevation gradient of our study system for three representative species with sites beyond their current thermal range limit (below range limit = solid points; beyond = hollow points; see Fig.S2) varies by species and life stage.Predictions are generated from model averaging results testing the effects of microhabitat on likelihood of recruitment (a-c), relative recruit counts (d-f), and 1-year seedling survival (g-i).Note that these predictions are not based on elevation, and that y-axis scales differ.Points are colored by the warmer, wetter West (green) or cooler, drier East (black) transects and size of points corresponds to observed recruitment (a-c), relative recruit counts (d-f), or proportion of seedlings that survived to one year (g-i) in a plot.Lighter colors for both transects indicate plots without recruitment and thus this size does not reflect an observed value.Full predictions are shown in Fig.S7.