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

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,


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 andSchleuning 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 andBenedetti-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 andMontesinos-Navarro 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 andLeibold 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 andRaes 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 andMontesinos-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 andLoreau 2008). Moreover, the responses of species associations to environmental changes 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 andLockwood 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.

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.
structure describes the structural complexity of the association network (Fig. 2b). We then analysed the relationships between the spatiotemporal 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. , 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.
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): (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 (SES i,j ) as follows (Eq. 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 SES i,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). 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). 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).
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).
with c max the maximum possible number of 3-to n-cliques, c obs the observed number of 3-to 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 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.

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 (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 (2003Wood ( , 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).
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.

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 generalism in time (but not in space Fig. 4e). Clique structure was positively related to species turnover in space and time (Fig. 4c, Table 1) and negatively to community generalism in space and time (Fig. 4f, Table 1).
In woodland, intensity was positively related to species turnover but not to community generalism (Table 1, Fig. 4). Attractiveness was not significantly related to species turnover and negatively related to community generalism. Clique structure was positively related to species turnover but not related to community generalism.
In grassland, intensity was not significantly related to species turnover and negatively related to community generalism. Attractiveness was negatively related to species turnover and positively related to community generalism. Clique structure was positively related to species turnover and negatively related to community generalism. Figure 4. Relationships between network indices versus spatial species turnover (β-diversity) and community generalism (community generalisation index (CGI)) for all habitats (black dots), woodland (green dots), grassland (yellow dots) and human settlements (red dots). First row: relationships between (a) intensity and β-diversity, (b) attractiveness and β-diversity, (c) clique structure and β-diversity. Second row: relationship between (d) intensity and CGI, (e) attractiveness and CGI, (f ) clique structure and CGI. Dots correspond to partial residuals of the regression models regressed over predictors and regression lines (solid lines) with confidence intervals (dashed lines) are shown when significant. Coefficients of partial determination (R 2 part ) are provided for each model.
In human settlements, intensity was not significantly related to species turnover and positively related to community generalism. Attractiveness was not significantly related to species turnover but positively related to community generalism. Clique structure was not significantly related to species turnover but negatively related to community generalism.

Discussion
Our study unravelled clear relationships between biotic homogenisation and changes in species associations. These relationships have been revealed thanks to the reconstruction of association networks from co-abundance data and to the ability of tracking modifications in the structure of those association networks. Biotic homogenisation (i.e. the replacement of a diversity of mainly specialist species by a few generalists, McKinney and Lockwood 1999) triggered by ongoing global change (Lockwood et al. 2000, Devictor et al. 2008, Godet et al. 2015 is considered as one of the most pervasive aspects of the biodiversity crisis (Olden et al. 2004). At the local scale, we measured the homogenisation of bird communities as a decrease in species turnover (McGill et al. 2015) and an increase in community generalism. We showed that biotic homogenisation was linked to weaker intensity and clique structure, and more positive attractiveness. In other words, more similar areas in terms of species composition sheltered weaker and relatively more positive associations but less structured association networks.
Networks of detailed interactions between birds remained limited to local communities (Orchan et al. 2013) and association networks to woodland assemblages (Lane et al. 2014) or mixed-species flocks (Mokross et al. 2014). Our results include communities from different habitats at large scale and provide relationships consistent in space and time. They Table 1. Result summary of the GAM models, coefficient estimates (species turnover (β-diversity) or community generalism (community generalisation index, CGI)), standard errors (SE), associated t value, significance level (p-value), coefficient of partial determination (partial r 2 ) and degree of freedom (df). p-value < 0.05 in bold.

All habitats
Species turnover (β-diversity) SE t-value p-value Partial r 2 df emphasize that biotic homogenisation and modifications in association networks are not independent processes. This brings about a new repercussion of environmental change and species community homogenisation (Li et al. 2018).
Overall, intensity decreases with biotic homogenisation (intensity declines with community generalism and increases with species turnover in space and time) implying that homogenised communities are mainly composed of habitat generalist species weakly associated with each other. The negative link between species turnover and community generalism (Supporting information) corresponds to an overall pattern visible at the European level (Le Viol et al. 2012) in which species contributing to increased community similarity are more likely to be habitat generalists. These two metrics are however not totally redundant, since the attractiveness decreases with species turnover but also with community generalism in time. It implies that negative associations are predominant in the differentiated communities (as opposed to homogenised communities), but also that negative associations become relatively predominant in communities where habitat generalists become more numerous. This could result, in particular, from competitive behaviour for instance for nest location or food as previously found in bird communities (Orchan et al. 2013, Lane et al. 2014. Finally, clique structure decreases with biotic homogenisation (clique structure increases with species turnover and declines with community generalism). It suggests that differentiated communities tend to have more complex network structures than homogenised communities. In other words, as the Eltonian filter is less and less visible, biotic homogenisation shapes sparse association networks in communities with weakly associated generalists. However, patterns found between network indices and biotic homogenisation were not similar among all habitats. In forest areas, spatial species turnover was high compared to other habitats and not related to community generalism (Supporting information). This indicates that forest communities were mostly composed of species with the same level of specialisation, as previously found in the French avifauna (Julliard et al. 2006). More particularly, specialists were together with other functionally close specialists, and generalists with other functionally close generalists (Supporting information). In those differentiated (either specialised or generalised) forest communities, associations were stronger and formed more complex networks. However, associations were more negative in communities with habitat generalists than in communities with specialists. This may imply that social information should be more important among habitat specialists and competition among habitat generalists.
Conversely, in human settlements, which are more perturbed than woodland, association intensity and species turnover were not related. Instead, intensity was negatively related to the amount of habitat specialists. The most differentiated communities were shaped by the strong environmental filter formed by the urban environment. This selected mostly specialist species, considered as urban winners (Guetté et al. 2017), forming differentiated communities of species able to deal with this environment. Such a filter selection for specialists in perturbed areas was previously shown on bird species (Gaüzère et al. 2020). But these specialists were functionally far from each other. Consequently, they may not strongly interact as they do not 'know' the other species (Mönkkönen et al. 2017). Social information that can be shared by those species is therefore limited which can explain the low attractiveness observed. This leads to a scenario in which differentiated communities with specialists are composed of weakly and negatively associated species, selected for their ability to subsist in an environment strongly modified by humans, forming sparse association networks. In other words, the Grinellian filter appears to strongly outweigh the Eltonian filter in those communities. Grassland communities had intermediate levels of diversity and specialisation compared to the two other types of habitats (Supporting information). Intensity of associations was not linked to species diversity but to the prevalence of habitat specialists. Specialists, in differentiated communities, were also more negatively associated and formed more complex association networks. Biotic homogenisation in grassland seems therefore to be at the expense of the competing habitat specialists while reshaping association networks toward weakly and less negatively associated generalist species forming sparse association networks.
Several studies have recently shown the difficulties of using species associations as reliable proxies for species interactions (Sander et al. 2017, Freilich et al. 2018, Blanchet et al. 2020. Species associations are indeed potentially affected by nonbiotic filters and some types of species interactions remain inaccessible from co-occurrence (e.g. amensalism (Morales-Castilla et al. 2015)). While our methodology takes into account non-biotic filters, it is still subject to remnant effects of those filters and additional processes linked to life-history traits (e.g. dispersal abilities). That means that some of the species associations we found are still likely to result from, for instance, fine grain habitat filtering. For example, the negative association between the short-toed treecreeper Certhia brachydactyla and the Eurasian robin Erithacus rubecula is probably due to the preference of the latter for young forest whereas the former is rather found in old stands (Laiolo et al. 2004). Another pitfall is the difficulty of estimating temporal variation in species associations from co-occurrence data. Currently, only state-space models allow to quantify species interrelations in varying environments (Deyle et al. 2016) and this approach requires long time-series generally not available across multiple sites and at large scales. This prevented us from estimating temporal variations in associations, although species interactions are known to vary in the short (Price et al. 2005, Olesen et al. 2008) and long term (Li andWaller 2016, Lyons et al. 2016) particularly in response to environmental changes (Rico-Gray et al. 2012, Bimler et al. 2018, Clark et al. 2018. In spite of these limitations, our approach was able to capture pairwise associations that could be related to existing knowledge on bird behaviour and interactions. For instance, the negative associations inside the Parus guild, in particular between the goldcrest Regulus regulus and titmice (Poecile montanus, Lophophanes cristatus) may be related to dominance behaviour of the last two species leading to the spatial exclusion of the goldcrest (Alatalo et al. 1985). Negative associations between the red-legged and the grey partridge are in line with the interspecific competition reported by Rinaud et al. (2020). The positive association between the Eurasian wryneck and the lesser spotted woodpecker can be related to the reuse of lesser spotted woodpecker's cavities by the wryneck (Pakkala et al. 2019). The calandra lark and the Eurasian skylark, positively associated, are known to reciprocally attract each other (Delgado et al. 2013). In addition to these examples, species associations also verified theoretical expectations based on species interactions (see the Supporting information for the relationship between species associations and functional distance). Furthermore, association network indices that we used correspond to aggregated indices that are likely to provide a useful proxy to explore the drivers of changes in ecological communities, even considering the gap between species associations and species interactions (Barner et al. 2018). We are therefore convinced that species associations encapsulate relevant information about the structure of these communities and its changes in space and time.
In conclusion, exploring the fate of species associations provides a new dimension to the biotic homogenisation process: in addition to the homogenisation of species composition, homogenised communities have, in general, weaker, less negative and more simple association networks. Using species associations could also help to discriminate between several forms of biotic homogenisation that took place in different habitats in which the role of Grinnellian and Eltonian filters varies. Accounting for species association has highlighted that differentiated communities, even when composed by habitat specialists, can hide very different processes resulting either in complex and strongly associated networks or very sparse association networks.
Acknowledgements -We warmly thank volunteers contributing to the FBBS. We particularly thank Alexandre Génin for his comments and help. Funding -This project was funded by the ANR project DEMOCOM.

Supporting information
The supporting information associated with this article is available from the online version.