DOI: 10.5061/dryad.rv15dv4gn
title: Variable species establishment in response to microhabitat indicates different likelihoods of climate-driven range shifts

description: Climate change is causing geographic range shifts globally, and
 understanding the factors that influence species’ range expansions is
 crucial for predicting future biodiversity changes. 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 effects on local population dynamics. We aim to better
 understand the role of microhabitat in range shifts in plants through its
 impacts on establishment by Q1) examining microhabitat variability along
 large macroclimatic (i.e., elevational) gradients, Q2) testing which of
 these microhabitat variables explain plant recruitment and seedling
 survival, and Q3) 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 five years, we recorded
 recruitment, survival, and microhabitat (i.e., high resolution soil, air,
 and light) 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 with
 elevation. We further observed extremely low recruitment and moderate
 seedling survival, 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-specific, as different species have different
 suitability and availability of microhabitat beyond their present ranges,
 thus calling into question low-resolution macroclimatic projections that
 will miss such complexities. We encourage further research on species
 responses to microhabitat and including microhabitat in range shift
 forecasts.
description: Study site We utilized the large
 macroclimatic and elevational 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). 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 graminoids. We selected 15 sites along each of the 2 transects
 (<em>n</em> = 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
 (<em>n</em> = 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. Within each block, we
 selected three replicate plots separated by ~2 m (<em>n</em> =
 180 total). Each plot contained three 50 cm x 50 cm quadrats ~ 20 cm from
 one another (<em>n</em> = 540 total). To ensure substrates had
 recruitment potential, we selected quadrats that had no saplings
 &gt;10 cm in diameter, no large rocks that covered &gt;10% of the
 quadrat, and ≥ 5% cover from seedlings or understory species. Of the three
 quadrats, we added seeds to two (separating confamilal pairs to separate
 plots to facilitate field identification of seedlings) and left the third
 quadrat unmanipulated to control for background recruitment. We grouped
 species by climate affinity (see next section) and sowed one quadrat per
 plot with seeds of species that generally occur in cooler temperatures and
 the other quadrat with seeds of species that generally occur in warmer
 temperatures. Species data We
 sowed 25 native species encompassing a variety of functional groups,
 climate affiliations, seed sizes, and regional prevalence to capture a
 variety of species responses in the community. Species included trees
 (<em>Abies grandis, A. lasiocarpa, Picea sitchensis, P. engelmannii,
 Pinus ponderosa, P. contorta</em>), shrubs (<em>Mahonia
 nervosa, M. aquifolium, Rubus ursinus, R. spectabilis, Sambucus cerulea,
 S. racemosa, Sorbus sitchensis, Vaccinium parvifolium, V.
 deliciosum</em>), forbs (<em>Eriophyllum lanatum, Anemone
 occidentalis, Erigeron perigrinus, Lupinus latifolius, Maianthemum
 dilatum, M. racemosum, Tolmiea menziesii, Tellima grandiflora</em>),
 and graminoids (<em>Carex stipata, Carex spectabilis</em>). 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. We also included unpaired species
 with large seeds (<em>L. latifolius</em>) or with high
 regional prevalence (<em>S. sitchensis, A. occidentalis</em>)
 to capture a variety of species in the community. We
 homogenized all seed sources for a given species with ~ 0.25 g of seed per
 species added to a quadrat. Therefore, more seeds of small seeded species
 were added than species with large seeds as per capita seedling emergence
 is inversely proportional to seed size (Bond et al. 1999). We sowed seeds
 in a sand mixture (~ 80 mL sand) scattered on quadrats with a lightly
 raked surface in September-October 2017. We added sand to aid seed
 dispersal within the quadrat and provide visualization for spreading
 seeds, rather than as a substrate for recruitment. We also added sand to
 control quadrats for consistency, and assume that the small amount of sand
 did not affect either soil or recruitment. Unfortunately, although seeds
 were collected regionally within Washington, we do not have information on
 the number of populations or individuals collected from. We also do not
 have information on collection and seed processing methods, although seeds
 of some species had obviously been cleaned. While seed provenance might
 cause different recruitment responses across microhabitat gradients,
 variation in germination rates did not significantly differ between seeds
 sourced locally versus from a nursery. In the following
 years, we tracked individual recruitment and seedling survival to measure
 success in early life history stages, and emphasize that success in early
 life stages does not necessarily correspond to successful transition to
 later life stages (e.g., Chang-Yang et al. 2021). We recorded recruitment
 (i.e., survival to end of first growing season) and yearly survival of
 each seedling during the growing season (May-September) in 2018, 2019, and
 2020, and recorded only surviving seedlings in 2022. We marked new
 recruits with colored toothpicks to track year of germination and years of
 survival. By recording new recruits over a 3-year period, we assume that
 we captured most delayed germination events, but note that seed
 germination timing is a vast area of research itself with life history
 traits and environmental conditions playing critical roles (e.g., Li et
 al. 1994, Donohue 2005, Soltani et al. 2015). 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 the sites in 2018 (wildfire closures) or any of the field
 sites in 2021 (closed USA border). Microhabitat
 data We measured a suite of abiotic and biotic
 parameters to quantify microhabitat and categorized these parameters as
 being directly influenced by climate or not to facilitate interpretation
 of the results and how this system may change in the future. We assume
 that microclimate (i.e., soil and above-ground moisture and temperature)
 is directly influenced by climate. We assume that soil microhabitat (i.e.,
 ratios of fungus to bacteria and carbon to nitrogen, water holding
 capacity) is indirectly influenced by climate, because it is partly driven
 by microbial activity, which might be affected by climate change. Unless
 distinguished, we refer to both as ‘microhabitat’ throughout. 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. Microhabitat data: climate loggers
 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
 placed one TMS-4 and one HOBO logger in a location representative of each
 block, for a total of 60 loggers. We calibrated volumetric water content
 (Kopecký 2021) in the TMS-4 sensors and cleaned the logged data (own
 functions and ‘myClim’ package from Man et al., 2023) for potential
 errors. We calculated seasonal variables from these data loggers, using
 ‘myClim’ functions 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. To measure winter (December, January,
 February) and summer (June, July, August) soil temperatures and calculate
 the duration of spring (March, April, May) snow cover, we deployed HOBO
 64K Pendant Temperature/Alarm data loggers (Onset, Bourne, Massachusetts
 USA) in each block (<em>n</em> = 2 loggers per site) from 2017
 through 2020. To measure summer soil temperature and moisture,
 soil-surface temperature, and plant-height temperature, we deployed
 Temperature-Moisture-Sensors (TMS-4) data loggers (TOMST, Prague, Czech
 Republic; Wild et al. 2019) in each block (<em>n</em> = 2
 loggers per site) from snow melt-out until late September 2022.
 We configured HOBO loggers to record temperature every 2 hours.
 We buried HOBO loggers approximately 1 cm below the soil surface (range of
 measurement: -20 to +70ºC, resolution: 0.14ºC, accuracy ± 0.53ºC). We
 configured TMS-4 loggers to record soil temperature and moisture every 15
 minutes. We installed TMS-4 loggers at approximate heights of -7.8, -0.2
 and +13.2 cm relative to the soil surface, corresponding to the three
 separate sensors on the logger (range of measurement: -60 to +85ºC,
 resolution: 0.0625ºC, accuracy ± 0.5ºC). We installed each logger so that
 the middle temperature sensor was covered with soil and affixed the
 provided sunshield to cover the upper sensor. We secured a wire cage
 around each logger to deter animal disturbance. We
 processed logger data by first removing data from sensors that had been
 compromised by animal disturbances, which we visually assessed when
 retrieving the loggers at the end of the field season. If a TOMST logger
 was disturbed, we only removed the sensors that was compromised (e.g.,
 entire logger removed out of the ground = all three sensors removed;
 sunshield removed from top sensor = top sensor removed). We only included
 the logger data in the species establishment analyses if loggers remained
 undisturbed at every block in which the species germinated. If a HOBO
 logger was disturbed, we removed that logger from the dataset entirely. We
 set the time series to both loggers corresponding to +1 day from a field
 deployment and -1 day from field retrieval. We processed the HOBO data by
 removing unlikely values (&lt; -20ºC or &gt; 60ºC) and 20ºC
 temperature spikes, and checked for alignment. We processed the TMS-4 data
 with the package ‘myClim’ (Man et al. 2023) as well as our own functions.
 Specifically, we cleaned data with ‘mc_prep_clean’, which removes
 duplicates, reorders misaligned measurements, and identifies missing
 measurements. We calibrated volumetric water content according to Kopecký
 et al. (2021), where volumetric Water Content = -0.0000000134*x^2 +
 0.000249622*x – 0.157889. We calculated all climate
 variables using the average of daily minimum or maximum of the season(s)
 corresponding to data logger deployment. We calculated summer variables
 corresponding to the approximate growing season at each site (from snow
 melt-out between May and July to plant senescence in late September), thus
 these variables encompass different temporal lengths.
 Microhabitat data: soil sampling 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 plot 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 <em>in situ</em> measurements 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. Upon each site snow melt-out in 2022, we took a 10
 cm soil core from a representative location between each set of replicates
 (<em>n</em> = 6 samples per site). We used a microBIOMETER®
 Test Kit (Prolific Earth Sciences, Montgomery, New York, USA) to conduct
 <em>in situ</em> measurements of live microbial biomass (ug
 biomass C/gram soil), fungal:bacterial ratio, and % fungal and % bacterial
 biomass. The microBIOMETER® assay separates live microbes from soil by
 suspending them in an extraction solution, which is then deposited onto a
 membrane surrounded by a grayscale color wheel on a provided test card and
 imaged with the microBIOMETER® mobile app. The color of the deposited
 sample is analyzed against the color wheel and used to calculate the
 different microbial parameters by virtue of each microbe type varying in
 their pigmentation and light absorption (https://microbiometer.com/). We stored remaining soil cores at ambient temperatures before transferring them to the University of Washington for analyses of carbon, nitrogen, and water holding capacity. We dried (~ 65ºC until mass decreased by less than 2% overnight), sieved, ground, and tested soil samples for inorganic carbon using 11.6 M HCl. We ran 30-50 mg of processed samples through a PerkinElmer® 2400 Series II CHNS/O Elemental Analyzer (2400 Series II) to quantify soil organic carbon to nitrogen (C:N) ratio. To quantify water holding capacity (WHC), we saturated sieved soils in cling film-covered filter funnels for 20 hours and calculated WHC based on values derived from dry weight versus wet weight (% WHC = [weight of water absorbed/weight of dry soil] * 100). Study design limitations Our experimental design is not without limitations. Since we sowed all seeds in the same year, we cannot test if we sowed in favorable versus unfavorable years, though we aimed to capture a large germination window by recording recruitment over three years. Additionally, we collected some of the data on species responses and microhabitat measurements asynchronously and therefore emphasize that only relative differences between sites, not absolute values, should be interpreted for their effects. Statistical analyses: microhabitat variability (Q1) To answer how microhabitat varies with elevation, a robust proxy for macroclimate in this system, we fit Linear Mixed Models (LMMs; package ‘lmerTest’; Kuznetsova et al. 2017) to test the effects of elevation (linear and quadratic terms), transect, and their interaction on each microhabitat parameter, with a random intercept of site. We included a quadratic term for elevation to account for potential optima in microhabitat variables within the elevational sampling gradient. To further assess microhabitat variability, we tested how variation within sites compares to among sites of each of the microhabitat parameters. To do this, we calculated the variance partition coefficient from intercept-only LMMs that included block nested within site as a random intercept. In cases where LMMs did not converge due to insufficient, or no, variability in the random intercept, we changed the random intercept to include only site. We also used a principal component analysis plot to identify if microhabitat variables directly or indirectly affected by climate cluster together. We conducted all data processing and analyses in R version 4.2.3 (R Core Team 2023). Statistical analyses: species establishment (Q2) To interpret the results of microhabitat associations to likelihood of recruitment, recruit numbers, and seedling survival, we first transformed 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<sub>species i, site-rep j </sub>/ maximum recruit counts<sub>species i</sub>. We used all non-zero recruitment data to calculate the proportion of surviving seedlings in all replicates, including control replicates, for 1-year, 2-year, and multi-year proportion of surviving seedlings. Sampling error (i.e., when we overlooked first-year seedlings but tagged them as second-year seedlings) yielded a calculated survival probability &gt; 1 in 0.14% of the data, and we set these values to 1. While we have general <em>a priori </em>hypotheses of how certain microhabitat conditions affect species establishment, 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 only analyzed species for which we observed a response in at least 8/180 plots for any given life stage to limit this study to species with sufficient recruitment to test for microhabitat effects on establishment. We fit separate binomial Generalized Linear Models (GLMs) for each species’ life stage (recruitment and 1- to 2- year seedling survival) to determine which of the microhabitat variables described above are important for establishment (Q2). While we do not consider 1- and 2- year seedling survival as separate life stages, we analyze them separately due to declining sample sizes over time. Each global GLM included all 8 microhabitat variables described above, plus any quadratic effects identified as important by AICc in a reduced model (‘AICcmodavg’ package; Mazerolle 2020). In fitting the global GLMs, we removed variables with a variance inflation factor &gt; 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-averaged coefficients (i.e., zero values entered for coefficients for terms not included in a particular model; ‘MuMIn’ package; Bartoń 2022), and otherwise we selected the model with the lowest AICc as the best model. Because of the link function in the GLMs, the results can be interpreted as increasing or decreasing the log odds of recruitment or seedling survival, corresponding to lower or higher likelihoods, respectively. 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 and interpret their effects on establishment, we further included canopy openness and transect as fixed effects. In 1-year survival models, we also included year as a fixed effect to test for differences among years in overall seedling survival and different frequencies of site visits. Almost all of the models testing the effects of microhabitat on species establishment met the GLM assumption of independence of residuals. However, many of the models did not meet other assumptions, including homogeneity of variance or a variance that equals the mean, or the link function was sometimes inappropriate (&gt; 0.5 difference from 1 for slope of link function). For these reasons, we are cautious in interpreting the results and do not interpret the exact quantitative values here, especially because candidate models generally explained a low proportion of variance. However, we did not find an alternative, more suitable approach, and believe that these models still provide important insights into the broad impacts of microhabitat on establishment. Statistical analyses: microhabitat suitability (Q3) To estimate which sites are located within or beyond each species elevational distribution, we calculated macroclimatic range position for each site per species. To start, we obtained scale-free downscaled mean annual temperature (MAT) from 1981-2020 gridded monthly climate normal data (800 x 800 m) (https://climatena.ca/). To calculate thermal range position, we first defined the lower and upper thermal ranges as the 2.5% and 97.5% quantiles, respectively, of the distribution of MAT occurrences for each species based on its 1980-2022 GBIF occurrences (https://www.gbif.org/) in the Pacific Northwest (North American temperate conifer forest biome; https://ecoregions.appspot.com/). Second, we calculated where each site fell within each species’ thermal range, standardizing this for each species: (site MAT - species range midpoint) / (upper species limit - lower species limit) * 0.5. The calculated range positions are relative to each species’ cold limit (-1), range center (0), and warm limit (1). Values &lt; -1 indicate that transplant sites were beyond the species current (1980 - 2022) observed upper elevational distribution limit. To answer how microhabitat suitability changes at and beyond species’ range edges (Q3), 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 suggesting that microhabitat could facilitate range expansions, in the sense that species establishment will likely occur in these suitable microhabitats for an overall higher likelihood of range expansion. Conversely, decreasing suitability suggests that microhabitat constrains range expansion for a lower likelihood for range expansion. For these analyses, we used only the species where sites extend beyond their regional upper elevational range limit (Supporting Information) and with an observed response in at least 8 plots. Out of the 25 species in this study, this yielded recruitment and recruit count suitability predictions for 8 species, 1-year seedling survival suitability predictions for 5 species, and 2-year seedling survival suitability predictions for 4 species. References Bartoń, Kamil. 2022. “MuMIn: Multi-Model Inference.” https://cran.r-project.org/web/packages/MuMIn/index.html. Bond, W. J., M. Honig, and K. E. Maze. 1999. “Seed Size and Seedling Emergence:  An Allometric Relationship and Some Ecological Implications.” <em>Oecologia</em> 120 (1): 132–36. https://doi.org/10.1007/s004420050841. Chang-Yang, Chia-Hao, Jessica Needham, Chia-Ling Lu, Chang-Fu Hsieh, I-Fang Sun, and Sean M. McMahon. 2021. “Closing the Life Cycle of Forest Trees: The Difficult Dynamics of Seedling-to-Sapling Transitions in a Subtropical Rainforest.” <em>Journal of Ecology</em> 109 (7): 2705–16. https://doi.org/10.1111/1365-2745.13677. Donohue, Kathleen. 2005. “Seeds and Seasons: Interpreting Germination Timing in the Field.” <em>Seed Science Research</em> 15 (3): 175–87. https://doi.org/10.1079/SSR2005208. Fox, John, Sanford Weisberg, Brad Price, Daniel Adler, Douglas Bates, Gabriel Baud-Bovy, Ben Bolker, et al. 2023. “Car: Companion to Applied Regression.” https://cran.r-project.org/web/packages/car/index.html. Kopecký, Martin, Martin Macek, and Jan Wild. 2021. “Topographic Wetness Index Calculation Guidelines Based on Measured Soil Moisture and Plant Species Composition.” <em>Science of The Total Environment</em> 757 (February): 143785. https://doi.org/10.1016/j.scitotenv.2020.143785. Kuznetsova A, Brockhoff PB, Christensen RHB. 2017. lmerTest Package: Tests in Linear Mixed Effects Models. <em>Journal of Statistical Software </em>82(13): 1-26. doi:10.18637/jss.v082.i13. Li, X. J., P. J. Burton, and C. L. Leadem. 1994. “Interactive Effects of Light and Stratification on the Germination of Some British Columbia Conifers.” <em>Canadian Journal of Botany</em> 72 (11): 1635–46. https://doi.org/10.1139/b94-201. Man, Matěj, Vojtěch Kalčík, Martin Macek, Josef Brůna, Lucia Hederová, Jan Wild, and Martin Kopecký. 2023. “myClim: Microclimate Data Handling and Standardised Analyses in R.” <em>Methods in Ecology and Evolution</em> 14 (9): 2308–20. https://doi.org/10.1111/2041-210X.14192. Mazerolle, Marc J. 2020. “AICcmodavg: Model Selection and Multimodel Inference Based on (Q)AIC(c).” https://cran.r-project.org/web/packages/AICcmodavg/index.html. Tredennick, Andrew T., Giles Hooker, Stephen P. Ellner, and Peter B. Adler. 2021. “A Practical Guide to Selecting Models for Exploration, Inference, and Prediction in Ecology.” <em>Ecology</em> 102 (6): e03336. https://doi.org/10.1002/ecy.3336.
description: # Variable species establishment in response to microhabitat indicates
 different likelihoods of climate-driven range shifts
 [https://doi.org/10.5061/dryad.rv15dv4gn](https://doi.org/10.5061/dryad.rv15dv4gn) ## Description of the data and file structure Made available here are compiled data on yearly germination/seedling census and microhabitat variables, and R scripts used to analyze the effects of microhabitat on germination and seedling survival. These data and R scripts were used in Chardon et al. 2024. Variable species establishment in response to microhabitat indicates different likelihoods of climate-driven range shifts. Ecography. Files are listed under their parent folder. **Compiled Data: these data files are used in the R scripts below, which load these data files as listed** data/abbreviations.xlsx: description of all variables and species abbreviations, with datasets each variable is (1) or is not (0) used in data/census_sdlg.RData: yearly seedling survival + microhabitat variables data/census_spp.RData: yearly census data + microhabitat variables at quadrat level (3 quadrats per rep) data/census_rep.RData: yearly census data + microhabitat variables at rep level (3 reps per block, with 2 block per site) data/thermal2elevRL.RData: thermal to elevational range limit per species data/seedmix_added: weight or number of seeds added to seedmix per species **Analyses: to be run in numeric order (uploaded to Zenodo)** scripts/01_data_expl.R: choose best random effects structure, variance partitioning, choose fixed effects scripts/02_param_corr.R: choose biologically meaningful explanatory parameters for models scripts/03_varsel_recruitprob_bin.R: Binomial GLM for recruitment probability scripts/04_varsel_recruit_bin.R: Binomial GLM for relative number of recruits scripts/05_varsel_sdlg_bin.R: Binomial GLM for 1- and 2-year seedling survival scripts/06_therm2elevRL.R: calculate elevational range limit corresponding to thermal RL scripts/07_microhab_suitab.R: predict binary/continuous recruitment and 1/2 year survival to each plot **Results: scripts used to generate figures and tables for manuscript (uploaded to Zenodo)** scripts/ms_figs.R: figures for manuscript scripts/ms_tables.R: tables for manuscript **Package versioning** The folder renv/ contains the packages versions used to run the above scripts. To use, download the entire repository and run renv::update(). ## Sharing/Access information Available as GitHub repository: [https://github.com/nchardon/spp_estab_microhab/](https://github.com/nchardon/spp_estab_microhab/) ## Code/Software We conducted all data processing and analyses in R version 4.2.3 (R Core Team 2023). ## References R Core Team, R. 2023. “R: A Language and Environment for Statistical Computing.” R Foundation for Statistical Computing, Vienna, Austria. [http://www.r-project.org.](http://www.r-project.org.)