Broad‐scale patterns of geographic avoidance between species emerge in the absence of fine‐scale mechanisms of coexistence

The need to forecast range shifts under future climate change has motivated an increasing interest in better understanding the role of biotic interactions in driving diversity patterns. The contribution of biotic interactions to shaping broad‐scale species distributions is, however, still debated, partly due to the difficulty of detecting their effects. We aim to test whether spatial exclusion between potentially competing species can be detected at the species range scale, and whether this pattern relates to fine‐scale mechanisms of coexistence.


| INTRODUC TI ON
The need to forecast shifts in species distributions under global climate change is driving an emerging interest in understanding the factors that shape species ranges (Pacifici et al., 2015). The presence of a species in a given location, and thus the species' range, depends on the abiotic environment (climate, topography and physical environment), biotic interactions and movement factors that relate to species dispersal ability under constraints of its evolutionary history (Poggiato et al., 2021;Soberón & Peterson, 2005). The relative contribution of these factors across spatial scales is, however, still not well understood. Climatic factors are commonly thought to shape the distribution of species at a broad spatial scale, whereas the impact of biotic factors is thought to be more pronounced at the local scale (Willis & Whittaker, 2002). Hence, the role of biotic factors in shaping broad-scale patterns of species ranges remains controversial (Early & Keith, 2019;Wiens, 2011;Wisz et al., 2013), for example, see the Eltonian Noise hypothesis (Soberón & Nakamura, 2009).
Among the different types of biotic interactions, interspecific competition can result in competitive spatial exclusion (Gause, 1932).
While the effect of competition on local-scale patterns such as ensemble composition is broadly accepted and supported by studies (e.g. Fraterrigo et al., 2014), substantial evidence shows that competition effects on species presence can scale up and drive broad-scale assemblage patterns, as seen in avifauna in Denmark (e.g. Gotelli et al., 2010), and also can shape species broad-scale range limits.
There are several documented examples of range patterns being shaped partially by competition (reviewed in Wisz et al., 2013). In their review of the literature, Sexton et al. (2009) found that 39 out of 51 studies supported or partially supported the role of competition in shaping broad-scale species range limits. Our mechanistic understanding of the role of interspecific competition in shaping species ranges is, however, limited, in part because we still lack tools to identify such patterns. Yet recent theoretical studies suggest that the scaling-up of competitive exclusion depends on the development of fine-scale coexistence mechanisms (Godsoe et al., 2015).
The ultimate consequence of interspecific competition in the absence of coexistence mechanisms, such as partitioning of the trophic, spatial or temporal ecological niche, is spatial exclusion (Gause principle: Gause, 1932;Hardin, 1960). Consequently, if competitive exclusion scales up and has an effect on species' geographic distributions, the predicted detectable signal on species' occurrences would be a tendency to be absent from their environmentally suitable area where the competitor is present. Indeed, parapatric ranges between morphologically similar or phylogenetically related species were traditionally interpreted as the result of competitive interactions, especially in the absence of geographical barriers to dispersal and when sharp edges do not match clear environmental gradients (Bull, 1991).
The main difficulty when aiming to infer competition effects using patterns of species presences is separating true avoidance resulting from competition from the effects of differential ecological preferences (Bar-Massada, 2015).
One of the current methodological approaches that consider biotic interactions at macroecological scales is including the geographic range of a potential competitor as an additional predictor layer in species distribution models (SDMs) (Anderson, 2017).
Although this approach often improves model performance (e.g. Palacio & Girini, 2018), it cannot be used to separate biotic from environmental effects (Dormann et al., 2018). More recent methodological developments aim to separate biotic from environmental effects at the local community level by using patterns of species co-occurrences (e.g. D' Amen et al., 2018). Among these approaches, joint species distribution models (JSDMs) (e.g. Harris, 2015;Pollock et al., 2014) model simultaneously species presences using environmental variables and identify patterns of residual co-occurrences between species that are not explained by environmental predictors and therefore might reveal a signal of biotic interactions. However, the validity and the interpretation of co-occurrence patterns as inferences on biotic interactions are unclear (see Blanchet et al., 2020;Dormann et al., 2018;Peterson et al., 2020;Poggiato et al., 2021;Zurell et al., 2018). Moreover, the assemblage data that these methods require reduce their applicability for understanding broad-scale species range patterns.
Despite the paucity of available methods to identify broad-scale effects of biotic interactions on species ranges, there is increasing awareness that considering the role of biotic interaction is important for better understanding how climate change will impact diversity patterns (Alexander et al., 2016). Biotic interactions can modify responses towards climate change, for example, through preventing species from being able to maintain or establish populations in predicted future suitable range due to increased overlap with competitors either in existing range or in new suitable areas (HilleRisLambers et al., 2013).
We aim to test whether spatial exclusion between potentially competing species can be detected at the species' range scale, and whether fine-scale mechanisms of coexistence affect broad-scale species distributions. To this end, we develop a measure of geographic avoidance between pairs of species that uses SDM outputs.
We apply the measure to four sets of cryptic Palearctic bat species (10 bat species) that show different degrees of ecological similarity and range overlap. Cryptic species are morphologically similar but genetically distinct species and therefore are likely to show high ecological similarity and compete for the same resources, offering an excellent model system for testing competitive interactions (Novella-Fernandez et al., 2020). We predict that broad-scale geographic avoidance can be detected in pairs of species with higher ecological similarity and in the absence of fine-scale coexistence mechanisms. We expect that species with high geographic avoidance will be less likely to colonize their future suitable range in areas of overlap with their competitors, which will reduce the size of the suitable range available to them under climate change. Therefore, it is important to understand current patterns of geographic avoidance to be able to incorporate the effects of future range overlap with competitors in climate change vulnerability assessments.

| Developing a measure of geographic avoidance
While SDMs model species' suitable geographic areas, this predicted environmentally suitable range (herein predicted range) is not entirely occupied by the species, which are restricted to a smaller portion of it (herein realised range). The difference between modelled species predicted and realised ranges can be attributed partially to processes not included in the correlative environmental SDM framework, such as biotic interactions or dispersal limitation (Guisan et al., 2017). A species whose range is negatively affected by the presence of another is expected to realise less of its predicted suitable range in those areas where the competitor is present. Based on this, we present a measure of species geographic avoidance that uses SDM outputs for pairs of potential competitor species i and j with partially overlapping ranges. The measure establishes a relationship between the realised proportion of the range of species i and j in overlapping (sympatric) and non-overlapping (allopatric) areas. To convert this to a measure of avoidance, we take 1 − the ratio between these: See Figure 1 for an overview of the implementation of this methodology.

• P i (predicted range): (P i ):
The binary geographic area that is environmentally suitable for species i. Binary suitable areas are generated from SDM outputs using the thresholding method that minimizes the difference between sensitivity and specificity. This thresholding method has been shown to outperform other approaches (Liu et al., 2013).
• PO ij (predicted range overlap): The geographic overlap between the predicted ranges of species i and j. This is the expected geographic overlap between species according to their binary model outputs (Gutiérrez et al., 2014).
• R i (realised range): The area of the predicted range where species i is present. A low proportion of predicted range realised suggests that processes not included in the abiotic model, such as the effects of biotic interactions and/or dispersal limitations (Soberón & Peterson, 2005), are preventing the species from occupying its entire environmentally suitable range. The realised range is estimated by clipping the predicted range of the species by its known range limits (e.g. Marcer et al., 2013). This can be based on detailed occurrence records, when such information is available, or expert-generated distribution maps, such as those available from the IUCN Red List (https://www.iucnr edlist.org). The reliability of the measure depends on accurately measuring species realised ranges, and therefore, it should not be applied to data-deficient species or cryptic species that have not been recently assessed if using IUCN red list range maps.
• RO ij (realised range overlap): The geographic overlap between the realised ranges of species i and j. It represents the overlap between the two species that occurs at the model resolution (cell size). High values denote species coexistence.
Inserting these parameters into Equation 1 and rearranging gets: Values of GA ij = 0 indicate that for a given pair of species i and j, an equal proportion of predicted range is realised in overlapping areas and in non-overlapping areas. Values of GA ij < 0 indicate that a higher proportion of predicted range is realised in overlapping areas.
Finally, values of GA ij > 0 indicate a lower proportion of predicted range is realised in overlapping areas, as expected if they occupy less of their suitable range in the presence of a competitor. The maximum value of GA ij will tend to 1 in the extreme case where species i and j do not realise any of their predicted range in their overlapping areas.
Because the measure is based on ratios, values based on very small proportion of realised range overlap (around less than 0.5%) are prone to high variability and are not reliable, and therefore, the measure is not suitable under these circumstances. Code to implement this measure in R is provided in dryad: https://doi.org/10.5061/ dryad.rbnzs 7hbq.

| Study system
Bats offer good case studies for assessing the role of interspecific competition in shaping species ranges to due their high number of cryptic species (e.g. Ibáñez et al., 2006). To test the performance of the measure with an empirical dataset, we used four sets of potentially competing cryptic bat species ( proportion of range realised in overlapping areas between i and j proportion of range realised in non − overlapping areas between i and j kolombatovici, Rhinolophus euryale and Rhinolophus mehelyi. All species within each genus share a very similar morphology but have different degrees of geographic and ecological overlap. Consequently, we expect different levels of potential competition among them.
See Table 1 for an outline of geographic, morphological and ecological overlap among these cryptic bat species and Tables S1.1-S1.3 for a more detailed overview of species ecological similarity. In ad- and 4,540 were unpublished records provided by seven bat researchers. Only records with validated species identification were used. In areas where ranges of similar species overlap, we used records that were confirmed genetically or morphologically by bat experts, considering identification year relative to when the species were split taxonomically. When original published records consisted only of a locality descriptor, geographic coordinates were obtained manually whenever possible (<1% of records). Spatial quality of records was checked, removing low-quality records in terms of spatial resolution and confirmed identification. We accounted for uneven sampling across the study area by removing clustered records in intensively F I G U R E 1 Flow chart representing the implementation of the approach to measure geographic patterns of species avoidance expected from macroecological effects of competition. 1-For a given pair of potentially competitor species, occurrence records and species environmental drivers are used to model species' predicted ranges. 2-Binary predicted ranges are cropped by species realised range limits to obtain species' realised ranges. 3-Ranges of both species are combined to estimate extent of overlap in predicted and realised ranges. 4-The measure of Geographic Avoidance (GA ij ) is applied. In this example, arrows describe a larger decrease from predicted to realised range in areas that overlap with the competitor (large arrows) than in non-overlapping areas (small arrows). Or, equivalently, a lower proportion of realised range occurs in overlapping areas compared to non-overlapping areas. This pattern is described by GA ij with a value >0 [Colour figure can be viewed at wileyonlinelibrary.com] sampled areas (Kramer-Schadt et al., 2013), using the ArcGIS toolbox "SDMtools" (Brown, 2014) with unclustering distances between 20 and 60 km based on species' degree of cluttering and range size. Most records used in the models (>95%) had a spatial precision <2 km. We included records with lower precision (up to 10 km) only in regions with very low sampling intensity (Eastern Europe and North Africa). This accounted for <5% of records. were selected based on Akaike information criterion (AIC) scores using enmtools v1.3 (Table S3.6 for final model features). To assess model performance, we used tenfold cross-validations, with 75% of records used for training and 25% for model testing. The ability of the models to discriminate between presence locations and background pseudo-absences was evaluated with area under the curve (AUC) of the receiver operator characteristics (ROC) and True Skill Statistic (TSS). The 10 model replicates were combined to obtain a final predictive map for each of the five modelling methods.
Ensemble models were obtained by using AUC values to proportionally weight each method according to its predictive power, excluding models with AUC < 0.75. Binary presence-absence maps were generated based on the thresholding criteria that minimizes the difference between sensitivity and specificity because these criteria has been shown to outperform others (Liu et al., 2013). However, we analysed whether using other four commonly used thresholding criteria resulted in changes in the measure. It is important to note that measures of accuracy of thresholded maps in general assume presence-absence data, while our models were generated using presence-only data and background points.
Species' realised ranges were calculated by clipping the SDM output binary maps using a convex hull polygon drawn around occurrence records (r package concaveman, Gombin et al., 2017) and adding a buffer of 30 km. This distance considers a species' home TA B L E 1 The four cryptic bat species groups included in the study with information on their overlap in distribution, morphology and ecology (see Tables S1.1-S1.3 for overview on species ecological similarity). Information taken from Dietz and Kiefer (2016) range (3-10 km depending on species; Boye & Dietz, 2005), and seasonal or reproductive movements that typically occur outside species home range (e.g. Robinson & Stebbings, 1997). Before using these maps to calculate the measure of GA ij , they were projected to the Gall-Peters equal-area projection to allow for area calculations.

| Assessing the performance of the measure
We calculated the measure of GA ij for the pairs of potentially competitor species and tested whether observed values were higher than distributions of null values obtained for pairs of virtual species.
For that we created for each species 15-30 sets of virtual ranges (r package virtualspecies, Leroy et al., 2015) that follow environmental gradients, using the PCA method. We selected from the resulting virtual species ranges, the areas with highest suitability scores totalling the equivalent size of the realised range of the real species. We distributed within that area the same number of random occurrence records as used in the models of the real species. With this procedure, for each species we obtained different sets of occurrence records with the same number of records and covering the same area as the real species, but following different, though ecologically meaningful, environmental gradients and having different spatial distributions. We used these randomly generated sets of occurrence record to model virtual species' predicted ranges, and we clipped them to obtain their realised ranges, following the same procedures outlined above for the empirical dataset. Finally, we calcu- We only used pairs of sets of virtual ranges whose realised range overlapped by more than 0.5% to reduce inaccuracies due to the calculation of ratios with very small values. We compared for each pair of potentially competing bat species, the observed value of GA ij with at least 100 null values. Observed values of GA ij higher than 95% of null values were considered as significantly higher than null expectations. We also compared GA ij values of potential competitor pairs with GA ij between all combination of non-competitor pairs (115), which includes species pairs from different cryptic groups as well as pairs formed with the additional seven species. GA ij values between potential competitor pairs higher than 95% of the values between non-competitor pairs were considered as significantly higher.
We analysed whether the measure was robust to missing data and spatial sampling biases. For each species, we randomly removed 10, 20 and 30% of the occurrence records across the entire study area.
In addition, we divided species' range extents to nine equal quadrats; we randomly selected five and removed from these 10%-30% of the occurrence records to simulate spatial bias in the data. Both thinning processes were repeated 10 times at each of the three data removal percentages. We used these sets of thinned records and the original complete set of records to run SDMs (Table S3.7). We assessed whether GA ij values obtained using the full datasets fell within 95% confidence intervals of GA ij values from thinned datasets.

| Comparing measure outputs to patterns of ecologically similarity
Ecological similarity of potential competitor bat pairs was ranked from 1 (minimum) to 5 (maximum) in terms of foraging habitat, diet and roosting ecology by five bat researchers based on the literature and expert opinion (Appendix S1 com). We included in the models the same climatic variables as in the main models (Tables S3.2-S3.5), excluding habitat variables because of high uncertainty in land cover change projections, and altitude due to collinearity with climate. We combined the outputs from the three GCMs and considered future suitable areas if identified by two or more GCMs. We estimated range losses through comparing the future models to models generated for present conditions using the same variables. Then, we calculated predicted future geographic overlap between pairs of cryptic bat species within a 500 km buffer around their current range. This distance was selected to represent the maximum distance that these bats are likely to be able to disperse in response to climate change by the end of the century. As none of the studied species are long-distance migrants, they are not likely to reach all suitable future areas projected by the models.
Species pairs with high geographic avoidance and an increased future geographic overlap may not realise their future range change predictions due to the effects of competition.

| Geographic avoidance between cryptic bat species
Ensemble SDMs had good discrimination ability (AUC test range:    (Figure 4).

| D ISCUSS I ON
Our study suggests, using a new approach based on SDM outputs and empirical data from sets of cryptic Palearctic bat species, that interspecific competition can limit the geographic ranges of morphologically similar species in the absence of fine-scale mechanisms of coexistence. This study provides additional evidence for the potential of competition to shape species' geographic ranges (e.g. Gotelli et al., 2010;Sexton et al., 2009;Staniczenko et al., 2018). Moreover, it supports predictions made by mathematical models that competitive exclusion scales up and impacts species ranges only when coexistence mechanisms, such as trophic or habitat partitioning, are not developed at finer spatial scales (Godsoe et al., 2015). Biotic interactions, and in particular interspecific competition, can slow down climate tracking and prevent species from colonizing new habitats (HilleRisLambers et al., 2013).
Yet despite the urgent need to better understand the effect of biotic interactions on species vulnerability to climate change, there is currently no adequate methodology to infer competitive effects at the regional or continental spatial scale.

| The link between broad-scale geographic avoidance and fine-scale mechanisms of coexistence
Biotic interactions are primarily thought to play a role in driving local assemblages, which are, as a result, often composed of morphologically and functionally different species (Bocher et al., 2014;Codron et al., 2015). However, there is increasing evidence that biotic in- Observed geographic avoidance (GA ij Obs). Null distribution of GA ij values using virtual species, including mean (GA ij null mean), and 95% confidence intervals (GA ij null CI). Extent of ecological similarity between species pairs (range 1-5) according to foraging habitat, diet and roost selection and evidence of resource partitioning when sympatric based on expert opinion and the literature (details in Appendix S1, Tables S1.1-S1.3) species with similar morphology and broad-scale range overlap that coexist through fine-scale habitat partitioning (Tables S1.1-S1.3).

Species pair
Broad-scale biogeographic patterns in bats are mainly driven by temperature and water availability (McCain, 2006;Ulrich et al., 2007) especially in areas where access to water can be limiting (Razgour et al., 2018). However, several studies point to biotic interactions and in particular interspecific competition as important processes in the group (see Salinas-Ramos et al., 2020). Competition has been suggested to limit morphological similarity among rhinolophid bats in Malaysia (Kingston et al., 2000) and bat assemblages in southern Africa, where patterns of body size are regularly spaced (Schoeman & Jacobs, 2008). Studies using exclusion experiments in both tropical (Kalka et al., 2008) and temperate forests (Böhm et al., 2011) show that bats can control arthropod abundance, suggesting that exploitative competition may occur between bats due to prey depletion. Competition among bats can also occur through interference when species are attaining prey resources. Large aggregations of bats are thought to forage less effectively as a consequence of echolocation interference (Amichai et al., 2015), and bats can actively use jamming calls to disrupt competitors and make them miss targets (Corcoran & Conner, 2014). Our study suggests that competition among bats can also scale up to affect their broad-scale geographic distributions in the absence of fine-scale mechanisms of resource partitioning.

| Performance of the measure of geographic avoidance
No currently available methodology enables direct inference of competitive interactions from analysing co-occurrence patterns at broad spatial scales given the commonly available density of occurrence data.
Our measure of geographic avoidance quantifies the expected geo-   Table 2 are because these models did not include land cover variables The applicability of these methods is however limited by their assumptions of identical environmental preferences, sampling effort and detectability for both species in the potential predicted overlapping area. The method of including a competitor's realised range in SDMs is thought to primarily suit simple systems (Anderson, 2017) and correlations between environment and competitors make difficult to distinguish the effect of competition from environmental effects (Godsoe et al., 2017). While JSDM approaches directly incorporate biotic interactions into the models, they are most suitable for spatially-discrete community level studies, and the common lack of available representative community data across large spatial extents limits their application in broad-scale studies. The use of JSDMs with insufficient sampling intensity would result in a large amount of false absences that the model may interpret as negative residual correlations.
The approach presented here is based on the assumption that SDMs are able to estimate representatively the predicted and realised ranges of the species, which, equally to any other study based on SDMs, is a function of the quality of the occurrence data used and the inclusion of all the relevant environmental variables in the models (Guisan et al., 2017). While we show that the measure is considerably robust to missing data and spatially biased data, the occurrence data used should offer a good representation of the realised range of both species. The realised ranges for our case study species have been extensively studied by the authors and therefore are best represented by our data. However, whether using either own data or IUCN range maps, uncertainties in realised range and extent of realised overlap among species should always be carefully taken into account. Hence, this approach is not suitable for data-deficient species with under-studied range extents.
Lower ratio of realised range in overlapping areas relative to nonoverlapping areas could occur for other reasons than competitive effects, such as dispersal limitations in that area (Guisan et al., 2017;Soberón & Peterson, 2005) or by chance due to inaccuracies in model predictions. Therefore, interpretations of the measure outputs should consider the potential effects of dispersal limitations. Dispersal limitation is less likely to bias the conclusions with volant animals like bats and birds, but more so with dispersal-limited taxa, such as other small mammals, reptiles and amphibians. The comparison between the observed patterns with null models based on virtual species allows to identify the cases where the observed pattern is greater than expected by chance and therefore most likely to be result of competition.
Altogether, this approach is most suitable for pairs of species with well and equally known ranges and little dispersal restrictions.
Inferences given by the measure should be taken as a geographic pattern consistent with a process of broad-scale competitive exclusion, and not as a direct inference of the effect of competition. A sensible interpretation of current biogeographical patterns also requires knowledge of evolutionary history to tease-apart biotic and historical effects, in particular for species with limited dispersal and colonization abilities (Dormann et al., 2018;Warren et al., 2014), and information on ecological similarity between species.

| Including biotic interactions in range shift projections under climate change
Biotic interactions can modify predicted responses of species to climate change (HilleRisLambers et al., 2013), as highlighted F I G U R E 4 Present and future predicted suitable range for E. serotinus (a), E. isabellinus (b), P. austriacus (c) and P. kolombatovici (d) and extent of range overlap with their potential competitors. Orange shows current range that will become unsuitable in future. Green shows current range that will still be suitable in future, and blue shows current unsuitable range that will become suitable. Black shows future suitable range that will overlap with future suitable range of a competitor [Colour figure can be viewed at wileyonlinelibrary.com] (a) (c) (d) (b) by theoretical models assessing species extinction risk (Norberg et al., 2012). Climate change can lead to changes in the intensity of interactions or the appearance of novel interactions (Alexander et al., 2016). While the outcome of the interaction between biotic interactions and range shift processes is difficult to predict, considering the patterns of species geographic avoidance presented in our study could be of relevance. In species pairs with high geographic avoidance, if future suitable range overlap is predicted to increase, model projections that ignore the impact of strong biotic interactions will overestimate future range suitability. In contrast, a decrease in predicted future suitable overlap could lead to an unexpected larger suitable area due to competitive release. In our study system, among the pairs of species whose current ranges are most likely shaped by competition, future models predict a decrease in range overlap between E. serotinus and E. isabellinus, which would lead to competitive release effect and larger realisation of their future suitable ranges than at present. Conversely, the ranges of P. austriacus and P. kolombatovici are predicted to overlap substantially more in the future, which may limit the ability, in particular of P. kolombatovici that already has a very restricted range, to shift its range to track future suitable conditions. We show that disregarding biotic interactions can affect our ability to accurately predict species future distributions and their vulnerability to climate change. We also stress the importance of generating future range suitability projections that indicate areas where interspecific competition may limit range shifts.

| CON CLUS IONS
Our study suggests that in absence of fine-scale mechanisms of resource partitioning, the effects of interspecific competition can scale up to impact species broad-scale geographic ranges. Therefore, this study highlights the importance of considering the effects of biotic interaction when predicting future range suitability under climate change. While a better comprehension of the operation of biogeographical processes across spatial scales requires better integration of community ecology process with larger-scale species distribution models, the measure of geographic avoidance that we present can be used to identify range patterns compatible to a broad-scale effect of interspecific competition. Its low data and computational requirements allow it to be widely applied as a screening tool to identify cases where biotic interactions could impact future climate change predictions.

ACK N OWLED G EM ENTS
We are grateful to Petr Benda and Helena Santos for providing location records and publications and to Jane Catford, Felix Eigenbrod and Kelvin Peh for comments on the work. RN-F was funded through a University of Southampton PhD studentship.
OR was funded through a NERC Independent Research Fellowship (NE/M018660/1). This work was supported by the Bavarian State

Ministry of Science and the Arts via the Bavarian Climate Research
Network bayklif (project "mintbio").

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