Intrinsic ecological dynamics drive biodiversity turnover in model metacommunities

Metacommunity model and asymptotic community assembly

We built a large set of model metacommunities (detailed in full in “Methods”) describing competitive dynamics within a single guild of species across a landscape. Each metacommunity consisted of a set of patches, or local communities, randomly placed in a square arena and linked by a spatial network. The dynamics of each population are governed by three processes: inter- and intraspecific interactions, heterogeneous responses to the environment and dispersal between adjacent patches (Fig. 1). Competition coefficients between species are drawn at random and the population dynamics within each patch are described by a Lotka-Volterra competition model. We control the level of environmental heterogeneity across the network directly by generating an intrinsic growth rate for each species at each patch from a random, spatially correlated distribution. To ensure any turnover is purely autonomous, we keep the environment fixed throughout simulations. Dispersal between neighbouring patches declines exponentially with distance between sites. This formulation allows precise and independent control of key properties of the metacommunity–the number of patches, the characteristic dispersal length and the heterogeneity of the environment.

To populate the model metacommunities, we iteratively introduced species with randomly generated intrinsic growth rates and interspecific interaction coefficients. Between successive regional invasions we simulated the model dynamics, and removed any species whose abundance fell below a threshold across the whole network. Through this assembly process and the eventual onset of ecological structural instability, both average local diversity, the number of species coexisting in a given patch, and regional diversity, the total number of species in the metacommunity, eventually saturate and then fluctuate around an equilibrium value—any introduction of a new species then leads on average to the extinction of one other species (Supplementary Fig. 1). In these intrinsically regulated metacommunities we then studied the phenomenology of autonomous community turnover in the absence of regional invasions or abiotic change.

In our metacommunity model, local community dynamics and therefore local limits on species richness depend on a combination of biotic and abiotic filtering (non-uniform responses of species to local conditions)^33,34,35 and immigration from adjacent patches, generating so called mass effects in the local community^36,37,38. Biotic filtering via interspecific competition is encoded in the interaction coefficients A_ij, while abiotic filtering occurs via the spatial variation of intrinsic growth rates R_ix. For simplicity, and since predator-prey dynamics are known to generate oscillations³⁹ through mechanisms distinct from those we report here, we restrict our analysis to competitive communities for which all ecological interactions are antagonistic. The off-diagonal elements of the interaction matrix A describe how one species i affects another species j. These are sampled independently from a discrete distribution, such that the interaction strength A_ij is set to a constant value in the range 0 to 1 (in most cases 0.5) with fixed probability (connectance, in most cases 0.5) and otherwise set to zero. Intraspecific competition coefficients A_ii are set to 1 for all species. This discrete distribution of the interaction terms was chosen for its relative efficiency. In the Supplementary discussion (and Supplementary Fig. 2) we show that outcomes remain unaffected when more complex distributions are modelled. Intrinsic growth rates R_ix are sampled from spatially correlated normal distributions with mean 1, autocorrelation length ϕ and variance σ² (Supplementary Fig. 3).

Dispersal is modelled via a spatial connectivity matrix with elements D_xy. The topology of the model metacommunity, expressed through D, is generated by sampling the spatial coordinates of N patches from a uniform distribution ({mathcal{U}}(0,sqrt{N})times {mathcal{U}}(0,sqrt{N})), i.e., an area of size N. Thus, under variation of the number of patches, the inter-patch distances remain fixed on average. Spatial connectivity is defined by linking these patches through a Gabriel graph⁴⁰, a planar graph generated by an algorithm that, on average, links each local community to four close neighbours⁴¹. Avoidance of direct long-distance dispersal and the sparsity of the resulting dispersal matrix permit the use of efficient numerical methods. The exponential dispersal kernel defining D_xy is tuned by the dispersal length ℓ, which is fixed for all species.

The dynamics of local population biomasses B_ix = B_ix(t) are modelled using a system of spatially coupled Lotka-Volterra (LV) equations that, in matrix notation, takes the form²³

with ∘ denoting element-wise multiplication. Hereafter this formalism is referred to as the Lotka-Volterra Metacommunity Model (LVMCM). Further technical details are provided in Methods and the Supplementary Discussion.

In order to numerically probe the impacts of ℓ, ϕ and σ² on the emergent temporal dynamics, we initially fixed N = 64 and varied each parameter through multiple orders of magnitude (Supplementary Fig. 4). In order to obtain a full characterisation of autonomous turnover in the computationally accessible spatial range (N ≤ 256), we then selected a parameter combination found to generate substantial fluctuations for further analysis. Thereafter we assembled metacommunities of 8–256 patches (Fig. 2a) until regional diversity limits were reached (with tenfold replication) and generated community time series of 10⁴ unit times from which the phenomenology of autonomous turnover could be explored in detail. We found no evidence to suggest that the phenomenology described below depends on this specific parameter combination. While future results may confirm or refute this, autonomous turnover arises over a wide range of parameters (Supplementary Fig. 4) and as such the phenomenon is robust.

Autonomous turnover in model metacommunities

For small (N ≤ 8) metacommunities assembled to regional diversity limits, populations attain equilibria, i.e., converge to fixed points, implying the absence of autonomous turnover²³. With increasing metacommunity size N, however, we observe the emergence of persistent population dynamics (Supplementary Fig. 5 and external video) that can produce substantial turnover in local community composition. This autonomous turnover can be represented through Bray-Curtis⁴² (BC) dissimilarity matrices comparing local community composition through time (Fig. 2b), and quantified by the number of compositional clusters detected in such matrices using hierarchical cluster analysis (Fig. 2a, c).

At intermediate spatial scales (Fig. 2, 16 ≤ N ≤ 32) we often find oscillatory dynamics, which can be perfectly periodic or slightly irregular. With increasing oscillation amplitude, these lead to persistent turnover dynamics where local communities repeatedly transition between a small number of distinct compositional clusters (represented in Fig. 2 by stripes of high pairwise BC dissimilarity spanning large temporal ranges). At even larger scales (N ≥ 64) this compositional coherence begins to break down, and for very large metacommunities (N ≥ 128) autonomous dynamics drive continuous acyclic change in community composition. The number of compositional clusters detected over time typically varies within a given metacommunity (Fig. 2a node colour), however we find a clear increase in the average number of compositional clusters, i.e., an increase in turnover, with increasing total metacommunity size (Fig. 2c).

Metacommunities in which the boundaries of species ranges along environmental gradients are clumped are termed Clementsian, while those for which range limits are independently distributed are  referred to as Gleasonian⁴³. We consider the block structure of the temporal dissimilarity matrix at intermediate N to represent a form of Clementsian temporal turnover, characterised by sudden significant shifts in community composition. Metacommunity models similar to ours have been found to generate such patterns along spatial gradients⁴⁴, potentially via an analogous mechanism⁴⁵. Large, diverse model metacommunities manifest Gleasonian temporal turnover. In such cases, species colonisations and extirpations are largely independent and temporal occupancies predominantly uncorrelated, such that compositional change is continuous, rarely, if ever, reverting to the same state.

Mechanistic explanation of autonomous turnover

Surprisingly, the onset and increasing complexity of autonomous turnover as system size N increases (Fig. 2) can be understood as a consequence of local community dynamics alone. To explain this, we first recall relevant theoretical results for isolated LV communities. Then we demonstrate that, in presence of weak propagule pressure, these results imply local community turnover dynamics, controlled by the richness of potential invaders, that closely mirror the dependence on system size seen in full LV metacommunities.

Application of methods from statistical mechanics to models of large isolated LV communities with random interactions revealed that such models exhibit qualitatively distinct phases^46,47,48. If the number of modelled species, S, interpreted as species pool size, lies below some threshold value determined by the distribution of interaction strengths (Supplementary Fig. 7), these models exhibit a unique linearly stable equilibrium (Unique Fixed Point phase, UFP). Some species may go extinct, but the majority persists⁴⁸. When pool size S exceeds this threshold, there appear to be no more linearly stable equilibrium configurations. Any community formed by a selection from the S species is either unfeasible (there is no equilibrium with all species present), intrinsically linearly unstable, or invadable by at least one of the excluded species. This has been called the multiple attractor (MA) phase⁴⁷. However, the implied notion that this part of the phase space is in fact characterised by multiple stable equilibria may be incorrect.

Population dynamical models with many species have been shown to easily exhibit attractors called stable heteroclinic networks⁴⁹, which are characterised by dynamics in which the system bounces around between several unstable equilibria, each corresponding to a different composition of the extant community, implying indefinite, autonomous community turnover (Fig. 3, red line). As these attractors are approached, models exhibit increasingly long intermittent phases of slow dynamics, which, when numerically simulated, can give the impression that the system eventually reaches one of several ‘stable’ equilibria, suggesting that turnover comes to a halt. We demonstrate in the Supplementary discussion that the MA phase of isolated LV models is in fact characterised by such stable heteroclinic networks (Supplementary Figs. 8 and 9). Note, we retain the MA terminology here because the underlying complete heteroclinic networks, interpreted as a directed graph^50,51 (Fig. 3, inset), might have multiple components that are mutually unreachable through dynamic transitions⁵², each representing a different attractor.

If one now adds to such isolated LV models terms representing weak propagule pressure for all S species (Supplementary Eq. (2)), dynamically equivalent to mass effects occurring in the full metacommunity model (Eq. (1)), then none of the S species can entirely go extinct. The weak influx of biomass drives community states away from the unstable equilibria representing coexistence of subsets of the S species and the heteroclinic network connecting them (blue line in Fig. 3). Typically, system dynamics then still follow trajectories closely tracking the original heteroclinic networks (Fig. 3), but now without requiring boundless time to transition from the vicinity of one equilibrium to the next.

The nature and complexity of the resulting population dynamics depend on the size and complexity of the underlying heteroclinic network, and both increase with pool size S. In simulations (Supplementary Fig. 10) we find that, as S increases, LV models with weak propagule pressure pass through the same sequence of states as we documented for LVMCM metacommunities in Fig. 2: equilibria, oscillatory population dynamics, Clementsian and finally Gleasonian temporal turnover.

Above we introduced the number of clusters detected in Bray-Curtis dissimilarity matrices of fixed time series length as a means of quantifying the approximate number of equilibria visited during local community turnover. As shown in Fig. 4a, b, this number increases in LV models with S in a manner strikingly similar to its increase in the LVMCM with the number of species present in the ecological neighbourhood of a given patch. Thus, dynamics within a patch are controlled not by N directly but rather by neighbourhood species richness. For a given neighbourhood, species richness depends on the number of connected patches, the total area and therefore total abiotic heterogeneity encompassed, and the connectivity, all of which can vary substantially within a metacommunity of a given size N. As illustrated in Fig. 4b, there is a tendency for neighbourhood richness to be larger in larger metacommunities, leading indirectly to the dependence of metacommunity dynamics on N seen in Fig. 2.

There is thus a close correspondence between dynamically isolated LV models and LVMCM metacommunities, in the sequence of dynamic states as propagule richness increases and in the resulting complexity of dynamics quantified by counting compositional clusters. This suggests that underlying heteroclinic networks, which are revealed by adding propagule pressure in isolated communities, explain the complex dynamics seen in LVMCM metacommunities.

For the isolated LV community, the threshold beyond which autonomous turnover is detected (>1 compositional cluster) occurs at a pools size of around S = 35 species, consistent with the theoretical prediction⁴⁷ of the transition between the UFP and MA phases (Supplementary Discussion). Close inspection of this threshold reveals an important and hitherto unreported relationship between the transition into the MA phase and local ecological limits set by the onset of ecological structural instability, which is known to regulate species richness in LV systems subject to external invasion pressure^23,24: in the Supplementary Discussion we show that the boundary between the UFP and MA phases⁴⁷ coincides precisely with the onset of structural instability²⁴ (Supplementary Eqs. (3)–(9)).

For LVMCM metacommunities, this relationship (demonstrated analytically in the Supplementary Discussion) is numerically confirmed in Fig. 5. During assembly, local species richness increases until it reaches the limit imposed by local structural instability. Further assembly occurs via the “regionalisation” of the biota⁵³—a collapse in average range sizes²³ and associated increase in spatial beta diversity—until regional diversity limits are reached²³. The emergence of autonomous turnover coincides with the onset of species saturation at the local scale. Autonomous turnover can therefore serve as an indirect indication of intrinsic biodiversity regulation via local structural instability in complex communities.

Thus, we have shown that propagule pressure perturbs local communities away from unstable equilibria and drives compositional change. In order to invade, however, species need to be capable of passing through biotic and abiotic filters^33,34,35. We would expect, therefore, that turnover would be suppressed in highly heterogeneous or poorly connected environments where mass effects are weak. Indeed, by manipulating the autocorrelation length ϕ and variance σ² of the abiotic filter represented by the matrix R and the characteristic dispersal length ℓ, we observe a sharp drop-off in temporal turnover in parameter regimes that maximise between-patch community dissimilarity (short environmental correlation or dispersal lengths, Supplementary Fig. 11). Thus, we conclude that it is not species richness or spatial dissimilarity per se that best predict temporal turnover, but the size of the pool of species with positive invasion fitness, i.e., those not repelled by the combined effects of biotic and abiotic filters.

The macroecology of autonomous turnover

We find good correspondence between temporal and spatio-temporal biodiversity patterns emerging in model metacommunities in the absence of external abiotic change and in empirical data (Fig. 6), with quantitative characteristics lying within the ranges observed in natural ecosystems.

Temporal occupancy

The proportion of time in which species occupy a community tends to have a bi-modal empirical distribution^54,55,56 (Fig. 6a). The distribution we found in simulations (Fig. 6e) closely matches the empirical pattern.

Community structure

Temporal turnover has been posited to play a stabilising role in the maintenance of community structure^57,58. In an estuarine fish community⁵⁹, for example, species richness (Fig. 6b) and the distribution of abundances were remarkably robust despite changes in population biomasses by multiple orders of magnitude. In model metacommunities with autonomous turnover we found, likewise, that local species richness exhibited only small fluctuations around the steady-state mean (Fig. 6f, three random local communities shown) and that the macroscopic structure of the community was largely time invariant (Supplementary Fig. 12). In the light of our results, we propose the absence of temporal change in community properties such as richness or the abundance distribution despite potentially large fluctuations in population abundances⁵⁹ as indicative of autonomous compositional turnover.

The species-time-area-relation, STAR

The species-time-relation (STR), typically fit by a power law of the form S ∝ T^w 12,60,61, describes how observed species richness increases with observation time T. The exponent w of the STR has been found to be consistent across taxonomic groups and ecosystems^12,13,62, indicative of some general population dynamical mechanism. However, the exponent of the STR decreases with increasing sampling area¹², and the exponent of the empirical Species Area Relation (SAR) (S ∝ A^z) consistently decreases with increasing sampling duration¹² (Fig. 6c, d). We tested for these patterns in a large simulated metacommunity with N = 256 patches by computing the species-time-area-relation (STAR) for nested subdomains and variable temporal sampling windows (see “Methods”). We observed exponents of the nested SAR in the range z = 0.02–0.44 and for the STR a range w = 0.01–0.44 (Supplementary Fig. 13). We also found a clear decrease in the rate of species accumulation in time as a function of sample area and vice-versa (Fig. 6g, h), consistent with the empirical observations. Meta-analyses of these patterns in nature have reported exponents which are remarkably consistent, with z typically in the range 0.1–0.3⁶³, and w typically in the range 0.2–0.4¹³, in both cases largely independent of location or taxonomic group¹³.

Thus, the distribution of temporal occupancy, the time invariance of key macroecological structures and the STAR in our model metacommunities match observed patterns. This evidence suggests that such autonomous dynamics cannot be ruled out as an important driver of temporal compositional change in natural ecosystems.

Turnover rate in simulated metacommunities

How do the turnover rates that we find in our model compare with those observed? Our current analytic understanding of autonomous turnover is insufficient for estimating the rates directly from parameters, but the simulation results provide some indication of the expected order of magnitude, that can be compared with observations. Key for such a comparison is the fact that, because the elements of R are 1 on average, the time required for an isolated single population to reach carrying capacity is ({mathcal{O}}(1)) unit times. Supplementary Fig. 12b suggests that transitions between community states occur at the scale of around 10–50 unit times. This gives a holistic, rule-of-thumb estimate for the expected rate of autonomous turnover, depending on the typical reproductive rates of the guild of interest. In the case of macroinvertebrates, for example, the time required for populations to saturate in population biomass could be of the order of a month or less. By our rule of thumb, this would mean that autonomous community turnover would occur on a timescale of years. In contrast, for slow growing species like trees, where monoculture stands can take decades to reach maximum population biomass, the predicted timescale for autonomous turnover would be on the order of centuries or more. Indeed, macroinvertebrate communities have been observed switching between community configurations with a period of a few years^64,65, while the proportional abundance of tree pollen and tree fern spores fluctuates in rain forest bog deposits with a period of the order of 10³ years⁶⁶—suggesting that the predicted autonomous turnover rates are biologically plausible.

Conclusions

Current understanding of the mechanisms driving temporal turnover in ecological communities is predominantly built upon phenomenological studies of observed patterns^2,67,68,69 and is unquestionably incomplete^10,59. That temporal turnover can be driven by external forces—e.g., seasonal or long term climate change, direct anthropogenic pressures—is indisputable. A vitally important question is, however, how much empirically observed compositional change is actually due to such forcing. Recent landmark analyses of temporal patterns in biodiversity have detected no systematic change in species richness or structure in natural communities, despite rates of compositional turnover greater than predicted by stochastic null models^1,70,71,72. Here we have shown that empirically realistic turnover in model metacommunities can occur via precisely the same mechanism as that responsible for regulating species richness at the local scale. While the processes regulating diversity in natural communities remain insufficiently understood, our theoretical work suggests local structural instability may explain these empirical observations in a unified and parsimonious way. Therefore, we advocate for the application of null models of metacommunity dynamics that account for natural turnover in ecological status assessments and predictions based on ancestral baselines. Future work will involve fitting the model described here to observations by estimating abiotic and biotic parameters from empirical datasets. In the Supplementary Discussion we show how different combinations of parameters lead to different quantitative outcomes (Supplementary Fig. 4), likely representing different types of empirical metacommunities. Understanding where in this parameter space natural systems exist may provide the foundation for a quantitative null model, a baseline expectation of turnover against which observations can be compared.

Our simulations revealed a qualitative transition from “small” metacommunities, where autonomous turnover is absent or minimal, to “large” metacommunities with pronounced autonomous turnover (Fig. 2). The precise location of the transition between these cases depends on details such as dispersal traits, the ecological interaction network, and environmental gradients (Supplementary Fig. 4). Taking, for simplicity, regional species richness as a measure of metacommunity size suggests that both ‘small’ and ‘large’ communities in this sense are realised in nature. In our simulations, the smallest metacommunities sustain 10s of species, while the largest have a regional diversity of the order 10³, which is not large comparable to the number of tree species in just 0.25 km² of tropical rainforest (1100–1200 in Borneo and Ecuador⁷³) or of macroinvertebrates in the UK (>32,000⁷⁴). Within the ‘small’ category, where autonomous turnover is absent, we would therefore expect to be, e.g., communities of marine mammals or large fish, where just a few species interact over ranges that can extend across entire climatic niches, implying that the effective number of independent “patches” is small and providing few opportunities for colonisation by species from neighbouring communities. Likely to belong to the ‘large’ category are communities of organisms that occur in high diversity with range sizes that are small compared to climatic niches, such as macroinvertebrates. For these, autonomous turnover of local communities can plausibly be expected based on our findings. Empirically distinguishing between these two cases for different guilds will be an important task for the future.

For metacommunities of intermediate spatial extent, autonomous turnover is characterised by sharp transitions between cohesive states at the local scale. To date, few empirical analyses have reported such coherence in temporal turnover, perhaps because the taxonomic and temporal resolution required to detect such patterns is not yet widely available. Developments in biomonitoring technologies⁷⁵ are likely to reveal a variety of previously undetected ecological dynamics, however and by combining high resolution temporal sampling and metagenetic analysis of community composition, a recent study demonstrated cohesive but short-lived community cohorts in coastal plankton⁷⁶. Such Clementsian temporal turnover may offer a useful signal of autonomous compositional change in real systems.

Thus, overcoming previous computational limits to the study of complex metacommunities^11,77, we have discovered the existence of two distinct phases of metacommunity ecology—one characterised by weak or absent autonomous turnover, the other by continuous compositional change even in the absence of external drivers. By synthesising a wide range of established ecological theory^{11,23,24,47,48,49}, we have heuristically explained these phases. Our explanation implies that autonomous turnover requires little more than a diverse neighbourhood of potential invaders, a weak immigration pressure, and a complex network of interactions between co-existing species.

Source: Ecology - nature.com