The speed of change in climate and land cover is associated with the speed of biodiversity changes in avian assemblages of the United States

PLOS One | https://doi.org/10.1371/journal.pone.0330153 August 29, 2025 4 / 20 selected climatic variables known to limit the distribution of birds physiologically [50–52]: mean precipitation (mm) minimum temperature(°C), and maximum temperature (°C). Each year, we calculated the mean values of these variables for the breeding season (i.e., May 1st to July 31st) using the methodology described in Supporting Information (S1 File). LULC data. To calculate the rate of LULC change between 1990 and 2019, we utilized the “Modeled Historical Land Use and Land Cover for the Conterminous US: 1938-1992” [53] and the scenario A1B of the “Conterminous United States Land Cover Projections - 1992 to 2100” [54] data sets. A detailed description of these data sets and the rationale for the selection of scenario A1B is provided in the Supporting Information (S1 File). We grouped the annual layers and retained for the analysis six combined layers containing the proportion of LULC per map cell for: barren land, crop/pasture, forest (deciduous forest + evergreen forest + mixed forest), grassland/shrubland, urban, and wetland (herbaceous wetland + woody wetland). S1 Fig in S1 File shows the methodology we used to obtain the grouped annual LULC layers. Elevation data. We included elevation in the analysis to account for its potential effect on biodiversity gradients. We used the North American Elevation map at 1 km resolution [55] to extract the elevation of each bird assemblage. Elevation

information is given in meters, with values smaller than/ at/ above zero indicating below/ at/ above sea level, respectively. To work with all data sets in the same spatial projection, we re-projected the environmental and biodiversity data to the Lambert Azimuthal Equal Area projection. We joined the annual climatic and land conversion layers and extracted the environmental information for each bird assemblage geographic coordinates. Rate of change in environmental conditions For each avian assemblage, we extracted the mean value of each annual climate and LULC map using three buffer sizes: 12.5 km, 25 km, and 50 km. We performed this step with the aim of accounting for potential effects that the scale of the analysis could have on the relationship between the rate of change in environmental variables and the rate of change of richness and assemblage dissimilarity [56]. At each buffer size, we calculated the rate of change of each environmental variable at each BBS route as the slope of the linear regression of the mean variable value with time (i.e., mean value ~Year) for the period 1990–2019 (Fig 1B II.). A positive slope indicated an increase in mean precipitation, minimum or maximum temperature, or proportion of LULC with time. A negative slope indicated a decrease in mean precipitation, mini- mum or maximum temperature, or proportion of LULC with time. Effect of the rate of environmental change on the biodiversity of species assemblages We built a generalized additive mixed-effect regression model (GAMM) for each assemblage dissimilarity metric (i.e., total

assemblage dissimilarity, species replacement, and species loss) and species richness in each ecoregion (gamm, “mgcv” R package, [57]. This type of model is frequently used in ecology to model non-linear and non-monotonic relationships [58]. GAMMs were fitted with a Gaussian distribution of errors and an identity link. We utilized GAMMs instead of Gener- alized Additive Models to account for spatial autocorrelation among the data points. The bird data were collected across space, which can result in biodiversity indices calculated from routes close together being more similar than those from routes further apart. To account for the possible spatial autocorrelation, we included a Gaussian spatial correlation struc- ture of the form ~ Latitude + Longitude in the models. We used the rates of change of the assemblage dissimilarity metrics (log-transformed) and species richness as the response variables and the rate of change of climatic and LULC variables as predictors (all as smooth terms in the GAMMs). Dissimilarity metric values were log-transformed to standardize the variance. For each GAMM in each ecore- gion, we reduced the number of predictor variables by excluding the LULC variable with the lowest mean proportion of LULC across the study period and the variables showing a correlation coefficient greater than 0.6. We included elevation