No country for small birds: Potential positive association among medium‐sized, aggressive species in Australian bird communities

Australian woodlands have been intensively cleared since European settlement and, in parallel, many species of birds inhabiting this habitat type have experienced a marked decline. Conversely, some species such as noisy miner Manorina melanocephala respond positively to habitat disturbance and due to their hyper‐aggressiveness can end up driving away specialized and fragmentation‐sensitive species. Recent studies have suggested that the negative impact of miners is exacerbated by means of synergistic interactions between these and other aggressive species, including nest predators. However, it is not clear if these positive associations arise through similar habitat requirements or due to potential mutual benefits (‘facilitation effects’), which should predominate in harsh environments as the stress‐gradient hypothesis (SGH) predicts.


| INTRODUC TI ON
Agriculture has fundamentally altered both the spatial configuration and function of ecosystems worldwide, often resulting in a mosaic of residual natural habitat patches surrounded by a farmland matrix.This phenomenon is particularly severe in some regions of the Australian continent where vast areas of temperate woodlands have been removed since European settlement (Lindenmayer & Fischer, 2006).Temperate grassy of southeastern Australia (the so-called wheat-sheep belt) represents a heavily modified natural ecosystem, with >85% cleared or altered to create grazing pastures and croplands for agriculture (Fischer et al., 2009).The present rate of land clearance in Australia is the sixth highest in the world, and the highest in the developed world (Australia Bureau of Rural Sciences, 2010;Evans, 2016;Reid, 1999).Although this agriculturedominated biome currently hosts high biodiversity value supporting over 150 bird species (Lindenmayer et al., 2016), there is increasing awareness that avian communities inhabiting the open woodlands of southeastern Australia are at risk due to habitat loss and degradation (Mac Nally et al., 2009).In addition to the direct impacts of clearance and its synergistic effect with fragmentation on bird populations, there are other indirect and interactive effects of land clearing on birds.Those include the increased potential for predation and competition by 'winner' species, that is, those that have benefited from the expansion of cleared areas and edge areas (e.g.Mac Nally et al., 2014).
Interspecific interactions like competition can play an important role in shaping community assembly (Wisz et al., 2013).The direction and strength of interactions between species are likely to be context-dependent (Tikhonov et al., 2017).For example, according to the 'stress-gradient hypothesis' (SGH), facilitative and competitive interactions will vary inversely across abiotic stress gradients (Bertness & Callaway, 1994; see also Barrio et al., 2013;García-Navas et al., 2021).Consequently, it has been suggested that facilitative interactions will be more common in conditions of high abiotic stress (e.g.low resource availability) relative to more benign abiotic conditions.On the contrary, under relatively benign environmental conditions, positive interactions should be rare and thus, competitive interactions should be the dominant structuring force (Callaway, 2007;Maestre et al., 2009).While positive indirect associations may be related to shared habitat preferences (species' convergence) and niche partitioning, negative indirect associations may reflect different habitat requirements (species' divergence) instead of the displacement of one species by another (Blanchet et al., 2020;D'Amen et al., 2018).Therefore, considering the effect of environmental variables is crucial to ascertain if the observed non-random associations arise as a result of abiotic factors (e.g.similar habitat preferences) or interspecific interactions.
Joint Species Distribution Models (jSDM) enable the simultaneous modelling of multiple species' responses to environmental variables, indicating whether species-to-species associations are negative, positive or random, while accounting for environmental constraints (Warton et al., 2015).This modelling framework has been recently expanded to incorporate phylogenetic and functional information to examine whether closely related species or species that share functional traits respond in a similar fashion to specific environmental variables (Hierarchical Modelling of Species Communities, HMSC; Elo, Jyrkänkallio-Mikkola, et al., 2021;Ovaskainen & Abrego, 2020;Ovaskainen et al., 2017).From a conservation perspective, this information is critical to detect those clades more vulnerable to the impacts of a specific threat.
Recent studies have shown that competition between generalist species thriving on agricultural landscapes ('winners') and species requiring a minimum proportion of natural vegetation may be an important driver of biotic homogenization of bird communities (Muñoz-Sáez et al., 2020).South-eastern Australia constitutes a region where the high prevalence of aggressive interactions among the bird species inhabiting this geographical zone makes it a wellsuited scenario for investigating co-occurrence networks.There is ample evidence that Manorina honeyeaters (miners) exert a strong negative influence on bird assemblages of this region due to their despotic behaviour (Mac Nally et al., 2014;Piper & Catterall, 2003;Thomson et al., 2015).Previous studies have revealed that synergistic interactions between miners and other hyper-aggressive species (butcherbirds Cracticus spp.) are likely to play an important role in suppressing threatened and vulnerable -mostly relationship between species trait dissimilarity and species-to-species association's strength.
Main Conclusions: Our results suggest that aggressive species like miners, kookaburras, lorikeets, magpies and butcherbirds can form potential synergies with each other that may intensify the direct negative influence of each species on small-sized songbirds.Since these species thrive in anthropized landscapes and have drastically increased their numbers in some regions, their associative potential should be considered in conservation actions.

K E Y W O R D S
Australia, bird assemblages, competition, HMSC, joint species distribution models, Manorina melanocephala, species associations, woodlands small-bodied bird species in Australia (Westgate et al., 2021).
In turn, these impacts can be buffered or intensified depending on habitat features (e.g.extent of midstorey cover or land-use intensity, Montague-Drake et al., 2011;Mortelliti et al., 2016;Robertson et al., 2013;Westgate et al., 2021).However, little is known about the nature and strength of these associations after controlling for species' environmental preferences, which is often done by disentangling co-occurrence patterns generated by ecological interactions from those generated by co-variation in the species responses to environmental variation.Thus, it is not clear if the existence of positive associations between miners and other large birds (e.g.butcherbirds, ravens) reported by other authors (Val et al., 2018;Westgate et al., 2021) are the result of potential mutual benefits that may arise from improved defence and food supplies (i.e.facilitation effects) or due to similar habitat requirements (preference for edge areas and open and degraded sites).
Similarly, it is not clear if the negative associations detected in previous studies between aggressive species and small-sized birds arise through indirect effects (e.g.shared predators).
In this study, we analysed the relative influence of different assembly processes in shaping local bird communities along a gradient of habitat degradation in south-eastern Australia.We modelled bird occurrences as a function of environmental covariates (climate and habitat variables), species traits and phylogeny using joint species distribution models (Figure 1).Using this approach, we addressed four objectives.First, we investigated how climate and landscape attributes influence species occurrence, and determined how much variation in species occurrence is due to environmental filtering and/or random processes.Second, we examined whether species-to-species associations are non-randomly distributed with respect to the species traits and whether, for instance, functionally more similar species tend to exhibit negative associations after controlling for species' habitat requirements, which could be indicative of competitive exclusion (Graves & Gotelli, 1993).Third, we focused on the noisy miner (Manorina melanocephala), one of the winners in human-altered landscapes and renowned for aggressively excluding other bird species from areas they occupy (e.g.Maron et al., 2013).
We tested for trait differences between those species with which noisy miners are negatively associated and those species with which they exhibit random associations.Lastly, we explored how residual species-to-species associations change along a productivity (stress) gradient.Negative associations could indicate a potential relationship with competition, whose confirmation would require mechanistic experiments (König et al., 2021;Zurell et al., 2018).
Hence, according to the SGH (García-Navas et al., 2021;Maestre et al., 2009), the incidence of negative associations (and thus, the overall degree of competitiveness) should be higher in communities inhabiting more benign (i.e. more productive and less anthropized) environments.

| Study area and bird data
We used bird surveys carried out between December and Febru- (Figure 2).Most of this area (mostly the Moreton Basin) has been cleared for grazing and intensive agriculture, whereas eucalyptus woodlands with tussock grass understorey and warm temperate rainforests constitute other representative habitats in this region (Keith, 2017).
We only kept surveys with count data for all species (surveys with presence/absence data were not used).In addition, we only kept sites with a minimum of three surveys, wherein surveys were carried out in at least two different seasons (some sites were surveyed up to seven times per season).Our dataset accounted for a total of 2476 detections (i.e.observed abundance) for 134 species across 268 surveys in 46 plots.Raptors like kites and eagles were not considered in this study as the point count method is not suited to their large home ranges.

| Phylogenetic data
We computed a Maximum Clade Credibility (MCC) tree from a sample of 1000 trees retrieved from BirdTree (www.birdtree.org;Ericson backbone).This MCC tree was obtained with the maxClad-eCred function of the R package 'phangorn' (Schliep et al., 2022) and pruned to incorporate only those species included in our analyses.
From this consensus tree, we computed a matrix of phylogenetic distances (P) using the cophenetic function of the 'ape' package (Paradis & Schliep, 2019).

| Trait data
For each species, we collated mean trait values for the following morphological variables: body mass (log-transformed), beak size (length) and wing aspect ratio (quantified using the Hand Wing index, HWI).Traits measurements were compiled from the literature and published datasets (Garnett et al., 2015;Tobias, 2021;Tobias et al., 2022).The HWI has been widely used as a proxy of dispersal ability and is computed as 100*D K /L w , where D K is Kipp's distance (the distance between the tip of the first secondary feather and the tip of the longest primary feather) and L w is wing length (Sheard et al., 2020).Since beak size and body mass (size) are highly correlated, we performed a phylogenetic size correction (Revell, 2009), and the residuals (size-corrected beak size) were used as trait values.We also collected information on diet preference (% of insects, fruits, seeds, plants and nectar in diet) and foraging strategies (% of time spent foraging on the ground, understorey, mid-high, canopy and in the air) from the literature F I G U R E 2 Sampling locations (purple dots) in eastern Australia.Bird surveys were conducted at each site following the '2-ha, 20 min' method (see Section 2 for details).The five species involved in the most associations are shown on the right.Pictures by C. Schultz, R. Smith, L. Kee, P. Kerrawn and F. Hort respectively.(Garnett et al., 2015;Wilman et al., 2014).Both subsets (diet preferences and foraging strategies) were reduced to two principal components by means of a principal component analysis (PCA).
The first two axes summarizing diet preferences accounted for 65.2% of the total variance.The first axis (PCd1) was negatively associated with the proportion of invertebrates in diet (−0.70) and positively with the percentage of fruits and plants (0.44 and 0.36 respectively).The second axis (PCd2) was negatively associated with the proportion of nectar in the diet (−0.68), and positively with the seeds proportion (0.56).The first two axes summarizing foraging strategies accounted for 67% of the total variation.The first axis (PCf1) was positively associated with the proportion of foraging on the ground (0.66) and negatively with the percentage of foraging at mid-high (−0.58) and at canopy level (−0.46).The second axis (PCf2) was explained primarily by variation in the proportion of understorey (−0.74) and aerial foraging (0.52).Thus, our final trait matrix (T) comprised seven variables: log-body mass, size-corrected beak size, HWI, PCd1, PCd2, PCf1 and PCf2.From this matrix, we computed pairwise functional dissimilarity using Gower distance.Since there is a tight link between beak attributes and resource acquisition and it has been shown that beak morphology constitutes the trait in which the footprint of interspecific competition is more likely to arise (Chira et al., 2020), we built a second matrix describing species dissimilarity in terms of beak size.that characterize the microclimatic conditions and degree of habitat heterogeneity of each sampling site and conformed our matrix of environmental covariates (X).GPP was used as a proxy for environmental stress since it has been suggested that, for vertebrates, resource gradients can be equated to inverse stress gradients (Barrio et al., 2013;Beaudrot et al., 2020; see also Fugère et al., 2012).

| Environmental variables
Consequently, environments with low GPP are more stressful to vertebrates than highly productive ones, allowing us to test the SGH (originally formulated by plan ecologists) in animal communities (Adams et al., 2022)

| Multi-species abundance model
In order to estimate abundance-based community diversity while accounting explicitly for the imperfect detection of individual birds during the censuses, we used a multi-season community N-mixture model implemented in a Bayesian framework (Royle, 2004;Yamaura et al., 2011).Multispecies abundance models (MSAM) combine and analyse the history of detections and non-detections of all individuals of each species encountered during replicated surveys at a set of sites (Iknayan et al., 2014).Our model follows the overall structure of the MSAM implemented by Montaño-Centellas et al. (2021) wherein the observed abundance y for species i at site j and season t is modelled using a binomial distribution parameterized using the latent true abundance N and the detection probability p.
Latent true abundance is estimated using a Poisson distribution with parameter being the expected abundance.
Detection probability p is defined using a species intercept-only model, such that each species has a constant detection probability estimated across sites, seasons and surveys.
In the season of the first survey f 1 for a given site, the expected abundance is defined using a species and site intercept-only model.
In the following seasons until the season of last survey f 2 for a given site, the expected abundance is defined using an auto-logistic formulation, with the same species-site-specific intercepts as during the first season, in addition to the product of the true abundance N in the previous season and the site-season-specific correlation term ϕ.
Parameters ϕ, and were drawn from community-wide normal distributions with community mean and precision parameters defined using, respectively, normal and gamma hyper-priors.
The analysis of this model was performed in a Bayesian framework using Markov chain Monte Carlo simulation with Gibbs sampling (Plummer, 2003) with the R package 'R2jags' (Su & Masanao, 2015) to call JAGS and export the result into R v.4.0.5.We used the simulation of 3 chains for 150,000 iterations, with the first 75,000 used as burn-in and sampling every 375 iterations.The thinning of MCMC chains was applied due to memory and computation limitations (Link & Eaton, 2012).Thus, we saved 200 iterations per Convergence was assessed for all underlying parameters using the Gelman-Rubin diagnostic (Gelman & Rubin, 1992), with potential scale reduction factors being below 1.1 for 99.97% of underlying parameters.Abundance from the posterior draws was capped at 300 individuals (the estimated value was initially above 300 for 0.16% of posterior draws).Three species (Lopholaimus antarcticus, Apus pacificus and Acanthiza reguloides) were removed all together from the model because their inclusion caused convergence issues for detection probability estimates.From our fitted model, we extracted 20 posterior estimates of N ij (or 'N-matrix', that is, the expected true cross-seasonal abundance of species), which were used as input in subsequent analyses.
The HMSC models were fitted using both abundance and occurrence data from the MSAM.Since for the abundance data, the model did not converge, we only report results yielded by the occupancy models (i.e.presence/absence matrix; Figure S2).

| Hierarchical Modelling of Species Communities
Original data included 133 bird species.For model stability, we focused on relatively common species and selected those species that were observed in a least 10 plots based on the raw data.This resulted in 75 bird species.We then applied Bayesian joint species distribution models (jSDM) within the Hierarchical Modelling of Species Communities (HMSC) framework (Ovaskainen & Abrego, 2020;Ovaskainen et al., 2017).HMSC is a flexible framework in which both species occupancies and abundances can be modelled in response to different environmental factors (climate, habitat).The residual variation in species occurrences is captured by latent variables.Speciesto-species residual association matrices (Ω) are then estimated as variance-covariance matrices of the loadings of the latent variables.Thus, all the variation not explained by the fixed effects (climate, habitat, etc.) is interpreted as being the result of biotic interactions (species-to-species associations) (Ovaskainen et al., 2016).HMSC includes a hierarchical layer which allows us to determine how species responses (Y) to these environmental variables (X) depend on species traits (T) and phylogenetic relationships (P) (Abrego et al., 2017).
In this way, it is possible to ascertain whether, for instance, taxonomically related species have more similar environmental responses than can be expected based on their traits.In the HMSC, the n y × n s response matrix Y consisted of the estimated presence-absence (the output of the MSAM) of the n s = 75 species observed in the n y = 46 sampling plots.We modelled Y with probit-regression, including the environmental variables (i.e.environmental predictors) (matrix X), the trait matrix (T) and the phylogenetic matrix (P) as covariates.
Specifically, we fitted a probit model including the environmental variables as fixed effects, and plot coordinates as a spatially structured random effect.
We fitted the models using the R package 'Hmsc' (Tikhonov et al., 2020) assuming the default prior distributions (Ovaskainen & Abrego, 2020).We sampled the posterior distribution with two MCMC chains, each of which was run for 5000 iterations, out of which the first 1000 were removed as burn-in.The iterations were thinned by 10 to yield 400 posterior samples per chain, and thus 800 posterior samples in total.We ran the model for 20 presenceabsence (O ij ) matrices randomly chosen and used the average and distribution of all model estimates.We pooled the 16,000 posterior samples (800 × 20 models) to obtain a single combined posterior distribution.From this posterior distribution we computed the support values of species' environmental responses ( estimates; i.e. how species in Y responds to covariates in X) and those of the species-to-species associations (Ω matrix) for two probability thresholds, 85 and 95%.The variance-covariance (Ω) matrix contains both direct and indirect associations whereas its inverse (precision matrix, R −1 ) represents the (residual) partial correlation between species while controlling for the effects of the other species (Otso Ovaskainen, personal communication; Poggiato et al., 2021; see Figure S1).
We examined the explanatory power (in-sample validation) of the models through species-specific AUC and Tjur's R 2 values (averaged across species), which measure how well the model discriminates those communities at which the species occurs from those where it does not occur (Pearce & Ferrier, 2000;Tjur, 2009).To compute predictive power (out-of-samples validation), we performed a two-fold cross-validation, in which the sampling units were assigned randomly to two folds, and predictions for each fold were based on model fitted to data on the remaining fold.
Second, we applied Mantel's tests to explore the relationship between species functional dissimilarity (estimated from all traits and only for beak size) and (i) residual species association scores (Ω values; i.e. direct and indirect associations) and (ii) partial correlation scores (R −1 values; i.e.only direct associations).We then focused on noisy miners' associations and assessed whether species-to-species associations are randomly or non-randomly distributed with respect to the species' traits.Lastly, we examined whether communities from a priori more stressful (i.e. less productive) environments are made up of species showing positive associations with other species in comparison with those conforming communities located in more benign environments.To that end, we obtained a precision species-to-species association matrix by averaging (residual partial) pairwise species association scores across our 20 pairwise matrices resulting for each model.From this final matrix, we computed a single value for each species (by averaging across species) indicative of their competitive or facilitative nature ('species association score', SAS).Species values were scaled up to the community level using a community-weighted mean (CWM) approach, in which values of species attributes are averaged across species in the community occurring at the site.We computed community-weighted means for each of the 20 bird abundance matrices.The overall degree of competitiveness (or facilitation) in each site (CWM of SAS averaged across matrices) was regressed against an index of gross primary production (GPP).All analyses were implemented using R 4.0.5.(R Core Team, 2022).

| Co-occurrence analyses
The number of species pairs for which the posterior probability (PP) of having positive or negative co-occurrence (residual covariance) was at least 0.85 was 42 (i.e.3.6% of the total of possible pair combinations; Table 1; Figure 3).About half of tentative species-to-species associations were negative (55%) and involved one of the following five aggressive species: the rainbow lorikeet, the noisy miner (M.melanocephala), the laughing kookaburra, the Australian magpie (Gymnorhina tibicen) and the grey butcherbird.The set of species that were the most sensitive to the presence of aggressive birds were the willie wagtail, the Lewin's honeyeater, the green catbird (Ailuroedus crassirostris) and the wonga pigeon.Thus, a dozen species brought together the bulk of (tentative) negative associations.Most (tentative) positive associations were between two of the five mentioned aggressive species, the strongest one being the association between the noisy miner and the rainbow lorikeet (Ω: 0.85, PP: .93;Table 1; Figure 3).None of the associations were significant for a .95probability threshold.
When controlling for indirect associations, we found that the top-ranked species in terms of positive (direct) associations (Table 2) and the strongest pairwise positive associations (Table 2) were almost the same as those obtained using the variance-covariance matrix (Figure 4).However, for the negative associations, the identities of both the top-ranked species and the strongest pairwise associations were less consistent between groups (indirect and direct effects and only direct effects; Tables 1 and 2).

| Species functional similarity and species associations
There was no positive relationship between species functional dissimilarity and species-to-species association's strength and sign regardless of the suite of traits employed to compute the dissimilarity matrix (all traits or a single one; beak size) and if we controlled for indirect effects or not (Mantel test: all p-values > .50).

| Environmental harshness and competitiveness
We found a marginally significant relationship between the average community-weighted mean of the 'species association score' (CWM of SAS) and an index of productivity (GPP) indicating that positive associations tend to be stronger (or more common) in more stressful, low-productivity environments (Pearson's correlation: ρ = −.26,p = .076;Figure 5).This relationship between CWM of SAS (computed separately for each matrix) and GPP was significant or marginally significant in 12 of the 20 analysed models (see Table S1).

| DISCUSS ION
The transformation and reconfiguration of natural landscapes and native vegetation via agriculture and urbanization can have a profound effect on avian communities (Martin et al., 2012).Beyond the direct consequences of the reduction in suitable habitat, there may be indirect impacts through changes in the competitive interactions between species (Mac Nally et al., 2012).Many generalist species ('winners') like the noisy miner, the rainbow lorikeet and TA B L E 1 Top-ranked negative and positive associations based on residual correlations (indirect and direct effects) and partial residual correlations (only direct effects) along the strength of the association (correlation scores averaged across 20 models) detected in bird communities of eastern Australia.

Species pair Strength (mean [range])
Omega (posterior probability) Note: The omega value (Ω, association score) and the posterior probability obtained after pooling the posterior samples of the 20 HMSC models are also given.
butcherbirds benefit from habitat degradation and thrive in anthropized landscapes.Some of them exert a negative influence on other species which end up being extirpated from the area.In turn, these winner species can interact positively as recent studies have reported (Hingee et al., 2022;Westgate et al., 2021).However, these positive associations may merely reflect shared ecological preferences and/or indirect effects.Here, we adopted a two-fold approach for modelling pairwise species correlations while accommodating imperfect detection.Using this methodology, we provide weak evidence that some aggressive species are positively associated, which may contribute to exacerbating their negative impacts on other species.

| Environmental predictors of species occurrence
The most important environmental predictors, as measured by average median loading size across all 75 species and across models, were landcover type (cleared or woodland), percentage of woody cover and mean temperature.At the species level, the set of species that accounted for the majority of positive associations were species negatively associated with habitat heterogeneity, woody cover and/or GPP.Thus, these species exhibit a high potential to thrive in cleared habitats and anthropized landscapes.This is of utmost relevancy in regions like Australia's eastern coast where agricultural F I G U R E 3 Species-to-species association matrix, estimated by pooling the posterior distribution of 20 HMSC models.Each matrix cell corresponds to a species pair and the gradation of blue-red colours depicts the species-to-species estimated association strength.Blue colours stand for the species that tend to co-occur (tentative co-occurrence; posterior probability >0.85%) more often than expected based on the individual species' occurrence probabilities.Contrary, red colours depict the species that tend not to co-occur.intensification and associated land use changes are leading to extensive biodiversity loss (Bradshaw, 2012;Lindenmayer et al., 2016).

| Species' associations
We found evidence of 'compartmentalism' (sensu Bascompte, 2010), with five medium-sized species more strongly associated with each other than with other species in the assemblage.The noisy miner, the rainbow lorikeet, the laughing kookaburra, the Australian magpie and the grey butcherbird accounted for 79% of the total number of tentative associations (posterior probability .85)detected in our models.The tendency for these swooping birds to be more common in sites with the presence of other taxa belonging to this group seems not to be associated with shared ecological requirements (preference for more-degraded sites) or indirect associations.

F I G U R E 4
Average residual correlation matrix for bird species after controlling for environmental variables.This matrix was obtained by averaging species-to-species associations (Ω) scores obtained for each of the 20 HMSC models.Associations vary between 1 (most positive; blue colour) and −1 (most negative; red colour).Each species is coded by including the first four letters from its genus and species name (see Table S2 for the full scientific names of the species).The five main aggressive species (Manorina melanocephala, Dacelo novaeguineae, Trichoglossus haematodus, Gymnorhina tibicen and Cracticus torquatus) are highlighted with an asterisk.
In line with previous studies, we found positive associations between noisy miners and other aggressive species like butcherbirds and allies while controlling for ecological preferences and the effect of the other species.In this vein, Grey et al. (1998) reported that grey butcherbirds abandoned the habitat from which noisy miners were experimentally removed, suggesting a synergistic relationship between these species.Miners may benefit from the capacity of the butcherbirds and Australian magpies to repulse larger nest-predators (snakes, gliders and small marsupials) since they place their nests in unusually conspicuous locations, and in return, butcherbirds and Australian magpies may benefit from abundant arthropod resources that are not exploited by miners and other advantages in the form of social communication (Fulton, 2008;Hingee et al., 2022).Hence, for instance, Dawson Pell et al. (2018) reported that Australian magpies can gain information on the type or location of predatory threat by eavesdropping on noisy miner alarm calls.The lack of aggression between aggressive species and parrots (rosellas, lorikeets) has also been previously documented (Bain et al., 2020;Hingee et al., 2022;Westgate et al., 2021).Rainbow lorikeets dominate and aggressively protect feeding and nesting resources, thereby excluding other species (Higgins, 1999).
Thus, the agonistic behaviour of noisy miners and rainbow lorikeets is quite similar.In this sense, Smith and Smith (2020) documented a rapid increase of both species to the detriment of others in a Sydney suburb over a 40-year period.Our results suggest the existence of positive associations between species that are not closely related nor exhibit similar morphological characteristics but may rely on complementary feeding or predator avoidance strategies, reinforcing the need to consider facilitative interactions along with competition when seeking to explain species assembly (Seppänen et al., 2007;Sridhar et al., 2012).
We found a smaller strength for negative associations, which accounted for more than half (55%) of the tentative associations.Smallsized (<30 g) insectivorous birds such as swallows and fantails were the species whose occurrence was most affected by the presence of noisy miners and rainbow lorikeets, the two species with a more negative impact on others.When controlling for indirect effects, we found that the strength of the negative associations involving the grey butcherbird or the Australian magpie and small passerines decreased (i.e.there was a weaker negative association).This result agrees with those reported by Westgate et al. (2021), who found that the presence of butcherbirds amplified the negative impacts of the noisy miner on many species, but they did not reduce the occurrence of small woodland birds in isolation.At this point, it must be stressed that in cases of asymmetrical competition (i.e. one species excludes the other one), the capacity of our methodology to infer negative associations is limited.It is due to the fact that the realized niches entirely explain the negative correlation between both species.Thus, no information on biotic interactions is left in the residuals, preventing jSDMs to suggest a competitive interaction from the residual correlation matrix (see Figure 1b in Poggiato et al., 2021).Consequently, we should have observed a larger number of negative residual correlations, mostly between aggressive species (miners, butcherbirds and lorikeets) and small-sized birds that are totally excluded by them ('Narcissus' effect; Colwell & Winkler, 1984).It suggests that even in systems where there is ample evidence that interspecific competition constitutes a key process, it is not straightforward to infer its signature in ecomorphological traits (Elo, Jyrkänkallio-Mikkola, et al., 2021, but see also Rabosky et al., 2007).The low incidence of tentative associations and the lack of significant associations agree with that reported in several studies of pairwise co-occurences, which concluded that non-random (significant) associations are rare in natural ecological communities (e.g.Elo, Ketola, & Komonen, 2021;Gotelli & Ulrich, 2010;Pitta et al., 2012).

| Species functional similarity and species associations
The inability of this approach to detect partly or fully hidden signals of biotic interactions on codistributions (see above) may explain why we did not detect more negative associations between functionally similar species as one would expect if competition drives the species-to-species associations.When the species pool is itself competitively structured it becomes difficult to detect overdispersion in subsets of species drawn from the species pool (Moulton & Lockwood, 1992).It is likely that the set of species with which miners and lorikeets are negatively associated constitute a non-random subset (biased by a 'Narcissus' effect) of the original one in terms of functional traits.We then tested for differences in trait similarity between species with which noisy miners, our focal species, were positively associated and those with which they were negatively associated.We found that noisy miners are more dissimilar in terms of beak morphology with respect to those species with which they are positively associated in comparison with those with which they are more likely to compete (negative associations).This finding agrees with a previous study by Howes et al. (2014), which pointed out that bird assemblages in sites with high noisy miner abundance were dominated by species with longer bills, although they did not control for body size.Our results suggest that this 'reverse keystone species' could act as a biotic filter shaping the ecological profile of bird communities.Also, this study suggests that, in general terms, trait overdispersion as a result of competition constitutes a phenomenon that is not easy to detect in natural communities.This agrees with recent studies that have suggested that, for common bird species, the sign of competition, if present, gives rise to patterns that are not discernible from those of environmental filtering (Elo et al., 2023).

| Environmental harshness and competitiveness
Albeit weak, we observed a relationship between the average level of 'competitiveness' estimated for each community across models and the level of Global Primary Production (GPP) for the plot housing it.A priori this trend supports the stress (productivity)-gradient hypothesis (SHG), which predicts the prevalence of positive interactions in low-productivity environments (Barrio et al., 2013;Beaudrot et al., 2020;García-Navas et al., 2021).However, in our study system, this pattern may reflect the tolerance/preference of butcherbirds, miners and kookaburras (i.e.those species that comprised most of the positive associations) for low-productivity environments since all these species displayed negative coefficients with GPP.This can explain why we observed a clear pattern in the strength and direction of the relationship between species-to-species associations and GPP across models despite the existence of limited variation in environmental conditions among plots.Thus, our results are not conclusive regarding the applicability of the SGH on vertebrates due to the peculiarities of these communities.

| Limitations and perspectives
Although the conditions of our study system (small community size and environmental heterogeneity) are optimal to get the most out of the jSDM, our results must be interpreted prudently.The power of jSDM to recover biotic interactions is limited as previous authors have shown (Zurell et al., 2018).Considering that jSDM are correlative, phenomenological models and thus should not be used to infer causal relationships, which would require mechanistic experiments.These kinds of experiments (e.g.species removal) are difficult to implement at the multispecies level (i.e. on a community scale).In addition, we cannot discard the influence of missing environmental variables, which could lead to spurious associations (e.g.Pollock et al., 2014).As mentioned above, it should also be noted that jSDM only retrieve the realized niche of species, while the fundamental ones remain uncovered.This circumstance hinders the capacity of this approach to infer the existence of competitive interactions among species.Lastly, computational limitations prevented us from using all posterior draws from the multispecies abundance model, which would have allowed us to incorporate the whole uncertainty from the MSAM output in the HMSC analysis.The implementation of models that simultaneously account for imperfect species detection, spatial autocorrelation in density for multiple species, as well as changes in spatial distributions over time (joint dynamic species distribution models; jDSDM), may contribute to overcoming some of the mentioned biases (Blanchet et al., 2020;Thorson et al., 2016).This methodology in combination with species trait information can provide ecological insights into species-to-species associations and contribute to a better understanding of species assembly rules.

| CON CLUS IONS AND CONS ERVATION IMPLIC ATIONS
As far as we know just a few studies have confronted the challenge of simultaneously modelling imperfect detection and species correlations (Tobler et al., 2019).Here, we combine a multispecies abundance model and joint species distribution models within the HMSC framework to obtain robust inferences on species co-distributions.
We found weak evidence of segregated/aggregated pairs in line with previous studies (see Korňan & Svitok, 2018 and references therein).
Only considering (tentative) non-random associations, our results suggest that aggressive species like miners, kookaburras, rainbow lorikeets and butcherbirds tend to be negatively associated with other species (mainly small songbirds) and, in turn, are positively associated with each other after controlling for habitat preferences and indirect effects.Some of these native species have drastically increased their numbers in some regions of mainland Australia and Tasmania being catalogued as a pest or invasive species (Low, 2009;Threatened Species Scientific Committee, 2013).Their capacity to form strong synergies should be considered by local managers to increase the effectiveness of conservation actions such as removal experiments (e.g.noisy miner culling; Beggs et al., 2020;Melton et al., 2021) and restoration/reintroduction programmes (e.g.Bennett et al., 2012).
ary across 17 seasons (from season 1998-1999 to season 2000-2001, and from season 2006-2007 to season 2019-2020) using the '2-ha, 20 min' method of the Australian Atlas (https://birda ta.birdl ife.org.au).This census method involves searching for birds in a two-hectare area for 20 min.During this period, the volunteer records the geographical reference at the centre of their site and all birds seen or heard in their survey area, including those flying overhead.The data are then validated by experts, ensuring that the geographical references are sensible and that the species are within their known range.Our original dataset included a total of 52 sites across the Clarence-Moreton bioregion, which covers large parts of the Great Dividing Range of south-east Queensland F I G U R E 1 Graphical summary of the workflow used to examine bird species co-occurrence patterns.This scheme includes the four main objectives of the present study (see further details in Figure S1).
Microclimate and landscape structure are major driving forces that influence avian community composition and abundance.We selected six environmental factors (annual precipitation [mm], mean annual temperature [°C], gross primary production [GPP], landcover type [binary variable: cleared habitat/woodland], woody vegetation cover [%], average vegetation height [dm] and habitat heterogeneity [estimated using Simpson diversity index, see Tuanmu & Jetz, 2015]) . Environmental variables were retrieved from global databases (annual precipitation and temperature: world clim.org;Fick & Hijmans, 2017; habitat heterogeneity: https:// www.earth env.org/texture) and different national repositories like the Terrestrial Ecosystem Research Network (TERN) and the Centre for Water and Landscape Dynamics of the Australian National University (woody vegetation cover: Gill, 2012; vegetation height: Scarth, 2013; GPP: Yebra et al., 2015).All environmental variables were summarized at a 2-ha level of resolution to match the grain from bird surveys.
The five main aggressive species (Manorina melanocephala, Dacelo novaeguineae, Trichoglossus haematodus, Gymnorhina tibicen and Cracticus torquatus) are highlighted with an asterisk.Top-ranked species in terms of negative and positive associations considering both indirect and direct effects (left) and only direct effects (right).For each species, the mean value resulting from averaging pairwise partial correlation scores is shown.Total average values obtained by averaging the mean value across the total number of species (n = 75) are given at the bottom.
Relationship between the average communityweighted mean of the 'species association score' (CWM of SAS) and an index of global primary production (GPP).The 'species association score' represents the level of competitiveness of each species in relation to the remaining ones.The average CWM was computed across 20 models (one for each abundance matrix).Each data point represents a single community (plot).Blue dots denote plots with the presence of noisy miners and green dots denote those plots in which this species was not registered.