Introduction

Primary productivity (hereafter “productivity”) and diversity (here defined as species richness) are closely related at global and regional scales. The most productive latitudes and elevations have the most species, and vice versa. Termed the productivity–diversity relationship, this correlation has been documented across plants (Currie et al. 2004, Šímová et al. 2011), invertebrates (Kaspari et al. 2000a, Beck et al. 2011, Classen et al. 2015), and vertebrates (Hawkins et al. 2003a, Hurlbert 2004, Evans et al. 2005, Jetz and Fine 2012, Ferger et al. 2014, Guan et al. 2016, Seoane et al. 2017, McCain et al. 2018), in many localities across the globe. Proposed mechanisms include physiological adaptations (e.g., Turner et al. 1987, Currie 1991, Hawkins et al. 2003b), diversification rates (e.g., Brown et al. 2004, Mittelbach et al. 2007), habitat heterogeneity (e.g., MacArthur and MacArthur 1961, Stein et al. 2014, 2015), the more-individuals hypothesis (e.g., Srivastava and Lawton 1998, McCain et al. 2018, Storch et al. 2018) among others. But despite numerous empirical studies, the support for various proposed mechanisms are mixed and inconclusive.

The most commonly invoked mechanism is the more-individuals hypothesis (MIH; Srivastava and Lawton 1998, McCain et al. 2018, Storch et al. 2018). The MIH proposes that productivity determines the quantity of food resources, which then regulates the number of individuals that can survive at a site. The number of species at a site is then determined by the number of individuals, as a species with low abundance is more likely to go locally extinct. In this way, more productive sites produce more food resources, which can support larger numbers of individuals and in turn can be divided into more species with stable populations (Wright 1983, Wright et al. 1993, Srivastava and Lawton 1998; Fig. 1a). Although the MIH mechanism can be criticized for lacking the potential for bidirectional causation (e.g., diversity influencing abundance) and the absence of interspecific effects (among others), it remains the most commonly tested mechanism of the productivity–diversity relationship.

Empirical evidence for the MIH remains equivocal even after numerous tests in plant (e.g., Currie et al. 2004, Pautasso and Chiarucci 2008, Šímová et al. 2011), invertebrate (e.g., Kaspari et al. 2000a, Currie et al. 2004, Yee et al. 2007), and vertebrate systems (e.g., Gaston and Evans 2004, Buckley and Jetz 2010, McCain et al. 2018). Studies include meta-analyses (e.g., Currie et al. 2004, Pautasso et al. 2011), manipulative experiments (e.g., Srivastava and Lawton 1998, McClain et al. 2016), field surveys (e.g., Kaspari et al. 2000b, Seoane et al. 2017), as well as a few mechanistic tests of multiple components of the MIH (Ferger et al. 2014, Classen et al. 2015, McCain et al. 2018). When taken as a whole, MIH results are a mixed bag of significant and nonsignificant relationships (e.g., Ferger et al. 2014, Classen et al. 2015, McCain et al. 2018, Storch et al. 2018, and references therein).

Improper formulation of the quantitative predictions of the MIH, ignoring potential feedback between levels of the MIH mechanism (e.g., diversity influencing community abundance), may explain the inconclusive nature of the literature (Storch et al. 2018). On the other hand, mixed results may also be explained by sampling issues. MIH studies are biased toward bird and insect systems, both of relatively small average body size (McCain et al. 2018, Storch et al. 2018). These systems typically exhibit relatively high interannual population variability (Ricklefs 1990, Begon et al. 1996, Morris and Doak 2002) therefore producing single-year abundance measurements not indicative of the long-term trend. Variation introduced by single-year abundance measurements could mask the relationship between long-term abundance and diversity, providing evidence against the MIH. Notably, all three studies testing multiple components of the MIH used organisms prone to high levels of population variability while sampling for only 1 or 2 yr (birds: Ferger et al. 2014, bees: Classen et al. 2015, small mammals: McCain et al. 2018). Multiyear sampling could compensate for population variability, but this is rarely done in MIH studies. There are a few long-term meta-analyses of the MIH (e.g., Carnicer and Díaz-Delgado 2008, Dobson et al. 2015). However, we are not aware of any long-term (>3 yr) MIH-specific field studies.

Thus, in empirical tests of the MIH, it is possible that natural population variability could be masking the “true” process—a positive abundance–diversity relationship—that exists on a longer time scale. Herein, using a simulation framework, we investigate whether natural population variability could be masking the MIH in empirical studies. These results will determine our confidence in the ability of previous studies to detect the MIH, and will provide insight to the sampling effort required to detect the MIH effectively, if it exists as a mechanism in nature.

Methods

To simulate the mechanism of the more-individuals hypothesis (MIH), linear relationships were fixed among the components of the mechanism (excluding food resources for simplicity). First, a one-to-one linear relationship was fixed between productivity (range: 500–6,450) and diversity (range: 3–20; herein diversity represents species richness) across simulated sites to mimic a productivity–diversity relationship (Fig. 1b). This range of diversity was chosen based on the three empirical studies examining multiple components of the MIH (mean = 21; Ferger et al. 2014, Classen et al. 2015, McCain et al. 2018). To examine this question in highly diverse systems, we also ran simulations with greater diversity (i.e., 100, 500 maximum species richness; Appendix S1). Additionally, a one-to-one linear relationship was fixed between productivity and long-term average community abundance (Fig. 1b). This produces a fixed linear relationship between long-term average community abundance and diversity, as the MIH predicts (Fig. 1b).

At each site, the long-term average community abundance was divided among the number of species present according to a log-normal distribution to mimic the abundance structure of a natural community (few abundant and many rare species; Preston 1948, Ulrich et al. 2010). This produces a long-term average abundance for each species. The log-normal species abundance distributions were simulated using the sim_sad function in the R package mobsim (May et al. 2018; Appendix S1). We simulated log-normal distributions across three levels of community evenness, as more even communities may show reduced effects of population variability (Appendix S1). Here, we present moderate community evenness as it is the most biologically relevant (details in Appendix S1).

To maintain stable average abundance over time, the initial population and the carrying capacity ( $N_0$ and $K$ , respectively) were set equal to the species’ 50-yr average abundance determined by the log-normal distribution, and the average log intrinsic growth rate was set equal to one for all simulations ( $\bar{r}=1$ ). Population variability values ranged from 0.01 to 2.0 across simulations, based on known values from the literature (McCain et al. 2016, and references therein). Species extinctions were prevented by setting a minimum population size of two individuals.

Species abundances from all species in each year were summed, producing 50 yr of community abundance for each site (Fig. 1c). The community abundance was then sampled for 1-, 2-, 3-, 4-, 5-, 10-, 20-, and 40-yr durations, starting at year 10 (Fig. 1d). We assumed perfect detection of all individuals in the community for each year in the sampling duration. Mean sampled community abundance was used in simulations with more than 1 yr sampled, so that each site had a single value of sampled community abundance regardless of the number of years sampled.

For each sampling duration across all sites, linear regressions of diversity on sampled community abundance were constructed. Prior to regression, feature scaling (normalization to range of 0–1) was used so that a “perfect” regression in the absence of population variability would yield a coefficient of one. Additionally, this means that the original fixed relationship between community abundance and diversity would have a slope of one. We then evaluated the linear models (here called detected relationships) according to three model attributes: detected slope (coefficient of the sampled community abundance regression term), significance of the slope (P < 0.05), and variation explained (R²). We evaluated these model attributes across population variability, sampling duration (number of years sampled), and number of simulated sites (n = 18, 36, 54), across 1,000 simulations for each combination of parameters.

In addition to the above procedure, two separate alterations were made for two additional sets of simulations. (1) We sampled abundances in nonconsecutive years to compare with consecutive sampling (e.g., sampling first and third years vs. first and second years) for sampling durations of less than 5 yr. A simulation study evaluating range shift detections in the face of population variability found that in some cases nonconsecutive sampling improved detection efficiency (McCain et al. 2016). We included nonconsecutive sampling to assess this possibility in these simulations. (2) We designed simulations with a random fixed relationship between productivity and 50-yr average community abundance. This allows us to assess how likely we are to detect a spurious relationship that resembles the MIH solely due to population variability. In each of these two instances, all other aspects were kept consistent with the above methodology; more details can be found in Appendix S1.

All simulations were coded and run in R, version 3.5.2 (R Development Core Team 2018). Code is available in Data S1, and on GitHub (http://doi.org/10.5281/zenodo.3686307).

Results

Higher interannual PV resulted in decreased detection of the more-individuals hypothesis (MIH) across all three metrics. As PV increased, mean detected slope decreased and the range of detected slopes widened (Fig. 2a). For example, for single-year sampling across 18 simulated sites, a PV of 0.1 yields a mean detected slope of 0.99, and a range of [0.77, 1.2], and a PV of 2.0 yields a mean detected slope of 0.51 and a range of [−0.25, 0.97]. With greater PV, the proportion of significant (P < 0.05) detected relationships decreased (e.g., PV = 0.1: 100% sig., PV = 2.0: 48% sig.; 18 sites, 1 yr; Fig. 2a, b). The amount of variation explained by the detected relationship (median R²) decreased with increasing PV. For example, the variation explained at highest PV (2.0) was roughly 1/5 of that at lowest PV (0.01) for 18 sites and single-year sampling (Fig. 2c).

With a longer sampling duration (i.e., more years sampled), mean detected slopes tended to approach the fixed slope of one (1 yr: 0.51, 10 yr: 0.88, 40 yr: 0.98; PV = 2.0, 18 sites), and the range of detected slopes shifted towards and narrowed around the fixed slope of one (1 yr: [−0.25, 0.97], 10 yr: [−0.02, 1.19], 40 yr: [0.49, 1.23]; PV = 2.0, 18 sites; Fig. 2a). The proportion of significant detected relationships increased with a longer sampling duration (1 yr: 48%, 2 yr: 69%, 5 yr: 93%; PV = 2.0, 18 sites), as did the variation explained (median R²; 1 yr: 0.21, 2 yr: 0.3, 5 yr: 0.5; PV = 2.0, 18 sites; Fig. 2b, c).

With an increase in the number of simulated sites, mean detected slopes grew closer to the fixed slope of one (18 sites: 0.51, 36 sites: 0.58, 54 sites: 0.63; PV = 2.0, 1 yr), and the ranges of detected slopes shifted upwards towards the fixed slope (18 sites: [−0.25, 0.97], 36 sites: [−0.26, 1.07], 54 sites: [−0.12, 1.14]; PV = 2.0, 1 yr; Fig. 3a). The proportion of significant detected relationships increased substantially across all PV values given an increase in the number of simulated sites (18 sites: 48%, 36 sites: 73%, 54 sites 82%; PV = 2.0, 1 yr; Fig. 3b). But an increase in the number of simulated sites led to a slight decrease in the amount of variation explained (median R² for 1 yr, 18 sites: 0.21, 36 sites: 0.17, 54 sites 0.16; Fig. 3c).

Nonconsecutive sampling did not yield substantial differences from consecutive sampling (Appendix S1: Fig. S2). Detection of a spurious relationship due to PV given a random underlying fixed relationship was rare. Mean detected slopes were near zero (though with a range of [−0.79, 0.87]; Appendix S1: Fig. S3), the proportion of significant detected relationships was near 5%, and median R² was less than 0.04 (details in Appendix S1). Results for varying levels of community evenness and for higher diversity simulations are shown in Appendix S1.

Discussion

Population variability (PV) can obscure the more-individuals hypothesis (MIH) relationship given low sampling effort in the form of a small number of years sampled across a small number of sites. This is especially true at the highest PV (2.0) and the lowest sampling effort (1 yr, 18 sites) where less than half of detected relationships were significant. Additionally, detected relationships remain poor in moderate–high PV (>0.5) simulations for sampling effort of 1–2 yr across any number of sites. With low PV (≤0.5), detected relationships fit the fixed relationship relatively well, even with low sampling effort. For reliable detection of the underlying fixed MIH relationship in these simulations, multiyear sampling (5+ yr) across many sites (36, 54) was required.

Most MIH studies occur in systems such as birds (e.g., Gaston and Evans 2004, Ferger et al. 2014), and invertebrates (e.g., Currie et al. 2004, Classen et al. 2015), which are prone to moderate–high levels of PV (Ricklefs 1990, Begon et al. 1996, Morris and Doak 2002). Additionally, most of the MIH field studies we found only sampled their communities for a single year, and none had more than 3 yr of data. Thus, negative MIH results in the literature could be due solely to population variability. We can be slightly more confident of supportive evidence for MIH results in the literature, because spurious positive relationships were very rare given a random underlying relationship (though there are still potential confounding factors).

To avoid additional assumptions and complexity, each species in a simulated community had the same population variability. This is unlikely in nature, where communities are made up of species with different degrees of population variability (e.g., Crowley and Johnson 1992, Brady and Slade 2004), and may have resulted in higher-than-expected variability in community abundance. However, in communities with a log-normal species abundance distribution, low-abundance species have less influence on the community variability (in simulation; but see Appendix S1 for highly even communities). Any variation added beyond what we would expect in nature due to constant population variability across species is likely swamped by the enormous variation inherent in field studies.

We designed these simulations to represent a conservative estimate of the degree to which population variability could mask the MIH. These simulations fixed perfect, one-to-one linear relationships, which are rare or nonexistent in nature. Imperfect (those with unexplained variation) underlying relationships would be more easily masked by population variability. These simulations assumed perfect sampling of abundance and diversity, though sampling bias or missed species detections due to small population size could further confound detection of the MIH in field studies (Link et al. 1994, Wilson et al. 1996, Magurran 2004). Additionally, food resources were excluded from these simulations, but in practice would have their own variability over time. Real-world considerations make it less likely to detect the MIH in the field reliably, if it exists, leaving these simulations as a best-case scenario for detection of the MIH.

Future MIH studies need to be designed with population variability in mind. Choosing study systems with low population variability (e.g., trees, ungulates, large birds; Harper 1977, Clutton-Brock et al. 1997, Sæther and Engen 2002) can compensate for a low sampling effort. However, these systems typically lack adequate variation in richness, productivity, and abundance to test the MIH effectively. In the systems commonly used for MIH studies (birds, insects, small mammals), 5+ yr of sampling across many sites (18+) is required to ensure high probability of detection of the MIH if it exists. To be confident that negative MIH results are not due to population variability, given sampling inaccuracies and additional food resource variability, even longer sampling durations are recommended (10+ yr).

The mixed evidence from empirical studies of the MIH leaves us uncertain about whether it is truly a mechanism of the productivity–diversity relationship. However, despite the large body of literature on the subject, few to no studies have the sampling effort required to be confident in their results from underlying interannual population variability. In particular, studies without support for the MIH may be suspect because of lack of adequate sampling. These simulations point to the gravity of considering population variability for past and future studies.