Effects of density, species interactions, and environmental stochasticity on the dynamics of British bird communities.

Our knowledge of the factors affecting species abundances is mainly based on time-series analyses of a few well-studied species at single or few localities, but we know little about whether results from such analyses can be extrapolated to the community level. We apply a joint species distribution model to long-term time-series data on British bird communities to examine the relative contribution of intra- and interspecific density dependence at different spatial scales, as well as the influence of environmental stochasticity, to spatiotemporal interspecific variation in abundance. Intraspecific density dependence has the major structuring effect on these bird communities. In addition, environmental fluctuations affect spatiotemporal differences in abundance. In contrast, species interactions had a minor impact on variation in abundance. Thus, important drivers of single-species dynamics are also strongly affecting dynamics of communities in time and space.

INTRODUCTION

While single‐species dynamics are reasonably well understood, much less is known about the potential importance of various deterministic and stochastic processes on the community level (e.g., Engelhardt et al., 2020; Kinzey & Punt, 2009).

From time‐series studies of single, or few, species, we know that density dependence, i.e., a negative effect of current, or previous year's population densities on the population growth rate, is an important process affecting population dynamics (e.g., Lande et al., 2003; Royama, 1992). This can occur, e.g., through intraspecific competition for limited resources (e.g., Hansen et al., 1999; Newton, 1998) or through interference competition that excludes individuals' access to space (e.g., Denac, 2006; Moksnes, 2004). Brook and Bradshaw (2006) found that about 75% of 1198 invertebrate and vertebrate species, including birds, showed intraspecific density dependence in their per‐capita growth rates. Intraspecific density dependence in birds has been found to vary in strength and form both between (Sæther et al., 2002) and within species (Sæther et al., 2008), influenced for example by species‐specific life‐history strategies such as generation time or by local adaptations within species to spatial variation in the biotic and abiotic environment (Sæther et al., 2002, 2005).

Interspecific competition is also a potentially important structuring force in natural communities (e.g., Cody et al., 1975; Elton, 1946; Gause, 1934; Hardin, 1960). Both intra‐ and interspecific interactions are fundamental to classical niche theory (Hutchinson, 1947, 1967, 1978), according to which stable coexistence of species is possible through resource partitioning (e.g., Lande et al., 2003; Mittelbach, 2012; Tilman, 1982). Interspecific density dependence, i.e., the effects of changes in one species' abundance on that of other species, is considered important in regulating species abundances (May et al., 1974; Rosenzweig, 1995) and ample experimental and observational quantification shows its importance also in birds (e.g., Alatalo et al., 1987; Dhondt, 2012; Gamelon et al., 2019). Central to coexistence theory is that density dependence within species must be stronger than that between species (Chesson, 2000) and recently, Adler et al. (2018) found evidence for this pattern in their review of plant studies.

Population fluctuations are additionally influenced by environmental stochasticity, i.e., between‐year variation in populations' expected survival and fecundity caused by random variation in the environment, such as, e.g., weather variables, pollution, or habitat loss, which affects all, or groups of, individuals exposed to it similarly (Lande, 1993; Lande et al., 2003). Environmental stochasticity tends to reduce the long‐run population growth rate (Lande et al., 2003; Lewontin & Cohen, 1969) and has been shown to strongly affect many species (Dobinson & Richards, 1964; Grant, 1986; Grøtan et al., 2009; Sæther et al., 2016).

Demographic stochasticity, i.e., within‐year variations between individuals in their survival and fecundity (Engen et al., 1998; Lande et al., 2003; MacArthur & Wilson, 1967) is also widely accepted as being an important process in population dynamics. It is known to have the strongest influence at low population sizes, because variation in individual fitness lowers the long‐term population growth rate, increasing risk of extinction (e.g., Lande, 1993; Melbourne & Hastings, 2008). Demographic stochasticity is the major form of stochasticity in neutral theories of community dynamics (Hubbell, 2001).

The scale, both in time and space, at which communities are observed, determines what patterns are perceived (Allen & Hoekstra, 1991; Levin, 1992). Species can perform differently in different environments (e.g., Levine & Rees, 2002; Wiens, 1989), so that their growth rates vary across space, and two species may thus co‐exist in areas that are not optimal for one/both of them, while they might be competitively exclusive in more favorable environments. Additionally, movement can greatly influence occurrence patterns (Martin et al., 2017; Wiens, 1989). In other words, interspecific interactions can affect occurrence patterns not only locally, but on a larger spatial scale as well (Araújo & Rozenfeld, 2014; Gotelli et al., 2010). Also, temporal variation in the environment allows coexistence of species even with exactly overlapping niches, due to temporal openings of local niches (Chesson, 2000; Hutchinson, 1951). To obtain a complete picture of the species community, studies should be conducted for long time periods and at large spatial scales (Ives et al., 2003; Magurran, 2007; Magurran et al., 2010), ideally retaining the species‐level information.

The aim of this study was to increase our understanding of the drivers of community changes over time occurring at different spatial scales using unique large‐scale, long‐term time‐series data existing on British bird communities. We use joint species distribution models (JSDMs) with latent variable structures, taking into account both abiotic and biotic processes in an integrated way (Warton et al., 2015) to examine the relative contribution of intra‐ and interspecific density dependence at different spatial scales, as well as the influence of environmental stochasticity, to spatiotemporal variation in abundance.

Recent studies have investigated various community drivers at different spatial scales (Karp et al., 2018, Frishkoff et al., 2019) based on information on the individual species of the community. However, studies based on time‐series data are still rare (but see Hendershot et al., 2020).

MATERIALS AND METHODS

Common Bird Census data

We utilized estimates of territories from the British Trust for Ornithology's (BTO) Common Bird Census (CBC) and the BTO/Joint Nature Conservation Committee/Royal Society for the Protection of Birds, Breeding Bird Survey (BBS). The CBC was a program monitoring populations of common bird species in farmland and woodland plots using territory mapping from 1962 to 2000 (Marchant et al., 1990). The CBC was replaced by the BBS (using line transects) in 2001, with a period of overlap where the established plots were continued, to ensure compatibility of the two methods (Freeman et al., 2007). We included only data from 1962 to 2002 (i.e., all data obtained from territory mapping) and used plots monitored continuously for ≥10 years (a few plots had two such stretches of observations; in which case, we included both), resulting in community time series of 10–36 years between 1965 and 2002, from a total of 121 plots (68 farmland and 53 woodland plots) distributed across large parts of the UK (excluding Northern Ireland), see Figure 1. The mean size in square kilometers for farmland plots was 0.8 (minimum, 0.19; maximum, 3.13) and 0.25 for woodland plots (minimum, 0.11; maximum, 0.6). We excluded species ecologically linked to aquatic communities as well as species observed in <30% of the plots. While one of the advantages of joint species distribution modeling is that parameter estimates for rare species are facilitated by information on responses of more common species (Ovaskainen & Abrego, 2020; Ovaskainen & Soininen, 2011), we nevertheless adopted a conservative approach by ensuring a minimum quality of the input data. For the full species list, see Appendix S1: Table S1.

FIGURE 1

Distribution of farmland (F, orange) and woodland (W, blue) plots monitored by the British Trust for Ornithology (BTO; Common Bird Census [CBC] and Breeding Bird Survey [BBS]). Time‐series length (in years) is indicated by circle size, with a minimum of 10 consecutive years in the period 1962–2002 (n = 121 plots)

Environmental data

The effects of environmental variables on population abundances of birds are often complex, including direct, indirect and time‐lagged effects (Pearce‐Higgins et al., 2015; Sæther & Engen, 2010). We applied a sliding time‐window approach (van de Pol et al., 2016) to detect the critical time window during which two key environmental factors, temperature and precipitation, influence the abundance of each species the most. We extracted monthly mean temperature (°C) and precipitation sum (mm) for each plot from the corresponding 5‐km2 resolution climate data grid obtained from the Met Office UKCP09 databases (13.09.2017; MetOffice, 2017) for the duration of the plot‐specific sampling period (>10 years, ranging from 1965 to 2002). There were only 11 out of our 121 plots for which the climate data from the same grid was shared with one (and in one single case, two) other plot(s). To identify the start and end of the time windows of the species‐specific variables that best predicted the observed population dynamics, we followed the regression based sliding time‐window approach by van de Pol et al. (2016). We transformed the estimates of territories n into an estimate of log‐transformed density (number of pairs of individuals per unit area) because of the well‐known influence of area on bird abundances (e.g., Ambuel & Temple, 1983; van Dorp & Opdam, 1987), by , where is the area of the plot, and we added 1 to n to assure positive values. For each species, we evaluated the best combination of time windows of the two environmental covariates based on the Akaike information criterion (AIC) value of the baseline model (a simple first‐order autoregressive model of the log population densities with plot as random effect) compared to a model including the aggregate statistics of a certain time‐window (either of one, both, or none of the covariates). We allowed for the selection of a window spanning up to 12 months prior to June at year t, the last month of the census (hence, the longest possible time‐window tested for extended from June at time t − 1 to May at time t). The best combination sometimes included only one, or even none, of the covariates; see Sandal et al. (2022) for the three top‐ranking combinations for each species.

Hierarchical joint species distribution modeling

We analyzed the data with hierarchical modeling of species communities (HMSC), a JSDM, using the implementation for time‐series data by Ovaskainen et al. (2017). This approach utilizes first‐order multivariate autoregressive [MAR(1)] models, a tool used to infer strengths of both between‐ and within‐species interactions (Ives et al., 2003; Novak et al., 2016; Wootton & Emmerson, 2005). MAR(1) models can be used as relatively simple approximations to complex natural systems, and essentially represent multispecies Gompertz dynamics (Certain et al., 2018; Ives et al., 2003). We followed the approach by Ovaskainen et al. (2017), using three different methods of including species interactions, compared to a variant excluding any interspecific interactions. See Appendix S2 for a detailed outline of the model structure. Common for all model variants, the response variable was the vector of log‐transformed species densities. We scaled these so that each species had zero mean and unit variance to make the species comparable among each other. Subsequently, we added plot‐specific constants (the same for all species) so that , i.e., all plots had the same average summed species density (y being the standardized, log‐transformed densities used as the response variable in the models). This step accounted for unknown variation in general habitat quality and sampling effort at the community level, while species‐specific variations are accounted for by the random effect described in eq. 1 below (see Appendix S3 for details on the linear transformations applied to the data before model fitting).

The linear predictor of the standardized log‐transformed density of species i at plot in year () was defined as where includes the fixed effects, the random effects, and the residual, which we assumed to be normally distributed for the log‐transformed densities. The fixed part was further decomposed as , where modeled the environmental effects (by the environmental covariates included in the model) and the species interactions (by the previous years' densities, see the following paragraphs).

As environmental covariates, we used the two species‐specific climatic variables identified before. Habitat type (farmland and woodland) was included as a factor, and log area of the plot was included as a continuous covariate, to account for possible nonlinear effects of plot area. To account for within‐species density dependence, we included the density of the focal species in the previous year as a species‐specific covariate. We calculated the density of the species in the previous year at four different spatial scales, as explained in the section titled Spatial scale of intraspecific density dependence and interspecific interactions.

All continuous covariates were mean centered and scaled to unit variance to make them compatible with the default priors of HMSC. To account for unmeasured covariates that create permanent spatial variation among the plots (i.e., variation in habitat quality among sites), we included plot (n = 121) as a spatially explicit random effect using latent variables with exponentially decaying spatial covariance structure (Ovaskainen et al., 2016). See Appendix S4: Table S1 for the estimates of the spatial scaling of the spatial random effect. Additionally, to be able to quantify the magnitude of temporal variability in the community abundance through time due to environmental conditions affecting all species simultaneously, we included year (n = 38) as a temporally explicit random effect, also with a exponentially decaying covariance structure (to account for the distance in time between sampling events).

We modeled species interactions as where is the number of species, models the effect of species j on species i, and is the standardized, log‐transformed density of species j at plot p the previous year. The interaction term was specified differently for the different variants of the model as follows (see Appendix S2 and Ovaskainen et al., 2017 for more details).

Model variant M1: No interspecific interactions

Only diagonal elements that model intraspecific density dependence were estimated, whereas the effects of interspecific interactions were excluded by fixing for .

Model variant M2: Full interactions

We included as predictors the densities of all the species in the previous year, thus estimating a full matrix of species‐to‐species interaction coefficients.

Model variant M3: Sparse interactions

We assigned for interaction coefficients a multiplicative prior in which the first factor modeled whether the interaction either was zero (with prior probability 0.9) or non‐zero (with prior probability 0.1), and the second factor had the same normally distributed prior as in model variant M2.

Model variant M4: Community‐level drivers

We modeled interspecific competition as Herein, species are assumed to contribute to one or more of “community‐level drivers”, defined as linear combinations of species densities, ecologically corresponding to e.g., total community biomass. As in Ovaskainen et al. (2017), we estimated both the contribution () of each species j to each driver k, as well as the influence of each driver to each species . Additional effects of intra‐specific interactions (beyond those captured by the community‐level drivers) are added to the interaction matrix via the Kronecker's delta , which is 1 for , and otherwise 0. We assumed the prior of Ovaskainen et al. (2017) that increases shrinkage with the index of the driver and fixed the number of drivers to two.

We acknowledge that such a model would be better if implemented as a state‐space model, i.e., accounting for a discrepancy between true and observed population sizes. In its current form, this is not fully possible in HMSC. Observation error is accounted for in the response, but not in the predictor, and it would be better to separately account for the process and observation error. However, our simulations confirm that even the current model does detect species interactions, if present (Appendix S7: Table S1 and Figure S1).

Spatial scale of intraspecific density dependence and interspecific interactions

We tested the performance of the four different model variants by calculating the log‐transformed densities of the species in the previous years at four spatial scales. Thus we applied spatial smoothing, using the exponentially decaying weighting function , where is the pairwise distance between plots in kilometers and models the spatial scale of interactions. We compared four contrasting spatial scales called henceforth local (denoted by L) scale (, thus applying no smoothing), proximal (P) scale (r 1 = 18 km equaling the average distance to the nearest plot), average (A) scale (r 2 = 204 km equaling the average distance between plots), and regional (R) scale (r 3 = 711 km equaling the maximum distance between plots). Note that, in addition to modeling the influence of large‐range interactions, spatial smoothing space may increase the predictive power of the models by reducing noise in the predictor.

Model fitting

The combination of four model variants and four spatial scales results in 16 models. We denote these by MA‐B, where A is the model variant and B is the spatial scale, so that, e.g., M3‐R is the model variant 3 (sparse interaction), with the species densities in the previous year measured at the regional scale. We fitted the 16 models used the R package Hmsc‐R (Tikhonov et al., 2020) assuming the default prior distributions (for details, see Ovaskainen & Abrego, 2020: Chapter 8). We sampled four MCMC chains with 37,500 iterations each, discarded the first one‐third (12,500 iterations) as burn‐in, and thinned the remaining by thin = 100 to yield 250 MCMC samples per chain and thus 1000 samples in total. We evaluated MCMC convergence visually by checking posterior trace plots (see Appendix S6: Figure S1) as well as quantitatively by examining the potential scale reduction factors of the model parameters (see Appendix S6: Figure S2; Tikhonov et al., 2020).

Deriving ecological inferences from the fitted models

To examine the predictive power of the different model variants at the different spatial scales, we applied for each of the fitted 16 model variants a five‐fold cross‐validation strategy, a technique that ensures the independence of model predictions from the data used to evaluate the accuracy of the predictions (Refaeilzadeh et al., 2009). We randomly split the years into five equal folds and fitted each model variant five times, each time masking the observations from one of the folds, and then combining the predictions from different folds to yield predictions for the entire data. For each species and each plot, we measured the predictive power as the correlation between these predicted and the observed species densities. For some cases (i.e., a time series of a particular species in a given plot) the calculation of correlation was not meaningful because of lack of sufficient variation, e.g., due to the species being mostly absent from the plot. To account for this, cases where the observed values differed from the modal value of that case for <5 years were excluded from the calculation of predictive and explanatory power.

We averaged the correlations over the plots to obtain a mean correlation for each species, and we then averaged over the species to obtain a mean correlation for each of the 16 model variants. We successfully validated the ability of the cross‐validation procedure to reveal the true structure of the interaction network with the help of simulated data (see Appendix S7 for details). We evaluated explanatory power by comparing the predictions from models fit with the full data (no years masked) to the observed data.

To estimate how much other variables than species interactions contributed to the explained variance, we utilized the function computeVariancePartitioning() from the Hmsc package (Tikhonov et al., 2020) to partition the explained variance of model variant M1 (i.e., without species interactions) among the fixed and random effects; see Appendix S5 for details. To evaluate additionally the percentage of unexplained variance, we first estimated the amount of explained variance for each species as the total variance (measured in the log‐transformed data, before normalization) multiplied by the mean R 2 value of the model, which in Hmsc corresponds to the proportion of explained variance. We subsequently found the amount of unexplained variance as the total variance minus the amount of explained variance. To facilitate comparison in the variance decomposition among species, we then scaled the sum of both explained and unexplained variance to 1 for each species and calculated the mean values for the percentage of variance allocated to each variable. To express intraspecific density dependence in terms of the effect of previous years' log population densities on growth rate (rather than on the current log population density), we subtracted 1 from the model estimates of the effect of intraspecific density dependence. We used the ggplot2 package (Wickham, 2016) to visualize our results.

RESULTS

Evaluation of alternative models for species interactions

The explanatory and predictive power of the models decreased substantially with the spatial scale of intraspecific density dependence and interspecific interactions, with the local scale L providing the highest explanatory and predictive power (mean explanatory and predictive power ranging from L, 0.37 to R, 0.27 and from L, 0.33 to R, 0.07, respectively). Using the model without interactions (M1) as baseline for comparison, neither explanatory or predictive power increased substantially by including species interactions into our models (M2‐4; Table 1) at any of the spatial scales (local [L], proximal [P], average [A], and regional [R]). If the bird communities were structured by species interactions to a large degree, one would expect M1 to have a low explanatory power, while the other model variants would perform better in comparison. At the local and proximal scale, models M2, M3, and M4 had a somewhat higher explanatory power than the simplest model M1, because these models are more complex and provide a tighter fit. The predictive power for the respective models was slightly decreased, further indicating that the slight increase in explanatory power was due to a higher number of parameters included in these models. Appendix S7: Table S1 and Figure S1 distinctly show the expected increase in both explanatory and predictive power in M2–M4 compared to M1 for a simulated data set with clearly defined species interactions.

TABLE 1

Explanatory and predictive power of the models M1–M4 at the four different spatial scales of species interactions

Scale	Explanatory power				Predictive power
	M1	M2	M3	M4	M1	M2	M3	M4
L	0.366	0.392	0.370	0.369	0.336	0.320	0.334	0.333
P	0.346	0.382	0.351	0.352	0.293	0.275	0.288	0.293
A	0.275	0.367	0.297	0.226	0.111	0.114	0.104	0.095
R	0.267	0.365	0.280	0.168	0.078	0.048	0.064	0.065

Notes: Explanatory power is measured as the mean of the correlation for each species in each plot between the observed time‐series and time‐series predicted by the different models fit to the observed data. Predictive power is measures as the mean of the correlation for each species in each plot between the observed time‐series and independently predicted time‐series obtained by fivefold cross‐validation. Scales are L (local scale), no spatial smoothing; P (proximal scale), average distance to the nearest neighboring plot (18 km); A (average scale), the average distance between plots (204 km); R (regional scale), maximum distance between plots (711 km).

This pattern was retained also when examining the comparison of explanatory and predictive power of models M2–M4 to the baseline model M1‐L at the level of individual species (Figure 2a,b). Especially for M2 we observed a somewhat increased explanatory power compared to M1 for most species, even at larger spatial scales of species interactions, but with a simultaneous decreased predictive power. For M3 and M4, locally, explanatory power was slightly better for most species, but not at any of the larger scales. Including species interactions did not improve the models' ability to predict the dynamics of any particular individual species, as M2‐L, M3‐L and M4‐L did not give substantially better predictions than M1‐L for any of the species (Figure 2b). The predictive power further decreased for all species as the local scale of species interactions was increased. Mean values of explanatory power and predictive power (Table 1) were similar between all model variants at the local and the proximal scale, but at the average, and especially at the regional scale, explanatory power was higher than predictive power, indicating that these models might be overfitted.

FIGURE 2

The difference for each species in (a) explanatory power and (b) predictive power between the focal model and the generally best model M1‐L. The horizontal lines indicate a difference of zero. The four subpanels in each panel, L–R, show results at each of the four spatial scales (L, local; P, proximal [18 km]; A, average [204 km]; R, regional [711 km]; increasingly light color), for the four different models, so that L1 = (M1‐L) − (M1‐L) = 0, L2 = (M2‐L) − (M1‐L), L3 = (M3‐L) − (M1‐L), L4 = (M4‐L) − (M1‐L), P1 = (M1‐P) − (M1‐L), …, R4 = (M4‐R) − (M1‐L). The mean difference is indicated by a dot. Values for M1‐L are exactly zero, since they are being compared with self

In summary, we do not find a strong signal of interspecific interactions. Our simulations show that the modeling approach would be able to capture biologically significant interactions, if present, and we thus conclude that our results are not a statistical artifact.

Relative contributions of intraspecific density dependence and environment

For M1, we partitioned the variance among fixed and random effects, as well as unexplained variance (Table 2). M1 was able to explain on average between 63.6% and 70.8% of the variation in the species' log‐densities across the four spatial scales. Note that the proportion of explained variance was higher than the average explanatory power (0.37). This is because the measure of explained variance includes, e.g., the variation between plots, while such variance is not included in the within‐plot correlations performed to measure explanatory power due to the normalization of the data. The partitioning revealed a strong, localized effect of intraspecific density dependence, as evidenced by a reduction from constituting 20.1% of the variance at the local scale to only 4.8% at the proximal scale. At an even larger average and regional scale, intraspecific density dependence explained only between 0.6% and 1%. The species‐specific temperature and precipitation time‐windows included into the model explained only a small proportion of the explained variation at the local scale. Habitat and log area constituted, relative to density dependence, a moderate amount of variation. Note that only density dependence was subject to spatial scaling, hence it is most meaningful to analyze the partitioning at the local scale. Random plot effects increased from L, 31.7% to R, 53.5%, while the random temporal effect decreased from L, 15.5% to R, 2.6%.

TABLE 2

Average percentage of variance partitioned among fixed and random effects, as well as unexplained variance, for model 1 (no species interactions) at the four different spatial scales of species interactions

Factor	Scale (%)
	Local	Proximal	Average	Regional
Fixed effects
Density dependence	20.1	4.8	0.6	1.0
Habitat	1.2	1.6	2.7	1.9
Log area	1.4	4.9	5.5	4.3
Precipitation sum	0.3	0.2	0.2	0.1
Mean temperature	0.6	0.4	0.4	0.3
Random effects
Plot	31.7	42.4	51.2	53.5
Year	15.5	14.3	4.0	2.6
Unexplained variance	29.2	31.4	35.4	36.3

Notes: Note that only density dependence is subject to the spatial smoothing. Local scale, no spatial smoothing; proximal scale, average distance to the nearest neighboring plot (18 km); average scale, the average distance between plots (204 km); regional scale, maximum distance between plots (711 km).

DISCUSSION

We studied the relative importance of intra‐ and interspecific density dependence and environmental effects on the composition of temperate bird communities. Intraspecific density dependence was of great importance locally (Table 2). The influence of intraspecific density dependence rapidly decreased with increasing spatial scale. Consequently, local population dynamics were largely unaffected by conspecifics at larger spatial scales. Environmental stochasticity, such as temperature and precipitation, showed only a comparably small influence on the community dynamics (Table 2), as will be discussed in more detail later. We did not find a strong structuring effect of interspecific interactions on the bird community, as including interspecific density dependence to the different model variants did not improve explanatory power nor predictive power at any spatial scale (Table 1 and Figure 2).

Inter‐ and intraspecific density dependence in the light of stable coexistence

Model performance was not enhanced by accounting for interspecific density dependence. Similarly, Mutshinda et al. (2009) did not find any strong signature of cross‐correlations in species abundances in their comparative analyses of long‐term studies of communities of moths, fishes, macrocrustaceans, birds, and rodents (data covering up to 40 years in time). In their review of studies in plants, Adler et al. (2018) report a weak influence of interspecific compared to intraspecific effects in the majority of cases.

Our simulations verified that interspecies interactions should be detectable in our data, if present (Appendix S7: Table S1 and Figure S1). Thus, the lack of an increase in model performance is not a statistical artifact and could be attributed either to (1) the nature of data collected in large‐scale, long‐term community ecology studies, (2) an actual low influence of species interactions on the contemporary community composition, or (3) a combination of the two. Naturally, the area of the plots sampled are limited, and there will be differences in species detectability as well as in between‐observer repeatability (Buckland et al., 2011; O'Connor & Marchant, 1982; Williamson & Homes, 1964).

On the other hand, a low signature of interspecies interactions supports the theory of stable coexistence, in which competitive exclusions can be expected to have already taken place (see, e.g., Brown & Wilson, 1956; Connell, 1980; Lack, 1971) or as an ongoing process confining each species to their current realized niches (e.g., Dhondt, 2012; Tuck et al., 2018), leading to the contemporary species assembly in which intraspecific competition dominates over interspecific competition (e.g., Barabás et al., 2016; Chesson, 2000; MacArthur & Levins, 1967). If, for instance, interspecific interactions primarily influence species' carrying capacities, and the species are currently fluctuating around these, we would not be able to detect the interactions that in the past led to the current community composition from the time‐series data. Our findings are in accordance with this theory of niche differentiation, as well as recent studies supporting it (Adler et al., 2018; Barraquand et al., 2018; Tuck et al., 2018): while including interspecific interactions to the model did not enhance explanatory power, intraspecific density dependence clearly stands out as the main variable driving the changes in population sizes across species at the local scale (Table 1 and Figure 2).

In contrast, the neutral theory of biodiversity (Hubbell, 2001) proposes that species can coexist (if only temporarily so) by assuming functional equivalence among species. However, this theory assumes population dynamics to be solely driven by demographic stochasticity and dispersal (Chesson, 2000; Mittelbach, 2012), which is not supported by our findings. Note that, within species, census errors cause autocorrelations in time‐series data that can artificially increase the observed strength of intraspecific density dependence (Freckleton et al., 2006), while across species, census errors inflate the variance of the data and can decrease the cross‐correlations among species. However, density regulation appears to be important in temperate bird species (Sæther et al., 2008), operating through both survival and reproduction (e.g., Grøtan et al. 2009; Perrins, 1979; Sæther et al., 1998), and recent studies have emphasized its importance for the structure and stability of communities (Adler et al., 2018; Kuang et al., 2017; Sæther et al., 2016; van Altena et al., 2016), thus further supporting our findings.

Observations of interspecific interactions at finer scales

While modeling interspecific interactions did not improve explanatory power, it does not mean that interspecific interactions do not take place or are not important. The influence of species interactions acts on the level of the individual and might not be captured by large‐scale data. Many studies investigating species interactions are carried out on a much finer scale, and in a setting where interactions may be more likely to occur. For instance, supplementary feeding is known to increase population sizes to (artificially) higher levels (e.g., Källander, 1981; Robb et al., 2008), likely intensifying competition at such easily accessible study sites (Francis et al., 2018). Nesting boxes may also be preferred over natural alternatives, further intensifying competition. Several studies have shown that e.g., the Blue Tit and the Great Tit compete over nest boxes (e.g., Minot & Perrins, 1986; Newton, 1994), but in a natural setting, their nest‐site preferences, while overlapping, differ (Minot & Perrins, 1986). Competition observed in these settings may thus be relatively localized and might not be representative for the avian community in general.

Fluctuating environments

Temperature and precipitation had a comparably small effect on the species' population fluctuations (Table 2), even if we included the best possible climate windows for each species in each plot separately. The contribution of habitat and log area to the total variance were also relatively small. Note that at larger spatial scales, the variance contributions simply changed in response to the changes made to the predictor variable density dependence.

Our models simultaneously include the effect of variables acting both within and across plots, and variables such as temperature and precipitation that might have a strong effect within single plots can appear to be less important as the focus is moved to a more general picture, as we do here. Previously, Mutshinda et al. (2009), detected a large effect of environmental factors on community composition across taxa, but environmental variance in their study includes both fixed weather variables and also, e.g., random and area effects.

The random plot effect constituted a relatively prominent proportion of the variation even at the local scale, indicating permanent variation in habitat suitability among plots. The importance of the random plot effect increased proportional to the decreased ability of density dependence to explain the growth rates of the local populations as the spatial scale of interactions is increased (Table 2), emphasizing the importance of the scale at which community dynamics are being studied.

The temporal random effect caused by environmental stochasticity (Lande et al., 2003) explains annual variation in community abundance across sites. Compared to the fixed environmental effects, the spatially structured temporal random effect explained more of the total variation, indicating that other environmental factors than those included in our model, have an additional effect on the single species dynamics. One such variable is the so‐called beechmast. In large parts of Europe, seeds from the beech (Fagus sylvatica) are an important, annually highly variable, food source for many birds during winter, and is known to strongly influence next years' densities for many species (Chamberlain et al., 2007).

It is also known that an incorrect specification of density regulation can lead to biased estimates of environmental effects (Lindén et al., 2013). Here, we observed that especially the estimations of the contributions from the two random effects were affected by the degree of spatial smoothing of the density predictor variables. We also note that, in spatial models, the environmental predictors and the spatial random effect can be partially confounded with each other, a phenomenon called spatial confounding (Hanks et al., 2015; Marques et al., 2021). However, we did not find signals of spatial confounding in our study (see Appendix S4).

It is very difficult to capture the patterns observed in nature in all its complexity by the use of simple models. Important factors may remain unaccounted for, and noise in the data generated by, e.g., demographic stochasticity and observational errors might be quite strong. However, the current implementation of HMSC does not allow for the estimation of demographic stochasticity, see Ovaskainen and Abrego (2020: 39).

In summary, utilizing a Joint Species Distribution Modeling approach to investigate the factors structuring a temperate bird community, we found intraspecific density dependence to be of great importance for community composition locally, while climate variables accounted for a small proportion of the variation in the composition of the bird communities. Including species interactions did not result in a better explanatory power or predictive power at any spatial scale. The observed pattern of weak inter‐ but strong intraspecific density dependence is in line with expectations from the theory of niche differentiation. Thus, it is advisable to conduct model selection prior to including interspecific interactions into community models, in terms of developing efficient models while avoiding overparametrization. At the same time, heterogeneity in the dynamics of single species, especially density dependence, is important to include to obtain a realistic understanding of community changes over time occurring at different spatial scales.

CONFLICT OF INTEREST

The authors declare no conflictof interest.

Appendix S1

Click here for additional data file.

Appendix S2

Click here for additional data file.

Appendix S3

Click here for additional data file.

Appendix S4

Click here for additional data file.

Appendix S5

Click here for additional data file.

Appendix S6

Click here for additional data file.

Appendix S7

Click here for additional data file.