Global changes alter the dynamics of biodiversity, and are forecasted to continue or worsen in the decades to come. Modelling approaches used to anticipate these impacts are mainly based on the equivalence between spatial and temporal response to environmental forcings, generally called space-for-time substitution. However, several processes are known to generate deviations between spatial and temporal responses, potentially undermining the prediction based on space-for-time substitution.

We here used high-resolution data from the french breeding bird survey to quantify and map the deviation between spatial and temporal patterns of bird abundances resulting from the dynamics of 124 species monitored in 2133 sites between 2001 and 2012. Using independent empirical data, we then tested specific predictions linked to the determinants (anthropogenic activities) and processes (lagged responses to environmental changes) potentially generating these deviations. We found that deviations between spatial and temporal patterns of abundances were particularly structured in space for bird communities. Following our predictions, these space–time deviations were positively correlated with the human influence on ecosystems, and linked with colonization–extinction ratios and community completeness, two markers of ongoing delayed responses to environmental changes. Our results suggest that the deviations between spatial and temporal patterns are related to recent anthropogenic environmental changes and disequilibrium responses to these changes. Investigating deviations between spatial and temporal patterns of biodiversity might open promising perspectives for a formal quantification of disequilibrium state of biodiversity at large spatial scale.

For details on the methodological workflow used to built the dataset, see Materials and Methods section of the paper. 

Raw bird data are not open but can be acessed via request to the French National Museum of Natural History (contact frederic.jiguet@mnhn.fr ). Data were taken from the French Breeding Bird Survey (FBBS) over the 2001–2012 period (Fig. 2a). In this survey, skilled volunteer ornithologists count birds at a given site, fol-lowing a standardized protocol, at the same site, year after year (Jiguet et al. 2012). Species abundances have been recorded for 2133 sites, each covering 4 km2. Each volunteer provides his/her home locality to the national coordinator, and a 2 ×2 km plot to be prospected is randomly selected from within a 10 km radius (i.e. among 80 possible plots) by the national coordinator. Each spring, volunteers carry out 10 point counts separated by at least 300 m within the selected 2 × 2 km site, for a fixed period of five minutes. Two sampling sessions are realized from 1 April to 8 May, then from 9 May to end of June (in order to detect both early and late breeders), with 4–6 weeks between both sessions. Counts must be repeated on approximately the same date each year (± 7 days), and at dawn (1–4 h after sunrise) by the same observer, in the same order. The maximum per-point count from these two sessions is retained as the measure of point-level species abundance. The French breeding bird survey have been designed to moni-tor population trends of the common bird species in France. Therefore, uncommon species are inadequately monitored by this scheme, and might inflate inter-annual variability at site scale and more generally lead to spurious estimations of population trends. We therefore excluded very rare birds and incorrectly monitored species (e.g. specific waterbirds in rare habitats) from our dataset by limiting the final data sample to the 116 bird species representing 99% of total abundance.Moving windowWhile the FBBS provides a large number of observations and a high density of sites across France, all sites do not have equal temporal coverage (average temporal coverage of FBBS sites: 4.73 years ± 3.20 [mean ± SD]). In order to provide reliable estimates of local population sizes and trends of birds in a

Moving window. While the FBBS provides a large number of observations and a high density of sites across France, all sites do not have equal temporal coverage (average temporal coverage of FBBS sites: 4.73 years ± 3.20 [mean ± SD]). In order to provide reliable estimates of local population sizes and trends of birds in a patial continuum covering the whole study area and period (among other reasons stated in the introduction), we used a spatial moving window approach (Fig. 1a). The principle is to calculate relevant metrics within a moving window delin-eated by a circle around a given focal site. This is a straight-forward way to summarize both local spatial and temporal trends emerging from regional dynamics (Gaucherel 2007, Gaucherel et al. 2008) and has been already valuably used with this dataset (Devictor et al. 2010, Godet et al. 2015, Gaüzère et al. 2017). Each window is a 80 km-radius circle centered on a focal site, and encompassing at least 20 sites. This resulted in 2094 windows of similar size. Note that the radius of 80 km was a compromise between spatial resolution, minimum number of sites within the window, and the result-ing number of windows to cover the study area (for more details, Gaüzère et al. 2015). Because the same site can be integrated in more than one window, the estimates produced by this approach are correlated in space. We therefore explicitly integrated this spatial structure in subsequent analysis.

Calculating deviation between spatial and temporal dynamics. For each window, space–time deviation was calculated as the difference between the relative population trend and the rela-tive population size of birds during the period 2001–2012. The relative population trend was estimated as thscaled change in bird abundance (all species together) per unit of time (year). First, we used linear mixed models performed with the lmer{lme4} function in R (<www.r-project.org>, Douglas Bates et al. 2015) to estimate temporal population trend. In these models (one for each window), the response variable was the log-transformed total abundance of birds (calculated for a given site at a given year), regressed over a fixed effect of year (2001–2012). To take into account the uncontrolled variability between sites (different observers, habitat) and species (different life-history traits such as mass) we included the species nested within the site as nested ran-dom effects on the intercept in our model. The fixed-effect year coefficient (i.e. the slope) indicated the temporal popula-tion trend over the 12 years of the survey. Temporal trends values were then scaled to unit standard deviation in order to compute relative population trends. We do not centered pop-ulation trend values to zero mean in order to keep the natural sign of population trend in subsequent analyses.The relative population size was calculated for each window as the scaled-centered average bird abundance over the study period and across all sites within the window. We first calculated the average population size of each site over the period of survey. Average population sizes were then centered to zero mean and scaled to unit standard deviation in order to compute the relative population size. Negative values of relative population size indicates a bird population is smaller than the average breeding bird population size in France, and conversely, positive values indicate populations larger than average. Note that we here investigated spatial and temporal patterns of abundance, but this method can be applied to other variables of interest. For example, we calculated space–time deviation based on community species richness, total bird biomass and species-by-species populations.

Human influence index. We used the global human influence index (HII) (ver. 2, Wildlife Conservation Society and Center for International Earth Science Information Network, Columbia University, 2005) from the NASA Socioeconomic Data and Applications Center as a proxy of a multi-dimensional anthropogenic pres-sure on ecosystems. This index estimates the direct human footprint on ecosystems (Sanderson et al. 2002) as a proxy of recent human pressure on biodiversity. It incorporates nine global data layers corresponding to: human population pres-sure (population density); human land use and infrastructure (built-up areas, nighttime lights, land use and land cover); and human access (coastlines, roads, railroads, navigable riv-ers). HII values at 1-km resolution were extracted for each of the monitored sites, and averaged for each window. The average HII value was 28.7 ± 5.23 [mean ± SD].Natural area extentWe used the percent of natural area as a proxy measuring the extent of natural surface in the window. The percent of natu-ral area in the window was based on landscape composition data from CORINE Land Cover 2006 raster data (European Environment Agency 2010) for each 80 km-radius window. In order to be able to consider large-scale biogeographical gradi-ents, CORINE Land Cover classes were aggregated as a func-tion of their top-level nomenclature, namely: artificial surfaces, agricultural areas, forest and semi-natural areas, wetlands and water bodies (European Environment Agency 2007). We then computed the percentage land cover (calculated as the habitat class area, in square meters, over total land area of the window) for the first three of these classes, in order to define landscape composition variables: % artificial area, % agricultural area and % natural area. In the CORINE nomenclature, natural areas encompasses broad-leaved forest, coniferous forest, mixed forest, natural grasslands, moors and heathland, sclerophyl-lous vegetation, transitional woodland-shrub, beaches, dunes, sands, bare rocks, sparsely vegetated areas, burnt areas, glaciers and perpetual snow. The average natural area extent was 31 ±18 [mean ± SD, in % of the total window area].

Colonization–extension ratio. We calculated colonization–extinction ratio by first selecting all sites monitored for at least three consecutive years (repre-senting 58% of the total dataset). For each site, we used the betadiver function in the R vegan package (<www.r-project.org>, Oksanen et al. 2015) to compute, for each pair of con-secutive years: the number of species present in both years (a); the number of colonizing species, i.e. species recorded in year t + 1 that were not present in year t (b); the number of extinct species (c); and the total number of species S (a + b + c). For each site, the colonization–extinction ratio was calculated as the average of b/c for pairs of consecutive years. These site-scale colonization–extinction ratio values were then averaged for each window. A value > 1 indicates that the site in the focal window experience more colonization than extinction, and conversely for values < 1. The average colonization–extinction ratio was 0.98 ± 0.030 [mean ± SD].Completeness of communityThe completeness of each community index was calculated as the log-ratio of dark and observed diversity (Pärtel et al. 2013): ln(observed richness/dark diversity). Dark diversity represents the number of species that would be environ-mentally and ecologically suited to conditions in the target local species pool, but which are absent from the observed community. It was estimated at point-count scale over all the study period using the co-occurrence method presented in (Lewis et al. 2016). An absent species is considered to be part of the local species pool if it typically co-occurs with the species that are observed in the focal c–ommunity. We selected this method because it is more accurate than other methods based on similarity in species’ ecological require-ments (Lewis et al. 2016), and because we lack information on point-scale environmental conditions. The probability of species co-occurrence was estimated by Beals index (Beals 1984, De Cáceres and Legendre 2008) calculated from the point-scale community matrices (2001–2012) using the beals{vegan} R function. Each species absent from the observed species pool was considered to be in the local pool (i.e. dark diversity) if the probability of co-occurrence was > 0.95 (5% threshold). This threshold is arbitrary, and the presence of outliers, which are common in co-occurrence probability dis-tributions can affect the outcome (Botta-Dukát 2012). We therefore modeled the first, fifth and tenth percentile thresh-olds in order to ensure the robustness of our estimations (Supporting information). Site-scale completeness values were then averaged for each window. The average coloniza-tion–extinction ratio was 2.22 ± 0.236 [mean ± SD].

All the aforementioned processed data were compiled into the ssdat_28062020.Rdata dataframe provided in the repository.

Statistical analyses. See Res_and_figures.R. We tested the non-random spatial distribution of the space–time deviations (prediction i) by modelling the spatial varia-tion of space–time deviation. To do so, we implemented a generalized additive mixed model (GAMM) where space–time deviation was regressed over a two-dimensional spline smoothing based on geographic coordinates (with a free degree of freedom). We included the factorial variable bio-geographic domain (Mediterranean, Continental, Alpine Atlantic) as a random effect to account for regional variability in climate and species pools. In short, this model tests for the presence of a significant effect of geographic coordinates on the space–time deviation. We also tested linear variation of space–time deviation along latitude and longitude by imple-menting two linear mixed models (LMM) where space–time deviation was regressed over latitude or longitude, with ran-dom effect of biogeographic domain. We tested predictions (2–5) linked the relationships between space–time deviation and (2) impact of anthropo-genic activities on ecosystems, (3) extent of natural areas, (4) ratio between local colonization and extinction and (5) completeness of community by implementing generalized additive mixed models (GAMM) testing the statistical rela-tionships between space–time deviation and relevant pre-dictors. We ran four models (one for each prediction) in which space–time deviation (calculated for each window, n = 2094) was the response variable regressed over each pre-dictor (2: human influence index, 3: natural area extent, 4: colonization–extension ratio, 5: completeness of com-munity). We scaled each predictor to zero mean and unit standard deviation prior to modeling in order to facilitate direct comparisons among them. We considered each model as an independent test as the predictors were uncorrelated (Supporting information). In order to explicitly account for the non-random spatial structure in the data and prevent issue linked to spatial autocorrelation in model residuals, we added a the geographical coordinates as co-variable with a two-dimensional spline-based smoothing (with a free degree of freedom) based on latitude and longitude, following the method of (Wood 2006). The biogeographic domain (facto-rial variable: Mediterranean, Continental, Alpine, Atlantic) was introduced as a random effect to account for regional variability in climate and species pools. For the model testing the relationship between space–time deviation and commu-nity completeness (prediction 4) we added the average site species richness as a fixed-effect co-variable in the model to account for its potential confounding effect (local richness is correlated to community completeness). Models residuals were formally checked for normality, homoscedasticity and absence of correlation with predicted values and spatial auto-correlation, and passed in all cases.

The code is also relevant to understand the different steps of the analyses. 

Please feel free to contact me if you need any complementary information. 

Data and code to replicate the analyses and results from the paper "Mismatches between birds' spatial and temporal dynamics reflect their delayed response to global changes" https://onlinelibrary.wiley.com/doi/10.1111/oik.08289

The data and code are distributed under the Licence creative Commons CC0 1.0 Universal (CC0 1.0) Public Domain Dedication. Which means you can reuse the code and data as your convenience without asking permission. However, we appreciate if you can cite the paper when using the data or the code.

ssdat_28062020.Rdata is a R data.frame containing all informations (in columns) for all window (in rows) allowing to recreate the figures and run the models presented in the paper. 

Res_and_figures.R is R code that allows to compute deviation between spatial and temporal dynamics, do the paper's  figures, and run models using the ssdat_28062020.Rdata dataset.