• Energy Research
• US close
• CN close
• ES close
• DRYAD close

# Land use change converts temperate dryland landscape into a net methane source Raw flux data for methane, carbon dioxide, and other species were measured using a paired Picarro-Licor trace gas analyzer from June – August 2021 (flux data is in `ghg_raw.csv`, data for statistical analysis in `ghg_stats.csv`). Net nitrogen mineralization data was collected through ion exchange resins (data is in `n_raw.csv`). Bulk density, soil texture, pH, and soil carbon and nitrogen were completed on soil samples analyzed at Yale University (data is in `soilcn.csv`). Carbon and nitrogen were measured on an elemental analyzer. Soil texture was estimated using particle size analysis. pH was measured using a bench top pH meter. Land cover classification was done using the NatureServe database for Sublette County, WY (data is in `landcover.xlsx` and in landcover.zip). In 'included' tab on landcover.xlsx, there are two sections of data on the left and right, not one rectangular table of data. ## Description of the data and file structure #### ghg\_raw\.csv: * date: date of sampling (mm/dd/yy) * week: week of sampling (Week 1 - Week 11) * id: soil collar ID * ranch: site name (site1-3) * location: one of the four land cover types measured [meadow (hay meadow); bog (introduced wetland); slsage (sloping sagebrush); upsage (upland sagebrush)] * sh/is: indicates whether sample was taken under shrub or in the interspace * stemp: soil temperature ('b0C) * pcent_sand: % sand content * ph: soil pH * chamber offset: height of soil collar from the surface of soil (measured in cm) * volume add: volume added through extension of chamber attachment (cm3) * start time: start time of sampling (hour:minute:second) * end time: end time of sampling (hour:minute:second) * time off: offset from the time read by the two instruments (s) * co2 flux: carbon dioxide flux measured over 8 minute period (umol co2 m-2 s-1) * co2 rsq: R-squared of exponential fit of values over 8 minute period * ch4 flux: methane flux measured over 8 minute period (umol ch4 m-2 s-1) * ch4 rsq: R-squared of exponential fit of values over 8 minute period * n2o flux: nitrous oxide flux measured over 8 minute period (umol n2o m-2 s-1) * n2o rsq: R-squared of exponential fit of values over 8 minute period * nh3 flux: ammonia flux measured over 8 minute period (umol nh3 m-2 s-1) * nh3 rsq: R-squared of exponential fit of values over 8 minute period * h2o flux: h2o flux (% water vapor) * h2o rsq: R-squared of exponential fit of values over 8 minute period * notes: anything important to note from data collection of that day #### ghgstats.csv: * date: date of sampling (mm/dd/yy) * week: week of sampling (Week 1 - Week 11) * id: soil collar ID * ranch: site name (site1-3) * location: one of the four land cover types measured [meadow (hay meadow); bog (introduced wetland); slsage (sloping sagebrush); upsage (upland sagebrush)] * stemp: soil temperature ('b0C) * h2o flux: h2o flux (% water vapor) * h2o rsq: R-squared of exponential fit of values over 8 minute period * sand: % sand content (0-5 cm) * clay: % clay content (0-5 cm) * pH: soil pH (0-5 cm) * carbon: % soil carbon (0-5 cm) * nitrogen: % soil nitrogen (0-5 cm) * co2 flux: carbon dioxide flux measured over 8 minute period (umol co2 m-2 s-1) * co2 rsq: R-squared of exponential fit of values over 8 minute period * ch4 flux: methane flux measured over 8 minute period (umol ch4 m-2 s-1) * ch4 rsq: R-squared of exponential fit of values over 8 minute period #### n\_raw\.csv: * ranch: site name (site1-3) * location: one of the four land cover types measured [meadow (hay meadow); bog (introduced wetland); slsage (sloping sagebrush); upsage (upland sagebrush)] * no3: nitrate concentrations (ug no3/10cm2/burial length) * nh4: ammonium concentrations (ug nh4/10cm2/burial length) * fe: iron concentrations (ug fe/10cm2/burial length) * s: sulfur concentrations (ug s/10cm2/burial length) #### soilcn.csv: * ranch: site name (site1-3) * location: one of the four land cover types measured [meadow (hay meadow); bog (introduced wetland); slsage (sloping sagebrush); upsage (upland sagebrush)] * depth: depth of soil sample (cm) * finebulkd: fine soil bulk density (g soil cm-3) * c: % soil carbon (0-5 cm) * n: % soil N (0-5 cm) * cpool: total soil carbon for 0-5 cm (g C m-2) * npool: total soil N for 0-5 cm (g N m-2) #### landcover.xlsx: Sheet1 (ecosystems): * gridcode: NatureServe grid code assignment * Shape Area: area covered by a given land cover type (m2) * Value: grid code assignment * Count: Number of land cover types classified under a given name within the county boundary * ECOLSYS_LU: Land cover type classification under NatureServe categories Sheet2 (included): * ecotype: Land cover type classification under NatureServe categories * shape area: area covered by a given land cover type (m2) * Total shape area: area covered by landscape of included cover types (m2) * % introduced wetland: % of area represented by introduced wetland cover types within landscape * % hay meadow: % of area represented by hay meadow cover types within landscape * % big sagebrush: % of area represented by big sagebrush cover types within landscape #### **landcover.zip** * landcover.dbf: database file, attribute data * landcover.prj: projection file, coordinate system * landcover.shp: shapefile, stores geometry of spatial data * landcover.shx: shape index file ## Code/Software #### figs.R: used to make the primary figures for the manuscript and supplementary figures Completed in R (version 4.4.2) Packages: `ggplot2`, `wesanderson`, `stringr`, `zoo`, `tidyverse`, `scales`, `ggpubr`, `cowplot`, `gt`, `dplyr`, 'sf' #### stats.R: -used to do the statistical analysis for the manuscript -Completed in R (version 4.4.2) -Packages: 'pracma', `tidyverse`, `ggpubr`, `rstatix`, `lme4`, 'lmerTest', `emmeans`, 'car', `MASS`, `glmm` ## **Supplementary ** scalefigure.pptx * used to make figure 6 for manuscript (weighted cumulative methane flux by cover type) ## Changes from previous version: * Created a few different figures (one land cover map and made some edits to existing figures). The script 'figs.r' reflects those changes. I added the land cover shapefile I used to make the land cover map (landcover.zip). I ran a few new analyses on cumulative gas fluxes, which is included in the script, 'stats.r'. All other files are the same.  Drylands cover approximately 40% of the global land surface and are thought to contribute significantly to the soil methane sink. However, large-scale methane budgets have not fully considered the influence of agricultural land use change in drylands, which often includes irrigation to create land cover types that support hay or grains for livestock production. These land cover types may represent a small proportion of the landscape but could disproportionately contribute to greenhouse gas exchange and are currently omitted in estimates of dryland methane fluxes. We measured greenhouse gas fluxes among big sagebrush, introduced wetlands, and hay meadows in a semi-arid temperate dryland in Wyoming, USA to investigate how these small-scale irrigated land cover types contributed to landscape-scale methane dynamics. Big sagebrush ecosystems dominated the landscape while the introduced wetlands and hay meadows represented 1% and 12%, respectively. Methane uptake was consistent in the big sagebrush ecosystems, emissions and uptake were variable in the hay meadows, and emissions were consistent in the introduced wetlands. Despite making up 1% of the total land area, methane production in the introduced wetlands overwhelmed consumption occurring throughout the rest of the landscape, making this region a net methane source. Our work suggests that introduced wetlands and other irrigated land cover types created for livestock production may represent a significant, previously overlooked source of anthropogenic methane in this region and perhaps in drylands globally. Raw flux data for methane, carbon dioxide, and other species were measured using a paired Picarro-Licor trace gas analyzer from June – August 2021 (flux data is in ghg_raw.csv, data for statistical analysis in ghg_stats.csv). Net nitrogen mineralization data was collected through ion exchange resins (data is in n_raw.csv). Bulk density, soil texture, pH, and soil carbon and nitrogen were completed on soil samples analyzed at Yale University (data is in soilcn.csv). Carbon and nitrogen were measured on an elemental analyzer. Soil texture was estimated using particle size analysis. pH was measured using a bench top pH meter. Land cover classification was done using the NatureServe database for Sublette County, WY (data used for scaling is in landcover.xlsx and shapefile is in landcover.zip). In 'included' tab on landcover.xlsx, there are two sections of data on the left and right, not one rectangular table of data.

Snow insulates the soil from air temperature, decreasing winter cold stress and altering energy use for organisms that overwinter in the soil. As climate change alters snowpack and air temperatures, it is critical to account for the role of snow in modulating vulnerability to winter climate change. Along elevational gradients in snowy mountains, snow cover increases but air temperature decreases, and it is unknown how these opposing gradients impact performance and fitness of organisms overwintering in the soil. We developed experimentally validated ecophysiological models of cold and energy stress over the past decade for the montane leaf beetle Chrysomela aeneicollis, along five replicated elevational transects in the Sierra Nevada mountains in California. Cold stress peaks at mid-elevations, while high elevations are buffered by persistent snow cover, even in dry years. While protective against cold, snow increases energy stress for overwintering beetles, particularly at low elevations, potentially leading to mortality or energetic trade-offs. Declining snowpack will predominantly impact mid-elevation populations by increasing cold exposure, while high elevation habitats may provide refugia as drier winters become more common. Climate Data This is the full data set that includes all temperature measures used in cold exposure pricipal component analysis and energy use model output for each winter at each site included in study. The second tab includes metadata for each column. Physiology Data This is the data for all beetles used in the field snow manipulation experiment and all biochemichal assays performed. It includes total protein, glycerol, sorbitol, glucose, and triacylglycerides. Respirometry This is the data set containing all respirometry data used in making the energy use model.

Field leaf trait and spectroscopy data – We used leaf trait and spectral data from an extensive field campaign along an elevation gradient (from 3500 m to 220 m elevation) in the Peruvian Amazon where leaf traits for 60-80% of basal area of trees >10cm DBH were measured within a well-studied 1 ha plot network from April – November 2013 (Enquist et al., 2017). In each one ha plot (N=10 plots), we sampled the most abundant species as determined through basal area weighting (enough species generally to cover ~80% of the plot’s basal area). For each species, we sampled the five (three in the lowlands) largest trees (based on diameter at breast height (DBH)) and sampled one sun and one shade branch. On each of these branches, leaf chemistry and leaf mass area (LMA) were measured with the methodology detailed in Asner et al. (2014). On five randomly selected leaves for each branch, we measured hemispherical reflectance with an ASD Fieldspec Handheld 2 with fiber optic cable, a contact probe that has its own calibrated light source, and a leaf clip (Analytical Spectral Devices High-Intensity Contact Probe and Leaf Clip, Boulder, Colorado, USA) following (Doughty et al., 2017). We measured leaf spectroscopy (400-1075 nm) on the same branches where the leaf traits were collected. Both LMA and Chlorophyll A had previously been shown with this dataset to have a correlation with leaf spectroscopy (Doughty et al., 2017). However, we had not previously tried to compare leaf spectral data with DBH directly. Plot data – Aboveground biomass - We used 2,102 of 19,160 total AGB field plots between +30° and -30° latitude classified as broadleaf evergreen trees by MODIS PFT using public data from Duncanson et al 2022 that was organized and publicly available through ORNL DAAC as an RDS (R data serialization) file. Distribution plots are shown in Fig S1 (AGB) and S2 (residuals). NPP and GPP - We also used 21, 1 ha plots where NPP and sometimes GPP were measured following the GEM protocol (Malhi et al., 2021). We focused on two regions: a Peruvian elevation transect with both NPP + GPP (n= 10, RAINFOR plot codes are ALP11, ALP30, SPD02, SPD01, TRU03, TRU08, TRU07, ESP01, WAY01, ACJ01(Malhi et al., 2017)) and a Bornean logging transect with only NPP (n= 11 RAINFOR plot codes are DAN-04, DAN-05, LAM-01, LAM-02, MLA-01, MLA-02, SAF-01, SAF-02, SAF-03, SAF-04, SAF-05 (Riutta et al., 2018). These plots were chosen because there are large changes in NPP/GPP across the elevation or logging gradient. GEDI data – We used the vertical forest structure (L2A and L2B, Version 2) and biomass (L4a) products from the GEDI instrument (R. Dubayah et al., 2020) between April 2019 to December 2022 for tropical forest regions (R. O. Dubayah et al., 2023). We used a quality filtering recipe developed in collaboration with GEDI Science Team members from the University of Maryland and NASA Goddard to identify the highest quality GEDI vegetation shots (R. Dubayah et al., 2022). A data layer that this iterative local outlier detection algorithm uses to exclude data is publicly available at R. O. Dubayah et al., (2023). For instance, some of the key data filters we applied were: included degrade flags of 0,3,8,10,13,18,20,23,28,30,33,38,40,43,48,60,63,68, L2A and L2B quality flags = 1 (only use highest quality data), sensitivity >= 0.98. With the GEDI data, we used canopy height, the height of median energy (HOME), and the number of canopy layers following Doughty et al 2023 (Doughty et al., 2023). Across all tropical forests, we created 300 by 300 m pixels containing all averaged (mean) GEDI data between 2019 and 2022. Using the centroid coordinates from each of the 2,102 plots, we found the 300 by 300 m averaged GEDI pixel that encompassed the plot. If the plot was not encompassed by the GEDI data, we searched a wider area by incrementally averaging a gradually increasing area of 1, 3, 5, and 10 pixels. In other words, if no 300 by 300 m pixel encompassed the plot, then we averaged all GEDI data an area one pixel out (4 by 4 = 1200 by 1200 m, 6 by 6 = 1800 by 1800m, 11 by 11 = 3300m by 3300m), gradually increasing the square until it encompassed an area with GEDI data. To compare with the NPP/GPP plots we compared RS trait and GEDI data for individual footprints within a 0.03 km radius of the plot coordinates. Remotely sensed leaf trait data – Based on a broader set of field campaigns, Aguirre-Gutiérrez et al., (2021) used Sentinel-2, climatic, topography, and soil data to create remotely sensed canopy trait maps for P=phosphorus % leaf concentration, WD = wood density g.cm-3, and LMA=Leaf mass area g m-2. Other data layers – We compared % one peak to several other climates, soils, leaf traits, and ecoregion maps listed below for the Amazon basin. Each dataset had its own resolution, which we standardized to 0.1 by 0.1 degrees. We used total cation exchange capacity (CEC) from soil grids (Batjes et al., 2020) from 0-5cm in units of mmol(c)/kg. We averaged TerraClimate (Abatzoglou et al., 2018) data between 2000 and 2018 for Vapor Pressure Deficit (VPD in kPa), Mean Monthly Precipitation (MMP) (mm/month), potential evapotranspiration (PET) and maximum and minimum temperature (°C). Statistical analysis – We used the Matlab (Matlab, MathWorks Inc., Natick, MA, USA) function “fitlm” to fit linear models to compare variables such as soil data, environmental data, leaf trait data (at 0.1° resolution) and GEDI structure data (300m and bigger resolution) to field biomass and NPP/GPP estimates. The P values listed are for the t-statistic of the two-sided hypothesis test. We used R to create a linear model to predict the best model ranked by AIC and parsimony using the dredge function from the MuMIn library (Bartoń, 2009). We also used the CAR package (Fox J & S, 2019) and the VIF command to test for multi-collinearity between variables. To account for spatial autocorrelation, we used Simultaneous Auto-Regressive (SARerr) models (F. Dormann et al., 2007) using the R library ‘spdep’ (Bivand, Hauke, & Kossowski, 2013). We tested different neighborhood distances from 10 km to 300 km and found that AIC was minimized at 80 km (Fig S3) and the corresponding correlogram showed reduced spatial autocorrelation (Fig S4). To predict leaf traits with the spectral information, we used the Partial Least Squares Regression (PLSR) (Geladi & Kowalski, 1986) using the PLSregress command in Matlab (Matlab, MathWorks Inc., Natick, MA, USA). To avoid over-fitting the number of latent factors, we minimized the mean square error with K-fold cross-validation. We use 70% of our data to calibrate our model and then the remaining 30% to test the accuracy of our model using r2. We use adjusted r2 which penalizes for small sample sizes throughout the manuscript. # Satellite-derived trait data slightly improves tropical forest biomass, NPP, and GPP predictions [https://doi.org/10.5061/dryad.ttdz08m5n](https://doi.org/10.5061/dryad.ttdz08m5n) The dataset contains leaf trait and spectral data to create Figures 1 and 2. It contains plot biomass data and satellite-derived leaf trait and structure data to create Figures 3-6. It contains plot NPP, GPP, and satellite-derived leaf trait and structure data to create Figures 7-8. ## Description of the data and file structure, including the associated Code/Software The Matlab code Finalcode_GEDIbiomass_Doughty2024.m contains all the code and data to create all the figures in the paper. The code has several sections that can be run independently. The first section starting on line 1 uses the dataset traitcompare.mat to create Figure 1. This dataset contains two tables of leaf trait data called carnegiechem and merged. Units and column names are contained within the table. The second section starting on line 62 uses the dataset traitgedidat.mat to create Figures 3-5 and Figures S1 and S2. This dataset contains a table called agball with plot biomass and coordinates. It also has the trait and GEDI data for these plots in nested structures. Units and descriptions are given in the code. The third section starting on line 443 uses the datasets traitgedidat.mat and soilclimdata.mat to create Figure 6. The dataset traitgedidat.mat is the same as described above and soilclimdata.mat contains 0.1 by 0.1 degree gridded data for climate variables like Tmax (C) or VPD (Pa) or soil chemistry like CEC. Units and description are given in the code. The forth section starting on line 536 uses the dataset Tamtreeheight.mat to create Figures 7 and 8. The dataset gedivsplot.mat contains table data with plot data, and nearby trait and GEDI data for several GEM plots. Units and description are given in the code and the tables. The fifth section starting on line 791 uses the dataset CombspecDBH.mat to create Figure 2. The dataset has variables specallz which is the leaf spectral data from 350-1075 nm for each leaf and dbhz1 with is the corresponding tree dbh (cm). It also has LMAz which is the LMA data with datz as the corresponding spectral data. To estimate spatial autocorrelation and the best model by AIC to create Table 1 and Figures S3 and 4, we used the R code processgedidata.r and the dataset biomass_trait_GEDI.xlsx. This dataset contains a table with latitude, longitude, field biomass and remote sensed biomass (Mg Ha-1), and traits LMA (g m2), Phosphorus (%), tree height (m), HOME (m) and % one peak (unitless). ## Sharing/Access information Original GEDI data are available from the USGS. Improving tropical forest biomass predictions can accurately value tropical forests for their ecosystem services. Recently, the Global Ecosystem Dynamics Investigation (GEDI) lidar was activated on the international space station (ISS) to improve biomass predictions by providing detailed 3D forest structure and height data. However, there is still debate on how best to predict tropical forest biomass using GEDI data. Here we compare GEDI predicted biomass to 2,102 tropical forest biomass plots and find that adding a remotely sensed (RS) trait map of LMA (Leaf Mass per Area) significantly (P<0.001) improves field biomass predictions, but by only a small amount (r2=0.01). However, it may also help reduce the bias of the residuals because, for instance, there was a negative relationship between both LMA (r2 of 0.34) and %P (r2=0.31) and residuals. This improvement in predictability corresponds with measurements from 523 individual trees where LMA predicts Diameter at Breast height (DBH) (the critical measurement underlying plot biomass) with an r2=0.04, and spectroscopy (400-1075 nm) predicts DBH with an r2=0.01. Adding environmental datasets may offer further improvements and max temperature (Tmax) predicts Amazonian biomass residuals with an r2 of 0.76 (N=66). Finally, for a network of net primary production (NPP) and gross primary production (GPP) plots (N=21), RS traits are better at predicting fluxes than structure variables like tree height or Height Of Median Energy (HOME). Overall, trait maps, especially future improved ones produced by surface biology geology (SBG), may improve biomass and carbon flux predictions by a small but significant amount.

<b>Abstract</b><br/>Leaf energy balance may influence plant performance and community composition. While biophysical theory can link leaf energy balance to many traits and environment variables, predicting leaf temperature and key driver traits with incomplete parameterizations remains challenging. Predicting thermal offsets (δ, Tleaf – Tair difference) or thermal coupling strengths (β, Tleaf vs. Tair slope) is challenging. We ask: 1) whether environmental gradients predict variation in energy balance traits (absorptance, leaf angle, stomatal distribution, maximum stomatal conductance, leaf area, leaf height); 2) whether commonly-measured leaf functional traits (dry matter content, mass per area, nitrogen fraction, δ13C, height above ground) predict energy balance traits; and 3) how traits and environmental variables predict δ and β among species. We address these questions with diurnal measurements of 41 species co-occurring along an 1100 m elevation gradient spanning desert to alpine biomes. We show that 1) energy balance traits are only weakly associated with environmental gradients, and 2) are not well predicted by common functional traits. We also show that 3) δ and β can be partially approximated using interactions among site environment and traits, with a much larger role for environment than traits. The heterogeneity in leaf temperature metrics and energy balance traits challenges larger-scale predictive models of plant performance under environmental change.

Available files to conduct a Bayesian integral projection model (IPM) and population viability analysis (PVA) for the California tiger salamander (CTS) include: -The preliminary frequentist CTS IPM and PVA script created by Christopher A. Searcy is: "CTS-Frequentist-IPM.R" -- Frequentist code to build the IPM and run the PVA. -The primary scripts to run the IPM and PVA are: "CTS_SOURCE.R" -- Source code for the vital rate functions. "CTS_IPM_PVA.R' -- Code to build the IPM, conduct sensitivity and elasticity analyses, and run the PVA. -The scripts used to build the vital rate functions that inform the IPM are: "CTS_Best_Survival.R" -- The best Cormack-Jolly-Seber model of metamorph and juvenile/adult CTS survival and recapture probabilities. "CTS_Growth.R" -- CTS metamorph and juvenile/adult growth functions. "CTS_Fertility.R" -- CTS fertility function. "CTS_Maturity.R" -- CTS maturity function. "CTS_Larval_Survival.R" -- CTS larval survival given egg density function. "CTS_Females_Precip.R" -- The proportion of CTS females breeding given annual December-January precipitation function. "CTS_Replacement.R" -- Adult-only replacement and reproductive success functions to construct piecewise environmental-dependency function. -All necessary data files to run the CTS_IPM_PVA.R script and support the findings of our study are: "adults-v2.txt" -- The adult CTS capture histories from the capture-mark-recapture study at Jepson prairie Preserve, CA. "covariates-v2.txt" -- The CTS capture histories from the capture-mark-recapture study at Jepson prairie Preserve, CA specifying individual body masses (ln-transformed; g). "metamorphs-v2.txt" -- The metamorph CTS capture histories from the capture-mark-recapture study at Jepson prairie Preserve, CA. "precip.csv" -- Study rain year-specific November-February and October-June precipitation values (mm). "density-distance.csv" -- Proportion of the post-metamorphic life stage individuals (from 0 to 1) found within 100-m distance radius increments from the pond shoreline. "Larval-Survival-Density.csv" -- Ln-transformed larval survival and Ln-transformed egg density (eggs/m^3) data. "meta-size-by-egg-density.csv" -- Study year-specific metamorph size and egg density (eggs/m^3) data. "Olcott20**.txt" -- Body size distributions of each cohort of CTS from the mark-recapture study across the 122 body size bins, where "**" indicates study year in the file name. "Stochastic_Climate_Pool.txt" -- Historic precipitation record from 1893-2012 (Vacaville and Nut Tree Airport station records). "Stochastic_Climate_Pool_Rev.txt" -- Historic precipitation record from 1893-2008 (Nut Tree Airport-only station records). "females-precip-posterior-samples-CENTERED.csv" -- The 500 random samples from the posterior distribution of the function of female CTS breeding given December-January precipitation (mm). "fertility-posterior-samples-CENTERED.csv" -- The 500 random samples from the posterior distribution of the function of clutch size (# eggs) given female CTS body mass (g). "growth-posterior-samples-CENTERED.csv" -- The 500 random samples from the posterior distributions of the life stage-specific growth functions. "larval-survival-posterior-samples.csv" -- The 500 random samples from the posterior distribution of the function of larval survival probability given egg density (eggs/m^3). "maturity-posterior-samples-CENTERED.csv" -- The 500 random samples from the posterior distribution of the function of maturity probabilitygiven CTS body mass (g). "replace-success-infection-samples-CENTERED.csv" -- The 500 random samples of the inflection point from the posterior distribution of the function of probability of metamorph recruitment being above the replacement rate given October-June precipitation (mm). "repro-success-infection-samples-CENTERED.csv" -- The 500 random samples of the inflection point from the posterior distribution of the function of probability of reproductive success given October-June precipitation (mm). "survival-posterior-samples-CENTERED.csv" -- The 500 random samples from the posterior distribution of the Cormack-Jolly-Seber model of life stage-specific survival probabilities given body mass (g). Scripts were developed and run using R version 4.0.0. Integral projection models (IPMs) can estimate the population dynamics of species for which both discrete life stages and continuous variables influence demographic rates. Stochastic IPMs for imperiled species, in turn, can facilitate population viability analyses (PVAs) to guide conservation decision-making. Biphasic amphibians are globally distributed, often highly imperiled, and ecologically well-suited to the IPM approach. Herein, we present the first stochastic size- and stage-structured IPM for a biphasic amphibian, the U.S. federally threatened California tiger salamander (Ambystoma californiense; CTS). This Bayesian model reveals that CTS population dynamics show the greatest elasticity to changes in juvenile and metamorph growth and that populations are likely to experience rapid growth at low density. We integrated this IPM with climatic drivers of CTS demography to develop a PVA and examined CTS extinction risk under the primary threats of habitat loss and climate change. The PVA indicates that long-term viability is possible with surprisingly high (20–50%) terrestrial mortality, but simultaneously identified likely minimum terrestrial buffer requirements of 600–1000 m while accounting for numerous parameter uncertainties through the Bayesian framework. These analyses underscore the value of stochastic and Bayesian IPMs for understanding both climate-dependent taxa and those with cryptic life histories (e.g., biphasic amphibians) in service of ecological discovery and biodiversity conservation. In addition to providing guidance for CTS recovery, the contributed IPM and PVA supply a framework for applying these tools to investigations of ecologically-similar species. Please see the associated manuscript for full methodological details.

# Vegetation growth responses to climate change: a cross-scale analysis of biological memory and time-lags using tree ring and satellite data The dataset includes tree-ring data for individual trees across three species, encompassing dimensionless tree-ring width (TRW) measurements, as well as data on the enhanced vegetation index (EVI), leaf area index (LAI), gross primary productivity (GPP), and various climate parameters. The TRW serves as an indicator of radial stem growth at the tree-species level. Remote sensing-based data of EVI, LAI and GPP were used to monitor ecosystem-scale canopy dynamics, leaf growth, and ecosystem carbon sequestration capacity, respectively. ## Description of the data and file structure 1. Climate_1956_2017.csv: The dataset includes the mean air temperature, mean maximum air temperature, mean minimum air temperature, mean sunshine duration, and total precipitation from 1956 to 2017 on a daily basis in the study area. *Notes*: Lat, Latitude; Lon, longitude; Elev, Elevation; MTEM, mean air temperature (ºC); MaxTEM, mean maximum air temperature (ºC); MinTEM, mean maximum air temperature (ºC); X20to20PRE, accumulated precipitation at 20-20 (mm); SSD, mean sunshine duration (h). 2. TRW_LF.csv: This dataset comprises data for each core of individual trees belonging to the Liquidambar formosana (LF), coded as LF_01A, where 'LF' denotes the tree species, '01' represents the tree number, and 'A' indicates the core sample number taken from each tree. The units for this tree-ring data are in 0.001mm. 3. TRW_CE.csv: This dataset comprises data for each core of individual trees belonging to the Castanopsis eyrei (CE), coded as CE_01A, where 'CE' denotes the tree species, '01' represents the tree number, and 'A' indicates the core sample number taken from each tree. The units for this tree-ring data are in 0.001mm. 4. TRW_CH.csv: This dataset comprises data for each core of individual trees belonging to the Castanea henryi (CH), coded as CH_01A, where 'CH' denotes the tree species, '01' represents the tree number, and 'A' indicates the core sample number taken from each tree. The units for this tree-ring data are in 0.001mm. 5. Dimensionless_TRW_data_of_the_three_tree_species.csv: Between October 2020 and July 2022, we sampled 25-29 mature and healthy trees per species, collecting one-to-two cores from each tree at 1.3 m above the ground using a 5.15 mm increment borer. The tree-ring cores were fixed, dried, polished, and visually cross-dated under a binocular microscope. We measured tree-ring width with the LINTAB™ 6 system to a 0.01-mm accuracy, covering data from 1957 to 2017. Standardization of tree-ring width data involved two phases. First, COFECHA software ensured the quality of cross-dating results by evaluating the synchronization of growth patterns across samples. Next, we used the detrend function from the dplR package in R to fit a modified negative exponential curve to each raw tree-ring series for detrending. Standardized indices were calculated by dividing the original ring widths by the fitted values and combining them into a single standardized chronology using a bi-weight robust mean to mitigate outlier influence. *Notes*: CE, Castanopsis eyrei; CH, Castanea henryi; LF, Liquidambar formosana. 6. EVI_MOD13Q1_16days.csv: The dataset consists of the enhanced vegetation index (EVI) for the study area, measured over 16-day periods. *Notes*: Start, date of start; End, date of start; EVI, enhanced vegetation index (unitless). 7. LAI_MCD15A2H_16days.csv: The dataset consists of the leaf area index (LAI) for the study area, measured over 16-day periods. To ensure a consistent time resolution for remote sensing-based vegetation indicators, the 8-day time periods of LAI was aligned with the 16-day time periods of EVI. This alignment was achieved by averaging LAI values from two consecutive 8-day periods. *Notes*: Start, date of start; End, date of start; LAI, leaf area index (m2/m2). 8. GPP_MOD17A2H_16days.csv: The dataset consists of the gross primary productivity (GPP) for the study area, measured over 16-day periods. To ensure a consistent time resolution for remote sensing-based vegetation indicators, the 8-day time periods of GPP was aligned with the 16-day time periods of EVI. This alignment was achieved by calculating GPP as the cumulative value of two consecutive 8-day periods. *Notes*: Start, date of start; End, date of start; GPP, gross primary productivity (kg C/m2). Vegetation growth is affected by past growth rates and climate variability. However, the impacts of vegetation growth carryover (VGC; biotic) and lagged climatic effects (LCE; abiotic) on tree stem radial growth may be decoupled from photosynthetic capacity, as higher photosynthesis does not always translate into greater growth. To assess the interaction of tree-species level VGC and LCE with ecosystem-scale photosynthetic processes, we utilized tree-ring width (TRW) data for three tree species: Castanopsis eyrei (CE), Castanea henryi (CH, Chinese chinquapin), and Liquidambar formosana (LF, Chinese sweet gum), along with satellite-based data on canopy greenness (EVI, enhanced vegetation index), leaf area index (LAI), and gross primary productivity (GPP). We used vector autoregressive models, impulse response functions, and forecast error variance decomposition to analyze the duration, intensity, and drivers of VGC and of LCE response to precipitation, temperature, and sunshine duration. The results showed that at the tree-species level, VGC in TRW was strongest in the first year, with an average 77% reduction in response intensity by the fourth year. VGC and LCE exhibited species-specific patterns; compared to CE and CH (diffuse-porous species), LF (ring-porous species) exhibited stronger VGC but weaker LCE. For photosynthetic capacity at the ecosystem scale (EVI, LAI, and GPP), VGC and LCE occurred within 96 days. Our study demonstrates that VGC effects play a dominant role in vegetation function and productivity, and that vegetation responses to previous growth states are decoupled from climatic variability. Additionally, we discovered the possibility for tree-ring growth to be decoupled from canopy condition. Investigating VGC and LCE of multiple indicators of vegetation growth at multiple scales has the potential to improve the accuracy of terrestrial global change models. The dataset includes tree-ring data for individual trees across three species, encompassing dimensionless tree-ring width (TRW) measurements, as well as data on the enhanced vegetation index (EVI), leaf area index (LAI), gross primary productivity (GPP), and various climate parameters. The TRW serves as an indicator of radial stem growth at the tree-species level. Remote sensing-based data of EVI, LAI and GPP were used to monitor ecosystem-scale canopy dynamics, leaf growth, and ecosystem carbon sequestration capacity, respectively. Dimensionless tree-ring width (TRW) measurements method: Between October 2020 and July 2022, we sampled 25-29 mature and healthy trees per species, collecting one-to-two cores from each tree at 1.3 m above the ground using a 5.15 mm increment borer. The tree-ring cores were fixed, dried, polished, and visually cross-dated under a binocular microscope. We measured tree-ring width with the LINTAB™ 6 system to a 0.01-mm accuracy, covering data from 1957 to 2017. Standardization of tree-ring width data involved two phases. First, COFECHA software ensured the quality of cross-dating results by evaluating the synchronization of growth patterns across samples. Next, we used the detrend function from the dplR package in R to fit a modified negative exponential curve to each raw tree-ring series for detrending. Standardized indices were calculated by dividing the original ring widths by the fitted values and combining them into a single standardized chronology using a bi-weight robust mean to mitigate outlier influence.

Field Data - We estimate canopy temperature at the km 83 eddy covariance tower in the Tapajos region of Brazil 1–3 using a pyrgeometer (Kipp and Zonen, Delft, Netherlands) mounted at 64 m to measure upwelling longwave radiation (L↑ in W m-2) with an estimated radiative-flux footprint of 8,000 m2 4. Data were collected every 2 seconds and averaged over 30-minute intervals between August 2001 and March 2004. We estimated canopy temperature with the following equation: Eq 1 – Canopy temperature (°C) = (L↑/(E*5.67e-8))0.25-273.15 We chose an emissivity value (E) of 0.98 for the tower data, as this was the most common value used in the ECOSTRESS data (SDS_Emis1-5 (ECO2LSTE.001) and the broader literature for tropical forests 5. We compared canopy temperature derived from the pyrgeometer to eddy covariance derived latent heat fluxes (flux footprint ~1 km2), air temperature at 40 m, which is the approximate canopy height (model 076B, Met One, Oregon, USA; and model 107, Campbell Scientific, Logan, Utah, USA) and soil moisture at depths of 40 cm (model CS615, Campbell Scientific, Logan, Utah, USA). Further details on instrumentation and eddy covariance processing can be found in 1,3. This site was selectively logged, which had a minor overall impact on the forest 6, but did not affect any trees near the tower. Leaf thermocouple data - We measured canopy leaf temperature at a 30 m canopy walk-up tower between July to December of 2004 and July to December of 2005 at the same site. We initially placed 50 thermocouples on canopy-exposed leaves of Sextonia rubra, Micropholis sp., Lecythis lurida) (originally published in Doughty and Goulden 2008). Fine wire thermocouples (copper constantan 0.005 Omega, Stamford, CT) were attached to the underside of leaves by threading the wire through the leaf and inserting the end of the thermocouple into the abaxial surface. The thermocouples were wired into a multiplexer attached to a data logger (models AM25T and 23X, Campbell Scientific, Logan, UT, USA) and the data were recorded at 1 Hz. Additional upper-canopy leaf thermocouple data from Brazil7, Puerto Rico8, Panama9, Atlantic forest Brazil10 and Australia 11, were generally collected in a similar manner. Satellite data - ECOSTRESS data (ECO2LSTE.001) – The ECOsystem Spaceborne Thermal Radiometer Experiment on Space Station (ECOSTRESS) mission is a thermal infrared (TIR) multispectral scanner with five spectral bands at 8.28, 8.63, 9.07, 10.6, and 12.05 µm. The sensor has a native spatial resolution of 38 m x 68 m, resampled to 70 m x 70 m, and a swath width of 402 km (53°). Data are collected from an average altitude of 400 ± 25 km on the International Space Station (ISS). ECOSTRESS is an improvement over other thermal sensors because no other sensors provide TIR data with sufficient spatial, temporal, and spectral resolution to reliably estimate LST at the local-to-global scale for a diurnal cycle 12. To ensure the highest quality data, we used ECOSTRESS quality flag 3520, which identifies the best quality pixels (no cloud detected), a minimum-maximum difference (MMD) indicative of vegetation or water (Kealy and Hook 1993), and nominal atmospheric opacity. We accessed ECOSTRESS LST data through the AppEEARS website (https://lpdaac.usgs.gov/tools/appeears/) for the following products and periods: SDS_LST (ECO2LSTE.001) from a long longitudinal swath of the Amazon for 25 December 2018 to 20 July 2020 (SI Fig 1a red box) and then a larger area of the western Amazon for 18 September to 29 September 2019 (SI Fig 1a green box), Central Africa for 1 August to 30 August 2019 (SI Fig 1b), and SE Asia for 15 January to 30 February 2020 (SI Fig. 1c). The dates were chosen as all ECOSTRESS data available at the start of the study for the smaller regions and for warm periods with low soil moisture for the larger areas. We calculated “peak median,” which is defined as the average of the highest three medians of each granule (i.e., for the Amazon SI Fig. 1a, there were 934 granules) for each hour period. Comparison of LST data – We compared ECOSTRESS LST to VIIRS LST (VNP21A1D.001) and MODIS LST (MYD11A1.006). A more detailed comparison and description of these sensors can be found in Hulley et al 202113. Details for the sensors and quality flags used are given in Table S1. Broadly, G1 for ECOSTRESS and VIIRS is classified as vegetation (using emissivity) and of medium quality. G2 is classified as vegetation, but of the highest quality. MODIS landcover classifies this region as almost entirely broadleaf evergreen vegetation, but using MMD (emissivity) only 18% (VIIRS) and 12% (ECOSTRESS) of the data are classified as vegetation, rather than as soils and rocks (Table S2). Therefore, we use the vegetation classification (from MMD) as a very conservative estimate of complete forest canopy cover and not farms, urban, or degraded forest where rocks or soils are more likely to appear to satellites. SMAP data – To estimate pantropical soil moisture, we use the Soil Moisture Active Passive (SMAP) sensor and the product Geophysical_Data_sm_rootzone (SPL4SMGP.005). SMAP measurements provide remote sensing of soil moisture in the top 5 cm of the soil 14 and the L4 products combine SMAP observations and complementary information from a variety of sources. We accessed SMAP data from the AppEEARS website for the following products and periods: Amazon for 25 December 2018 to 20 July 2020 (SI Fig 1a), Central Africa for 25 December 2019 to 20 July 2020 (SI Fig 1b), and Borneo for 25 December 2018 to 20 July 2020 (SI Fig 1c). Warming experiments – For model validation, we used the results of three upper-canopy leaf and branch warming experiments of 2°C (Brazil), 3°C (Puerto Rico), and 4°C (Australia). The first experiment (Brazil), was 4 individual leaf-resistant heaters on each of 6 different upper-canopy species at the Floresta National (FLONA) do Tapajos as part of the Large-Scale Biosphere–Atmosphere Ecology Program (LBA-ECO) in Santarem, Brazil14. On the same six species, black plastic passively heated branches by an average ~2°C. Initially, heat balance sap flow sensors and the passive heaters were added to 40 branches, but we had confidence in the data from 9 heated and 4 control in the final analysis. The second experiment (Puerto Rico) had two species (Ocotea sintenisii (Mez) Alain and Guarea guidonia (L.) Sleumer where leaves were heated by 3 °C at the Tropical Responses to Altered Climate Experiment (TRACE) canopy tower site at Sabana Field Research Station, Luquillo, Puerto Rico8. The final experiment (Australia), which increased leaf temperatures by 4 °C, was conducted at Daintree Rainforest Observatory (DRO) in Cape Tribulation, Far North Queensland, Australia. Leaf heaters were installed using a pair of 30-gauge copper-constantan thermocouples, one reference leaf, and one heated with a target temperature differential of 4°C. There were two pairs in the upper canopy of each tree crown installed in 2–3 individuals across four species with the thermocouples installed on the underside of the leaves. Two absolute 36-gauge copper-constantan thermocouples were installed in each species to measure the leaf temperatures of the reference leaves. Thermocouple wires connected into an AM25T multiplexer from Campbell Scientific connected to a CR1000 Campbell datalogger. More details about the experiment and sensors can be found in 11. Model – We created a model of individual leaves on a tree (100 by 100 grid where each leaf is a pixel) to estimate the upper limit of tropical canopy temperatures with projected changes in climate. At the start of the simulation, we randomly applied the measured distribution (ambient Fig 1c) of canopy leaf temperatures >31.2 °C (chosen to give a mean canopy temperature of 33.2 ± 0.4 °C, matching the canopy average Fig 1b) to the entire grid. Each year we increased the mean air temperatures by 0.03°C to simulate a warming planet. As air temperatures reached +2, 3, and 4°C, we applied the leaf temperature distributions (but subtracted out the air temperature increases) from the different warming experiments (+2°C (Brazil), +3°C (Puerto Rico), and +4°C (Australia), respectively (Fig S7)). We ran the model at a daily time step with leaves flushing once a year (all dead leaves reset to living each year). In addition, to take into account the effect of climate inter-annual variation - specifically drought, these mean canopy temperatures were further increased or decreased by deviations from mean maximum air temperatures at 40 m pulled each day from the Tapajos eddy covariance tower1–3 and soil moisture at 40 cm depth (m3 m-3) which controlled canopy temperatures following equation 2 (Fig S6). Eq 2 – Canopy temperature (°C) = 46.5-33.6*soil moisture (m3 m-3) For example, in a non-drought year, on a day when max air temperatures were 0.1 °C higher than average and soil moisture was 0.01 m3 m-3 lower than average (which would add 0.3 °C to canopy temperatures (Eq 2)), we would add 0.4 °C to the grid canopy temperature that day. Every year, there was a 10% random probability of either a minor (80% probability) drought which reduced soil moisture by 0.1 m3 m-3 and increased air temperatures by 0.5 °C or severe drought (20% probability), which reduced soil moisture by 0.2 m3 m-3 and increased air temperatures by 1 °C. This is similar to the Amazon-wide temperature increases during the last El Niño 15. If an individual leaf temperature increases to above 46.7 °C (Tcrit) the leaf died, following Slot et al. (2021). Prior research has suggested that irreversible damage could begin at 45 °C 16 and T50 for tropical species is 49.9 °C 17, and we use these values in a sensitivity study. We further explore the impact of duration of Tcrit on mortality in a sensitivity study (ranging between needing a single exposure to four exposures to Tcrit to die). Over the season, if a leaf died, then it did not contribute towards canopy evapotranspiration. We ran simulations as a 3D canopy with an LAI of 5 where if the top leaf died, then it was replaced by a shade-adapted leaf with a Tcrit 1 °C lower 18. If each of the 5 LAIs died, then all leaves in that grid cell were dead and canopy evaporative cooling decreased by that percentage. Several lines of evidence suggest that under normal hydraulic conditions, when radiation load increases from ~350 to 1100 W m-2 (e.g. between shady and sunny conditions) average canopy temperature increases by ~3 °C and therefore, evaporative cooling for a full 1100 W m-2 is ~4.4°C4,19 (we vary this in a sensitivity study between 3.7 and 5.1°C). For example, if, over a year, 1000 leaves (10% of all leaves) surpass Tcrit and die, evaporative cooling for all leaves in the grid will be reduced by 10% (1000/(100 by 100 grid)) or 0.44 °C and 0.44 °C will be added to mean canopy temperature. Therefore, mean canopy temperature could heat up by a maximum of 4.4°C either due to a reduction of soil moisture or from an increase in dead leaves. We ran each simulation until the point where all leaves were dead and repeated this 30 times. We assumed loss of tree function following the death of all leaves, but we discuss this further in the discussion. We then ran sensitivity studies for several of the key variables (bold indicates the standard model parameter) including: drought (0.05, 0.1, to 0.2 m3 m-3 decrease in soil moisture), change in Tcrit (Tcrit: 45, 46.7, 49.9 °C), Tcrit range (100 by 100 grid =random distribution of 46.7±2, 100 by 100 grid =46.7±0), Max evaporative cooling (3.7, 4.4°C), (Tcrit duration (exceed Tcrit once, exceed Tcrit more than 3 times) and soil moisture coefficient (-33.6 -38.2; i.e. change the slope from Fig S6 by ± 1 sd). Methods References Miller, S. D. et al. Biometric and micrometeorological measurements of tropical forest carbon balance. Ecol. Appl. 14, 114–126 (2004). da Rocha, H. R. et al. Seasonality of water and heat fluxes over a tropical forest in eastern Amazonia. Ecol. Appl. 14, 22–32 (2004). Goulden, M. L. et al. Diel and seasonal patterns of tropical forest co2 exchange. Ecol. Appl. 14, 42–54 (2004). Kivalov, S. N. & Fitzjarrald, D. R. Observing the Whole-Canopy Short-Term Dynamic Response to Natural Step Changes in Incident Light: Characteristics of Tropical and Temperate Forests. Boundary-Layer Meteorol. 173, 1–52 (2019). Jin, M. & Liang, S. An Improved Land Surface Emissivity Parameter for Land Surface Models Using Global Remote Sensing Observations. J. Clim. 19, (2006). Miller, S. D. et al. Reduced impact logging minimally alters tropical rainforest carbon and energy exchange. Proc. Natl. Acad. Sci. 108, 19431 LP – 19435 (2011). Doughty, C. E. An In Situ Leaf and Branch Warming Experiment in the Amazon. Biotropica 43, 658–665 (2011). Carter, K. R., Wood, T. E., Reed, S. C., Butts, K. M. & Cavaleri, M. A. Experimental warming across a tropical forest canopy height gradient reveals minimal photosynthetic and respiratory acclimation. Plant. Cell Environ. 44, 2879–2897 (2021). Rey-Sanchez, A. C., Slot, M., Posada, J. & Kitajima, K. Spatial and seasonal variation of leaf temperature within the canopy of a tropical forest. Clim. Res. 71, 75–89 (2016). Fauset, S. et al. Differences in leaf thermoregulation and water use strategies between three co-occurring Atlantic forest tree species. Plant. Cell Environ. 41, 1618–1631 (2018). Crous K Y, A W Cheesman, K Middleby, Rogers Eie, A Wujeska-Klause, A Y M Bouet, D S Ellsworth, M J Liddell, L A Cernusak, C V M Barton, Similar patterns of leaf temperatures and thermal acclimation to warming in temperate and tropical tree canopies., Tree Physiology, 2023;, tpad054, https://doi.org/10.1093/treephys/tpad054. Xiao, J., Fisher, J. B., Hashimoto, H., Ichii, K. & Parazoo, N. C. Emerging satellite observations for diurnal cycling of ecosystem processes. Nat. Plants 7, 877–887 (2021). Hulley, G. C. et al. Validation and Quality Assessment of the ECOSTRESS Level-2 Land Surface Temperature and Emissivity Product. IEEE Trans. Geosci. Remote Sens. 60, 1–23 (2022). Reichle, R., Lannoy, G. De, Koster, R. D., Crow, W. T. & 2017., J. S. K. SMAP L4 9 km EASE-Grid Surface and Root Zone Soil Moisture Geophysical Data, Version 3. Boulder, Color. USA. NASA Natl. Snow Ice Data Cent. Distrib. Act. Arch. Center. doi https//doi.org/10.5067/B59DT1D5UMB4. (2017). Jiménez-Muñoz, J. C. et al. Record-breaking warming and extreme drought in the Amazon rainforest during the course of El Niño 2015–2016. Sci. Rep. 6, 33130 (2016). Berry, J. & Bjorkman, O. Photosynthetic Response and Adaptation to Temperature in Higher Plants. Annu. Rev. Plant Physiol. 31, 491–543 (1980). Slot, M. et al. Leaf heat tolerance of 147 tropical forest species varies with elevation and leaf functional traits, but not with phylogeny. Plant. Cell Environ. 44, (2021). Slot, M., Krause, G. H., Krause, B., Hernández, G. G. & Winter, K. Photosynthetic heat tolerance of shade and sun leaves of three tropical tree species. Photosynth. Res. 141, 119–130 (2019). Doughty, C. E. & Goulden, M. L. Are tropical forests near a high temperature threshold? J. Geophys. Res. Biogeosciences (2009) doi:10.1029/2007JG000632. The critical temperature beyond which photosynthetic machinery in tropical trees begins to fail averages ~46.7°C (Tcrit) 1. However, it remains unclear whether leaf temperatures experienced by tropical vegetation approach this threshold or soon will under climate change. We found that pantropical canopy temperatures independently triangulated from individual leaf thermocouples, pyrgeometers, and remote sensing (ECOSTRESS) have midday-peak temperatures of ~34°C during dry periods, with a long high-temperature tail that can exceed 40°C. Leaf thermocouple data from multiple sites across the tropics suggest that even within pixels of moderate temperatures, upper-canopy leaves exceed Tcrit 0.01% of the time. Further, upper-canopy leaf warming experiments (+2, 3, and 4°C in Brazil, Puerto Rico, and Australia) increased leaf temperatures non-linearly with peak leaf temperatures exceeding Tcrit 1.3% of the time (11% >43.5°C, 0.3% >49.9°C). Using an empirical model incorporating these dynamics (validated with warming experiment data), we found that tropical forests can withstand up to a 3.9 ± 0.5 °C increase in air temperatures before a potential collapse in metabolic function, but the remaining uncertainty in our understanding of Tcrit could reduce this to 2.6 ± 0.6°C. The 4.0°C estimate is within the “worst case scenario” (RCP-8.5) of climate change predictions2 for tropical forests and therefore it is still within our power to decide (e.g., by not taking the RCP 8.5 route) the fate of these critical realms of carbon, water, and biodiversity 3,4.

# Data for: Recent changes in thermal niche position and breadth of bird assemblages in Spain in relation to increasing temperatures Name: David Ramón-Martínez ORCID:0000-0001-7537-6254 Institution: Doñana Biological Station (EBD-CSIC) Address: Amrico Vespucio 26, Sevilla 41092, Spain Email: Name: Javier Seoane ORCID:0000-0001-9975-4846 Institution: Centro de Investigacion en Biodiversidad y Cambio Global, Universidad Autonoma de Madrid (CIBC-UAM); Terrestrial Ecology Group, Department of Ecology, Universidad Autonoma de Madrid(TEG-UAM). Address: Darwin, 2. Madrid 28049, Spain Email: **Aim:** Animal communities around the world are responding to climate change by altering their taxonomic composition, mainly through an increase in the colonisation rate of warm-dwelling species and the local extinction of cold-dwelling ones. We assessed whether the taxonomic composition of bird assemblages in peninsular Spain has changed in accordance with the recent increase in temperature. We also evaluated the role of species' thermal affinities and population dynamics on these changes. **Location:** Peninsular Spain. **Taxon:** Birds. **Methods:** We compared assemblages reported in the last Spanish breeding bird atlases (1998-2002 vs 2014-2019) in 10x10 km squares. We described species’ thermal niches by overlaying global species breeding distributions and world temperature metrics (based on mean, minimum, maximum and range), and then aggregated them to obtain a set of community thermal indices for each assemblage (CTIs, and CTR for ranges). Long-term average temperatures and local current temperatures were related to changes in CTIs using spatial GLMMs, which considered habitat change. We identified the species most responsible for variation in assemblages and regressed species’ influence on thermal affinities and population dynamics. **Results:** CTIs increased with temperature and warm-dwelling species became more prevalent to the detriment of cold-dwelling ones. However, we found a counteracting effect of temperature and habitat. Cold-dwelling forest species were among the most influential species, mainly through colonisation, while warm-dwelling farmland species contributed through local extinctions (both attenuated local increases in CTI). The mean thermal breadth of assemblages (CTR) decreased with temperatures. **Main conclusions:** The taxonomic composition of bird assemblages shifted in line with the main expectations due to global change (thermophilisation), mainly due to local colonisation of warm-dwelling species, although it did not show the pattern of thermal homogenization suggested elsewhere. Our results add further evidence of the interplay between climate warming and land-use change in the ongoing adjustment of animal communities. ## Description of the Data and file structure The dataset is a dataframe that comprises the Community Thermal Indices (response variable) and the standardized and unstandardized environmental and geographic variables employed as predictors of the spatial GLMM. This model related the temperatures to the changes in CTI, considering the habitat (forest) change. The Community Thermal Indices were computed from the Species Thermal Indices (Devictor et al., 2008). We obtained four thermal indices for each species (Species Thermal Index STI) by combining the global breeding species distribution and the climate information. The STI1 (i) shows the mean temperature of the breeding season (April-July) throughout the species breeding distribution range. Similarly, the STI2 (ii) is the average of the maximum temperatures above the percentile 95 in July and the STI3 (iii) is the average minimum temperature below the percentile 05 in April in the species’ breeding distribution range. These three indices represent a species' thermal affinity. On the other hand, the fourth index (iv) (Species Thermal Range - STR) represents the average thermal range (April-July) throughout the breeding distribution area and can be understood as species thermal breadth. We calculated a set of community thermal indices (CTI) for the assemblage of bird species in each of the 10x10km UTM grid squares of each of the breeding bird atlases. We obtained four different CTIs: CTI1, CTI2, CTI3, and CTR. The first three were calculated as the average of the STI1, STI2, and STI3 of the species present in the assemblage, respectively. The CTR (Community Thermal Range) is based on the average temperature range of the species (STR) that make up the assemblage and thus informs on the average niche breadth (Gaget et al., 2020). We calculated CTIs for each of the four-year periods covered by the atlases. The dataset also includes the standardized and unstandardized local temperature and forest cover for each grid square and for each breeding bird atlas. It also includes the standardized and unstandardized coordinates of each grid square. Local temperatures were obtained from Chelsa (v.2.1., Karger et al., 2017), averaging data for each five-year sampling period in each square. We used the CORINE Land Cover Accounting Layers built for the years 2000 and 2018, to link forest cover with the community indices for the first and second sampling periods, respectively. The variables included in the dataset are the following: * **UTM10**: The identity of each 10x10 km square grid from the Spanish Breeding Bird Atlases. * **fperiod**: Each of the sampling periods considered (1998-2002; 2014-2019). * **longitude**: Longitude of the grid square centroid (CRS: WGS84; EPSG=4326 ). * **latitude**: Latitude of the grid square centroid (CRS: WGS84; EPSG=4326). * **sd_longitude**: Standardized longitude of the grid square centroid. * **sd_latitude**: Standardized latitude of the grid square centroid. * **forest_cover**: Forest landcover (ha) in each square in 2000 and 2018 CORINE LandCover Accounting Layers versions. The forest landcover in 2018 is assigned to the second period observations (2014-2019), whereas the forest landcover in 2000 is assigned to the first period observations (1998-2002). We considered as forest landcover the CORINE/Landcover categories 311 “Broad leaf forest”; 312 “Coniferous Forest” and 313 “Mixed forest”. * **sd_forest_cover**: The standardized forest landcover in each square in 2000 and 2018 CORINE LandCover Accounting Layers versions. The forest landcover in 2018 is assigned to the second period observations (2014-2019), whereas the forest landcover in 2000 is assigned to the first period observations (1998-2002). We considered as forest landcover the CORINE/Landcover categories 311 “Broad leaf forest”; 312 “Coniferous Forest” and 313 “Mixed forest”. * **temperature**: The mean annual temperature (ºC) of each square grid in each period obtained from Chelsa v.2.1 (Karger et al., 2017). This dataset is based on downscaled air temperature two meters above the ground modelized from the data collected from many sources (mainly weather stations, weather balloons, aircraft, ships and satellites). The mean annual temperature of the period 1998-2002 is assigned to the observations from the first period (1998-2002). The mean annual temperature of the period 2014-2018 is assigned to the observations from the second period (2014-2019). Temperature was downloaded in Kelvin*10, and then converted to ºC previous to the analysis. * **sd_temperature:** The standardized mean annual temperature of each square grid in each period obtained from Chelsa v.2.1 (Karger et al., 2017). This dataset is based on downscaled air temperature two meters above the ground modelized from the data collected from many sources (mainly weather stations, weather balloons, aircraft, ships and satellites). The mean annual temperature of the period 1998-2002 is assigned to the observations from the first period (1998-2002). The mean annual temperature of the period 2014-2018 is assigned to the observations from the second period (2014-2019). * **CTI1**: Community Thermal Index 1. Average of the STI1 (thermal optimum) of the species present in a square grid. The STI1 is computed as the mean temperature (ºC) of the breeding season (April-July) along the global distribution range of a species during the breeding season. Wordclim monthly average temperatures for 1970-2000 (Worldclim 2.0: (Fick & Hijmans, 2017)) were used for this purpose. * **CTI2**: Community Thermal Index 2. Average of the STI2 (thermal maximum) of the species present in a square grid. The STI2 is computed as the average of the maximum temperatures (ºC) above the percentile 95 in July along the global distribution range of a species during the breeding season. Wordclim maximum temperatures of July for 1970-2000 (Worldclim 2.0: (Fick & Hijmans, 2017)) were used for this purpose. * **CTI3**: Community Thermal Index 3. Average of the STI3 (thermal minimum) of the species present in a square grid. The STI3 is computed as the average minimum temperature (ºC) below the percentile 05 in April along the global distribution range of a species during the breeding season. Wordclim minimum temperatures of April for 1970-2000 (Worldclim 2.0: (Fick & Hijmans, 2017)) were used for this purpose. * **CTR**: Community Thermal Range. Average of the STR (thermal range) of the species present in a square grid. The STR is computed as the difference between STI3 and STI2. ## Sharing/access Information Temperature for obtaining STI and CTI was obtained from Wordclim 2.0 (Fick & Hijmans, 2017). Local temperatures of square grids were computed from Chelsa v.2.1. (Karger et al., 2017) Grid square forest cover was obtained from CORINE LandCover Accounting Layers (EEA, 2019). Species global distribution maps were facilitated by Birdlife-International () REFERENCES 1. Devictor, V., Julliard, R., Couvet, D., & Jiguet, F. (2008). Birds are tracking climate warming, but not fast enough. Proceedings of the Royal Society B: Biological Sciences, 275(1652), 27432748. 2. EEA. (2019). Corine Land Cover Accounting Layers. 3. Fick, S. E., & Hijmans, R. J. (2017). WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology, 37(12), 43024315. 4. Gaget, E., Galewski, T., Jiguet, F., Guelmami, A., Perennou, C., Beltrame, C., & Le Viol, I. (2020). Antagonistic effect of natural habitat conversion on community adjustment to climate warming in nonbreeding waterbirds. Conservation Biology, 34(4), 966976. 5. Karger, D. N., Conrad, O., Bhner, J., Kawohl, T., Kreft, H., Soria-Auza, R. W., Zimmermann, N. E., Linder, H. P., & Kessler, M. (2017). Climatologies at high resolution for the earths land surface areas. Scientific Data, 4(1), 120. The dataset is a dataframe that comprises the Community Thermal Indices (response variable) and the environmental and geographic variables employed as predictors of the spatial GLMM. This model related the temperatures to the changes of CTI, considering the habitat (forest) change. The Community Thermal Indices were computed from the Species Thermal Indices. We obtained four thermal indices for each species (Species Thermal Index – STI) by combining the global species’ distribution and the climate information. The STI1 (i) shows the mean temperature of the breeding season (April-July) throughout the species’ distribution range. Similarly, the STI2 (ii) is the average of the maximum temperatures above the percentile 95 in July, and the STI3 (iii) is the average minimum temperature below the percentile 05 in April in the species’ breeding distribution range. These three indices represent a species’ thermal affinity. On the other hand, the fourth index (iv) (Species Thermal Range - STR) represents the average thermal range (April-July) throughout the distribution area and can be understood as species thermal breadth. It is computed as STI3-STI2. We calculated a set of community thermal indices (CTI) for the assemblage of bird species in each of the 10x10km UTM grid squares of each of the breeding bird atlases. We obtained four different CTIs: CTI1, CTI2, CTI3, and CTR. The first three were calculated as the average of the STI, STI2, and STI3 of the species present in the assemblage, respectively. The CTR (Community Thermal Range) is based on the average temperature range of the species (STR) that make up the assemblage and thus informs on the average niche breadth (Gaget et al., 2020). We calculated CTIs for each of the four-year periods covered by the atlases. The dataset also includes the standardized and unstandardized local temperature (ºC) and forest cover (ha) for each grid square and for each breeding bird atlas. It also includes the standardized and unstandardized coordinates of each grid square in decimal degrees (WGS84). Local temperatures were obtained from Chelsa (v.2.1., Karger et al., 2017), averaging data for each five-year sampling period in each square. We used the CORINE Land Cover Accounting Layers built for the years 2000 and 2018, to link forest cover with the community indices for the first and second sampling periods, respectively  Aim: Animal communities around the world are responding to climate change by altering their taxonomic composition, mainly through an increase in the colonisation rate of warm-dwelling species and the local extinction of cold-dwelling ones. We assessed whether the taxonomic composition of bird assemblages in peninsular Spain has changed in accordance with the recent increase in temperature. We also evaluated the role of species' thermal affinities and population dynamics in these changes. Location: Peninsular Spain. Taxon: Birds. Methods: We compared assemblages reported in the last Spanish breeding bird atlases (1998–2002 vs 2014–2019) in 10x10 km squares. We described species’ thermal niches by overlaying global species breeding distributions and world temperature metrics (based on mean, minimum, maximum and range), and then aggregated them to obtain a set of community thermal indices for each assemblage (CTIs, and CTR for ranges). Long-term average temperatures and local current temperatures were related to changes in CTIs using spatial GLMMs, which considered habitat change. We identified the species most responsible for variation in assemblages and regressed species’ influence on thermal affinities and population dynamics. Results: CTIs increased with temperature and warm-dwelling species became more prevalent to the detriment of cold-dwelling ones. However, we found a counteracting effect of temperature and habitat. Cold-dwelling forest species were among the most influential species, mainly through colonisation, while warm-dwelling farmland species contributed through local extinctions (both attenuated local increases in CTI). The mean thermal breadth of assemblages (CTR) decreased with temperatures. Main conclusions: The taxonomic composition of bird assemblages shifted in line with the main expectations due to global change (thermophilisation), mainly due to local colonisation of warm-dwelling species, although it did not show the pattern of thermal homogenization suggested elsewhere. Our results add further evidence of the interplay between climate warming and land-use change in the ongoing adjustment of animal communities.

# Biodegradable microplastics can cause more serious loss of soil organic carbon by priming effect than conventional microplastics in farmland shelterbelts ## Description of the data and file structure ### dataset Raw data.csv Soil data for the variables tested in the paper. Variables are as follows: * Group = The number of the microplastics addition and control group * First emission of soil CO2 (mg g-1 SOC) = CO2 release rate of soil organic carbon at first sampling * δ13 C of CO2 in first samping(‰) = The δ13 C of CO2 at first sampling * DOC(mg kg-1) = Dissolved organic carbon content of soil samples after incubation * MBC(mg kg-1) = Microbial biomass carbon content of soil samples after incubation * DTN(mg kg-1) = Dissolved total nitrogen content of soil samples after incubation * MBN(mg kg-1) = Microbial biomass nitrogen content in soil samples after incubation * NH4+-N(mg kg-1) = Content of ammonium nitrogen in soil samples after incubation * NO3--N(mg kg-1) = Nitrate nitrogen content in soil samples after incubation * pH = pH value of soil sample after incubation * null = No data were obtained. In this study, conventional microplastics exhibited no degradation during the short-term incubation; therefore, 13C isotopes were not employed to differentiate soil CO2 emissions in the treatment group involving conventional microplastics. One-way analysis of variance (ANOVA) examined the differences in soil-derived CO2 emission, SOC loss, hydroxyl index (HI), soil physiochemical properties and microbial characteristics of different soils and MPs groups (P < 0.05). The above experimental data were conducted using SPSS 27.0. Structural equation model was analyzed using Amos 26.0 software to explore the pathways of MPs addition on cumulative soil-derived CO2 emissions. Globally, the widespread utilization of plastic products has resulted in the accumulation of microplastics (MPs) in the soil. MPs have the potential to impact the loss of soil organic carbon (SOC). Nevertheless, the influence of different types of MPs on SOC loss remains uncertain. In this study, a 38 d’ incubation experiment with two kinds of conventional MPs (polyethylene (PE), polypropylene (PP)) as well as two kinds of biodegradable MPs (polyhydroxyalkanoate (PHA), polylactic acid (PLA)) were added into three types of soil (loam, sandy loam, and sandy soil) in farmland shelterbelts, and the sources of CO2 emissions was distinguished by the difference in 13C isotope abundance between the biodegradable MPs (PHA and PLA) (-10.02 ~ -9.92 ‰) and the soil (-24.39 ~ -22.86 ‰) (>10‰). In conjunction with the structural characterization of MPs, as well as soil physicochemical properties and microbial characteristics, we observed that the conventional MPs did not degrade in short term incubation, but significantly enhance soil-derived CO2 emissions by altering the dissolved N content (NH4+-N and DTN) and reducing microbial biomass carbon (MBC) content only in sandy loam soil (P<0.05). Biodegradable MPs degraded significantly, and enhanced soil-derived CO2 emissions by reducing soil dissolved total N (DTN) and NO3--N contents in loam, sandy loam and sandy soil (P<0.05). Overall, the input of biodegradable MPs causes a more serious loss of SOC than conventional MPs as the soil sand content increased in short term incubation, which needs to be considered in predicting the global impact of increasing biodegradable MPs pollution.