Introduction

Ongoing human alteration of the Earth surface causes changes in biodiversity across scales (Gibson et al. 2011; Murphy and Romanuk 2014; Newbold et al. 2015). Globally, about 32% of all known vertebrate species show decreasing population sizes or range contractions (Ceballos et al. 2017; WWF 2018) with reported species extinction rates being several times higher than expected naturally (Brooks et al. 2002; Pimm et al. 2014). Yet, any change in biodiversity is scale and measure dependent (Sax and Gaines 2003; Chase and Knight 2013) and, perhaps surprisingly, whether local—opposed to global—biodiversity is declining is still debated (Thomas 2013; McGill et al. 2014).

Several global meta-analyses have demonstrated that some biodiversity measures, notably species richness, have not declined at local scales, e.g. the scale of biodiversity sampling (Vellend et al. 2013, 2017; Dornelas et al. 2014). However, these results have been questioned, particularly whether the data are spatially and temporally biased (Gonzalez et al. 2016) or whether sites with and without land change were differentiated (Cardinale et al. 2018). Changes in land use and land cover (LULC) have been identified as one of the dominant driver of terrestrial biodiversity loss (Díaz et al. 2019). This raises the question whether such changes on land can explain changes in local biodiversity measures across space and time.

Present differences in LULC influence local biodiversity globally. Studies have shown local biodiversity to be consistently reduced at sites with more intensively used land (Murphy and Romanuk 2014; Newbold et al. 2015; Alroy 2017), where on average 13.6% fewer species and 10.7% fewer individuals were observed compared to undisturbed primary vegetation (Newbold et al. 2015). However, these analyses relied on spatial comparisons of local biodiversity across land uses, therefore did not capture temporal biodiversity change. In addition, these studies ignored the influence of past changes in LULC at local (Midolo et al. 2018; Jung et al. 2019b) and landscape scales (Cousins et al. 2015), which can influence local biodiversity (Tscharntke et al. 2012; Turner and Gardner 2015; Miguet et al. 2016).

At the landscape scale, defined as the extent at which spatio-temporal dynamics affect ecological processes (Turner 1989; Pickett and Cadenasso 1995), local biodiversity is influenced by the variability of resources, such as food or nesting material, or by ecological processes, such as migration or fear of predation (Hanski and Ovaskainen 2000; Chase 2003; Turner and Gardner 2015; Fernández et al. 2016). However, these influences are dynamic as landscapes are constantly changing because of natural and anthropogenic factors (Pickett and White 1985; Manning et al. 2009; Turner and Gardner 2015). Previous studies have shown that landscape-wide changes in LULC may have a lasting influence on local biodiversity through ‘biotic lag’ effects (Metzger et al. 2009; Ewers et al. 2013). Yet, most studies focussed on small geographic regions and changes in forest cover (Rittenhouse et al. 2010) and did not investigate general impacts of landscape-wide changes in LULC on local biodiversity across spatio-temporal scales. A lack of data on local biodiversity change and landscape-wide LULC change has so far prevented comparative assessments (De Palma et al. 2018).

Increasing availability of satellite imagery enables the assessment of changes in LULC at broad spatial and temporal scales (Kennedy et al. 2014; Pasquarella et al. 2016). Long-term satellite missions, such as NASA’s Landsat, provide one of the best sources of time series to monitor land surface conditions (Kennedy et al. 2014; Vogelmann et al. 2016; Hermosilla et al. 2018; Song et al. 2018). Time series of land surface conditions, such as photosynthetic activity, can measure intra- and inter-annual vegetation dynamics (Pettorelli et al. 2005; Fisher et al. 2006). Changes in photosynthetic activity are directly linked to changes in LULC (DeVries et al. 2015; Jung et al. 2019b) and can be differentiated by attributes (Watson et al. 2014). For example, a loss or gain of vegetation cover causes abrupt changes in the magnitude of photosynthetic activity (Kennedy et al. 2010; DeVries et al. 2015), while greening or browning lead to increasing or decreasing trends in photosynthetic activity (de Jong et al. 2013; Mueller et al. 2014). These attributes can be robustly quantified at the landscape scale and related to changes in local biodiversity.

Birds are one of the best-surveyed taxonomic groups globally. Local biodiversity change quantified from repeated breeding bird surveys (BBS) has been widely studied (Harrison et al. 2014; Pardieck et al. 2018). Previous studies have shown that changes in bird diversity are dependent on the specific biodiversity measure considered (Schipper et al. 2016; Jarzyna and Jetz 2017) and are often non-linear (Gutzwiller et al. 2015; Barnagaud et al. 2017). Change in bird diversity varied across ecoregions and for birds of particular functional traits (Harrison et al. 2014; Schipper et al. 2016; Jarzyna and Jetz 2017), such as migratory or grassland dependent species, which declined notably in developed countries (Fewster et al. 2000; Sanderson et al. 2006; Stanton et al. 2018). Changes in LULC are most likely a driving factor of such declines, although most studies investigated only spatial correlations in remotely-sensed photosynthetic activity and local bird diversity (Rowhani et al. 2008; Goetz et al. 2014; Hobi et al. 2017). Notably Rittenhouse et al. (2010) found bird assemblage composition to be altered in landscapes with more “disturbed forests”, which they assessed using remotely-sensed time series. However, to our knowledge, no previous study has investigated whether landscape-wide changes in LULC in general correlate with and explain changes in local bird diversity.

The aim of this paper is to investigate whether changes in preceding and concomitant LULC, quantified as changes in remotely-sensed photosynthetic activity, influence bird diversity change. Previous studies have found differences in biodiversity between sites to be linked to changes in LULC at the local (Jung et al. 2019a, b) and landscape scale (Ewers et al. 2013; Cousins et al. 2015). Yet, it is unknown to what extent preceding and concomitant changes in LULC are linked to biodiversity change. The notion of both biodiversity and landscapes being dynamic can assist in a better understanding of the drivers of biodiversity change (Manning et al. 2009), and help inform conservation actions (Essl et al. 2015a) and landscape management practices (Lindenmayer et al. 2008; Reed et al. 2016). Consequently, this study hypothesizes that (i) changes in local bird diversity are driven by landscape-wide changes in LULC, (ii) incorporating lagged changes in LULC that occurred before bird diversity sampling increases explanatory power, and that (iii) the explanatory power of landscape-wide changes in LULC on local bird diversity change varies across ecoregions and functional groups of bird species. We combine 34 years (1984–2017) of annual BBS records collected at sites across the conterminous United States of America with time series of medium–high resolution (nominal ~ 30 m) satellite imagery from the Landsat missions. Using Breaks for Additive Seasonal and Trend (BFAST), a generic change detection algorithm, we detect changes in LULC as trends (e.g., shifts in greening or browning) in, or abrupt shifts in magnitude (e.g., immediate gain or loss in photosynthetic activity) of photosynthetic activity in the landscape surrounding each BBS route. Non-linear spatio-temporal models were used to correlate these landscape-wide changes in photosynthetic activity with changes in local bird diversity.

Methods

Bird diversity time series

Time series of local bird count records (1984–2017) were obtained for the conterminous USA (excluding Alaska, Hawaii and Puerto Rico) from the North American Breeding Bird Survey dataset (BBS, available from https://www.pwrc.usgs.gov/bbs/, Pardieck et al. 2018). Bird counts were conducted annually during the breeding season (April to August with > 83.3% routes surveyed in June) along approximately 39.4 km long roadside survey routes and usually follow a standard protocol that involves fifty 3 min stops at evenly spaced intervals (approximately 0.8 km) (Ralph et al. 1995). At each 3 min stop, volunteer observers record the number and identity of every bird species seen or heard within approximately 400 m distance from the route. For our analyses we included routes that followed the standard BBS protocol of fifty randomly selected stops (94.4% of all routes) and had at least ten years of sampling between 1984 and 2017, as many BBS routes were not sampled every year (mean proportion of missing years = 19.7%). The period from 1984 to 2017 was chosen to align with the availability of satellite data (but see “Time series of annual photosynthetic activity at the landscape scale”). We removed routes from the analyses with non-acceptable weather conditions according to BBS standards (Ralph et al. 1995) and excluded all nocturnal, crepuscular and aquatic species from the analyses as they are not well sampled by BBS methods (Gutzwiller et al. 2015; Jarzyna and Jetz 2017). All partially identified species (e.g. “sp.”), hybrids and species with unclear taxonomy (e.g. “A × B”) were removed from further analyses. In total, time-series from 2745 routes (out of 5248 in the BBS in conterminous USA) had suitable data for further analyses.

We calculated two different biodiversity measures commonly applied to BBS data. First, we calculated the geometric mean of relative abundance (GM), a composite indicator, which quantifies relative changes in both abundance and evenness, with the latter being affected even if overall abundance is not changing (Buckland et al. 2011, 2017; Harrison et al. 2014). The GM for the year y is defined as \(GM_{y} = \exp \left( {\frac{1}{S} \mathop \sum \limits_{i = 1}^{S} \log \left( {\frac{{A_{iy} + 1}}{{A_{{iy_{0} }} + 1}}} \right)} \right)\), where S quantifies the total number of species with i being an individual species, \(A_{iy}\) the abundance of species i in year y. The first four years of BBS data (1984–1987) were used to define the baseline years y0 (calculated from the median number of individuals for each observed species) and to align the analyses with the baseline years used in the land change detection (but see “Detection of landscape-wide changes in annual photosynthetic activity”). Whenever no BBS was conducted in the years between 1984 and 1987 on a given route, we used the first year of available BBS data to define the baseline year y0. The GM is unaffected by species detectability as it is based on within-species abundance trends, however it cannot be quantified for absent species and is unable to reflect changes in assemblage composition (Buckland et al. 2011). We added a constant (1) to all abundance values before calculating the annual GM to account for the species being absent in some years.

Second, as measure of change in assemblage composition, we calculated the progressive Bray–Curtis index (pBC,Bray and Curtis 1957; Rittenhouse et al. 2010) as the difference in composition between a baseline and all following years of sampled bird assemblages (Rittenhouse et al. 2010). The pBC is defined as \(1 - \mathop \sum \limits_{{{\text{i}} = 1}}^{{\text{S}}} \frac{{\left| {{\text{A}}_{{{\text{iy}}}} - {\text{A}}_{{{\text{iy}}_{0} }} } \right|}}{{\left( {{\text{A}}_{{{\text{iy}}}} + {\text{A}}_{{{\text{iy}}_{0} }} } \right)}}\), where S, Aiy and \({\text{A}}_{{{\text{iy}}_{0} }}\) are defined as above for GM. The pBC ranges between 1 and 0, with increasing values closer to 0 indicating a less similar bird assemblage composition relative to the baseline years y0.

Time series of annual photosynthetic activity at the landscape scale

Following previous studies, we define the “landscape” as the buffer with 19.7 km radius around the centroid of each BBS route which fully encompasses the majority of BBS routes and approximates the median natal dispersal distance of North American bird species (Sutherland et al. 2000; Pidgeon et al. 2007; Albright et al. 2011). To quantify changes in land use and land cover (LULC) in the landscape surrounding each BBS route, we calculated changes in photosynthetic activity using imagery from the Landsat 4, 5, 7 and 8 satellites (1984 to 2017, ~ 30 m nominal resolution) supplied by the United States Geological Service (USGS) available through Google Earth Engine (Gorelick et al. 2017). All Landsat images were radiometrically (Chander et al. 2009) and atmospherically calibrated to surface reflectances (Masek et al. 2006).

For each surface reflectance image, we masked out clouds and cloud shadows as identified by the ‘cFMask’ algorithm (Zhu and Woodcock 2012) and areas permanently covered with water, snow or ice according to a mask derived from the 2011 National Land Cover Database (NLCD) land-cover map at ~ 30 m resolution (Homer et al. 2015). BBS routes with less than 50% land area (N = 18) within the surrounding landscape were excluded from further analyses, assuming that breeding birds are less influenced by terrestrial landscape-wide changes in LULC at these routes. Furthermore, a lack of clear-sky images in certain years can lead to missing data for parts of the buffered BBS routes. Routes with more than 50% missing data (N = 2) from 1984 to 2017 were excluded from further analyses assuming that the Landsat satellites have missed most changes in LULC (median proportion of missing data = 1.06% ± 1.54 median absolute deviation [MAD], Fig. S2).

A spectral index of photosynthetic activity (the two-band enhanced vegetation index [EVI], Jiang et al. 2008) was calculated for each surface reflectance image. We composited all EVI data up to three months before the summer solstice (20th of March to 20th June of each year) into a single annual image (1984, 1985,…, 2017) retaining the greenest (95% percentile) EVI value. A 95% percentile instead of maximum was used to reduce the influence of extreme outliers. We used three months of EVI data to capture the greening onset in annual vegetation dynamics (Fig. S1), a period that can assist in distinguishing between land cover types (Pettorelli et al. 2005; Fisher et al. 2006; Zhang et al. 2006), and that matches the sampling period and conditions during which most BBS were conducted (March to June). All data pre-processing and compositing was done using the Google Earth Engine platform (Gorelick et al. 2017).

Detection of landscape-wide changes in annual photosynthetic activity

We quantify the proportion of grid cells within each 19.7 km buffer showing a trend in or an abrupt change in magnitude of photosynthetic activity as measured by EVI (Fig. 1d). Among all algorithms proposed to detect changes in remotely-sensed time series (Zhu 2017), we relied on the generalized fluctuation framework originally developed for econometrics (Bai and Perron 2003; Zeileis 2005), later adapted for remote sensing as the Breaks for Additive Seasonal and Trend (BFAST) algorithm (Verbesselt et al. 2010). For each annual EVI time series, we tested for single or multiple structural breaks in linear trend using a recursive Moving Sum of Residuals (Rec-MOSUM) test over each four year window period (Zeileis 2005). A statistically significant (p < 0.05) structural change test indicates whether at least a single structural break exists, in which case we iteratively fitted segmented linear regression models over the entire time series. The optimal number and position of all structural breaks were detected by minimizing both the Bayesian Information Criterion (BIC) and residual sum of squares (RSS) of the segmented regression models (Zeileis 2005; Verbesselt et al. 2010). The framework requires a gap-free time series (“strucchange” package in R, ver. 1.5-1) and similar to previous studies we filled missing data using linear interpolation between adjacent years (Verbesselt et al. 2010).

Fig. 1
figure 1

a Location of breeding bird survey (BBS) routes (grey lines) across the United States of America. b Hypothetic changes in LULC within a single grid cell in the landscape (buffered circle around the BBS route). For each grid cell within the landscape, time series of annual summarised March–June EVI were tested for a single or multiple changes in LULC (see “Detection of landscape-wide changes in annual photosynthetic activity”), defined as either a trend (greening [“dark green”] or browning [“brown”]) in or an abrupt change in the magnitude (abrupt loss [red line] or gain [“blue”]) of photosynthetic activity as measured by the EVI. For each route we summarize c changes in local bird diversity (as quantified by the geometric mean of relative abundance [GM] and by progressive Bray–Curtis index [pBC]) relative to a baseline year y0 (highlighted in grey) for an example BBS route; and d the proportion of all grid cells within the landscape (buffer with 19.7 km radius) with either a trend in or an abrupt change in the magnitude of EVI (colours as in a) per year. Map shown in Albers equal area conic projection (NAD83). Silhouette from phylopic.org released as public domain

Full size image

Per grid cell and year, we differentiated all the detected change events in photosynthetic activity as either a trend in or an abrupt change in the magnitude of EVI (Fig. 1). Abrupt changes in magnitude were quantified using the predicted EVI data (from the segmented linear regression model) before and after the detected change date \(({\text{EVI}}_{{{\text{After}}}} - {\text{EVI}}_{{\text{Before }}} )\) and categorized as either immediate loss or gain in photosynthetic activity in a given year if negative or positive, respectively. For trends, we assessed for each year whether the linear trend in annual EVI was significantly (p < 0.05) increasing (‘greening’), decreasing (‘browning’) or flat (‘stable’). Similarly, for time series with non-significant structural change tests, we fitted simple linear regression models to test whether the overall trend in EVI (across all 34 years) significantly increased or decreased.

For each landscape around a BBS route and year (Fig. 1b), we summarized the amount of land area that had either a trend (greening or browning) in or an abrupt change in magnitude (loss or gain) of EVI. Because the total land area differed among BBS routes, we calculated proportions relative to the total land area (see “Time series of annual photosynthetic activity at the landscape scale”). The change detection algorithm relies on a moving window (four years) and thus no change events could be detected in the first (1984–1987) and last four (2014–2017) years of each EVI time series. If a change event occurred within these years, the algorithm would set the date to the latest or earliest year possible (i.e. 1987 and 2014, respectively) causing an inflated number of incorrectly dated change events at the start and end of each time series. We therefore considered the first four years as ‘baseline’ (year0) and the last four as ‘overhang’ and removed them from further analyses.

Additional predictors and bird trait data

At continental scales, bird diversity at BBS routes has been shown to be influenced by a number of environmental variables (Rowhani et al. 2008; Goetz et al. 2014; Hobi et al. 2017; Barnagaud et al. 2017). For a coarse measure of overall vegetation activity (Rowhani et al. 2008; Hobi et al. 2017), we calculated the mean EVI across all 34 years of annual Landsat composites per buffered BBS route (see “Time series of annual photosynthetic activity at the landscape scale”). Previous studies have shown that the number of bird species varies with elevation (Jarzyna and Jetz 2017), therefore we included the mean elevation of the buffered BBS route from the global GMTED (~ 1 km resolution) product (Danielson and Gesch 2011). As precipitation-driven anomalies have been shown to affect the number and abundance of bird species (Barnagaud et al. 2017), we used the Standardized Precipitation-Evapotranspiration Index (SPEI), which quantifies anomalies relative to the conditions observed in a moving window before a given month (Vicente-Serrano et al. 2010, 2012). For each BBS route we extracted the average monthly SPEI from SPEIbase (ver. 2.5, https://spei.csic.es, (Vicente-Serrano et al. 2010)) calculated on a climatology from 1901 to 2015 and over a moving window of three months from January to March of each year (Vicente-Serrano et al. 2010), thus capturing precipitation anomalies in the winter months. We used the winter months in order to avoid potential collinearity with the EVI, while capturing winter conditions that might affect bird survival.

Similar to previous studies we used four functional trait groups—nesting status, migratory behaviour, habitat guild and body mass—to differentiate all bird species (Schipper et al. 2016; Barnagaud et al. 2017). Data on nesting (ground or canopy) and migratory behaviour (resident, short-distance and neotropical migrants) were obtained from Albright et al. (2011), while data on bird species habitat guilds (e.g. woodland, shrubland, grassland and urban birds) were extracted from the USGS website (https://www.mbr-pwrc.usgs.gov/bbs/guild/guildlst.html, accessed on 28/08/2018). The mean body mass (bm, measured in g) for all bird species was extracted from the Amniote database (Myhrvold et al. 2015) and grouped into terciles of all estimates, e.g. small, medium and large birds. For species without trait estimates, we filled the missing data with the most common (mode) trait within the same bird genus, provided more than 50% of all species within that genus had existing body mass estimates or identical categorical trait. For each BBS route and trait group we calculated separate GM estimates, where routes with at least 10 years of data and at least three different species within a trait group were available.

Spatio-temporal models

The aim of the statistical analyses was to investigate whether changes in local bird diversity (measured by GM and pBC) and landscape-wide changes in LULC, quantified as changes in remotely-sensed photosynthetic activity, are correlated. To do so we relied on generalized additive regression models (GAMs), which are commonly used to model species population trends (Fewster et al. 2000) and can handle complex non-linear, spatio-temporal and hierarchical datasets (Kneib et al. 2009; Wood 2011). All considered variables were included as thin-plate smooth (fixed to 4 residual degrees of freedom to prevent overfitting) in the GAMs and we applied a smoothing penalization for variable selection (Wood 2008, mgcv parameter: select = TRUE). The approximate significance of non-linear model terms was assessed using an approach by Wood (2013). All GAMs were fitted using the ‘mgcv’ package (Wood 2011, ver. 1.8-24) in R (R Core Team 2018, ver. 3.5.0).

We distinguished between four groups of variables to be included as thin-plate smooths in the full GAM. (1) As “environmental predictors” factors (fenvironment) we considered the mean EVI, elevation and, for each year, the SPEI. (2) For landscape-wide changes in LULC (flandscape), we included for each year the proportion of grid cells within the landscape (arcsine square root transformed) showing a trend (browning or greening) in or an abrupt change in the magnitude (immediate loss and gain) of EVI (Fig. 1d). (3) Incorporating spatial autocorrelation into regression models can improve predictive power (Kneib et al. 2009; Dornelas et al. 2012), especially when local biodiversity was surveyed over large scales such as the conterminous United States of America. We followed an approach by Kneib et al. (2009) and included the spatial coordinates (fspatial) of each BBS route using a non-linear smooth surface function \({\text{g}}\left( {{\text{x}}_{{{\text{Northing}}}} ,\;{\text{x}}_{{{\text{Easting}}}} } \right)\) with a tensor product P-spline. Northing and easting coordinates were obtained by projecting the centroid of each buffered BBS route to an Albers equal area conic projection (NAD83). (4) We included the BBS route ID (fID) as random intercept in all models, therefore estimating the effect of fenvironment, flandscape and fspatial on local biodiversity measures (GM and pBC) across all BBS routes, to account for varying species detectability and misidentifications (Sauer et al. 1994; Harris et al. 2018). We acknowledge that using the route ID does not fully account for differences in observer abilities (there can be multiple observers for a single route), but previous studies found limited influence of varying observers over large scales (Jarzyna and Jetz 2017; Barnagaud et al. 2017). All biodiversity time series were detrended by including time (year) as a linear predictor to avoid spurious correlations. To account for temporal autocorrelation, we included an autoregressive error structure (AR1), which we parametrized by visually assessing the autocorrelation function of the full model residuals at lag 1 (ρ = 0.5).

We tested if preceding (e.g. the years before a BBS) landscape-wide changes in LULC continued to influence bird diversity change in subsequent years. A ‘lagged’ correlation between two time series is commonly known as “Granger causality”, where one “time series xt contains information in past terms that helps the prediction of yt” (Granger 1969). We followed an approach by Papagiannopoulou et al. (2017) and assessed the relative improvement in explanatory power of models including preceding instead of concomitant changes in LULC. Preceding abrupt changes in the magnitude in EVI (loss or gain) of up to five years were included either individually, thus adding estimates for the preceding year i = 1,…, 5 only; or cumulatively, where aggregated estimates for the preceding years 1:i were included in the model (Jung et al. 2019a). The relative improvement in explanatory power was assessed using out-of-bag (OOB) coefficients of determination (R2). To do so we split all time series into training and test datasets (50/50) 100 times at random. All models included the fenvironment and fspatial variables to account for variation not directly attributable to landscape-wide changes in LULC.

Lastly, we assessed the explanatory power of each group of variables (fenvironment, fspatial, flandscape) spatio-temporally and for birds grouped by functional traits. To do so we fitted several GAMs using the GM (log-transformed) or pBC as response variable with a gaussian log-link distribution. We first fitted a “full” GAM including all variables, followed by separate GAMs where groups of variables (fenvironment, fspatial, flandscape) were explicitly excluded from the model. Models for both GM and pBC converged well (Fig. S3–4), although the largest changes in assemblage composition were generally poorly predicted by the models (Fig. S4). The explanatory power of all models was assessed by calculating the R2 of each model. The group of variables (fenvironment, fspatial, flandscape) explaining the most variation was then identified from the largest reduction (partial R2, relative to the full model) in R2 (Papagiannopoulou et al. 2017). We assessed patterns of the most important group of variables spatially and in relation to robust linear trends in biodiversity measures (fitted using the MASS package, ver. 7.3-49, (Venables and Ripley 2002)). Lastly, we investigated if the explanatory power of landscape-wide changes in LULC (flandscape) varied either with functional trait groups (see “Additional predictors and bird trait data”) or with BBS routes grouped by U.S. ecoregions (Level 1, Omernik 1987). For each functional trait group, we fitted two separate GAMs either including or excluding all flandscape variables before calculating the difference in R2 attributable to flandscape variables. For U.S. ecoregions we assessed the contribution of flandscape variables to the total R2 on overall GM change.

Results

Both local bird diversity and landscapes have changed across the conterminous USA. Across all BBS routes the geometric mean of relative bird abundances (GM) increased by 0.01% ± 0.002 standard error (SE) per year (mean first derivative) in the first two decades from 1984 to 2005, after which annual decreases of 0.01% ± 0.003 SE were observed (Fig. S5a). The compositional similarity of bird assemblages (pBC) decreased by 0.006% ± 0.001 SE per year (Fig. S5b). Landscapes surrounding each BBS route had on average 6% ± 6.42 SD (range 0.02–78.96%) of land area experiencing at least one land change in the period 1984 to 2017 (Fig. S6). Over the same period a decrease in landscape-wide changes in LULC were observed (mean robust linear trend = − 0.00015 ± 0.0198 SD, range -0.545 to 0.112) but with large spatial variability (Fig. S7). Across all BBS routes the mean proportion of land with an abrupt change in the magnitude of EVI (loss or gain) fluctuated strongly (Fig. S8a), while the mean proportion of grid cells showing a trend in EVI (significant increase or decrease) showed an inverse hump-shaped pattern for greening and a continuous decrease for browning (Fig. S8b). The proportion of landscape-wide abrupt changes in magnitude or trend were not correlated among BBS routes (Fig. S9) and across ecoregions (Fig. S10).

Bird diversity change is correlated with landscape-wide changes in LULC. The GM significantly decreased in years with a large proportion of landscape-wide abrupt gains of EVI (F4 = 10.8, p < 0.001; Fig. 2a, blue line). More landscape-wide abrupt losses of EVI led to a significant decrease in GM, but only after ~ 10% of the landscape had abrupt losses in a given year (F4 = 6.44, p = 0.001, Fig. 2a, red line). The GM also decreased with more land in the landscape browning (F4 = 37.89, p = 0.057, Fig. 2a), while greening land in the landscape had no significant effect on changes in GM (F4 = 0, p = 0.529, Fig. 2a). The pBC significantly increased with a large proportion of landscape-wide abrupt losses (F4 = 8.25, p < 0.001, Fig. 2b) or gains in EVI (F4 = 0.614, p = 0.1, Fig. 2b). The pBC also increased with a large proportion of browning (F4 = 13.81, p = 0.038, Fig. 2b) or greening land in the wider landscape (F4 = 74.25, p = 0.005, Fig. 2b).

Fig. 2
figure 2

Partial effects of landscape-wide changes in LULC, quantified as changes in photosynthetic activity (proportion of landscape) per unit of change in a the geometric mean of relative abundance (GM) and b the progressive Bray–Curtis index (pBC). Colours indicate either abrupt shifts in magnitude with losses (red lines) or gains (blue) in EVI or trend with greening (green) or browning (brown) land. Error margins show the estimated standard error of the partial effect (grey shading). Flat lines without error margins indicate that the term was penalized out during model fitting and thus had no effect on the biodiversity measure

Full size image

Other environmental factors strongly influenced local bird biodiversity change. The GM significantly increased (F4 = 1789.06, p < 0.001, Fig. S11a) and the pBC significantly decreased (F4 = 1923.71, p < 0.001, Fig. S11b) in landscapes with high mean elevation. GM significantly increased (F4 = 291.05, p < 0.001) in landscapes with overall low photosynthetic activity (EVI < 0.4) but decreased in landscapes with high photosynthetic activity; a pattern that was reversed for pBC (Fig. S11). Anomalies of precipitation between January and March as measured by the SPEI had no significant effect on GM or pBC change (Fig. S10).

Changes in LULC in one year continued to influence local bird diversity in subsequent years. The mean explanatory power (out-of-bag [OOB] R2) of abrupt shifts in magnitude in concomitant years (Lag 0, Fig. 3) was 0.03 (0.047 cumulatively) for GM and 0.126 (0.122) for pBC. A consideration of abrupt shifts in magnitude in preceding years explained modestly more variation than those in concomitant years (Fig. 3). Including one to five preceding years of abrupt shifts in magnitude individually explained similar amounts of variation (mean OOB R2 = 0.031) in GM, whereas for pBC only preceding abrupt shifts in magnitude more than three years increased the explanatory power (mean OOB R2 = 0.129, Fig. 3). Considering cumulatively preceding abrupt shifts in magnitude increased the mean explanatory power for both GM (Fig. 3), while for pBC the relative improvement in explanatory power was highest at three cumulatively included preceding years (mean OOB R2 of year three = 0.128).

Fig. 3
figure 3

Preceding landscape-wide changes in LULC, quantified as changes in photosynthetic activity, of one to five years improve predictions of mean relative bird abundance (GM) and bird assemblage composition (pBC). Bird diversity time series for all BBS routes were randomly split (100 times) into training and test datasets and the explanatory power (R2) was assessed relative to a model that only included concomitant abrupt shifts in magnitude (gain or loss of EVI) averaged across all random subsets. Past changes in LULC were either added to models individually (circles) or aggregated cumulatively (triangles). Error bars show the standard deviation of the out-of-bag (OOB) R2 values

Full size image

We assessed whether the explanatory power of all variables varied in space (Fig. 4a, c) and for linear trends of bird diversity change (Fig. 4b, d). The full model including all variables explained 64.7% of the total variation of changes in GM (69.3% for pBC), with most of the variation explained by unknown differences among BBS routes (partial R2 of fobs = 58.5%. for GM and 39.8% for pBC). Of all variables considered, landscape-wide changes in LULC were the most important predictor of GM change for 34.8% of BBS routes (partial R2 range 0–54%, Fig. 4a) and for pBC in 46.6% of BBS routes (partial R2 range 0–7%, Fig. 4c). Incidentally, landscape-wide changes in LULC were the best predictor for some of the greatest changes (increase/decrease per year) in local bird diversity measures (Fig. 4b, d).

Fig. 4
figure 4

Most important group of variables explaining changes in the a GM and c pBC at 2745 BBS routes across the conterminousUnited States of America. Point sizes in (a, c) indicate larger partial R2 of the most important variable group. Colours indicate which of the considered variable groups, fenvironment (red), flandscape (blue) or fspatial (green), explained most of the variation (greatest partial R2) in the full model. b, d Partial R2 of the most important variable group averaged per increase or decrease (robust linear trend per year) in GM or pBC

Full size image

The explanatory power of landscape-wide changes in LULC on changes in GM differed among bird species of varying functional traits and across ecoregions (Fig. 5, Fig. S12). On average landscape-wide changes in LULC did not explain (mean partial R2 = -0.02 ± 0.09 SD) changes in GM for birds of varying trait groups (Fig. 5a). Similar to spatial patterns of the most important group of variables (Fig. 4), landscape-wide changes in LULC were only important for a subset of BBS routes (Fig. 5a, blue outliers) in which they explained up to 71.9% of the total R2. For many BBS routes the inclusion of landscape-wide changes in LULC did not increase but decreased the R2 for explaining changes in GM (Fig. 5a, red outliers). A visual exploration could not identify any spatial patterns in these outlier BBS routes and there were no visually distinguishable differences between ecoregions (Fig. 5b) and especially in Southern Semi-Arid Highlands landscape-wide changes in LULC did not increase the explained variation in GM change, despite the on average large proportion of browning land (Fig. S10).

Fig. 5
figure 5

a Partial R2 of landscape-wide changes in LULC (difference in explanatory power after excluding flandscape variables) for explaining changes in GM grouped by functional trait group. Shown is the distribution (grey), median and 50% quantile (black) for each response variable. Red and blue points indicate outliers (1% smallest/biggest partial R2 values). b Absolute partial R2 of landscape-wide changes in LULC explaining trends in GM grouped by U.S. ecoregions. Coloured depending on whether landscape-wide changes in LULC increased (blue) or decreased (red) overall R2. Black points and error bars show mean and standard error of the Partial R2

Full size image

Discussion

This study investigated whether changes in land use and land cover (LULC), approximated as changes in trend in or abrupt changes in magnitude of photosynthetic activity, in the landscape surrounding the breeding bird survey (BBS) routes in the conterminous USA are correlated and explain changes in local bird diversity. We found that landscape-wide abrupt shifts in magnitude were correlated with a decrease of the geometric mean of relative abundances (GM, Buckland et al. 2011) and an increase in the progressive bird assemblage composition (pBC, Rittenhouse et al. 2010). A great proportion of browning land was correlated with a decrease in relative bird abundance and changes in assemblage composition, while more greening land affected assemblage composition only (Fig. 2). Confirming previous studies, some environmental factors (e.g. mean elevation and photosynthetic activity) influenced average bird diversity change. Changes in both relative abundance and assemblage composition were not only influenced by concomitant abrupt shifts in magnitude, but also by individual and cumulative effects of preceding changes in LULC (Fig. 3). On average, landscape-wide changes in LULC had some explanatory power (R2 > 0.1) for only a few BBS routes without any clear pattern in space (Fig. 4), across trait groups (Fig. 5a) or among ecoregions (Fig. 5b). We discuss how these results link to previous studies of local biodiversity change and landscape ecology in general.

Landscape-wide changes in land use and land cover as drivers of biodiversity change

Changes in land use and land cover (LULC) have previously been linked to local biodiversity change (Brooks et al. 2002; Ewers et al. 2013; Cousins et al. 2015; Jung et al. 2019b). In line with these studies, we found local biodiversity measures to be more affected by larger abrupt shifts in magnitude at the landscape scale (Fig. 2). A great proportion of abrupt shifts in magnitude and trend (for ‘browning’) in the wider landscape were associated with a significant decline in relative bird abundance (Fig. 2a), potentially indicating local bird population decline as fewer individuals were observed (Loh et al. 2005; Buckland et al. 2011). Meanwhile more abrupt shifts in magnitude and trend in the wider landscape increased the similarity in assemblage composition relative to y0 (Fig. 2b). Because we found the relative mean abundance of bird species (GM) to decline with a greater proportion of landscape-wide changes in LULC, it is likely that the changes in assemblage composition (pBC) are caused by an increase in species richness, a pattern shown before for the BBS data (Schipper et al. 2016). The greening of a landscape may be related to agricultural intensification (Zhu et al. 2016), which is known to cause declines in bird abundance (Stanton et al. 2018). Previous studies found compositional changes in bird assemblages to be particularly associated with changes in the occurrence of rare and specialist species, leading to a “homogenization” of assemblages (McKinney and Lockwood 1999; Olden 2006; Newbold et al. 2018). It could be that landscape-wide changes in LULC increase the heterogeneity of resources and bird habitats available, thus allowing a greater number of bird species, but fewer individuals overall, to thrive (Holt 2009; Stein et al. 2014), for instance through increased competition (Randall Hughes et al. 2007).

Changes in relative bird abundance and assemblage composition differed with environmental gradients (Fig. S11). Consistent with previous studies (Lomolino 2001; Jarzyna and Jetz 2017), relative bird abundance increased at BBS routes of high elevation (Fig. S11a), indicating that bird species increasingly utilize high elevation regions. These species appear to be different from the species previously inhabiting BBS routes at high elevations, given the strong influence of elevation on assemblage composition (Fig. S11b). Furthermore, we found changes in bird abundance to decrease and “flatten” in BBS routes with high average photosynthetic activity (EVI > 0.4, Fig. S11a), in contrast to changes in bird assemblage composition, where changes occurred in BBS routes of high photosynthetic activity (Fig. S11b). This is consistent with previous studies that demonstrated that high photosynthetic activity is negatively correlated with bird abundance (Barnagaud et al. 2017) but positively with richness (Rowhani et al. 2008; Goetz et al. 2014), which could drive changes in assemblage composition. Similar to previous studies (Barnagaud et al. 2017), we found no strong effect of precipitation anomalies prior to a BBS on bird abundance or assemblage composition (Fig. S11).

Lag effects of preceding changes in land use and land cover

Changes in land use and land cover can have immediate and delayed impacts on local biodiversity (Kuussaari et al. 2009; Hylander and Ehrlén 2013; Jung et al. 2019b). Theory suggests that—single and cumulative—preceding changes in LULC are correlated with larger changes in local biodiversity (Scheffer et al. 2001; Andersen et al. 2009; Dornelas 2010; Watson et al. 2014; Ratajczak et al. 2018). We demonstrated that considering preceding landscape-wide changes in LULC, assessed as changes in trend in or magnitude of photosynthetic activity, helps explain local changes in bird abundance and assemblage composition (Fig. 3). Increasing explanatory power of individual preceding years could be linked to an average “ecological memory” effect for birds, particularly for the 4th and 5th year prior and for changes in assemblage composition (Fig. 3), similar to what has been shown for plant species (Ogle et al. 2015). The impacts of cumulative preceding changes in LULC depends on their duration (Essl et al. 2015b) and frequency (Watson et al. 2014; Ratajczak et al. 2018) and a recent study found that considering cumulative preceding—relative to concomitant—changes explained more differences in local (bird) biodiversity (Jung et al. 2019a). Similarly, we found that a consideration of cumulative periods of preceding changes in LULC explained some differences in local biodiversity (Fig. 3). Preceding changes may have affected the resources available to birds thus directly influencing their fitness and persistence in subsequent years (Holt 2009; Harrison et al. 2011; Ogle et al. 2015).

Our understanding of “lagged” effects of changes in LULC on biodiversity change are still in their infancy. The majority of previous studies investigated climatic influences on richness and abundance change (Albright et al. 2011; Lindström et al. 2013; Valtonen et al. 2013; Martay et al. 2017), but little is known about the influence of past changes in LULC as assessed from remote sensing. Rittenhouse et al. (2012) investigated differences in the proportion of landscape-wide land cover on bird diversity, but only used bi-annual, thematically non-consistent estimates of land cover. Other studies investigated the link between preceding change in LULC and local biodiversity (Jung et al. 2019b), but only for spatial differences in local biodiversity rather than biodiversity change per se, which might mask lasting impacts (França et al. 2016; De Palma et al. 2018). It could also be that concomitant impacts of landscape-wide changes in LULC on bird diversity are not apparent yet, with some previous studies having found lags of up to several decades (Findlay and Bourdages 2000; Watts et al. 2020). Future studies could benefit from analysing impacts of preceding changes in LULC on both mean and variance of biodiversity change (Leung et al. 2017; Christensen et al. 2018), considering longer past time frames as well as considering varying trajectories of remotely-sensed change in LULC (Watson et al. 2014) and identifying the drivers of LULC change (e.g. deforestation, agricultural expansion).

Variability in explanatory power in space and functional traits

Quantifying local biodiversity change and identifying drivers of these changes is not trivial (Dornelas et al. 2012; Cardinale et al. 2018). Drivers of local biodiversity change are often unknown or cannot be reliably quantified (Hallmann et al. 2017). In an attempt to forecast local bird richness change, Harris et al. (2018) parametrized models with and without (‘naïve’) remotely-sensed photosynthetic activity and climatic data. Surprisingly, they found naïve models to predict bird richness change better than those including vegetation and climate variables, which they attributed to a lack of abrupt biodiversity changes. Compared to bird richness, which has been found to increase in the BBS data (Schipper et al. 2016), we found abundance and assemblage composition to decline (Fig. S5), but it remains unclear what has driven these changes.

Bird diversity can be constrained by “thresholds” of land-surface conditions—such as vegetation availability—in the wider landscape (Andersen et al. 2009; Gutzwiller et al. 2015). A global review of threshold responses to landscape-wide changes in LULC suggests that bird diversity is most affected if more than 27.9% of the landscape has changed (Melo et al. 2018). With exception of a few BBS routes (Figs. 4 and S6), the average proportion of change in LULC, assessed as changes in trend in or magnitude of photosynthetic activity, within landscapes was only 6% (Figs. S6, S7 and S10), which could explain why landscape-wide changes in LULC in our models explained on average little variation in bird diversity change and were important for a few BBS routes only (Fig. 4). However, it could also be that impacts of landscape-wide changes in trend in or magnitude of photosynthetic activity on bird diversity are poorly generalizable and depend on local context and functional traits of bird species. In this study we quantified changes in LULC as changes in trend in or magnitude of photosynthetic activity without identifying the underlying drivers of these changes, such as for example wild fires, vegetation defoliation by insects, or specific land use change. Future efforts should investigate whether distinguishing by such drivers could explain a greater proportion of variance of bird diversity change.

Changes in local bird diversity differ by functional trait groups (Fig. S12, Jarzyna and Jetz 2017; Barnagaud et al. 2017). Yet, the explanatory power of landscape-wide changes in LULC on bird diversity change did not vary by functional groups (Fig. 5a). Many American birds are migratory and as such are affected by human persecution and climatic anomalies on their migration paths (Sanderson et al. 2006; Tottrup et al. 2012). Although we did not find any difference in explanatory power between migratory and non-migratory species (Fig. 5a), our analysis only considered changes in LULC in bird breeding grounds as the location of wintering grounds are unknown. In contrast to previous studies (Schipper et al. 2016), a breakdown into habitat guilds also did not assist in identifying differences in explanatory power (Fig. 5a), which is surprising given the difference in trend between for instance woodland and grassland birds (Fig. S12). Potentially changes in LULC specific to certain bird habitats, e.g. changes in vegetation height (Goetz et al. 2014), are better predictors of bird diversity change in such cases.

Conclusion

In this study we investigated the influence of landscape-wide changes in LULC, quantified as changes in trend in or magnitude of photosynthetic activity, on biodiversity trends. Landscapes surrounding the BBS routes are constantly changing (Figs. S6 and S7) and such changes are expected to influence local biodiversity (Manning et al. 2009; Turner and Gardner 2015; Seppelt et al. 2016). However, the processes influencing local biodiversity at the landscape scale are difficult to quantify (Chase 2003), dependent on spatial scale (Miguet et al. 2016) and other environmental predictors (elevation, terrain, climatology). Overall our results indicate that landscape-wide changes in LULC are correlated with (Fig. 2) but on average did not explain bird diversity change across spatial scales (Fig. 4), functional groups (Fig. 5a) or ecoregions (Fig. 5b). Preceding changes in LULC assisted in explaining changes in bird diversity (Fig. 3), highlighting the importance of biotic lag effects. Overall, for most BBS routes, the drivers explaining local bird diversity change remain unknown (Figs. 4 and 5) and we suggest future studies to consider alternative attributes of remotely-sensed changes in LULC at the landscape-scale (Watson et al. 2014) or other spatio-temporal variables not quantifiable from optical remote sensing (e.g. pesticide use, human persecution). We furthermore suggest that more research is needed on scale-dependent effects (local vs landscape changes) of biodiversity change.