Mismatches between birds' spatial and temporal dynamics reflect their delayed response to global changes - Gaüzère, Pierre; Devictor, Vincent (Dryad, 2021)
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.
Keywords: FOS: Natural sciences