Accounting for demographic uncertainty increases predictions for species coexistence: A case study with annual plants.

Natural systems contain more complexity than is accounted for in models of modern coexistence theory. Coexistence modelling often disregards variation arising from stochasticity in biological processes, heterogeneity among individuals and plasticity in trait values. However, these unaccounted-for sources of uncertainty are likely to be ecologically important and have the potential to impact estimates of coexistence. We applied a Bayesian modelling framework to data from an annual plant community in Western Australia to propagate uncertainty in coexistence outcomes using the invasion criterion and ratio of niche to fitness differences. We found accounting for this uncertainty altered predictions of coexistence versus competitive exclusion for 3 out of 14 species pairs and yielded a probability of priority effects for an additional species pair. The propagation of uncertainty arising from sources of biological complexity improves our ability to predict coexistence more accurately in natural systems.

INTRODUCTION

One of the most fundamental goals in ecology is to understand the processes that maintain biodiversity and allow species to coexist (May, 1972). To this end, theoretical ecologists have focused on building models that can describe observed patterns of diversity and species co‐occurrence (Chesson, 2000; MacArthur, 1958). Modern coexistence theory (MCT) defined the necessary components for species to stably coexist, accounting for species’ niche overlaps and fitness differences (Chesson, 2000, 2020). Experiments testing coexistence mechanisms using individual fitness models have further investigated the determinants of local diversity (Hallett et al., 2019; Kraft et al., 2015; Levine & HilleRisLambers, 2009); often failing, however, to capture the true diversity observed in empirical systems and instead predicting the competitive exclusion of species found to commonly co‐occur over long time periods (Kraft et al., 2015; Wainwright et al., 2019). One possible reason for this discrepancy is that coexistence is often assessed using models that rely on single point estimates and disregard associated uncertainty in demographic and competition parameter values. We aimed to improve the realism of coexistence models by propagating the model estimated uncertainty using a Bayesian framework that yields a distribution of expected pairwise coexistence potential rather than a traditional point estimate of coexistence. We applied our framework using a controlled experiment in an annual wildflower community.

Modern coexistence theory is currently the most utilised framework for investigating species coexistence and the maintenance of diversity (Chesson, 2000, 2018). In this framework, species coexist if their population can recover from low density in the presence of other species in the community (Chesson, 2000; Turelli, 1980). This condition is central in MCT‐based methods for modelling coexistence such as the ‘invasion criterion’ (Chesson, 2000; Grainger et al., 2019; MacArthur & Levins, 1967; Turelli, 1980) and the ratio of niche differences to average fitness differences (we will refer to this as the ‘ratio approach’ (Chesson, 2000; Godoy & Levine, 2014)). Under the invasion criterion, two species can coexist if their average long‐term per capita growth rate is positive when introduced at a low density into a population of the other species at equilibrium. In the ratio approach, niche differences allow a species to increase from low density, and then regulate population growth at high densities (negative density dependence). Relative fitness differences represent the advantage that species have over each other when evaluated without niche differences (Chesson, 2018), leading to competitive dominance in the absence of niche differences. Hence, two species can coexist if the strength of their niche differences outweigh that of their average fitness differences.

Despite the developing sophistication of coexistence modelling, empirical tests of modern coexistence theory often produce models that predict far less coexistence than observations of long‐term co‐occurrences in natural systems suggest (Barabás et al., 2018; Kraft et al., 2015). For example, Godoy and Levine (2014) developed a framework derived from Lotka–Volterra competition models to estimate niche and fitness differences to predict coexistence between 153 species pairs. These authors found that only 12 out of 153 pairs were predicted to coexist, while for all other species pairs average fitness differences overwhelmed the niche differences, predicting competitive exclusion. In a recent study in the York gum–jam woodlands in Western Australia (the same system studied here), Wainwright et al. (2019) investigated six species pairs of annual plants and similarly found that fitness differences were almost always much stronger than niche differences, despite these species being observed to regularly co‐occur at local scales in the community. Other studies have found more convincing evidence of widespread coexistence. For instance, Adler et al. (2010) found that in a shrub and perennial grass system, there was an excess of niche differences, sufficiently strong to maintain diversity and stronger than necessary given small observed fitness differences. Similarly, Adler et al. (2018) synthesised across 21 studies and found, on average, stronger intraspecific competition than interspecific, suggesting possible stable coexistence.

We propose that the discrepancy between observation and MCT tests could also be a consequence of the tendency of coexistence models to exclude uncertainty in estimated demographic rates arising from unaccounted for sources of biological complexity. Uncertainty can be attributed to many sources including, but not limited to, demographic heterogeneity (genetic differences among individuals in their demographic rates; (Liancourt & Tielbörger, 2011), demographic stochasticity (variation in realised individual demographic rates despite having the same underlying rate; (Melbourne & Hastings, 2008)) and plasticity in individual trait values. These biological sources of variation can create uncertainty in parameter estimates, such as species’ intrinsic fecundity and their pairwise competition coefficients, even in the absence of observational errors. Therefore, by fully propagating uncertainty in parameter estimates and underlying variation—a long‐standing goal of coexistence theory (Chu & Adler, 2015; Stump et al., 2021) that is made relatively straightforward using Bayesian inference (Ellison, 2004; Gelman et al., 1995)—we can generate more realistic estimates for coexistence in complex biological systems.

Intuitive arguments suggest that the uncertainty generated by many sources of biological complexity may promote coexistence by ‘blurring’ competitive differences and subsequently yielding uncertainty in competitive hierarchies between species or substituting for niche differences. In a theoretical investigation, Hart et al. (2016) counterintuitively found that intraspecific variation promoted coexistence only under a specific set of circumstances, for instance when the variation was higher in otherwise inferior competitors. Rather, variation between individuals tended to either have no effect on coexistence or inhibited coexistence via three main mechanisms. First, individual variation tended to reinforce competitive dominance of superior competitors via non‐linear averaging. Second, niche variation between individuals tended to weaken coexistence by strongly reducing species‐level niche differentiation (Roughgarden, 1972; Slatkin, 1980). Third, variation in the intrinsic fecundity of individuals carries a demographic cost, particularly in small populations (Lande, 1993), reducing the likelihood of a species being able to recover from low density. Regardless of the underlying mechanisms (and even if we cannot separate the various sources of variation), incorporating uncertainty in model estimates arising from the multiple sources of biological variation has the potential to improve understanding of coexistence.

We used a Bayesian modelling framework to investigate coexistence in annual plant communities from the diverse York gum–Jam woodlands in Western Australia. In this framework, the probability of a coexistence outcome is a quantification of the uncertainty in that outcome given experimental bounds. We incorporated parameter uncertainty in the form of posterior intervals for density‐independent fecundity and intra‐ and interspecific interactions between co‐occurring species pairs from a single growing season and after removing the effects of spatial abiotic variation. We calculated probability distributions for both low‐density growth rates (LDGRs) and niche and fitness differences, providing estimates of the probability of coexistence outcomes we expect to see in natural systems (see Figure 1).

FIGURE 1

Panels (a) and (b) depict traditional point‐estimate tests of coexistence for a single species pair. Panel (a) represents the use of the invasion criterion and calculation of low‐density growth rate (LDGR). Here, the mean LDGR is calculated for two species, and a positive LDGR for each species when invading into the other at its equilibrium represents coexistence. If a species has a negative growth rate, we expect it to be competitively excluded. Likewise, panel (b) depicts the typical representation of coexistence predicted by the ratio of fitness differences to stabilising niche differences. The black dot represents a species pair that falls within the coexistence space (ρ < κ/κ green shaded region). In panel (c) (flipping the axes from panel (a) so that LDGR is on the x axis and different coloured lines represent the probability distributions of the two species involved in pairwise competition), we show that incorporating uncertainty into calculations of LDGR generates probability distributions of coexistence. In this example, species 2's distribution can overlap zero despite the mean LDGR (without propagating uncertainty; the point at the base of the distribution) falling in the positive region. In panel (d), we depict how, when incorporating uncertainty, the point estimate in panel b may yield a scatter of possible outcomes, with uncertainty about whether the pair falls within ‘coexistence space’. Similarly, this case can be represented as a joint probability distribution of coexistence

METHODOLOGY

Data collection

We collected data on an annual plant community from July to October 2018 in the York gum–jam woodlands in West Perenjori Nature Reserve (29°47'S, 116°20'E), south‐west Western Australia. This region experiences a Mediterranean climate, where winter rainfall triggers the germination of a diverse array of annual forb species.

We chose a mixture of common annual native (5) and exotic (3) species as study species for our field study, demographic modelling and coexistence analyses. The five native species were: Hyalospesrma glutinosum, Plantago debilis, Podolepis canescens, Trachymene cyanopetala and Goodenia rosea. The three exotic species were Arctotheca calendula, Medicago minima and Pentameris airoides. All these species are annual forbs, except for Pentameris airoides, which is an annual grass. To investigate the inter‐ and intraspecific interactions between each study species and their neighbours, and to disentangle estimates of the species’ density‐independent seed production from intraspecific density dependence, we manipulated the local interaction neighbourhood around individuals and recorded total seed production per individual as described below.

Using a spatially nested design, we established 16 (1 × 1 m) plots for each of our eight study species. We arbitrarily positioned plots across a study area of approximately 12 ha within West Perenjori Reserve. Within each plot, we placed six non‐overlapping 15‐cm diameter ‘neighbourhood rings’ centred on an individual or ‘phytometer’ of the study species for that plot (plots were species‐specific). Half of the rings within a plot were left un‐manipulated, with the abundance and identity of all individuals around the phytometer recorded, while the other half had all germinants except the phytometer carefully removed by hand at the beginning of the growing season. In the latter set of rings, we also removed delayed germinants during periodic checks of each neighbourhood throughout the growing season. We tracked the survival to reproductive maturity and seed production of each phytometer. To assess seed production, we placed thin mesh bags over immature fruiting bodies as these fruiting bodies were produced to minimise disturbance to each plant and minimise any loss of seeds. For most species, the phytometer's total seed production was counted by hand, excluding seeds that were unfilled. In the case of Pentameris airoides, the number of florets was recorded and then multiplied by two as a measure of fecundity as florets contain two seeds on average. We recorded all phytometers that died prior to seed production as having zero seed production (2% of total study plants). We conducted a test of the seed counting process itself (using the same seeds collected for this study) to demonstrate the very small amount of observation error (<6%) in our data set (Appendix 1, Table S1).

We chose a 15‐cm diameter ‘neighbourhood ring’ based on past studies in this system. This ring size has been shown to capture the majority of effects of local interactions on phytometer plants in this system (Martyn, 2020; Mayfield & Stouffer, 2017). Our plot size was sufficiently small that plants experienced near homogeneous abiotic conditions within each plot (Dwyer et al., 2015).

Statistical analyses

Annual plant fecundity model

To generate parameter estimates necessary to calculate predictions of coexistence, we fit Bayesian models for each study species’ fecundity within a given year following the commonly used Beverton–Holt annual plant model (Adler et al., 2007; Beverton & Holt, 1957). In this model, the expected fecundity of an individual of species i (F‐hat) is given by where λ is the intrinsic fecundity of species i (i.e. the average number of seeds produced by an individual in the absence of intra‐ or interspecific competition), α is the direct effect of species j on species i (accounting for intraspecific competition when j ≠ i and interspecific competition when j ≠ i) and N is the abundance of species j (here defined as the number of individuals). We fit the model for each of our eight study species. For each, we estimated α terms for the effect of each study species on the focal species. We additionally included a single α term for individuals of all other species not included in our list of eight study species. Thus, this term represented the average effect of individuals from non‐study (‘Other’) species on each focal species. If less than three individuals of one of the other study species co‐occurred with the focal species across all replicates, these individuals were added to the ‘Other’ category and no α term was estimated for that study species, assisting in model convergence. By discounting low‐abundance neighbours from our study species and grouping all other species into a single ‘Other’ category, we balanced tractability and accuracy. We acknowledge that some ecologically valuable information may be lost by grouping rare species in this study, but we felt the trade‐off with tractability was sufficient to merit this decision given the specific focus on this study.

To account for small‐scale environmental sources of variation in the data, we included a plot‐level random effect (ε) in all models except the model for A. calendula which would not converge with the inclusion of the environmental random effect. For all other species, this plot‐level random effect, ε, defined the impact of plot p on the expected fecundity of an individual of species i, allowing us to account for confounding effects of underlying environmental heterogeneity across plots. Observed fecundity (F) was thus modelled as F ~ poisson(ε) for all species except A. calendula, for which it was modelled as F ~ poisson(F‐hat).

In the Beverton–Holt model (Eq. 1), interaction coefficients that are facilitative (α < 0) may result in negative population densities. To prevent such occurrences, competition coefficients were defined as strictly positive (α ≥ 0) in all model fits and given standard, half‐normal priors: α ~ N +(0,1). Intrinsic fecundity (λ) was also strictly positive and given an uninformative gamma prior for most model fits: λ ~ gamma(0.001, 0.001). However, for Plantago debilis, this formulation of the model struggled to converge, so we instead defined λ as λ = exp(λ′) with λ given an unconstrained, uninformative prior: λ′ ~ N(0, 1000). The plot‐level random effects were drawn from a gamma distribution defined by a strictly positive hyperparameter σ (i.e. ε ~ gamma(σ, σ) ; (Lee et al., 2020)). In most model fits, the hyperparameter σ was given the same uninformative gamma prior as λ : σ ~ gamma (0.001, 0.001). However, similar to the issue arising from bounding λ in the model fit for Plantago debilis, this constrained model would not converge for T. cyanopetala, Podolepis canescens and M. minima. We used a similar approach for these species and defined σ as a transformation of the unrestricted parameter σ′ such that σ = exp(σ′) with σ′ given an unconstrained, uninformative prior: σ′ ~ N(0, 1000). For all model fits, we assessed convergence of chains using the Gelman–Rubin convergence diagnostic, checking trace plots for chain mixing and thinning as needed to remove autocorrelation. We ran three chains, sampling 6000 iterations and thinned by two iterations for all species except Plantago debilis for which we ran 9000 iterations and thinned by three iterations. We fit all models in R (Version 3.5.3) using the package Rstan (Stan Development Team, 2020).

Germination and seed survival rates

To model coexistence through both the invasion criterion and ratio approach, we estimated rates of germination (g) and seed survival (s) for each species from prior field and growth chamber experiments that focused on these two parameters (T. Martyn, unpublished data). Thirty mesh bags of 20 seeds of each species were placed in the field early in March 2017 and collected at the end of the growing season (October 2017). The number of germinated seeds were counted, and the ungerminated seeds were placed on germination paper and kept moist in growth chambers set to 19°C for 12 h of light and 9°C for 12 h of dark. After 2 weeks, the number of germinants were counted and added to field germinants to give the total number germinated. The remaining seeds were stained with tetrazolium dye and dissected to determine viability. The number of viable seeds was added to the total germinated to give the number that survived.

With these data, we fit the following models: G ~ Binomial(N,g) and S ~ Binomial(N,s), where G is the observed number of seeds that germinated, S is the observed number of seeds that survived and N is the number of seeds added to each replicate of the experiment.

For both g and s, we used Bayesian estimation with default uniform (0,1) priors.

Predicting coexistence with the invasion criterion

To predict coexistence using the invasion criterion, we first simulated each species as the ‘resident’ species to its steady‐state population size (i.e. distribution of equilibrium estimates) in the absence of heterospecific competitors (using Equation 2 below but setting N= 0). We propagated uncertainty in this simulation by performing 4500 simulation runs; for each run, we systematically drew each value of the posterior of intraspecific competition coefficients, germination and survival rates (which had a total of 4500 values), giving a distribution of resident equilibrium estimates. Overall seed count, N , (abundance) of species i, at the beginning of the growing season at time t, was modelled as:

The first term describes the carryover of seeds in the seed bank, where s is the annual survival of seeds in the seed bank and g is the fraction of seeds that germinate in a given year. The second term describes population growth and annual seed production. λ is the number of viable seeds produced in the absence of competition at time t and α describes the per‐capita intraspecific competitive effect between above‐ground stems (N , ). We simulated population growth from a starting point of two seeds, for 200 years so that each species’ population reached a steady‐state distribution (see Appendix 1, Figure S1).

Under the invasion criterion, stable coexistence occurs when each species as the invader, i, can invade the resident species, j, at the resident's steady‐state equilibrium, calculated using Equation 2 above, when i ≠ j and N is set to the resident's steady‐state population size. We again propagated uncertainty by performing multiple runs of the simulation to ‘invade’ species i into the population of species j, systematically drawing each posterior value for the invader's parameterisation and each estimate of the steady‐state population size of species j (4500 total). The same parameter values were kept at each time step within a run.

The ‘invader’ species’ growth rate is termed its ‘low‐density growth rate’ (LDGR), which determines how quickly a species can recover from low density. LDGR at time t is given by ln(N + /N. We calculated distributions of the LDGR for 14 species pairs (representing the pairs most commonly found to co‐occur in this natural system and that had adequate data to calculate each a and a pairwise combination). We started with an invader abundance of one seed, and the resident species at its steady‐state distribution (sampling a different value from the distribution each run) (see Appendix 1, Figure S2).

Predicting coexistence with the ratio approach

We also calculated distributions for both niche differences and average fitness differences given posterior values for α , α , g and s values.

Niche overlap, ρ, between any pair of species is defined as:

Niche differences (1 − ρ) are given by Chesson (2000), and Godoy and Levine (2014):

Thus, niche differences reflect the degree to which intraspecific competition outweighs interspecific competition. The average fitness difference between species j and i is defined as κ/κ, and is expressed as follows (Godoy & Levine, 2014):

The greater the ratio, κ/κ, the greater the fitness advantage of species j over i, and the faster species j excludes i in the absence of niche differences. The square root term describes the degree to which species i is more sensitive to intra‐ and interspecific competition than species j (i.e. ρ above). η describes the seeds produced per seed lost from the seed bank, and is given by:

Using this ratio approach, the condition for coexistence is ρ < κ/κ, where species j is the better competitor as per Godoy and Levine (2014). As above, we calculated niche and fitness differences for each of our posterior draws of intraspecific and interspecific competition coefficients, germination and survival rates (out of a total of 4500 values). The code for all models and analyses presented here are available at: https://github.com/cathbowler/Probability‐of‐coexistence‐

RESULTS

We found that accounting for uncertainty in model parameters yielded a probability of coexistence for 3 out of 14 species pairs based on both the invasion approach and the ratio approach. Additionally, accounting for uncertainty showed that priority effects (i.e. the order of species arrival) were likely to determine which species persists and which is excluded for one out of our 14 species pairs, where traditional approaches would instead have predicted competitive exclusion of one species by the other, without showing variation in the winner.

Invasion criterion

Traditional mean invasion criterion analyses predicted coexistence for only one out of 14 species pairs (Pentameris airoides with Podolepis canescens). For all other species pairs, with simulations based on the mean of each demographic parameter, one species was predicted to be competitively excluded by the other, as one species’ mean LDGR was positive while the other was negative (Figure 2—Case 3). Results were the same when using posterior median rather than mean values (see Appendix 1, Figure S3).

FIGURE 2

A framework for examining how incorporating uncertainty can alter species coexistence. When expanding MCT to consider underlying probability distributions of low‐density growth rate (LDGR), there are six categories of possible outcomes, or ‘cases’. In case 1, both LDGR distributions could be entirely negative, in which case there is a 100% probability that neither species can invade when the other species is present at equilibrium abundances. In this case, coexistence cannot occur, and rather priority effects (i.e. the order of species arrival) determine community composition. In case 2, both LDGR distributions could fall entirely in the positive range of LDGR values, in which case coexistence has a 100% probability of occurring, as both species can invade when the other is present at equilibrium. In case 3, one species’ distribution falls exclusively in the negative range, while the other falls exclusively in the positive range; the species in the positive range (species 2 in our example) is expected to outcompete the other with 100% probability. In cases 1–3, considering underlying probability distributions of LDGR does not alter our expectations of coexistence versus competitive exclusion. However, in cases 4–6, coexistence depends on LDGR probabilities. In cases 4 and 5, one species’ LDGR distribution overlaps zero, while the other does not. In case 4, we can calculate a probability of either coexistence or competitive exclusion, based on the density of the distribution that falls above zero for species 1 (in this example). In contrast, in case 5, there is a probability that either completive exclusion or priority effects occur, impacting coexistence outcomes. Finally, in case 6, there is a probability that competitive exclusion, priority effects or coexistence may occur

When considering the LDGR probability distribution for the same 14 species pairs, 3 additional species pairs were predicted to have some probability of coexistence. Of these distribution‐dependent cases, one pair was consistent with Figure 2—Case 4, where H. glutinosum, had a positive LDGR distributions, but the second species, Plantago debilis, overlapped zero, yielding a 47% probability of coexistence (Figure 4b). Another pair was consistent with Figure 2—Case 5, where one species, Pentameris airoides, had a negative LDGR distribution, but the second species, A. calendula, had a LDGR distribution overlapping zero, giving an 11% probability that priority effects would determine which species is competitively excluded and which persisted (Figure 4c). Two pairs fell into Figure 2—Case 6, with both species exhibiting distributions overlapping zero (H. glutinosum with Pentameris airoides and Pentameris airoides with Plantago debilis) (Figure 4a & d). Pentameris airoides with Plantago debilis had a 24% probability of coexistence (area under the curve where LDGR distributions overlapped in the positive region) and H. glutinosum with Pentameris airoides had a 9% probability of coexistence.

FIGURE 4

Examples of the low‐density growth rate (LDGR) distributions observed for mutually invading species pairs. The points correspond to each species’ LDGR calculated without the propagation of uncertainty (see Appendix 1, Figure S3). Coexistence is predicted when LDGR values for each species are above zero (vertical dashed line). The number in the top right corner of each plot is the probability of coexistence

Ratio approach

Using traditional mean‐only calculations, we again found that only one species pair (Pentameris airoides with Podolepis canescens) was predicted to coexist (Figure 5). For all remaining species pairs, fitness differences outweighed niche differences, leading to competitive exclusion of one of the two species. As with the invasion criterion, results were identical when using the posterior median versus mean values (data not shown).

FIGURE 5

Empirically observed coexistence patterns and correlations between stabilising niche differences and fitness differences across species pairs. Coloured dots in the main plot represent average fitness and stabilising niche differences for each pair of species, calculated using the mean values for all parameters, hence without propagated uncertainty. The green shaded region represents the region where the condition for coexistence is met (‘coexistence space’; ρ < κ/κ). Subplots show examples of the niche and fitness difference calculations where uncertainty has been propagated through the coexistence simulations, where black shapes in the subplots correspond to the average values (see Appendix 1, Figure S4)

When propagating uncertainty in demographic parameters, 3 out of 14 species pairs had substantial probabilities of coexistence (Figure 5). For the species pair that was predicted to coexist based on average demographic values, distributions gave consistent results; the probability of coexistence was 98% (Figure 5). For H. glutinosum with Plantago debilis and Pentameris airoides with Plantago debilis, we saw a 47% and 24% probability of coexistence, respectively. For most other species pairs, the distributions generally stayed outside coexistence space (less than 1% probability of coexistence), often producing a notable covariation trend, indicating strong correlation between log fitness differences and niche differences across posteriors.

DISCUSSION

We found that propagating uncertainty in coexistence modelling revealed the potential for different predicted outcomes of coexistence for 3 out of 14 species pairs and yielded a probability of priority effects determining which species will persist for an additional species pair.

Modern coexistence theory is built on the idea of variation (focusing on environmental) yielding windows of opportunities for coexistence. Despite the acknowledged importance of many sources of variation to coexistence, previous applications of this framework have tended to use demographically‐deterministic frameworks, where coexistence is calculated as a yes/no answer. These past applications minimise uncertainty from demographic sources, including demographic heterogeneity and stochasticity, and instead use either point‐estimate average demographic parameters (e.g. density‐independent growth rates and competition coefficients) or large‐scale manipulations of the environment (e.g. rainfall; Hallett et al., 2019). Our approach allowed for the calculation of coexistence probability distributions incorporating uncertainty in our demographic parameter estimates. This uncertainty, likely generated in large part by demographic heterogeneity and stochasticity, was directly included in our approach through the resampling of posterior values across different runs.

This approach is likely to be particularly useful in cases where point estimates fall close to the coexistence threshold (see Figures 3, 4, and 5). By excluding much of the uncertainty from environmental variation and observational error, this approach also points to the potential importance of underlying biological sources of complexity contributing to coexistence uncertainty.

FIGURE 3

The fraction of species pairs out of the 14 considered for whom coexistence is distribution dependent. In panel (a) mean only, we show the fraction of species pairs for which coexistence is predicted, versus where one species is competitively excluded, versus neither species being able to invade. In panel (b), based on distributions informed by a direct quantification of uncertainty in parameter estimates, we show the fraction of species pairs that fall into each category, but with the added category of the outcomes being distribution dependent because one or both of the low‐density growth rate distributions overlap zero; this occurs for four of our 14 species pairs

While uncertainty in demographic estimates is certainly driven in part by temporal and spatial variation in environmental conditions, the magnitude and impact of within‐population species’ performance uncertainty has historically been disregarded in the coexistence literature (Shoemaker et al., 2020). Using the traditional average value approach, our results were consistent with Wainwright et al. (2019), Godoy and Levine (2014), and Kraft et al. (2015) who all predicted low incidences of pairwise coexistence. Accounting for uncertainty yielded multiple potential outcomes of coexistence in our system for three of the fourteen species pairs. This could be due to unaccounted for sources of biological variation (including, but not limited to, variation among individuals and demographic stochasticity) structuring our model estimated uncertainty. If some individuals experience higher intrinsic fecundity or lower competition than average, it can lead to a higher LDGR, weaker fitness differences or greater niche differences. Hart et al. (2016) demonstrated mathematically that individual variation could promote coexistence under a specific set of circumstances, for instance, when there is a higher level of individual variation in otherwise inferior competitors. We found that species with the widest posterior probability (i.e. with the most uncertainty in the parameter estimate) for intrinsic fecundity were those that consistently outcompeted other species (see Appendix Figure S5). Additionally, species pairs at the edge of the coexistence region (Figure 5) or with one species whose LDGR was near 0 (Figure 4) often were most effected when accounting for uncertainty, as expected. However, importantly, it was not just location of the species pair means (i.e. close to the coexistence region or not), but also the amount of variability in LDGR distributions that influenced the probability of coexistence. For instance, M. minima and Pentameris airoides lay fairly far from the region of coexistence based on stabilising niche differences and fitness differences (Figure 5), but this species pair had a larger probability of coexistence than A. calendula and Pentameris airoides, whose point estimate fell much closer to the region of coexistence (Figure 5; Appendix 1, Figures S3 and S4).

Using a Bayesian approach to consider a posterior distribution of the abundance of the resident species and treating species as being composed of discrete individuals is important for more accurately predicting coexistence in natural systems (Schreiber et al., 2018). Future work should aim to disentangle multiple demographic sources of uncertainty, for example by directly testing the effects of demographic stochasticity versus heterogeneity, as the effect of stochasticity is expected to be stronger in small populations and hence may be most relevant in rare populations or those recovering from low density (Lande, 1993). Our empirical findings and the mathematical findings of Hart et al. (2016) demonstrate that the effects of this uncertainty (particularly between inferior and superior competitors) may be case and system dependent.

In addition to considering the potential influence of demographic sources of uncertainty, we also considered the potential impacts of observational error on our findings by collecting and counting individual seeds produced by each plant rather than extrapolating values from subsamples of seeds. We further rechecked the accuracy of our direct seed counts and found that observational error had a minimal influence on our estimated uncertainty and subsequent coexistence outcomes. Accounting for genotypes would have allowed for an even more robust study of the effects of demographic variation on the outcomes of coexistence but was beyond the scope of this study.

Improving real‐world applicability of coexistence modelling

Increasing real‐world applicability of coexistence models is of particular importance when considering the assembly of communities composed of native and exotic species. By incorporating uncertainty in the estimates of individual responses to neighbours, we found that one exotic species (Pentameris airoides) had a probability of coexistence with two native species (H. glutinosum and Podolepis debilis), when it would not have been predicted to coexist with these species using traditional models (see Appendix 1, Figure S3 & S4).

The invasion criterion itself has recently been critiqued by Pande et al. (2019) for its reliance on LDGR calculations as these require observed or simulated population growth through time. Pande et al. (2019) state that this approach fails to capture the effect of temporal random abundance variations (i.e. stochasticity) on the parameters of population persistence. For instance, an increase in random temporal or environmental variation could inherently lead to an increase in LDGR. We demonstrated here that uncertainty in parameter estimates can, in fact, impact LDGRs dramatically. We also demonstrated a simple way such uncertainty can be incorporated into coexistence modelling and empirical tests of coexistence theory. Determining the exact sources of uncertainty (e.g. temporal stochasticity) is an important next step in operationalising coexistence models for real‐world applications. Additionally, while here we focus on the effects of uncertainty in species’ density‐independent growth rates and interaction coefficients, future work could extend this framework to quantify probabilities of coexistence with temporal trends in uncertainty and environmental variability. We encourage further extensions, both theoretical and experimental, that examine the effects of variability using frameworks that incorporate both facilitation and competition between species. These extensions will require more targeted and long‐term experimental designs than this study used and modifications to the underlying model of species’ fecundity (e.g. Bimler et al. (2018)).

Despite accounting for uncertainty in our model estimates and predictions of coexistence, we still saw relatively lower occurrences of predicted coexistence for species pairs than is reflected in co‐occurrence patterns (in line with Kraft et al. (2015) and Wainwright et al. (2019)). This is likely due to the inherent issues with calculating coexistence in a single year and at the local neighbourhood scale, which are recognised to miss important stabilising factors like temporal and spatial heterogeneity (Hallett et al., 2019). Like most models of coexistence, we forced interactions to be competitive, but advancements allowing for facilitation may also further increase predictions of coexistence to match co‐occurrence patterns we see in nature.

Propagating uncertainty in coexistence analysis may reveal where priority effects are likely to determine which species will be competitively excluded. For instance, when both species’ LDGR distributions overlap zero, priority effects could determine which species will be competitively excluded (i.e. when the first species to arrive at a site either positively or negatively affect establishment, growth, or reproduction of species that arrive later (Grainger et al., 2019)) or which species is outcompeted and which persists. Determining if priority effects have come into play (i.e. determining which species will be competitively excluded) would require a mechanistic understanding of relevant species interactions, and likely experimental manipulations of the species’ environment; with our approach, however, we were able to identify scenarios where they are likely to occur (refer Figures 2 and 4).

CONCLUSION

Our Bayesian approach demonstrates a straightforward method for incorporating uncertainty into calculations of coexistence. Bayesian models have been extensively used in population ecology (Clark et al., 2010; Ellison, 2004; Vincenzi et al., 2014) but much less so in community ecology. Theory development in ecology involves a constant tension between effectively representing nature and keeping data requirements and complexity under control. The inclusion of uncertainty in model parameters increases the biological realism in our models and improves our ability to capture a more accurate prediction of coexistence based on quantifiable variability in nature.

AUTHORS CONTRIBUTION

CHB, MMM and LGS conceived this project; CHB and IRT designed and conducted the field study with guidance from MMM. CHB, LGS; and CWL created and analysed the models. CHB wrote the paper with input from LGS, CWL and MMM. MMM and LGS are co‐senior authors.

PEER REVIEW

The peer review history for this article is available at https://publons.com/publon/10.1111/ele.14011.

Supplementary Material

Click here for additional data file.