Unravelling species co‐occurrence in a steppe bird community of Inner Mongolia: Insights for the conservation of the endangered Jankowski’s Bunting

To evaluate the patterns of bird assemblage and distribution in an endangered grassland system, taking into accounts both environmental and biotic effects. To further focus on an endangered songbird and associated steppe birds.


| INTRODUC TI ON
Understanding ecological processes that drive species distribution and associations is a central question in ecology, evolution and conservation biology, especially for evaluating spatial protections for threatened species, tracking the spread of invasive species and interpreting the impacts of environmental change across species and communities (Thorson et al., 2016). Species distributions and co-occurrences result from environmental requirements and filtering, as well as biotic interactions such as competition, predation, parasitism and mutualism with other species (Tulloch, Chadès, & Lindenmayer, 2018). As an illustration, within bird communities, habitat specialist species co-occur more often with other species sharing a similar specialization degree (Julliard, Clavel, Devictor, Jiguet, & Couvet, 2006). However, strong competitive interactions decrease the rate of species co-occurrence, so of local species richness, while shared habitat preferences and the lack of strong antagonistic interactions lead to positive co-occurrence patterns. In the past decades, habitat loss, climate change and other forms of human-caused mortality (Calvert et al., 2013;Loss, Will, & Marra, 2012) have contributed to the global loss in abundance of animal populations and changes in the number and magnitude of species associations (Rosenberg et al., 2019).
Consequently, it appears essential to quantify species associations and their responses to environmental change in real communities to understand how to preserve community structure and ecosystems services.
However, these methods usually do not disentangle co-occurrence patterns generated by ecological interactions with those generated by covariation in species-specific responses to abiotic factors . Other methods incorporate species as predictors into the classical framework of species distribution models (SDM), allowing assessing whether and how biotic interactions affect species co-occurrence (e.g. Barbet-Massin & Jiguet, 2011;Sanín & Anderson, 2018). This approach can be inefficient for the estimation of too-numerous species associations in rich communities  and thus has been typically applied for dominant species only (le Roux, Pellissier, Wisz, & Luoto, 2014). Recently developed joint species distribution models (JSDMs) are an effective tool to capture the effects of both biotic and abiotic interactions in communities (Dormann et al., 2018;Hui, 2016;Inoue, Stoeckl, & Geist, 2017;Ovaskainen et al., 2017;Pollock et al., 2014;Tulloch et al., 2018;Warton et al., 2015). JSDM usually fits a multivariate hierarchical model that uses generalized linear models, ordination techniques and latent variables to simultaneously explore interactions across taxa and the response of abundance to environmental variables (Warton et al., 2015). Typically, the model parameters are estimated using maximum likelihood or Bayesian methods implemented by Laplace Approximation and Markov Chain Monte Carlo (MCMC) simulations, respectively. The JSDMs have been widely applied to explore fish (Inoue et al., 2017), bird (Royan et al., 2016), butterfly (Ovaskainen, Roy, Fox, & Anderson, 2016b) and fungus communities (Abrego, Dunson, Halme, Salcedo, & Ovaskainen, 2017) in order to model the joint distribution or abundance of species in communities and further estimate the response of each co-occurring species to abiotic environmental factors.
Within birds, those associated with grassland habitats are highly threatened (Traba & Morales, 2019). Across breeding biomes, grassland birds have the largest magnitude of total population loss since 1970 in North America (Rosenberg et al., 2019) and Europe (Jiguet et al., 2010). Despite the impacts of environmental covariates (e.g. land use or climate) on bird communities have been well studied, few studies have evaluated the effect of biotic interactions on species distribution and abundance. However, it is known that bird species interact with each other directly and indirectly in a number of mutualistic and antagonistic ways, including competition for food and space, facilitation, predation and parasitism. Consequently, a population loss in any species (rare or common, widespread or endemic) may have disproportionate influence on the community structure, modifying components of food webs and ecosystem functions within communities (Ehrlén & Morris, 2015).
In this study, we applied JSDMs to evaluate the patterns of species assemblage and distribution in an endangered grassland system, taking into accounts both environmental and biotic effects.
We focused on the grassland bird community of steppes in Inner Mongolia, where anthropogenic activities destroyed and fragmented this habitat, because of crop development, overgrazing, coal mining and urbanization (John, Chen, Lu, & Wilske, 2009;Tao et al., 2015). We were particularly interested in focusing on one songbird within this community, the Jankowski's Bunting Emberiza jankowskii. This passerine is an endangered species that has a restricted breeding range and suffers a steep population decline since decades (BirdLife International, 2019a;Gao, 2002), with a relict population size currently estimated at c.10,000 breeding individuals (Han et al., 2018). It shares its range with the closely related Meadow Bunting Emberiza cioides, which has a large distribution range, is a common breeding species in Asia with up to 400,000 pairs and presumably stable populations (BirdLife International, 2019b). Since the two buntings share a similar morphology, key aspects of life history and resource use in the study area, we were particularly interested in learning how their co-occurrences are linked to the sharing of similar niche requirements and potential competition. Such a conflict could lead to a potential progressive replacement of the highly specialized

K E Y W O R D S
Emberiza cioides, Emberiza jankowskii, joint species distribution models, niche sharing, species interactions and localized Jankowski's by the more generalist and widespread Meadow Bunting, fitting classical mechanisms operating in biotic homogenization (Le Viol et al., 2012). We further expected potential interactions between other species present in the community, such as potential predators (corvids, raptors) and a potential parasite (cuckoo), while other ground-nesting steppe birds should co-occur with the buntings because they might share similar ecological requirements (wheatears, larks). Unlike Jankowski's Bunting, the other steppe or farmland species have a favourable conservation status or are not considered as priority conservation targets. As conservation efforts focusing on a single species may be rather cost-ineffective within a broader conservation strategy, we also wanted to evaluate whether Jankowski's Bunting could represent a good indicator species as a surrogate for overall steppe bird conservation. Therefore, we addressed three specific objectives in this study. First, we wanted to quantify the species-specific responses to environmental variables, including climate, landscape and habitat predictors. Second, we wanted to examine patterns of significant positive or negative co-occurrence in this bird community and interpret some in the context of niche sharing and species interactions. Third, we wanted to test whether we could use the presence of Jankowski's Bunting as a surrogate of the abundance of other steppe birds for global conservation purpose.

| Study area and bird abundance data
The study area is located in the steppe zone of northeast Inner Mongolia ( Figure 1). The area is dominated by natural grassland, but cropland, small lakes, ponds and small villages are also characteristics of the landscape. Average annual temperature is 6-7°C, and annual precipitation averages 300-400 mm, with most of the rainfall occurring from June to September. We estimated bird species abundances by realizing point counts, a fixed observer recording all individual birds seen or heard within a 50 m radius during 10 min.
The sampling was conducted under good weather conditions during May and June 2018, during the peak of the bird breeding season. A total of 560 count points was surveyed in different landscape types, with a minimum distance of 500 m between two points. We restricted our analyses to those 30 species known to use grassland or farmland as their main breeding habitat and excluded seven species with low occurrence, detected at less than five sites (Chukar

| Defining independent environmental predictors
We ran a principal component analysis to reduce the original set of potentially correlated habitat and bioclimatic variables to a twodimensional space, respectively. The first two principal components of habitat variables (habitat PC1 and habitat PC2) captured 30% and 26% of the total variability, respectively. Habitat PC1 was negatively related to edge density, the number of land cover patches and Simpson diversity, indicating the gradient of a rather intact habitat at the landscape scale. Habitat PC2 corresponded to a gradient in local habitat parameters with an increase in bare ground percentage and a decrease in plant canopy and vegetation height (see Appendix Table S1 and Figure S1a). Climate PC1 and climate PC2 captured 59% and 24% of the total variability in climate variables, respectively. Climate PC1 was positively related to winter precipitation (Precipitation of Driest and Coldest Quarter), while negatively related to winter temperature (mean temperature of driest and coldest quarter). Climate PC2 was mainly negatively related to summer precipitation (precipitation of wettest and warmest quarter; see Appendix Table S2 and Figure S1b). The PCA was computed by "prcomp" function built in R.3.6.1 (R Core Team, 2019), and the visualization of its outputs was made by "factoextra" package (Kassambara & Mundt, 2017).

| Statistical analysis
We used the BORAL package (Hui, 2016) to investigate how different abiotic and residual components relate to the processes of bird community assembly. Raw bird data in the models are the matrix of the detected abundances across the 560 count points for the 23 bird species. A key feature of the BORAL package is the ability to fit species distribution to environmental variables using generalized linear models and incorporate latent variables as a parsimonious method to account for residual correlation among species. The latent variables here can be interpreted as a device to account for any residual covariation not explained by the measured factors, such as missing predictors or species interaction (Warton et al., 2015). If these "missing parts" are informative for the species co-occurrence, then they will induce a residual correlation between the species (Ovaskainen, Abrego, Halme, & Dunson, 2016a).
Following the method of Hui (Hui, 2016), we first fitted a pure latent variable model (LVM), in which species are regressed against a set of unknown covariates, leading to an unconstrained ordination biplot for visualizing site and species patterns (Appendix Figure S2a).
To further study how well the bird assemblage as a whole (both species abundance and composition) is explained by the environment, we fitted correlated response model (CRM), by including abiotic covariates, to separate the correlations between species due to environmental response and the residual correlation due to other sources of covariation. We used the default priors for MCMC parameters and excluded random site effects in our model process.
To obtain an overall measure of the model's predictive power (i.e. how well the predictor variables describe the species assemblage), we calculated a proportional difference in the trace of the estimated residual covariate matrix between the LVM and CRM (Warton et al., 2015). Among the outputs of correlated response model (CRM), we identified the relative importance of abiotic covariates for each species and constructed a horizontal line plot showing 95% highest posterior density (HPD) intervals for the column-specific regression coefficients by "ggboral" package. We further plotted the correlations between species due to environmental response and residual covariates, respectively. We also constructed another ordination plot analogous to LVM, as an alternative way to visualize the residual correlation between species (see Appendix Figure S2b).

| RE SULTS
The

| D ISCUSS I ON
Fine-scale species distribution patterns result from processes linked to environmental filtering and biotic interactions (Sanín & Anderson, 2018 Jankowski's Bunting. All responded positively to habitat PC1 and negatively to habitat PC2, indicating that they are more abundant in F I G U R E 3 Plots of the correlations between species due to the abiotic response (a) and residual correlations (b). Significant associations are displayed where species co-occurred more (less) frequently than by chance, with an alpha threshold of 0.05. The sign of the correlations is represented by colours (red and blue for negative and positive correlations, respectively), while the strength of correlations is represented by the size of the circles

(b) Correlations due to residual covariates
intact natural grassland while they avoid fragmented habitats with too much bare ground. In Inner Mongolia, many studies showed that the dominant grassland habitat was getting more fragmented, characterized by an increasing number of patches at the regional and biome scales (e.g. John et al., 2009). Furthermore, the increase in portions of barren cover caused by overgrazing suggests that the landscape in this arid and semi-arid ecosystem is becoming more homogeneous and water stressed (Wu & Ci, 2002). Such habitat changes should impact the occurrence and abundance of those sensible steppe songbirds.
Besides, we found that Jankowski's and Meadow Buntings share similar responses to all abiotic variables, demonstrating that the two species have a large niche overlap, and revealed a strong positive co-occurrence pattern, while the opposite would be expected in case of a strong output of competition for similar habitats leading to a local mutual exclusion. Indeed, no two species with the same niche should stably coexist due to competitive exclusion (HilleRisLambers, Adler, Harpole, Levine, & Mayfield, 2012). To mediate the competition, Meadow Bunting has developed spatial and temporal strategies to enlarge niche differences. For example, it usually inhabits more diverse habitats and breeds earlier in the season, sometimes nesting in small trees hence reducing nest destruction from grazing. As a consequence, Meadow Bunting appears as more adaptive and competitive for using limited resources than the sympatric Jankowski's Bunting. Changes in environmental conditions may accelerate the displacement from weaker competitors to stronger competitors, which may render Jankowski's Bunting more vulnerable due to its more restricted habitat requirements and potentially lower ability to exploit new opportunities. To disentangle the co-occurrence of the two bunting species, we adopted here a multiple species approach, and not the regular pairwise approach (Veech, 2014), giving a more general overview of species interactions beyond the two sibling interacting species. The two studied buntings did not appear as strong competitors here despite sharing very similar niche requirements.
The joint species distribution modelling revealed the degree of displacement (negative links) or coexistence (positive links) between species in the studied steppe bird community. Correlations due to abiotic covariates identified species with similar habitat requirements (positive interactions) and species with contrasting environmental requirements (negative interactions). The Jankowski's Bunting displays several positive co-occurrences with other steppe birds, including Meadow Bunting, Skylark and Pied Wheatear. All being open steppe ground-nesting species, we suggest that they likely respond in an analogous way to habitat variable, so should be sensitive to environmental changes in a similar way, as they share multiple dimensions of their environmental niche. Looking at the contribution of environmental variables to the prediction of their abundances (Figure 2), they share a negative response with climate PC1 (more precipitation and lower temperature in winter) and habitat PC2 (more bare ground and less plant canopy), indicating that a warmer and drier winter could benefit to the abundance of these species. This also informed us on how many, and which, potential species might benefit from recovery actions (Borthagaray, Arim, & Marquet, 2014). For example, the four ground-nesting birds that co-occurred were all likely to increase when livestock grazing was reduced, probably because of a lower risk of nest destruction by trampling (Bas, Renard, & Jiguet, 2009), but also as a more complex vegetation cover is generally associated with greater richness of small birds, particularly nectarivores and insectivores (Val, Eldridge, Travers, & Oliver, 2018). On the other hand, the negative correlations due to abiotic variables of species co-occurrences are less straightforward to explain, but might be due to non-overlapping niches, such as Barn Swallow and Common Cuckoo, or Isabelline Wheatear and Collared Dove.
Within the species interactions potentially captured by the JSDM, we considered a parasite (common cuckoo) and many possible host species (open-nesting passerines such as larks, wheatears, buntings). As a widespread species, the common cuckoo is occurring in different habitat types, such as natural forests, secondary growth forests, open wooded areas, steppe, scrub, meadows and reed beds.
Furthermore, as a main brood parasite, it has numerous host species including many insectivorous songbirds such as flycatchers, chats, warblers, pipits, wagtails and buntings. The positive correlations between cuckoo and five open-nest steppe songbirds here probably illustrate that the distribution of parasitic birds is not only driven by their own ecological requirements (climate and trophic availability), but also largely by the presence of their host species (Ducatez, 2014;Lee et al., 2014). Our results are also consistent with former research showing that cuckoo local occurrence is indicative of high overall bird species richness (Morelli et al., 2017). Here, we found positive correlations for the co-occurrence of common cuckoo and three of its identified hosts in Inner Mongolian steppes (the two buntings and the skylark).
The latent variable approach used here helped to isolate environmental-driven co-occurrence patterns from species associations. Positive co-occurrence patterns often tend to reflect mutualistic interactions such as facilitation or complementarity (Morales-Castilla, Matias, Gravel, & Araujo, 2015;Ottosson et al., 2014;Ovaskainen & Soininen, 2011). For instance, Eurasian Magpie exhibited a significant positive interaction with Tree Sparrow, suggesting that they may benefit from more efficient foraging in mixed-species feeding flocks (Fitzgibbon, 1990). Furthermore, after accounting for abiotic variables, the very common Skylark and two other songbirds (Chinese Greenfinch Chloris sinica and Mongolian Lark) exhibited strong negative co-occurrence patterns with Amur Falcon, a major passerine predator in Inner Mongolia.
A logical explanation is that prey species could adopt fine-scale spatial-temporal anti-predator behaviours and try to avoid co-occurring with predator (Cunningham, Scoleri, Johnson, Barmuta, & Jones, 2019). However, not all negative co-occurrences indicate negative biotic interactions, they could also be driven by missing environmental drivers, migration history or phylogenetic history (Dormann et al., 2018), especially when species have a special habitat preference and occupy a unique niche space. These could be an alternate explanation for the negative association pattern observed with Jankowski's Bunting or Mongolian Lark and other birds, especially for species usually considered as living close to human settlements, such as corvids (crows and magpies), Barn Swallow or Collared Dove.
Global changes had significant impacts not just on individual species but on the composition of ecological communities (Sax & Gaines, 2003). By using JSDM, we were able to confirm that bird assemblage depends not only on their individual response to the abiotic environment, but also on species interactions, with a predominant contribution of the later. Previous studies challenged the importance of pairwise species interactions for contributing to broad-scale distribution patterns, for example reporting an overriding influence of broad-scale climate on co-occurrence patterns (Copenhaver-Parry & Bell, 2018), while the importance and direction of interactions also depend on environmental stress (He, Bertness, & Altieri, 2013). By considering an overall community beyond the classical species pairwise approach (Veech, 2014), our results thus contributed to the growing evidence that species pools are filtered by environmental and biotic variables at different scales . Given that species exhibit such varied responses to environmental variables, we may identify species groups of relevance to habitat management and restoration. In our study, several species predominantly occur in natural grassland of high-quality and thus could benefit from the implementation of the National Key Public Forest Protection Project in China, which fences protected areas and excludes these habitats from grazing and other agricultural activities (China Forestry Sustainable Development Strategy Research Group, 2002). Other species are primarily found in fragmented and degraded grasslands and may thus be considered as disturbance specialists of little concern for conservation. Furthermore, in conservation biology, indicator, umbrella and keystone species are frequently used as surrogates for the presence of other species, of overall community or of specific environmental conditions (Caro, 2010;Caro & O'Doherty, 1999). However, the surrogate species approach has been criticized and might produce expensive mistakes for conservation planning (Andelman & Fagan, 2000). One solution is to consider surrogate groups (Wiens, Hayward, Holthausen, & Wisdom, 2008).
Another approach is to develop niche models to identify such surrogates (Nekaris, Arnell, & Svensson, 2015), while the development of further JSDM to identify surrogate species or communities might contribute to a solution. Since the current population of Jankowski's Bunting is mainly confined to Inner Mongolia and only two of the known sites have any form of official protection (Han et al., 2018), projects that focus on the conservation and restoration of its remnant grassland habitat are necessary. Besides, Jankowski's Bunting displayed strong positive co-occurrences with other ground-nesting birds and exhibited significant responses to all measured habitat and climate variables, indicating that this endangered bird has a high niche specialization and wide association with other sympatric steppe bird species. According to the coexistence theory, species that positively co-occurred with a declining or lost species in a threatened landscape are more likely to decline or go extinct than other species, due to shared environmental requirements and threats (HilleRisLambers et al., 2012). Consequently, we would argue that local losses or declines of the vulnerable Jankowski's Bunting can signal negative outcomes for the entire community, and conservation management on Jankowski's Bunting, such as reducing grazing intensity or setting up new protected areas (Jiang et al., 2008), should not only be important to this single engendered species but also beneficial for other co-occurring open-nest grassland birds.
In conclusion, JSDM allowed the combination of traditional species-level knowledge of responses to abiotic variables with co-occurrence analyses, permitted the systematic identification of indicator species and the examination of community structures. By thinking not only about individual species but about how they share space and resources with others, we can imply management actions that benefit the most vulnerable species and avoid actions that might lead to unexpected declines (Tulloch et al., 2018). Our findings are thus helpful for understanding how abiotic and biotic processes interactively alter bird communities and making effective management decisions to mitigate multiple threats in Inner Mongolian grasslands.
This flexible approach could also be adapted to suit conservation applications elsewhere. Retrieved from http://www.birdl ife.org