Biotic homogenisation in bird communities leads to large‐scale changes in species associations

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Biotic homogenisation in bird communities leads to large‐scale changes in species associations © 2021 the Authors Published version Rigal, Stanislas; Devictor, Vincent; Gaüzère, Pierre; Kéfi, Sonia; Forsman, Jukka T.; Kajanus, Mira H.; Mönkkönen, Mikko; Dakos, Vasilis Rigal, S., Devictor, V., Gaüzère, P., Kéfi, S., Forsman, J. T., Kajanus, M. H., Mönkkönen, M., & Dakos, V. (2022). Biotic homogenisation in bird communities leads to large‐scale changes in species associations. Oikos, 2022(3), Article e08756. https://doi.org/10.1111/oik.08756 2022 www.oikosjournal.org OIKOS Oikos 1 © 2021 The Authors. Oikos published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Jiri Reif Editor-in-Chief: Dries Bonte Accepted 2 November 2021 00: 1–14, 2021 doi: 10.1111/oik.08756 00 1–14 The impact of global change on biodiversity is commonly assessed in terms of changes in species distributions, community richness and community composition. Whether and how much associations between species are also changing is much less documented. In this study, we quantify changes in large-scale patterns of species associations in bird communities in relation to changes in species composition. We use network approaches to build three community-aggregated indices reflecting complementary aspects of species association networks. We characterise the spatio–temporal dynamics of these indices using a large-scale and high-resolution dataset of bird co-abundances of 109 species monitored for 17 years (2001–2017) from 1969 sites across France. We finally test whether spatial and temporal changes in species association networks are related to species homogenisation estimated as the spatio–temporal dynamics of species turnover (β-diversity) and community generalism (community generalisation index). The consistency of these relationships is tested across three main habitats, namely woodland, grassland and human settlements. We document a directional change in association-based indices in response to modifications in species turnover and community generalism in space and time. Weaker associations and sparser networks were related to lower spatial species turnover and higher community generalism, suggesting an overlooked aspect of biotic homogenisation affecting species associations and may also have an impact on species interactions. We report that this overall pattern is not constant across habitats, with opposite relationships between biotic homogenisation and change in species association networks in urban versus forest communities suggesting distinct homogenisation processes. Although species associations contain only partial signatures of species interactions, our study highlights that biotic homogenisation translates to finer changes in community structure by affecting the number, strength and type of species associations. Keywords: Avifauna, community, homogenisation, interaction network, species association, β-diversity Biotic homogenisation in bird communities leads to large-scale changes in species associations Stanislas Rigal, Vincent Devictor, Pierre Gaüzère, Sonia Kéfi, Jukka T. Forsman, Mira H. Kajanus, Mikko Mönkkönen and Vasilis Dakos S. Rigal (https://orcid.org/0000-0002-3180-4592) ✉ ([email protected]), V. Devictor, S. Kéfi and V. Dakos, ISEM, Univ. de Montpellier, CNRS, IRD, EPHE, Montpellier, France. – P. Gaüzère (https://orcid.org/0000-0003-1259-6131), Univ. Grenoble Alpes, CNRS, Univ. of Savoie Mont Blanc, LECA, Laboratoire d’Écologie Alpine, Grenoble, France. SK also at: Santa Fe Inst., Santa Fe, NM, USA. – J. T. Forsman and M. H. Kajanus, Dept of Ecology and Genetics, Univ. of Oulu, Oulu, Finland. JTF also at: Natural Resources Inst. Finland, Oulu, Finland. – M. Mönkkönen (https://orcid. org/0000-0001-8897-3314), Dept of Biological and Environmental Science, Univ. of Jyväskylä, Jyväskylä, Finland. Research 2 Introduction Among the major effects of global change on biological diversity, the modification or even the extinction of species interactions has early on been identified as being pervasive, but is still poorly understood (Janzen 1974, Diamond 1989). A perturbation in species interactions may be decoupled from changes in species richness or community composition because there are many more interactions than species (Poisot et al. 2015, Gravel et al. 2019). In particular, modifications affecting species interactions can be stronger (Valiente-Banuet et al. 2015) or weaker (Li et al. 2018) than those affecting species richness. The structure and dynamics of species interactions are among the main drivers of community dynamics (Davis et al. 1998, Barabás et al. 2016), and therefore represent a critical subject of study for ecology and biodiversity conservation (Kissling and Schleuning 2015, García-Girón et al. 2020). Despite the importance of integrating species interactions into conservation biology, we still have a limited understanding of the drivers and consequences of changes in the strength and the structure of species interactions. In the last decades, there has been an increasing use of network approaches to study species interactions in empirical and theoretical communities (Bascompte et al. 2003, Ings et al. 2009, Kéfi et al. 2015, Trøjelsgaard and Olesen 2016). Ecological communities can thus be depicted as interaction networks by defining nodes as individuals or species, and links between the nodes as species interactions (Newman et al. 2006). The estimation of species interactions is, however, subject to a conceptual question (how are the strength and type of interactions defined?) and a technical challenge (how to estimate an interaction?). In some cases (e.g. in simple trophic networks with few taxa), observations or experiments can address both issues as the existence and type of species interactions are clearly identified. However, these cases provide inference of species interactions in local and specific systems, which makes it difficult to derive general rules for interactions in larger communities (Whittaker et al. 2005, Denny and Benedetti-Cecchi 2012). The empirical identification and measure of interactions in species-rich communities in particular, is challenging because of the high number of potential interactions to be estimated (proportional to the square of species number) (Barner et al. 2018). An alternative approach is to assume that species associations (inferred from their spatial aggregation) are shaped, at least to some extent, by the combination of true interactions (i.e. clear ecological relationships such as competition or predation). In this case, studying communities with a large number of species and broad spatial coverage should be a good framework for estimating species associations, although the ability of spatial co-occurrence patterns to infer pairwise species interactions is still controversial (Blanchet et al. 2020). Nonetheless, species co-occurrence might be an information-rich proxy of the outcome of direct and indirect biotic interactions in communities (Delalandre and MontesinosNavarro 2018, Freilich et al. 2018). Indeed, the composition of a local community results from interspecific interactions as well as multiple intertwined processes generating patterns of spatial aggregation between species (Fig. 1) (Wisz et al. 2013). These factors include neutral processes (regional dispersion and local stochasticity (Hubbell 2001)), historical processes (phylogeography (Kraft et al. 2007)) and niche processes (HilleRisLambers et al. 2012, Letten et al. 2017). Niche processes combine what are sometimes referred to as Grinnellian and Eltonian processes (Chase and Leibold 2003, Devictor et al. 2010a). Grinnellian processes (Grinnell (1917); later extended by Hutchinson (1957)) consider the niche as the species’ response to environmental conditions acting as an environmental filter for the community. Eltonian processes (Elton 1927) consider the niche as the species’ impact on its environment and refers to the mutual dependency of species with each other, including the limiting similarity hypothesis (i.e. the niche overlap between two species that limits their coexistence) (MacArthur and Levins 1967, Abrams 1975, Martin and Bonier 2018) and facilitation between species (e.g. cooperation, exchange of social information (Seppänen et al. 2007, Gil et al. 2019, Tu et al. 2019)). Refining the Eltonian component (i.e. the part of the cooccurrences due to the biotic filter) can be done through multiple approaches (Kissling et al. 2012). Based on null models, one can control for the species associations that are simply expected by chance rather than grounded in ecological processes by testing whether species are found together more or less frequently than expected by chance (Gotelli 2000, Ulrich and Gotelli 2010, Kohli et al. 2018). Indirect effects between species (i.e. the effect of a third species on the association between two other species) can be evaluated using partial correlations (Faust and Raes 2012, Harris 2016). Recent progress with joint species distribution models have also provided ecologists with new tools for estimating species associations by studying residual co‐occurrence patterns after accounting for environmental niches from large datasets (Tikhonov et al. 2017, Zurell et al. 2018). Overall, recent methods removing non-Eltonian components from co-occurrences (Azaele et al. 2010, Faisal et al. 2010, Ovaskainen et al. 2010, Lindenmayer et al. 2015) are promising for uncovering species association networks (Araújo et al. 2011, Morueta-Holme et al. 2016). Co-occurrences are thus information-rich proxies of direct and indirect biotic interactions in communities (Delalandre and Montesinos-Navarro 2018, Freilich et al. 2018). Association networks based on species co-occurrences are useful for capturing community organisation through aggregated community indices, i.e. statistics summarising an aspect of the network at the community level (Barner et al. 2018) even if interaction networks may remain out of reach (Sander et al. 2017, Thurman et al. 2019). Tracking large-scale changes in species associations might represent a significant advance for macro-ecology and conservation biogeography. Indeed, ecological processes are ultimately influenced by which and how species interact (Cardinale et al. 2002, Goudard and Loreau 2008). Moreover, the responses of species associations to environmental changes 3 are not necessarily proportional to the responses of individual species (Valiente-Banuet et al. 2015). Therefore, measuring community changes through the change in species diversity within local communities or, at a larger scale, between communities (for instance using β-diversity) may mask important modifications of the structure and properties of those communities (Poisot et al. 2017). For instance, the replacement of a set of diverse and mainly specialist species by a few generalists (McKinney and Lockwood 1999, Olden et al. 2004) is a well-documented form of biotic homogenisation (Clavel et al. 2011). In those communities, the composition tends to be closer to random expectations (Barnagaud et al. 2017), i.e. with less and less visible niche processes. Yet, anthropogenic perturbations can also act as a strong filter selecting for more specialist species (Gaüzère et al. 2020), in which case the Grinellian filter may outweigh the Eltonian one. While many studies have evidenced the impact of global change on biotic homogenisation (Newbold et al. 2018), whether homogenisation in species composition (as a decrease in species turnover and an increase in community generalism) is related to a directional change in species associations remains to be explored (Li et al. 2018). In this study, we conducted a large-scale spatio–temporal analysis of bird species association networks. Birds form a relevant group to estimate ecologically meaningful species associations (Sridhar et al. 2012). Indeed, the Eltonian filter is likely to impact bird co-occurrences since competitive exclusion as well as social information exchanges have been frequently shown to occur between birds among other interacting processes (Thomson et al. 2003, Forsman and Thomson 2008, Magrath et al. 2015). Furthermore, the fact that bird species have been widely monitored for decades allows to track changes in species association patterns in space and time; an opportunity that is not possible for most of other groups of organisms. Using data from the French Breeding Bird Survey, we addressed the following two objectives: 1) reconstruct species association networks in communities from co-abundance data, 2) test whether biotic homogenisation was linked to directional changes in association networks. Material and methods Overview We first inferred species associations from co-abundance (co-occurrence with abundance) data (Fig. 2a) corrected for non-Eltonian co-occurrence processes. We then quantified different structural properties of the species association networks using three complementary network indices: intensity, attractiveness and clique structure of the network. Intensity corresponds to the mean association strength, attractiveness is the ratio of positive to negative associations, and clique Figure 1. Community assembly processes and species co-abundance. Species interactions that influence species spatial aggregation (or segregation) and temporal change in abundance are referred to as the Eltonian component of species co-abundance. In addition to the Eltonian component, co-abundances are also the result of habitat filtering (Grinellian), random processes due to neutral dispersal, as well as historical processes related to the species’ phylogeography. The result of all these processes leads to the observed species co-abundances. Each letter stands for a different species. Species U and Z share a common biogeographic region and random processes have not prevented them from co-occurring. As they live in a similar habitat and interact in a way that enables their coexistence, they can be observed together in the same location at the same time. 4 structure describes the structural complexity of the association network (Fig. 2b). We then analysed the relationships between the spatio– temporal dynamics in bird network indices and the spatio–temporal dynamics in biotic homogenisation (through species turnover and community generalism) at a large scale and separately in three main types of habitats: woodland, grassland and human settlements. All the analyses were done with the R software (ver. 3.4.4, <www.r-project.org>) and the R scripts and data are available on Dryad (Rigal et al. 2021). Bird data Bird data were extracted from the French Breeding Bird Survey (FBBS) (Jiguet et al. 2012). In this scheme, volunteer ornithologists monitored common bird species on 2514 sites (Fig. 2a) from 2001 to 2017, following a standardised protocol. Sites are 2 × 2 km squares in which abundances of breeding bird species were monitored on 10 homogeneously distributed sampling points across habitats in the landscape. Each sampling point was monitored for 5 min, 1–4 h after sunrise, twice a year, the first between 1 April and 8 May and the second after an interval of 4–6 weeks, between 9 May and 15 June, to account for early and late breeding birds. For each sampling point and each year, in addition to bird abundance and identity, the geographical coordinates, weather conditions, altitude, distance of the contacts from the observer and main habitat (44 types) were recorded. In order to avoid habitat classes with too few observations, we grouped the 44 types of habitat described in the field into 19 classes for the estimation of species associations. We also distinguished the three main types of habitats of our dataset (woodland, grassland and human settlements) for network analyses (for habitat details see the Supporting information). Among the 242 species recorded in the dataset, we selected the 109 most abundant species (representing 99% of the total abundance) to avoid any over-representation of rare species (that are therefore more difficult to monitor). Moreover, we used the abundances corrected for the detectability of the species (Supporting information) and kept only the maximum abundance of the two passages for each species. After removing rare species and the sites only monitored once, our dataset comprised 19 580 sampling points in 1969 sites and 109 species (species listed in the Supporting information). Hereafter, a community corresponds to birds monitored in a sampling point and for which a species association network can be calculated. Figure 2. From bird monitoring to association network indices. (a) Spatial distribution of the 1969 selected (out of 2514) sites from the 2001 to 2017 FBBS (STOC-EPS). On each site (2 × 2 km square), bird observations were recorded on 10 sampling points. Mainland France was split into 4 biogeographic regions (Alpine, Atlantic, Continental and Mediterranean). (b) Association network indices: intensity (i.e. mean strength of species associations), attractiveness (predominance of positive (light orange) or negative (dark purple) species associations) and clique structure (level of structuration of the species association network). (c) Geographically weighted regressions using data from moving windows were used to assess β-diversity and species association metric values. For each index, examples corresponding to low and high values are displayed from left to right, respectively and lines represent species associations. For β-diversity, site colours represent site composition, the more diverse the colours in the window, the higher β-diversity. For intensity, the thicker the line, the higher the absolute value of the association. For attractiveness, the high value example is 0.5 (three positive associations and one negative association out of 4 existing associations) and the low value is −0.5. For clique structure, the high value is 0.125 (two realised cliques out of 16 possible cliques) and the low value is 0. 5 Species associations We estimated associations between pairs of species from bird co-abundance data (Morueta-Holme et al. 2016) for each year (2001–2017), for each of the four biogeographic regions and for each of the 19 habitats using the five following steps (Fig. 3). Step 1. In order to limit the influence of phylogeography and habitat features on species associations, we first grouped the data by biogeographic region (Continental, Atlantic, Mediterranean, Alpine), by habitat (each sampling point comes with information about its habitat (among the 44 types) from the field observer and this habitat type was actualised according to the grouping made in this study (19 habitat classes, Supporting information)) and by year to estimate an association for each pair of bird species, for each year, for each of the four biogeographic regions (EEA 2016) and for each of the 19 habitat classes. Bird assemblages are indeed different in the four biogeographical areas as some species are specific to one biogeographical region (e.g. the Sardinian warbler Curruca melanocephala in the Mediterranean area and see species repartition in the Supporting information). Note that we used the most detailed habitat information available but, as finer habitat grain is out of reach, we admit that not all the influence of the Grinnellian filter has been removed by this first grouping (which captures 5.36% of the co-abundance variance, Supporting information). Step 2. In each biogeographic region and habitat, we used the log-transformed co-abundance data (to obtain normally distributed data) to calculate observed associations as partial correlations between each pair of species (Schäfer and Strimmer 2005) as follows (Eq. 1): Pc Oi,j ij ii jj () = - () () ´ () - -- å åå O OO , ,, 1 11 (1) with O the matrix of observed abundance (species × sites), Pc(O)i,j the partial correlation between species i and j, and Σi,j−1 the value for species i and j of the inverse of the covariance matrix. This approach partially removes the indirect effects of other co-occurring species on the estimated association between the two considered species by focusing on the conditional association (Harris 2016). Step 3. Partial correlations can be affected by species commonness, since common species have higher probabilities to co-occur than less abundant species because of a higher representativeness in the data (Blüthgen et al. 2008). To correct this bias, we computed partial correlations on 1000 random co-abundance datasets obtained by keeping constant the total number of individuals in a given sampling point, and assuming that the probability for a species to occur in a given sampling point was proportional to its frequency in the dataset. We then calculated standardised effect sizes of partial correlations between species i and j (SESi,j) as follows (Eq. 2): SESPc OPcN Pc N i,j i, ji ,j i,j = () - () () () () m s (2) where Pc(O)i,j is the observed partial correlation between species i and species j, μ(Pc(N))i,j and σ(Pc(N))i,j the mean and standard deviation of partial correlations from the 1000 randomly sampled datasets. Step 4. In order to identify significant associations, we calculated a two tail p-value for each pairwise association using the rank of the observed association in the Gaussian distribution of null associations obtained from step 3. That is, we determined the number of replicates for which the absolute value of the observed partial correlation is greater than the absolute null partial correlation (p-values were corrected for multiple comparisons following Benjamini and Hochberg (1995)). Significant associations therefore corresponded to SESi,j for which adjusted p-values were below 0.05. Step 5. For each species pair, for each biogeographic region and for each habitat, we averaged the significant associations over the 17 annual associations (one for each year, i.e. if an association was not significant for some years, these years were removed and the averaged association was calculated using association values from the remaining years). If there was no significant association across the 17 years, the association for this given species pair in this given biogeographic region and this given habitat was set to zero. This step was carried out for all 260 191 combinations (5886 pairs of species × four biogeographic regions × 19 habitats, but note that some species pairs were not observed in all biogeographical areas and all habitats). Association network indices We considered three mathematically independent indices that describe different aspects of the network built upon pairwise association estimates (Fig. 2b and see examples in the Supporting information). The three indices were calculated for 121 172 species association networks corresponding to the communities monitored in the 19 580 sampling points between 2001 and 2017 (17 years × 19 580 sampling points, but note that each sampling point was not necessarily monitored on each year). Intensity I quantifies the strength of associations in the species association network of a community. It reflects the average intensity of the associations in the network. It is weighted to account for the differences between species’ abundances (Eq. 3). 6 Figure 3. Workflow for estimating species association networks. For each biogeographical region and each habitat (step 1), an observed association matrix was obtained by partial correlations (step 2) (positive correlations in orange, negative correlations in purple, colour intensity proportional to correlation strength) from co-abundance data (species D, K, G, U as an example). Random association matrices were calculated using partial correlations on permuted datasets (1000 times) and were used to calculate standard effect sizes (step 3) of observed associations as well as their adjusted p-values to obtain significant SES of the observed association matrix (step 4). Step 1–4 were repeated for each year providing annual associations, which were then averaged over years for each species pair. species associations were finally added to the spatial co-abundance data to obtain a species association network for each of the sampling points (step 5). 7 I n’ n’ ij ij j n i n ij j n i n = ´ == == åå åå a 11 11 (3) with n the number of species, n′ij the number of pairs of species i and j in the community pool and αij the association (as defined in step 5) between species i and j (with i ≠ j; when i = j, αij = 0). I varies between 0 and |α|max. High values of I are reached in communities including mainly strong associations. Attractiveness A quantifies the prevalent sign of the associations as the number of positive associations minus the number of negative associations standardised by the total number of associations (Eq. 4). Attractiveness is analogous to the association ratio in plant networks (Saiz et al. 2014). However, we choose to label this metric attractiveness rather than association ratio because it also stands for methods estimating associations (Chiyo et al. 2011). A=- + - - pp pp + + (4) with π+ the number of positive associations and π− the number of negative associations. It varies between −1 (if all the associations are negative) and 1 (if all the associations are positive). Clique structure C quantifies the level of structuring of the species association network. It is calculated using the number of existing cliques (i.e. fully connected groups of species (Luce and Perry 1949)) with three or more nodes (obtained with the R package igraph, Csardi 2013), standardised by the number of potential cliques in a given network (Eq. 5). Cobs max obs 3 obs = = = c c c n i c n i n i i= n i= n i= æ è çö ø ÷ æ è çö ø ÷-æ è çö ø ÷ å å 00 22 å - () - =c nn+ n obs 21 2 1 (5) with cmax the maximum possible number of 3to n-cliques, cobs the observed number of 3to n-cliques, n the number of species in the network. C quantifies the complexity of the network architecture resulting from the interweaving of associated species (Supporting information). Networks with high C values have a complex structure, with multiple imbricated groups of interconnected species. Networks with low C only have a few small sized interconnected groups of species. Spatial estimates and temporal trends in association network indices Spatial estimates of association network indices For spatial analyses, we averaged the annual values of each index (I, A and C) for each sampling point, resulting in one value for each index for each sampling point. We then computed the spatial window values of each association network index, for each site and for each year, using an 80-km radius window. We determined the window size as a compromise between a large spatial coverage and a fine spatial resolution and we conducted the analyses for various radii to assess the robustness of our results to changes in the window size (Supporting information). Spatial window values were computed to analyse, on a similar spatial scale, the relationships between community indices and β-diversity which is an inter-site measure based on species data from multiple sites. It also provided more complete data when sampling points or sites were not monitored every year, in particular for calculating temporal trends. We estimated spatial window values using geographically weighted regression (GWR, using the R package spgwr, (Brunsdon et al. 1996, Gollini et al. 2015). In this approach, the centre of each site was consecutively considered as the centre of a fixed radius window. Each index was calculated using data from all sampling points encapsulated within the spatial window. A weight was attributed to each sampling point, which decreased with the distance to the central selected site following a bisquare kernel function. Temporal trends of association network indices We estimated the spatial window trends as the temporal trend of each association network index (I, A and C) following the same framework as for spatial window values. The trend of each index corresponds to the coefficient of a linear regression calculated using annual index values in the selected sampling points, weighted according to their proximity to the central site. In addition to being calculated on the whole dataset, spatial and temporal window values were also computed for each of the three main groups of habitat (woodland, grassland and human settlements). In this case, we selected sampling points that 1) belonged to the type of habitat considered, and 2) were located within sites dominated by this type of habitat. We recalculated spatial averages of network indices (too few trends were significant to conduct temporal analyses) using the subset of sampling points for each of the three main groups of habitat. Species turnover and community generalism Spatial and temporal variation We assessed biotic homogenisation using the spatial and temporal variations of species turnover (using β-diversity) and community generalism (using community generalisation index, CGI) following the same framework as for spatial window values of I, A and C. β-diversity corresponds to species diversity between a set of sites and is the result of two components, species turnover sensu stricto and nestedness 8 (Baselga 2010). A raw observed species turnover between two communities, κ1 and κ2, can result from a simple difference in community size if κ2 includes fewer species. In this case, part of the observed turnover is due to the nested composition of κ2 in the composition of κ1. That is, a simple species loss from κ1 to κ2 generates a turnover due to a variation in community size. In contrast, the species turnover s.s. corresponds to the species not shared by the two communities, i.e. resulting from the replacement of one species by another. We first aggregated species data from sampling point level to site level to obtain species data for each site. We then randomly selected 10 sites in each spatial window (Devictor et al. 2010b) and computed the species turnover s.s. (hereafter referred to as species turnover) of the set of sites in that window using the betapart R package (Baselga and Orme 2012). We repeated this selection step 10 times, and we took the mean of β-diversity. We used the species generalisation index (SGI) from Godet et al. (2015) to calculate for each community (i.e. at each sampling point) a community generalisation index (CGI). SGI corresponds to the habitat specialisation value of a species independently from its abundance and range. For the 109 bird species, SGI ranges from 0 to 46.06 (mean = 28.73, SD = 9.74). CGI represents the average habitat specialisation of species in a given assemblage, weighted by their local abundances. CGI by sampling point ranges from 0 to 40.26 (mean = 28.26, SD = 3.86) and by site (after GWR) from 23.34 to 32.05 (mean = 28.19, SD = 1.61). Association network indices versus species turnover and community generalism We analysed the spatial relationships between species turnover (β-diversity) and the three association network indices as well as the spatial relationships between community generalism (CGI) and network indices by performing generalized additive models (GAM, using the R package mgcv, Wood and Wood 2015) to assess the linear relationship between each of the three network indices and species turnover or community generalism, while explicitly modelling the spatial autocorrelation (Eq. 6). That is, each index (I, A and C) of association networks i was successively considered as the response variable R (Gaussian family, link identity) regressed over β-diversity or CGI (explanatory response E). We explicitly modelled the spatial autocorrelation using a tensor product of a thin plate regression spline based on geographic coordinates (longitude lon and latitude lat) of sites following Wood (2003, 2017). As species richness may influence network indices, we added species richness SR as a covariable to disentangle the effect of species turnover or community generalism from the effect of species richness on network indices. In addition, as part of the results could be driven, at least to some extent, by unchecked structural relationships between association indices and species turnover, we tested whether the observed relationships could be due to the intrinsic redundancy between network indices and β-diversity using simulated and permuted data (Supporting information). RESR f iiii ii =+ ++ () +ab ge lonlat, (6) with α the intercept, β the effect of the explanatory variable (β-diversity or CGI), γ the effect of species richness and εi ~ N(0,σ). Using a similar model, we also tested the temporal relationship between the trends in the three association network indices and the trend in species turnover or community generalism (controlling for the trend in species richness) using their spatial window trends. Limits of relying on space-fortime substitution (i.e. relying only on spatial gradient to infer temporal relationships) are well documented (Damgaard 2019) and required to complement spatial relationships with temporal ones when possible. In this study, this temporal analysis was possible at the national scale but not in each of the main types of habitats. Diagnostics for all models are in the Supporting information. Results Species associations from co-abundance We found 8.1% of positive associations, 38.3% of negative associations, whereas 53.6% of associations were non-significant. 40% of the species pairs showed qualitatively constant associations (i.e. significant associations that were positive or negative in more than 90% of cases) across habitat/biogeographic region combinations. On average, each species was associated with 1–93 other species (mean = 41, SD = 26) with variations between habitats and biogeographic regions (associations available in the Supporting information). In particular, in woodland, we found 4.3% of positive associations (e.g. between the Eurasian wryneck Jynx torquilla and the lesser spotted woodpecker Dryobates minor), 28.7% of negative associations (e.g. between the European robin Erithacus rubecula and the common cuckoo Cuculus canorus) and 67.0% of nonsignificant associations. In grassland, we found 7.8% of positive associations (e.g. between the calandra lark Melanocorypha calandra and the Eurasian skylark Alauda arvensis), 29.1% of negative associations (e.g. between the red-legged partridge Alectoris rufa and the grey partridge Perdix perdix) and 63.1% of non-significant associations. In human settlements, we found 4.6% of positive associations (e.g. between the house sparrow Passer domesticus and the Eurasian collared dove Streptopelia decaocto), 30.1% of negative associations (e.g. between the house sparrow Passer domesticus and the carrion crow Corvus corone) and 64.7% of non-significant associations. Relationships between association network indices, species turnover and community generalism Intensity was positively related to species turnover in space and time (Fig. 4a, Table 1) and negatively to community generalism in space and time (Fig. 4d, Table 1). Attractiveness was negatively related to species turnover in space and time (Fig. 4b, Table 1) but negatively related to community