Description
| Abstract
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. (2024-05-24)
Methods
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 (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. Within each block, we selected three replicate plots separated by ~2 m (n = 180 total). Each plot contained three 50 cm x 50 cm quadrats ~ 20 cm from one another (n = 540 total). To ensure substrates had recruitment potential, we selected quadrats that had no saplings >10 cm in diameter, no large rocks that covered >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 (Abies grandis, A. lasiocarpa, Picea sitchensis, P. engelmannii, Pinus ponderosa, P. contorta), shrubs (Mahonia nervosa, M. aquifolium, Rubus ursinus, R. spectabilis, Sambucus cerulea, S. racemosa, Sorbus sitchensis, Vaccinium parvifolium, V. deliciosum), forbs (Eriophyllum lanatum, Anemone occidentalis, Erigeron perigrinus, Lupinus latifolius, Maianthemum dilatum, M. racemosum, Tolmiea menziesii, Tellima grandiflora), and graminoids (Carex stipata, Carex spectabilis). 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 (L. latifolius) or with high regional prevalence (S. sitchensis, A. occidentalis) 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 (n = 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 (n = 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 (< -20ºC or > 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 in situ 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 (n = 6 samples per site). We used a microBIOMETER® Test Kit (Prolific Earth Sciences, Montgomery, New York, USA) to conduct in situ 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 countsspecies i, site-rep j / maximum recruit countsspecies i. 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 > 1 in 0.14% of the data, and we set these values to 1.
While we have general a priori 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 > 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 (> 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 < -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.” Oecologia 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.” Journal of Ecology 109 (7): 2705–16. https://doi.org/10.1111/1365-2745.13677.
Donohue, Kathleen. 2005. “Seeds and Seasons: Interpreting Germination Timing in the Field.” Seed Science Research 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.” Science of The Total Environment 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. Journal of Statistical Software 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.” Canadian Journal of Botany 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.” Methods in Ecology and Evolution 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.” Ecology 102 (6): e03336. https://doi.org/10.1002/ecy.3336. (2024-05-24) |