The spatial configuration of biotic interactions shapes coexistence-area relationships in an annual plant community
Study systemWe conducted our study in Caracoles Ranch, located in Doñana National Park (SW Spain 37° 04′ N, 6° 18′ W). The study area has a Mediterranean climate with mild winters and an average 50-year annual rainfall of 550–570 mm. Vegetation is dominated by annual grassland species, with no perennial species present. A subtle topographic gradient (slope 0.16%) is enough to generate vernal pools at the lower border of the ranch from winter (November–January) to spring (March–May), while upper parts do not get flooded except in exceptionally wet seasons. In our study, an extreme flooding event occurred during the growing season of 2018. A strong soil salinity–humidity gradient is structured along this topographic gradient.In September 2014, we established nine plots of 8.5 m × 8.5 m along a 1 km × 200 m area. Three of these nine plots were located in the upper part of the topographic gradient, three at the middle, and the last three at the lower part. The average distance between these three locations was 300 m and the average distance between plots within each location was 30 m (minimum distance 20 m). In addition, each plot was divided into 36 subplots of 1 m × 1 m with aisles of 0.5 m in between to allow access to subplots where measurements were taken (total of 324 subplots). This spatial design was established to parameterize population models including an intrinsic fecundity component and the effect of intra- and interspecific pairwise interactions. Specifically, the core of the observations involved measuring, for each focal individual, per germinant viable seed production as a function of the number and identity of neighbors within a radius of 7.5 cm including individuals of the same species. This radius is a standard distance used in previous studies to measure competitive interactions among annual plant species29,34, and has been validated to capture the outcome of competition interactions at larger scales (1 m²) under locally homogeneous environmental conditions35. From November 2014 to September 2019, we sampled 19 species present in the study area each year. We sampled one individual per subplot for widespread species and several individuals per subplot when species were rare (max. 324 individuals/species). This sampling design ensured that all species are balanced in terms of number of observations, and that we capture the full range of observed spatial interactions among species across the study area. Furthermore, we obtained independent estimates of seed survival and seed germination rates in 2016 (see17 for details on obtaining these rates). These 19 species belong to disparate taxonomic families and exhibit contrasted functional profiles along the growing season (Supplementary Table 1). The earliest species with small size and open flowers, such as C. fuscatum (Asteraceae), peak at beginning of the growing season (February), while late species with succulent leaves, such as S. soda (Amaranthaceae) and S. splendens (Amaranthaceae), grow during summer and peak at the end of the growing season (September-October). All these species represent up to 99% of plant cover in the study area.Estimating species interaction networks and intrinsic growth ratesWe estimated the effect of nearby individuals on individual fecundity via a Ricker model of population dynamics, which allowed us to estimate the strength of positive or negative interactions among pair of species, and therefore, to build a matrix of interactions among species. This approach has been previously applied to study annual plant systems under Mediterranean-type climates36, and it has also recently been shown to have several advantages compared to other formulations34. For example, this model implemented using a negative-binomial distribution for individual fecundities is more flexible in terms of modeling over-dispersion than a Poisson model, while maintaining predictions as positive integers. The model is of the form$${F}_{i,t}={lambda }_{i}{e}^{-({sum }^{}{alpha }_{i,j}{N}_{j,t})}$$
(1)
where ({lambda }_{i}) is the number of seeds produced by species i in the absence of interactions, ({alpha }_{i,j}) is the per capita effect of species j over species i (which can be positive or negative, thus allowing both competitive and facilitative effects), and ({N}_{j,t}) is the number of individuals of species j within 7.5 cm of the focal individual at timestep t. We fitted this model to the empirical data using Bayesian multilevel models with a negative-binomial distribution34. For model fitting, we used non-informative priors with MCMC settings of 5000 iterations (of which 2500 were warm-up) and 6 chains. The model was implemented using the brms R package37. The effect of changes in environmental conditions on species persistence can be phenomenologically evaluated by allowing models to vary in their estimates of species’ intrinsic growth rates and the reorganization of species interactions38. In our case, to evaluate the effect of environmental heterogeneity on species persistence (Question 1), we developed two complementary models. In both cases, we modeled the observed viable seed production per individual as a function of the identity and abundance of neighboring species. For the model assuming that plant species interact within a homogeneous environment across plots, we pooled together observations from the whole study area, and allowed the intercept and slope of the relationships to vary across years by including year as a random effect. Thus, the ({lambda }_{i}) and ({alpha }_{i,j}) values in Eq. 1 vary across years, but are homogeneous for the whole study area. We used the means from the obtained posterior distributions as estimates in the subsequent analyses. For the model that assumes that heterogeneous environments across space and time impact species population dynamics, we included an additional crossed random effect “plot”, thus obtaining spatially and temporally differentiated seed production in the absence of neighbors (({lambda }_{i})) and interaction coefficients (({alpha }_{i,j})). Importantly, our modeling approach does not evaluate the magnitude per se of the spatiotemporal variability in our system. It rather evaluates the response of plant species to changes in environmental conditions through their effects on vital rates and interaction coefficients (see39,40,41 for similar approaches). Likewise, this approach does not model the spatial dynamics of the community or spatially explicit mechanisms such as dispersal, but rather uses observed spatially explicit associations of individuals to infer their vital rates and interaction coefficients. In the following, we refer to the two developed models as “homogeneous parameterization” and “heterogeneous parameterization”, respectively (Fig. 1).The statistical methodology generates a posterior distribution of estimates for each parameter inferred, i.e., for each intrinsic fecundity rate (({lambda }_{i})) and interaction coefficient (({alpha }_{i,j})). These means, by definition, do not capture the full variability obtained with the statistical model, and may potentially be biased, especially for species pairs that have comparatively few observations. To ensure that our results were not biased by using the posterior mean as a fixed value in subsequent analyses, we replicated our analyses using random samples from the posterior distributions instead of the mean values. We generated 100 random draws from each parameterization and compared the obtained curves to the ones derived from the posterior means (Supplementary Note 1 and Supplementary Fig. 3).Finally, we assume that the study system presents a rich soil seed bank but we do not explicitly model its direct influence on driving the spatial pattern of species interactions or intrinsic vital rates: rather we use fixed field estimates of seed survival and germination rates in our modeling framework (see section “Analyzing species persistence”). This assumption implies that we cannot evaluate the contribution of a spatially or temporally varying seed bank to the shape of CARs and SARs, which remains an open question for future studies.Analyzing species persistenceTo analyze which species are predicted to persist and coexist with others in our system, we built communities based on the species’ spatial location. At the smallest spatial scale, given a community of S species observed in the field in a given plot and a given year, we calculated the persistence of each species within every community combination, from 2 species to S. Thus, we obtained for each species, plot, and year, two estimates of persistence, one from the homogeneous and another from the heterogeneous parameterization. To scale-up our predictions of species persistence at increasingly large areas, we aggregated species composition and persistence patterns from increasing numbers of plots. We consistently evaluated species persistence using a structuralist approach because prior work has shown it is compatible with the model used to estimate interaction coefficients (Eq. 1)14. Specifically, for a given community we first used the strength of sign of intra- and interspecific interactions to compute its feasibility domain (note that the structuralist approach can accommodate different signs in the interaction coefficients). Broadly speaking the feasibility domain is the structural analog of niche differences, and it represents the possible range of intrinsic species growth rates compatible with the persistence of individual species and of the entire community14. Indeed, the larger the feasibility domain, the larger the likelihood of species to persist. Yet, computing the feasibility domain does not tell us which species can persist. To obtain such information, we need to check whether the vector containing the observed differences in intrinsic growth rates between species fits within the limits of the feasibility domain. If so, then all species are predicted to coexist. If not, then one or more particular species is predicted to be excluded (see14 for a graphical representation).In order to quantify the feasibility of ecological communities, the intrinsic growth rates and interaction coefficients must be formulated according to a linear Lotka-Volterra model, or an equivalent formulation14. We transformed the parameters obtained from Eq. 1 to an equivalent Lotka-Volterra formulation with the following expression (Supplementary Note 2):$${r}_{i}={log}left(frac{1-(1-{g}_{i}){s}_{i}}{{g}_{i}}right)+{lambda }_{i}$$
(2)
where ({g}_{i}) is the seed germination rate of species i and ({s}_{i}) is its seed survival rate. Thus, we quantified the feasibility of our communities using the ri intrinsic growth rates from Eq. 2 and the ({alpha }_{i,j}) coefficients, which are not modified. For our main analyses, we used empirical estimates of seed survival and germination rates. We further explored the influence of these vital rates in the transformed intrinsic growth rates in Supplementary Note 2 (see also Supplementary Fig. 4 and Supplementary Table 4).The structuralist methodology further allowed us to dissect which specific configuration of species interactions is behind species persistence in our system (Question 2), among three possibilities: first, a given species may be able to persist by itself, and hinder the long-term persistence of neighboring species (category dominant). Second, pairs of species may be able to coexist through direct interactions (category coexistence of species pairs). The classic example of two-species coexistence is when the stabilizing effect of niche differences that arise because intraspecific competition exceeds interspecific competition overcome fitness differences41. Lastly, species may only be able to coexist in more complex communities (category multispecies coexistence)23, thanks to the effect of indirect interactions on increasing the feasible domain of the community14. A classic example of multispecies coexistence is a rock–paper–scissors configuration in which the three species coexist because no species is best at competing for all resources24,42. Because species may be predicted to persist under different configurations in a given community, we assigned their persistence category to the simplest community configuration. For instance, if we predicted that a three-species combination coexists as well as each of the three pairs separately, we assigned these species to the coexistence of pairs category26,43. Finally, if a species is not predicted to persist but it is observed in the system, we classify it as naturally transient, that is, it will tend to become locally extinct no matter what its surrounding community. In order to ascertain our classification of species as transient, we further analyzed whether these species shared ecological traits known to be common to transient species. In particular, a pervasive characteristic of transient species is their comparatively small population sizes. We explored the relationship between our classification as transient and species abundance through a logistic regression with logit link (supplementary Table 3).In addition to our main analyses, based on the structuralist approach, we explored the local stability44 of the observed communities, which evaluates their asymptotic response to infinitesimal perturbations, and thus provides a complementary view to the potential coexistence of the system (Supplementary Note 3).Species–area and coexistence–area relationshipsTo answer Question 3, we obtained standard SARs for each year, by calculating the average diversity observed when moving from 1 plot (72 m²) to 9 plots (650 m²) of our system. In the classification of SAR types proposed by Scheiner et al.45, the curves from our system are thus of type III-B, i.e., plots in a non-contiguous grid, with diversity values obtained using averages from all possible combinations of plots. Likewise, the yearly CARs were built taking the average number of coexisting species in each combination from 1 to 9 plots. In this case, a species was taken to persist in a given area if it was persisting alone or if it was part of at least one coexisting community within that area. We obtained CARs for the two parameterizations, i.e., assuming homogeneous interaction coefficients and individual fecundity throughout the study area, or explicitly including spatial variability in these terms. We fitted the CARs from Fig. 2 to power-law functions and obtained their associated parameters (Supplementary Table 2) using the mmSAR v1.0 package in R46. In the last step of the analyses, to evaluate the role of species identity in driving these empirical fits of CARs, we compared them to two complementary null models that reshuffle the strengths of per capita interactions between species pairs across the interaction matrix. In particular, as baseline we took the CARs from the homogeneous parameterization, in order to have a unique interaction strength value per species pair. In the first null model, and taking the inferred interaction matrix from a given plot and year, we redistributed the pairwise interaction coefficients randomly. That is, we fixed the number of species observed in a certain plot and year, as well as the structure of the interaction matrix, but randomized the magnitude of observed pairwise interactions (both intra and interspecific interactions) in that community. The second null model is similar, but keeping the diagonal coefficients of the interaction matrix, i.e., the intraspecific terms, fixed. While the first null model accounted for the effect of interspecific competitive responses, as well as self-limiting processes on driving CARs, the second null maintained self-limiting processes fixed by avoiding changes in the diagonal elements of our interaction matrices. We ran 100 replicates of each model for each year, and obtained the average CARs across replicates. All analyses were carried out in R v3.6.3, using packages tidyverse47 v.1.3.1 for data manipulation and visualization, and foreach48 v1.5.1 and doParallel49 v1.0.16 for parallelizing computationally intensive calculations.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More
