DOI: 10.5061/dryad.tdz08kq05
title: Mismatches between birds' spatial and temporal dynamics reflect their delayed response to global changes

description: 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.
description: For details on the methodological workflow used to built the
 dataset, see Materials and Methods section of the paper. 
 <strong>Raw bird data are not open but can be acessed via
 request to the French National Museum of Natural History</strong>
 (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 <strong>Moving window</strong>. 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. <strong>Calculating deviation between
 spatial and temporal dynamics. </strong>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
 (&lt;www.r-project.org&gt;, 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.
 <strong>Human influence index. </strong>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]. <strong>Colonization–extension
 ratio. </strong>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
 (&lt;www.r-project.org&gt;, 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 &gt; 1 indicates
 that the site in the focal window experience more colonization than
 extinction, and conversely for values &lt; 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
 &gt; 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]. <strong>All the aforementioned processed
 data were compiled into the ssdat_28062020.Rdata dataframe provided in the
 repository.</strong> <strong>Statistical
 analyses. See Res_and_figures.R. </strong>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. 
description: 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.
 <strong>ssdat_28062020.Rdata </strong>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.  <strong>Res_and_figures.R</strong>
 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.