Main

A positive interspecific relationship between abundance and distribution —abundance-occupancy relationships (AORs) — is considered one of the most general and robust patterns in ecology ^1–4. Sometimes referred to as a macroecological law ^5,6, the AOR asserts that empirically locally abundant species tend to be widely distributed, and conversely, locally rare species tend to be geographically restricted in their range. The mechanism driving this relationship was never proven, and it remains unresolved why species distribution should affect per-unit-area abundance (or vice versa). Nonetheless, the existence of a pervasive AOR has underpinned many practical applications in ecology and conservation ⁷, for example, setting harvest rates for fisheries ⁸, managing invasive species by restricting expansion rather than local elimination, and identifying species at high risk of extinction in biodiversity inventories such as the IUCN Red List Criteria ⁹. Given the increasing human-induced land-use changes in the Anthropocene ¹⁰, concomitantly with increasing debate about global biodiversity change [cf. ¹¹], fully understanding the relationship between abundance and range size is of increasing importance.

Many plausible biological mechanisms have been proposed for AORs, yet none of them has unequivocal support ^3,4,12,13. Among all mechanisms, it is noteworthy that a spatially explicit neutral model of biodiversity and biogeography can generate AORs ^14,15. Specifically, this macroecological ‘null’ model can produce a positive correlation between species range (or occupancy) and their per-unit-area local as well as total global abundance. This observation, in turn, supports the utility of neutral theory as a null model of community and macroecology ¹⁴. Although neutral theory may provide a biological null model, an additional null hypothesis is that AOR does not exist. Indeed, sampling bias can create AORs because locally rare species are more likely to be missed, resulting in an underestimation of range size or occupancy, thereby generating a positive relationship ^13,16. Yet, this sampling explanation has long been discarded as a plausible mechanism leading to observed patterns ^2,3,17. This is because of substantial empirical evidence for positive interspecific relationships, including a meta-analysis of 279 effect sizes with an overall effect of r = 0.58 (or its Fisher’s transformation: Zr = 0.66) in 2006 ¹. It does not seem that sampling bias alone could explain this remarkably strong relationship.

Nonetheless, a large amount of variation does exist in empirical patterns of AORs, including strikingly negative relationships ^12,18-20. Some of the observed heterogeneity is likely to be due to different aspects of sampling, such as the number of species and spatial and temporal coverage ^3,4,12. Also, other types of bias could generate artefactual AORs: namely ‘confirmation bias’, where sampling is prejudiced to support one’s hypothesis and ‘publication bias’, where statistically significant relationships are preferentially reported and published. Although both biases are widespread, including in ecological studies ^21–23, no studies so far systematically considered or quantified both biases in the context of AORs [cf. ¹]. Furthermore, there has, until recently, been a lack of large and methodologically consistent data resources, therefore leaving a traditional meta-analytic approach as the best available option for testing the validity and generality of the AOR.

A citizen-science dataset to test AORs

Here, we use data from eBird — a global citizen science dataset aimed at counting birds — to quantify the relationship between abundance and occupancy (AOR). Large citizen science datasets collected for non-hypothesis-driven purposes are not random samples [see ²⁴], but they have the advantage of avoiding biases such as confirmation and publication bias. Also, using the eBird dataset allows us to estimate heterogeneity due to sampling intensity (e.g., the duration of a sampling event directly influences the number of species recorded). Specifically, we can quantify how AOR will change in relation to increases in species richness and sampling duration, both of which are predicted to reduce the magnitude of AORs ^3,4,19.

For occupancy, we use global range size not only because global range size should be relatively stable — ‘local’ range sizes for one species could vary dramatically — but also because different types of occupancy measures were deemed to contribute less to the observation heterogeneity ^25,26. Fortunately, for birds, a large database of global range sizes has already been compiled ²⁷. For abundance, we use two different measurements: local species counts and local mean density, as follows. First, we carry out the largest known meta-analysis by synthesizing correlations between global range sizes of 7,635 species and local species counts collected across 16,562,995 eBird checklists (resulting in 16,562,995 Zr values and corresponding sampling variances; Fig. 1). These checklists all included counts of each species present and the duration of observation (hereafter, effort time).

Figure 1

Second, we conduct a phylogenetically controlled comparative analysis, regressing species range sizes on 7,464 estimates of globally derived species’ mean density, equivalent to mean local density (per 5-degree grid cell), estimated in earlier work ²⁸ (see Methods for more details). Given the different potential biases mentioned before, we expected a more modest relationship in relation to that of the previous meta-analysis (r = 0.58; note that this relationship included many different taxa; if restricted to bird species, it was even stronger r = 0.74 or Zr = 0.95) ¹⁷. Also, although no such empirical evidence appears to exist, it seems feasible that in filling in an eBird checklist, some people may undercount common and widespread species while they may overcount rare and geographically restricted species. If this is the case, the relationship (AORs) could be further weakened. Yet, if such overcounting and undercounting were present, we expect it would introduce large heterogeneity into our dataset because that type of behaviour would not be consistent across all contributors, and they would sometimes result in negative AORs, increasing variability among the 16,562,995 Zr values.

Overwhelming support against AOR

Surprisingly, the overall (aggregated) relationship between local abundance and global occupancy was near-zero (r = 0.015), although this relationship was statistically significant due to our extremely large sample size (p = 0.0005, z = 2.805, Zr = b_{[overall mean]} = 0.015, 95% confidence interval, CI = [0.004, 0.025]; Fig. 2, Table S1). However, this significant relationship disappeared (r = 0.0009) once we controlled for species number and effort time (p = 0.863, z = 0.173, Zr = b_{[overall mean]} = 0.0009, 95% CI = [-0.0092, 0.0111]); both variables were statistically significant predictors of the effect. As expected, the increase in species number (modelled as the inverse of species number - 3, which is equivalent to sampling error for Zr) and effort time on the natural log scale, decreased the strength of the relationship (sampling variance: p < 0.0001, z = 140.29; b_{[sampling variance]} = 0.0147, 95% CI = [-0.0149, -0.0145]); ln(effort time): p < 0.0001, z = -183.45, b_{[ln(effort time)]} = 0.230, 95% CI = [0.226, 0.233]; marginal R² = 5.1% for the model with these two predictors; ²⁹; Table S2-4). These observations are consistent with the explanation that sampling protocols can create positive artifactual relationships between range and abundance.

Figure 2

Even more remarkably, our meta-analysis suggested that the AOR is likely indistinguishable from zero even with a larger dataset because the observed heterogeneity among effect sizes was very small (i.e., most effect sizes were effectively zero after accounting for sample size). A measure of relative heterogeneity I²_[total] was 13.5%, meaning that 86.5% of all the observed variation (in Fig. 2) is due to sampling error, therefore, is neither biological nor ecological (country level; I² = 5.5%, σ² = 0.005; state level: I² = 6.3%, σ² = 0.005; effect-size level; I² = 1.5%, σ² = 0.001); in contrast, the average I²_[total] across 86 ecological meta-analyses was approximately 92% ³⁰, making our observed heterogeneity unusually low. Low relative heterogeneity, however, does not necessarily mean absolute heterogeneity is also low ³¹. We found the absolute heterogeneity, σ²_[total] = 0.011, approximately one-thirtieth of the heterogeneity (σ²_[total] = 0.323) found in the previous meta-analysis ¹. Also, this is around one-tenth of the average heterogeneity (σ²_[total] = 0.125; median = 0.105), found in 31 meta-analyses in ecology and evolution ²³. Overall, low relative and absolute heterogeneities indicate that our dataset of 16,562,995 effect sizes does not have much variability left to be explained despite our observations coming from many different locations across the globe and contributed by tens of thousands of individual birdwatchers. Importantly, we emphasize that this combination of zero effect and very small heterogeneity is only expected when a particular phenomenon is not real.

Moreover, our phylogenetic comparative analysis, which accounted for phylogenetic uncertainty ³², corroborated our meta-analytic results [cf. ³³]. The global range sizes had little predictive power on mean species density (both on log10; p = 0.808, t_99.6 = 0.0227, b_[slope] = 0.0928, 95% CI = [-0.1615, 0.2068]; Fig. 3, Table S5). Taken together, our results provide overwhelming evidence against the fundamental relationship between species range and local abundance, while the results are consistent with this relationship as a sampling artifact. Nevertheless, our results are also consistent with previously published empirical evidence. This is because we have shown that relationships between global species ranges and local counts can be null, strongly negative, or strongly positive, which can be generated primarily by sampling (error) variance (shown in Fig. 2).

Figure 3

Lawless macroecology and non-neutral theory: implications

Our results demonstrate clearly that the AOR is not observed in a very large global dataset, with both applied and theoretical ramifications. First, we must reconsider fishing quotas, conservation priorities, and invasive species control strategies based on AORs [cf. ⁷]. Second, it is probably futile to pursue all ecological or biological mechanisms proposed for AORs [see ^3,13]. We cannot exclude, however, the possibility of AORs occasionally emerging in some restricted areas because there was a small unexplained variance in the meta-analytic dataset. Most notably, our near-zero results with small heterogeneity suggest that contrary to earlier suggestions ^14,15, the spatial neutral model is not a suitable null model of macroecology. Within the neutral theoretical framework, AORs can be broken by local adaptation (an alternative hypothesis) ¹⁴. If local adaptation were to disrupt a predicted positive relationship (AOR), we would have observed substantial heterogeneity and a reduced relationship, not a near-zero relationship. This is because it is extremely unlikely that local adaptation and neutral processes are in a perfect balance, resulting in an exact-zero relationship with little heterogeneity. In other words, local adaptions are expected to create local specificities and global variability/heterogeneity; our meta-analysis did not find such variability.

We point out that any model of a positive AOR posits a mechanism that connects species range and local abundance, whereas Rabinowitz effectively decoupled these two variables ³⁴. Our results are consistent with this decoupling. By adding habit specificity to species range and abundance, she suggests seven forms of rarity, which could reflect different underlying macroecological mechanisms; for example, two forms of rarity are geographically restricted but locally abundant with narrow or broad habitat types. Unlike the world created by the spatial neutral model, our results support many such ‘rare’ species ³⁵. In this regard, it is no longer surprising that AORs do not exist. Therefore, we believe Rabinowitz’s framework, rather than the AOR, has more empirical support from global-scale patterns of species abundance and provides a useful conceptual structure for future theoretical and applied work, although this line of work is still limited.

A credibility revolution in ecology beyond biases and crises

Whilst we provided an explanation for the non-existence of AOR, it still feels hard to comprehend the extent of overestimation in the previous meta-analysis (r = 0.58 for all taxa; r = 0.74 for birds) ^1,17. We have shown that some bias may be due to sampling bias. However, we speculate that much of the overestimation originates from publication bias and confirmation bias, which is supported by mounting evidence from meta-research studies ^21–23. Although we do not have direct evidence, our eBird datasets are free from these two types of biases (i.e., birdwatchers generally do not think of the macroecological patterns that would later be tested with the data they submit), while the literature-based meta-analyses are not. Regarding publication bias, the original meta-analysis of AOR states, “the fail-safe number indicates that more than half a million unpublished null results would be required to nullify an effect of this magnitude” ¹. Indeed, we provided much more than a half million null effects to reach our null conclusion (Fig. 2), although we should note that large datasets like eBird have other biases than publication or confirmation biases ²⁴. For example, it is possible that by excluding checklists with a single ‘X’ (see Methods), we are preferentially removing abundant species as birdwatchers may report ‘X’ for more common species with high abundances.

A recent study reexamining 86 ecological and evolutionary meta-analyses demonstrated a 23% reduction in overall effects due to publication bias, turning 33 of 50 statistically significant meta-analytic conclusions (66%) into non-significant ²³. Similarly, a study examining 83 topics in life sciences showed that the effect size of non-blind studies, which are at risk of confirmation bias, was twice as large as blinded counterparts protected against confirmation bias ²¹ [see also ²²].

We finish with an intriguing parallel topic to AORs in psychology, where the current replication crisis started ^36,37. There have been over 100 studies, and many theoretical models support the hypothesis of ‘ego depletion’ where self-control is a finite resource, so self-control will decrease once it is exerted ³⁸. The first meta-analysis of ego depletion, like AOR, suggested a very strong support for it (standardized mean difference, or d = 0.62). Yet a subsequent multi-lab replication found that ego-depletion does not exist is so weak as to be negligible (d = 0.04) ³⁹. Indeed, a series of multi-lab replications has indicated that several psychological phenomena, which were once believed to be real beyond a reasonable doubt, are too weak to be useful or are nonexistent⁴⁰. In ecology, recently collated large datasets collected for non-hypothesis-driven purposes offer a unique opportunity to revisit and retest longstanding ideas.

Taken together, we call for reexamining all ecological laws, rules, and patterns, as very few topics are free from sampling, confirmation, and publication biases [cf. ⁴¹]. To counter such biases, we urgently require a ‘credibility revolution’, a more optimistic name for a replication crisis, turning this crisis into an opportunity to improve science. A credibility revolution in ecology, like in psychology, needs to embrace non-traditional to avoid confirmation and publication bias, such as pre-registration ⁴², registered reports ⁴³, prospective and living meta-analyses, open synthesis communities ⁴⁴, and big-team-science collaborations ⁴⁵ involving community (citizen) scientists ⁴⁶.

Methods

Quantifying abundance-occupancy relationship at the local-scale

We used the eBird dataset ^47,48 to assess the relationship between local-scale abundance and occupancy (i.e., global range size). eBird, launched in 2002 by the Cornell Lab of Ornithology, is a global citizen science project which enlists volunteer birdwatchers to submit ‘checklists’ of birds seen and/or heard while birdwatching. Data undergo a semi-automated filtering process before being entered into the dataset, and expert reviewers additionally review species (or counts of species) that surpass pre-set filters before being accepted into the dataset ⁴⁹. Importantly, birdwatchers must indicate whether they are submitting a ‘complete’ checklist representing all birds that an individual birdwatcher was able to identify during their birdwatching outing. Further, birdwatchers can either submit the count of a species during their birding, or they can submit an ‘X’ to signify that a species was present but not estimate the number of individuals present during their birdwatching outing.

We downloaded the eBird basic dataset (version ebd_rel-May2020) and considered all eBird checklists between January 1st, 2005 and May 31st 2020. We then performed some quality assurance, applying an additional set of filters to the data, potentially removing any ‘outliers’ that could produce undue leverage on our results. The following filtering was completed [sensu ^24,50]. We only included: (1) complete checklists; (2) checklists that were <240 minutes and >5 minutes; (3) checklists that travelled < 5km; and (4) checklists that travelled < 500 ha. Because birdwatchers will sometimes use an ‘X’ to signify presence, and this is most likely to happen for more abundant species, we excluded any checklist which had at least any “X” on it, as this could potentially influence the correlation between the abundance of a species and range size by disproportionately removing the most abundant species from the correlation. We further only considered checklists that had at least ten species recorded on them, and a correlation test was performed only if we had range size data (see below) for a minimum of 4 species on the checklist.

We used range size maps from BirdLife International ²⁷. When an eBird checklist met the aforementioned criteria, we performed a correlation test using Pearson’s correlation coefficient from the cor.test function in R ⁵¹. Both the counts of every species and the range size were log-transformed before estimating a correlation (for a workflow, see Fig. 1). We obtained 16,562,995 correlations based on 3,005,668,285 individual bird observations, including 7,635 species. We note that we conducted all computational and statistical work using R and we created plots using the R package, ggplot2 ⁵², patchwork ⁵³ and their dependencies.

Meta-analysis of Big Data

We transformed correlations between species abundance and range into Fisher’s Z o Zr to unbound and calculated sampling variance for each Zr value; note that the inverse of the sampling variance of Zr is N (the number of species) – 3 (see Fig 1). We used the R package, asreml ⁵⁴ to run a multilevel random-effects model ⁵⁵; note the asreml is a commercial package, so it is not free. Our large meta-analyses with ~17 million effect sizes were only able to run with asreml given the computational time required for such a large dataset. We had ‘country’ (245 levels) and state code (2,871 levels) as random factors in the model to control for non-independence. In addition, to quantify the variance component for these two clustering factors and also at the level of effect sizes (16,562,995 levels), we modelled ‘units’ (the effect size level random effect or residuals) in the asreml function with ‘the number of species - 3’ as the ‘weights’ argument and asr_gaussian(dispersion = 1) as the ‘family’ argument. We also obtained the multilevel versions of I^{2 30,56} to obtain relative heterogeneity for our meta-analytic model (Table S1).

To gauge the impacts of potential biases, we fitted two moderators: 1) the z-transformed version of ln(checklist duration) as a surrogate for the amount of effort for observation and 2) sampling variance, which is usually used to detect publication bias, more specifically, small study bias where effect sizes from small studies can create ‘funnel asymmetry’, creating bias in meta-analytic overall mean ⁵⁷ (cf. Fig. 2). We ran two uni-moderator models and one multi-moderator model with both moderators (three meta-regression models in total; Table S2-4). We estimated the multilevel model versions of R^{2 29}.

Quantifying the abundance-occupancy relationship at the macro-scale

To corroborate our local-level analysis described above, we quantified an additional macro-scale analysis of the relationship between abundance (i.e., density) and occupancy (global range size). For this, we used data from a recently published analysis of global abundances for birds within 5-degree grid cells ²⁸. This dataset was derived by integrating expert-derived abundance measures with a large, less structured global citizen science dataset using a multiple-imputation technique to estimate density within 5-degree grids for 9,700 bird species (575 girds). We used these predicted density estimates from each grid cell and, for each species, took the mean of all density estimates in the grid cells for which a species was found. This mean density was then our measure of macro-scale abundance. For our measure of occupancy, we used a summation of all range sizes for the grids a species was found in, calculated by using range maps from BirdLife International focusing on the entire extent of a species’ extant range, ignoring the effect of transient species. Our analysis incorporated a total of 7,464 species of bird species, corresponding to the species for which we had both range maps and estimated density, along with phylogenetic information included in ⁵⁸.

Phylogenetic comparative analysis

To statistically test whether there was an effect of abundance and occupancy at the macro-scale, we used phylogenetic comparative analysis. We used avian phylogeny from ⁵⁸ and analysed 100 phylogenetic trees using the R function phylolm ⁵⁹. Resulting estimates from the 100 models were merged using Rubin’s rules, as described in ³², to obtain current estimates and errors that accounted for phylogenetic uncertainty (Table S5); we implemented this procedure using the R function miInference from the norm2 package ⁶⁰.

Acknowledgements

We thank Szymek Drobniak for helping run large meta-analytic models on a cluster computer system. We are also grateful to Malgorzata Lagisz, Tim Parker, Daniel Noble, Tim Blackburn, and Diana Bowler, and Luís Borda de Agua for their comments on earlier versions of this manuscript.

Funding

Australian Research Council Discovery Project Grant, DP210100812 (SN)

Author contributions

Conceptualization: SN, WKC & CTC

Methodology: SN, WKC & CTC

Investigation: SN, WKC & CTC

Visualization: SN & CTC

Project administration: SN

Writing – original draft: SN, WKC & CTC

Writing – review & editing: SN, WKC & CTC

Data and materials availability

All data, code, and materials are available online unless they are too large to be archived (https://github.com/itchyshin/AORs), and they will be archived in a public repository after the manuscript has been accepted for publication.

Supplementary Materials

Table S1

Table S2

Table S3

Table S4

Table S5