Do large‐scale associations in birds imply biotic interactions or environmental filtering?

There has been a wide interest in the effect of biotic interactions on species' occurrences and abundances at large spatial scales, coupled with a vast development of the statistical methods to study them. Still, evidence for whether the effects of within‐trophic‐level biotic interactions (e.g. competition and heterospecific attraction) are discernible beyond local scales remains inconsistent. Here, we present a novel hypothesis‐testing framework based on joint dynamic species distribution models and functional trait similarity to dissect between environmental filtering and biotic interactions.


| INTRODUC TI ON
Biotic interactions are an integral part of communities and have a profound effect in shaping diversity (Gaüzère et al., 2022;Ratzke et al., 2020). These interactions take place at the level of individuals: negative interactions, such as competition (Cadotte & Tucker, 2017), result in decreased individual's fitness because of individuals of other species, whereas positive interactions, such as facilitation (Bruno et al., 2003), result in increased fitness. These local effects have been shown to scale up, and shape species' occurrences at large spatial and temporal scales (Blois et al., 2013;Bruno et al., 2003;Cavieres et al., 2014;King et al., 2021;Wisz et al., 2013). Still, studying biotic interactions at these macroecological scales has remained challenging as interactions are most often inferred from observational data on species occurrences, which are also shaped by other mechanisms (environmental filtering, dispersal and drift; Mutshinda et al., 2009). Here, we present an approach based on time series of species abundances and functional traits for disentangling biotic interactions from environmental filtering at a community level to explore the prevalence of biotic interactions at large spatial extent.
An increasingly popular method for inferring biotic interactions from observational data is joint species distribution modelling (JSDMs; Ovaskainen et al., 2017;Pollock et al., 2014;Thorson et al., 2015;Warton et al., 2015). JSDMs enable simultaneous modelling of occurrences for multiple species, and produce 'speciesto-species associations' that are excess or deficits in spatial or spatiotemporal co-occurrence, relative to a baseline occurrence rate set by the environmental variables or a random occurrence (Dormann et al., 2018). Hence, they capture the component of variation in occurrence not explained by the environmental variables used as covariates and indicate whether an association between each species pair is negative (segregation), positive (aggregation) or random. These species-to-species associations are sometimes interpreted as indicating the direction and magnitude of biotic interactions (D'Amen et al., 2018).
However, it is not straightforward to interpret species-to-species associations from JSDMs, even if the effect of environmental variables has been taken into account. First, it is widely known that exclusion of ecologically important abiotic variables may lead to obscure interpretation of species-to-species associations (Kissling et al., 2012;Poggiato et al., 2021;Warton et al., 2015). Hence, negative associations may be due to different preferences for environmental characteristics rather than negative interactions, such as competition, and positive associations may still be due to similar preferences for environmental characteristics rather than positive interactions, such as facilitation. Second, the use of static snap-shot data from dynamic interactions may not represent biotic interactions reliably (Dormann et al., 2018;Kilpatrick & Ives, 2003). JSDMs only allow modelling of spatial autocorrelation between species' occurrences, ignoring the temporal correlation structure (Dormann et al., 2018). Third, and perhaps most fundamentally, JSDMs model the biotic interactions solely in the residuals, and thus only retrieve the realized niche of species, while the fundamental niche remains uncovered (Poggiato et al., 2021). In summary, the usefulness of JSDMs for inferring biotic interactions has been widely questioned (Barner et al., 2018;Blanchet et al., 2020;Poggiato et al., 2021;Zurell et al., 2018).
A part of the solution for overcoming these issues is to use abundance instead of occurrence data, and studies of JSDMs have pointed out several times that abundance data should indeed provide more reliable inference on biotic interactions (Blanchet et al., 2020;Dorazio et al., 2015;Poggiato et al., 2021; see also Ulrich & Gotelli, 2010). This is because, for example, competition would be shown as decreased abundances of species with abundance data, which is a more realistic outcome of competition than total competitive exclusion (i.e. species not co-occurring). Second part of the solution is to use time-series data on abundance. Indeed, temporal variation in species' abundances has been found to be an effective way to study the relative strengths of environmental variation and biotic interactions (Houlahan et al., 2007;Mutshinda et al., 2009). This is because temporal (co)variation in species' abundances reveals whether two species showing similar average abundances show high abundances at the same time, possibly reflecting positive interactions, or whether one has lower abundance when the other has a high abundance, possibly reflecting negative interactions. Time series of abundances can be studied via joint dynamic species distribution models (JDSDMs) (Thorson et al., 2016). In JDSDMs, the correlation structure between species within different time steps can be modified. For inferring biotic interactions, the most optimal structure would be that the abundance of each species at time t directly depends on the abundance of other species at time t − 1 (Sebastian-Gonzalez et al., 2010;Thorson et al., 2016), thus enabling a causal link in abundance variation between years (e.g. Barraquand et al., 2021). In practice, this is often computationally too demanding for communities with tens of species but even the more simplistic JDSDMs on abundance should yield more reliable inferences on biotic interactions than static models. and functional trait similarity is beneficial in ecological interpretation of species-tospecies associations from data covering several decades and biogeographical regions.

K E Y W O R D S
competition, functional traits, heterospecific attraction, joint dynamic species distribution models, macroecology, VAST Together with JDSDMs, species' functional and ecological traits may be used in building a supplementary modelling framework with a priori predictions (Kohli et al., 2018;König et al., 2021;Mönkkönen et al., 2017;Snell Taylor et al., 2020). However, identifying the traits that reflect biotic interactions requires careful consideration. For instance, as species compete for food, traits associated with feeding ecology and behaviour may more directly represent the biotic interactions than some other traits (e.g. body size) where the association with interactions is only indirect. Therefore, methods where species with more similar or dissimilar functional or ecological traits are expected to interact more or less, respectively, than others (Elo et al., 2021;Schleuning et al., 2015) may improve interpreting whether associations from JDSDMs result from biotic interactions.
Here, we present a novel hypothesis-testing framework based on JDSDMs and functional trait similarity to dissect between biotic interactions and environmental filtering. We demonstrate our framework by studying within-trophic-level interactions at large spatial grain and extent with long-term breeding bird data from France and Finland. In bird communities, species show both competition (Cody, 1974;Martin & Martin, 2010;Robinson & Terborgh, 1995) as well as positive interactions, such as facilitation (Martin & Eadie, 1999) within a trophic level. Also, birds may base their breeding site choice on information provided by an individual of another species (Parejo & Avilés, 2016;Seppänen et al., 2007).
This so-called heterospecific attraction may result in positive spatiotemporal associations between potential competitors (Kivelä et al., 2014;Mönkkönen et al., 1990;Thomson et al., 2003). Both negative and positive local-scale interactions among birds can scale up to and be visible at large spatial scales (Belmaker et al., 2015;Gotelli et al., 2010;Heikkinen et al., 2007;Mönkkönen et al., 2017). By contrast, some other studies have found little support for biotic interactions at large spatial scales (Dorazio et al., 2015;König et al., 2021;Royan et al., 2016;Sandal et al., 2022), at least when compared to environmental drivers and averaged over species (Snell Taylor et al., 2020). Hence, while environmental filtering clearly affects bird community assembly, the importance of biotic interactions among bird species at large spatial scales is still unclear.

| Analytical framework
We build our framework on abundance-based JDSDMs and information on species functional traits. From JDSDMs, we derive spatial and spatiotemporal associations among species' abundances ( Figure 1).
Spatial associations estimate covariance in species' abundances, averaged across years, that are not explained by the environmental variables included in the model. They roughly correspond to the classical residual covariances of JSDMs. Spatiotemporal associations represent the covariation in abundances in space and time not captured by the spatial covariance. Thus, while spatial associations are likely to capture the effect of shared responses to missing environmental variables, the spatiotemporal associations may carry the sign of biotic interactions.
Using a phylogeny and functional traits to separate environmental filtering and biotic interactions is not straightforward as the two mechanisms may produce similar patterns in relation to phylogeny and functional traits, especially when studying co-occurrence (Cadotte & Tucker, 2017;Mayfield & Levine, 2010). In our framework, we develop hypotheses for environmental filtering (H1) and biotic interactions (H2) which predict unique combinations of (i) the presence of negative spatiotemporal associations, (ii) the slope and (iii) heteroscedasticity in the relationship between spatiotemporal association and species functional dissimilarity ( Figure 2). Environmental filtering leads functionally similar species to have the strongest positive spatiotemporal associations and dissimilar species to have the strongest negative spatiotemporal associations (H1). Thus, the slope of the relationship between spatiotemporal associations and functional dissimilarity is predicted to be negative ( Figure 2). Biotic interactions, on the other hand, lead to different predictions depending on whether we consider negative interactions such as competition, positive interactions, such as heterospecific attraction, or both.
Competition (H2a) predicts that species' spatiotemporal associations are negative and given that similarity of the functional traits triggers competition, the negative spatiotemporal associations are strongest between functionally similar species (Beaudrot et al., 2019;König et al., 2021;MacArthur & Levins, 1967). This results in a positive slope between spatiotemporal association and functional dissimilarity ( Figure 2) which is a unique prediction for competition and not predicted by environmental filtering, although otherwise their effects may be hard to discern (Mayfield & Levine, 2010). By contrast, heterospecific attraction predicts similar species to have the strongest and positive spatiotemporal associations and dissimilar species to have random spatiotemporal associations (H2b). This is because heterospecific attraction is based on information (e.g. suitable nesting site) acquired from other species, and the information value decreases with increasing ecological distance (Seppänen et al., 2007).
We note that the resulting negative slope between spatiotemporal association and functional dissimilarity is also predicted by environmental filtering (H1a, Figure 2). Together competition and heterospecific attraction result in similar species to have the strongest spatiotemporal associations, either positive or negative, whereas dissimilar species have random spatiotemporal associations (H2c).
Consequently, the heteroscedasticity of the spatiotemporal associations decreases with increasing functional dissimilarity ( Figure 2). Finally, the relationship between spatiotemporal association and functional dissimilarity described above applies to species with positive spatial associations but not for species with random or negative spatial associations, and this holds for both environmental filtering and biotic interactions ( Figure 2). If species' spatial associations are positive, they are more likely to be present in the same habitats, and thus the potential for direct interaction is high or the species are more likely to use the same resources. By contrast, with random or negative spatial associations, probability of direct interactions and using same resources decreases because the probability of cooccurrence decreases, and the pattern dilutes.

| Field data
We used annual census datasets for terrestrial birds collected dur-

| Data preparation
Our goal was to dissect biotic interactions from environmental filtering. Hence, we aimed for removing the most obvious environmental factors causing variation in abundance among species.
First, we divided the French data into four biogeographical regions (Atlantic, Continental, Alpine and Mediterranean) according to European Environment Agency (www.eea.europa.eu/) and similar to Barnagaud et al. (2012). These regions coincide closely with previous studies (e.g. De Heer et al., 2005) of breeding birds at biogeographical scale. We further divided the two largest regions into Northern F I G U R E 1 A simplified figure showing the difference between spatial and spatiotemporal associations between species. The data include abundances of three species (A, B, C) in 3 years and three sites (panels a, e; b, f; c, g). Spatiotemporal association reflects the spatiotemporal correlations among log e -transformed abundances of species (panel d). Here, species A and species B have strong positive spatiotemporal association, whereas species C is negatively associated with the other two species. However, when averaged across years (panels e-g) the log e -transformed abundances of the three species are highly positively correlated, and consequently, the spatial associations of all species pairs are positive (panel h). Here, the spatial and spatiotemporal associations are simply calculated as correlations between log e -transformed abundances, whereas our multivariate spatiotemporal delta model (VAST) uses a latent variable structure on the explicit spatiotemporal data to estimate the associations, while taking into account the environmental variables used as covariates. For both habitat types, we only included the terrestrial species breeding and/or foraging in farmlands and forests, separately, according to delHoyo et al. (2014). We excluded birds of prey and grouse, as the point count method used is not reliable for counting them (Andersen, 2007;Pakkala et al., 1983). Therefore, we consider the data to represent within-trophic-level associations. In the French data, the censuses were done twice per year, and the maximum number of individuals for each species in each point and year detected between the two censuses is considered as a proxy of the local abundance for that species on that point. In Finland, only one census per year was conducted because breeding is more synchronized due to the short breeding season at high latitudes. Each observed bird, pair or brood is considered as a breeding pair in the Finnish survey (Koskimies & Väisänen, 1991). Hence, we multiplied the observations by two to obtain the number of breeding adult individuals.

| Environmental variables
We obtained mean temperature (°C) and total precipitation ( (H2a), heterospecific attraction (H2b), and combination of competition and heterospecific attraction (H2c) for the relationship between species' spatiotemporal associations ('STA') and functional dissimilarity ('FD'). The pattern is shown as a small figure (each dot represents a species pair) and the slope, heteroscedasticity, and requirement of the presence of negative spatiotemporal associations ('Some STAs are negative') is presented. For instance, the hypothesis (H2c) predicts that the slope in the relationship is null, the heteroscedasticity of the data decreases with decreasing functional dissimilarity, and the data include negative spatiotemporal associations. For all hypotheses, it is predicted that the pattern is the strongest for species pairs having positive spatial associations. Sign '−' denotes that the hypothesis does not have a particular prediction for the pattern feature.

Hypothesis
The relationship between species spatio-temporal associations (STA) and functional dissimilarity (

| Vector autoregressive spatiotemporal model
For estimating species-to-species associations, we used a multivariate spatiotemporal delta model (R package 'VAST', version 3.7.11; Thorson, 2019;Thorson & Barnett, 2017). VAST uses the latent factor structure to jointly model abundances of multiple species in space and time, while accounting for spatial and spatiotemporal we specified a gamma probability distribution. In this study, we focused only on the first component of the delta model, that is, numbers density, because the variation in average individual biomass covariance matrices was very low. We included two dynamic environmental covariates to the models (year-and site-specific mean breeding-season temperature and total precipitation). Furthermore, we assigned all the sampling points within the 20 km × 20 km cell to the same location; this equals the grain size of the study. Previous studies Heikkinen et al., 2007) have shown that at least a grain size of 10 km × 10 km is robust for detecting biotic interactions at a macroecological scale.
We focused on the most common species in each dataset as the focal community consists mainly of these species (we included the most abundant species representing 80% of the observed individuals). Their species-to-species associations can be modelled reliably, whereas the abundance vector for rare species consisting mainly of zeros, contains only little information. This resulted in altogether 59 species and different number of species in each dataset (Appendix 3, and forest), the models did not converge, most likely due to the complexity of the models. We ran a similar VAST model as described above but without the environmental covariates, for the most abundant species representing 90% of the observed individuals, and with no autocorrelation structure. However, these two models produced different spatial and spatiotemporal associations (Appendix 4, Figure S4.1), and we do not consider the latter models further.
We converted the spatial and spatiotemporal covariances into correlations, referred to as spatial and spatiotemporal associations from now on, which ranged from −1 (strongest negative association) to 1 (strongest positive association). We classified both spatial and spatiotemporal associations as negative (−) or positive (+), according to the sign of the association, if the 95% confidence intervals of the correlation did not encompass zero. Otherwise, we classified spatial and spatiotemporal associations as random (0), meaning that statistically the association did not differ from zero.

| Species' functional dissimilarity
For measuring species' functional dissimilarity, we used three approaches: (i) morphological dissimilarity, (ii) diet dissimilarity and (iii) their combination. First, bird morphology is strongly correlated with species' feeding, locomotion and habitat use (Miles & Ricklefs, 1984). We included species' log 10 -transformed body weight as an overall indicator of the size of the species, as similar body weight in general (Kohli et al., 2018)

| Linking species' spatiotemporal associations and dissimilarity
To test our hypotheses, we used generalized least squares (GLS) linear models (function 'gls' from the R package 'nlme' version 3.1-153; Pinheiro et al., 2021) fitted with the maximum likelihood method.
GLS models enable us to model the slope as well as the heteroscedasticity predicted by our hypotheses (Figure 2). For each dataset, we fitted a GLS model with spatiotemporal association of a species pair as a response variable and spatial association group (factor: SA+ = positive, SA0 = random, SA− = negative), functional dissimilarity and their interaction as explanatory variables. We also added habitat dissimilarity as an explanatory variable to control for the differences between habitat preferences. Hence, we were able to test whether the slope in the relationship between spatiotemporal association strength and functional dissimilarity is negative (H1, H2b), positive (H2a) or null (H2c), and whether the relationship is strongest to species with positive spatial association (H1, H2a-c; Figure 2). To test whether the heteroscedasticity in the spatiotemporal associations decreased with increasing functional dissimilarity (as predicted by H2c; Figure 2), we ran (i) a model without any modelling of heteroscedasticity and (ii) a model where heteroscedasticity was modelled with the 'varExp' variance function (Pinheiro et al., 2021)

| Spatial and spatiotemporal associations
Spatial associations were most often positive (45%) or random (37%), and least often negative (18%). Spatiotemporal associations were mostly positive (80%), followed by random (15%), and only 5% were negative (Figure 3; Appendix 5, Table S5.1; Appendix 6, Figure S6.1-4). Majority of the species pairs had positive spatiotemporal associations coupled with a positive (39%) or random (29%) spatial association but also birds with negative spatial association showed positive spatiotemporal associations (3%-25%, depending on the dataset) (Appendix 5, Table S5.1). Negative spatiotemporal associations were due to a few species in those datasets where they occurred (e.g. Figure 3b,d,j) and non-existing in other datasets ( Figure 3h,l). Particularly, negative spatiotemporal associations for species with positive spatial association represented 0%-5% of all species pairs, depending on the dataset. Although there was some variation between the datasets (Appendix 5, Table S5.1), these general trends were clear.

| Linking spatiotemporal associations to species' dissimilarities
For the relationship between spatiotemporal associations and functional dissimilarity, the best GLS model in terms of AIC included a combination of morphology and diet as a measure for functional dissimilarity in five datasets, morphology in two and diet in one (Table 1).
We found five negative slopes for the relationship between spatiotemporal associations and functional dissimilarity among species with positive spatial associations (Table 1, Figure 4). This fits the predictions of environmental filtering (H1) and heterospecific attraction (H2b) (Figure 2). For the other three datasets, the mean slope was zero, as predicted by joint effects of competition and heterospecific attraction (H2c). For Continental South farmlands, heteroscedasticity tended to increase with increasing dissimilarity in diet (Table 1).
However, residual variance was not greatest in the least dissimilar species, as predicted by competition and heterospecific attraction (H2c), but for species with average dissimilarity (Figure 4d). For other datasets, the improvement in model fit due to the inclusion of the variance function was small or even reduced (Appendix 7, Table S7.1).
Thus, we did not find evidence for joint effects of competition and heterospecific attraction (H2c). Also, the very few combinations of positive spatial and negative spatiotemporal associations refute the joint effects of competition and heterospecific attraction (H2c) as F I G U R E 3 Spatial and spatiotemporal associations ('Correlation') between bird species' abundances in different biogeographical regions and habitats in France: Alpine forest (a, b), Atlantic North farmland (c, d), Continental North farmland (e, f) and forest (g, h), Continental South farmland (i, j) and forest (k, l); and in Finland: farmland (m, n) and forest (o, p). Spatial and spatiotemporal associations were estimated with a multivariate spatiotemporal model VAST including two environmental covariates. The estimated parameter values for which 95% confidence intervals include zero are shown with a black cross. Please see Appendix 6, Figure S6.1-4 for the corresponding figures with species' names. Correlation well as competition (H2a). As expected, species pairs had decreasing spatiotemporal association with increasing habitat dissimilarity in several datasets (Table 1).

| DISCUSS ION
Despite the broad interest in biotic interactions at large spatial scales and the rapid development of the statistical methods to obtain them from observational data, it is still unclear whether within-trophiclevel biotic interactions are discernible beyond local scales. Here, we presented a novel hypothesis-testing framework based on JDSDMs and functional trait similarity to dissect between competition and environmental filtering. We demonstrated our framework with interactions within a trophic level at large spatial grain and extent with long-term breeding bird data from France and Finland. We found that a vast majority of the species had positive spatiotemporal associations, and these were particularly strong for functionally similar species. Our results thus refute the unique prediction by competition (i.e. a positive slope between spatiotemporal association and functional dissimilarity) and support the view that environmental filtering and positive interactions, such as heterospecific attraction, dominate assembly rules within the most common species in the communities.
In general, approximately 80% of species pairs had a positive spatiotemporal association: they tended to be abundant at the same place and at the same time. Although this aggregation was somewhat expected as we focused on the species with a specific habitat type (i.e. either in forests and farmlands), we expected that controlling for the main environmental variables (temperature and precipitation) and biogeographical region would lead to more variation in spatiotemporal associations. For instance, it could have been expected that competition is more intense with areas of low environmental stress and decreases in areas with high abiotic stress (Cavieres et al., 2014), that is, alpine and boreal regions in our case. However, the prevalence of positive spatiotemporal associations was independent of the biogeographical area or habitat type. Moreover, positive spatiotemporal associations were also found for birds with negative spatial associations. This suggests that the common forest birds have an overall synchrony in their abundances, which may result from the variation in the general habitat quality, not restricted to specific resources, of forests and farmlands.
TA B L E 1 Results from linear models using generalized least squares for the spatiotemporal association and functional dissimilarity between bird species pairs in different biogeographical regions and habitats in France and Finland ('Dataset').

Morphology & Diet
We found not only the dominance of positive spatiotemporal associations, but also that the positive spatiotemporal associations were the stronger, the more functionally similar the species were.
Thus, the aggregations were most likely due to different environmental resources which species can efficiently track. However, we must bear in mind that the pattern is predicted by both environmental filtering and heterospecific attraction as our approach is not able to discern the two mechanisms. Both mechanisms rely on the occurrence of the resources. In environmental filtering, individuals seek preferred resources irrespective of others, while in heterospecific attraction, individuals are used to locate the shared resources more efficiently (Seppänen et al., 2007). Indeed, species for which heterospecific attraction has been empirically demonstrated (e.g. Parus major and Ficedula hypoleuca; Forsman et al., 2002) showed both positive spatial and positive spatiotemporal association (PARMAJ and FICHYP in Appendix 6, Figure S6.4c,d). This implies that heterospecific attraction may affect abundances of these cavity-breeding species also at large spatial and temporal scales.
In general, the high mobility of birds and factors affecting breeding habitat selection (Doligez et al., 2002;Morinay et al., 2020) may scale up interspecific interactions far beyond the scale of the breeding territories.
We did not find evidence for competition: only a few species pairs showed negative spatiotemporal association coupled with positive spatial associations, and the unique prediction from competition (i.e. the positive relationship between spatiotemporal association and functional dissimilarity) was not fulfilled. Hence, although it has been shown that bird species' interactions may scale up to regional extent (Belmaker et al., 2015;Gotelli et al., 2010;Heikkinen et al., 2007;Mönkkönen et al., 2017;Snell Taylor et al., 2020), our results suggest that, for common bird species, the signs of competition, if present, produce patterns that are not discernible from those of environmental filtering (Cadotte & Tucker, 2017;Mayfield & Levine, 2010). Moreover, the fact that habitat dissimilarity had a negative relationship with spatiotemporal associations suggests that species sharing habitat preferences tend to be also positively spatiotemporally associated. This supports the view that environmental factors are more important than competition for covariation in F I G U R E 4 The relationship between bird species' spatiotemporal associations and functional dissimilarity in different biogeographical regions and habitats in France: Alpine forest (a), Atlantic North farmland (b), Continental North farmland (c) and forest (d), Continental South farmland (e) and forest (f); and in Finland: farmland (g) and forest (h). The lines represent fitted linear regressions between spatiotemporal association strength and functional dissimilarity, separately for species having positive (red), random (black) and negative (blue) spatial associations. If the association is significant, the line is solid, otherwise it is dashed. The horizontal dashed grey line shows where spatiotemporal associations are zero, and the dotted grey line shows the average trend within all datasets together. Please note that the measure for the functional dissimilarity (morphology, diet, morphology & diet) varies according to which measure gave best fit to the data in terms of Akaike information criterion (AIC) (Appendix 7, Table S7.1) and the scale of the x-axes differ.
Curiously, we found a pattern which was not predicted a priori by either biotic interactions or environmental filtering. This was the positive relationship with functional dissimilarity and spatiotemporal associations for negative spatial associations, and less strikingly for random spatial associations in Alpine forest (Figure 4a). The pattern was driven by strong negative spatiotemporal associations between functionally similar species (which is expected for these species having negative spatial association, and can thus be expected to use different resources) coupled with very strong positive spatiotemporal associations with functionally dissimilar species, namely the great spotted woodpecker (Dendrocopos major) with the Eurasian blue tit (Cyanistes caeruleus), the great tit (Parus major) and the Eurasian blackcap (Sylvia attricapilla; DENMAJ, PARCAE, PARMAJ and SYLATR in Appendix 6, Figure S6.1a,b). This result can partly stem from the complex relationship between the species: as cavity nesters, blue tit and great tit benefit from the woodpecker producing the cavities, as has been also shown for owl species (Heikkinen et al., 2007). On the other hand, the great spotted woodpecker also feeds on other species' eggs and nestlings. It is possible that the nest predation overturns the benefits of creating cavities, resulting in negative spatial associations between the woodpecker and the other species in time t. By contrast, the benefits of cavities made by woodpeckers are gained in time t + 1 (or even later) resulting in the strong positive spatiotemporal association. The relationship with blackcap as a cup nesting species may partly be explained by its habitat preference for older deciduous-tree-dominated forests that include nesting trees for the great spotted woodpecker and consequently cavities for tits.
Our results only apply to the most abundant species in the community, these species including the pied flycatcher and great tit, for which both competition and heterospecific attraction has been experimentally shown . The focus on common species has several advantages. First, in terms of individuals, they make up the majority of the community and consequently also the majority of the interactions at the level of individuals. Second, their species-to-species associations can be more reliably estimated as the abundance vectors for the rare species contain little independent information (mostly consisted of zeros) about their spatial and spatiotemporal associations. Third, focusing on common species automatically excludes the species observed in atypical habitats.
A disadvantage of our approach is that we modelled only a part of the species' communities. To detect the possible role of biotic interactions involving rarer species, selecting just a handful of species that would be most likely to show interactions (e.g. Heikkinen et al., 2007) would provide a fruitful approach. The smaller number of species would allow for full use of JDSDMs: estimating the effect of each species at time t − 1 directly on other species at time t (Barraquand et al., 2021;Sebastian-Gonzalez et al., 2010) as well as asymmetric species-to-species associations, which are far more realistic than symmetric associations .
Dynamic abundance data coupled with species' functional traits have been suggested to improve the reliability of inferences concerning biotic interactions when using large-scale data on observed abundances (Blanchet et al., 2020;Dorazio et al., 2015;Snell Taylor et al., 2020;Thorson et al., 2016;Ulrich & Gotelli, 2010). Here, we show that our approach of combining species trait information with JDSDMs facilitates ecological interpretation of species-tospecies associations based on long-term abundance data. Instead of finding evidence for competition, we showed that processes leading to species aggregation (mixture between filtering and true social attraction) seem to dominate assembly rules within these species. We encourage using this approach for assessing assembly rules in other study systems for deriving ecological generalizations.

ACK N OWLED G EM ENTS
We thank the volunteers in Finland and France for the bird count data collection, and CSC-IT Center for Science, Finland, for computational resources. This study was funded by Kone Foundation

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data used in the analyses and the main codes are openly available in Dryad at https://doi.org/10.5061/dryad.ht76h drjh.