Heterogeneity in the isolation of patches may be essential for the action of metacommunity mechanisms

Borthagaray, Ana I.; Cunillera-Montcusí, David; Bou, Jordi; Tornero, Irene; Boix, Dani; Anton-Pardo, Maria; Ortiz, Esteban; Mehner, Thomas; Quintana, Xavier D.; Gascón, Stéphanie; Arim, Matías

doi:10.3389/fevo.2023.1125607

ORIGINAL RESEARCH article

Front. Ecol. Evol., 31 March 2023

Sec. Population, Community, and Ecosystem Dynamics

Volume 11 - 2023 | https://doi.org/10.3389/fevo.2023.1125607

This article is part of the Research Topic Patterns and Processes in Ecological Networks over Space View all 5 articles

Heterogeneity in the isolation of patches may be essential for the action of metacommunity mechanisms

$\nAna I. Borthagaray$ Ana I. Borthagaray¹

David Cunillera-Montcusí^1,2,3

Jordi Bou⁴

Irene Tornero^3,5

Dani Boix³

Maria Anton-Pardo^3,6^†

Esteban Ortiz¹

Thomas Mehner⁷

Xavier D. Quintana³

Stéphanie Gascón³

Matías Arim^1,8^*

¹Departamento de Ecología y Gestión Ambiental, Centro Universitario Regional del Este (CURE), Universidad de la República, Maldonado, Uruguay
²FEHM-Lab, Section of Ecology, Department of Evolutionary Biology, Ecology and Environmental Sciences, University of Barcelona, Barcelona, Spain
³GRECO, Institute of Aquatic Ecology, University of Girona, Girona, Spain
⁴LAGP-Flora and Vegetation, Institute of the Environment, University of Girona, Girona, Spain
⁵Institute of Aquatic Ecology, Centre for Ecological Research, Budapest, Hungary
⁶Cavanilles Institute of Biodiversity and Evolutionary Biology, University of Valencia, Valencia, Spain
⁷Leibniz Institute of Freshwater Ecology and Inland Fisheries (IGB), Berlin, Germany
⁸Centro Interdisciplinario de Ciencia de Datos y Aprendizaje Automático (CICADA), Universidad de la República, Montevideo, Uruguay

The spatial isolation gradient of communities and the gradient in the species dispersal ability are recognized as determinants of biodiversity in metacommunities. In spite of this, mean field models, spatially explicit models, and experiments were mainly focused on idealized spatial arrangements of communities leaving aside the combining role of dispersal and isolation gradients in metacommunity processes. Consequently, we have an incipient understanding of the role of the real spatial arrangement of communities on biodiversity patterns. We focus on six metacommunities for which confident information about the spatial arrangement of water bodies is available. Using coalescent metacommunity models and null models that randomize the location of water bodies, we estimated the potential effect of the landscape on biodiversity and its dependence on species dispersal ability. At extremely low or high dispersal abilities, the location of ponds does not influence diversity because different communities are equally affected by the low or high incoming dispersal. At intermediate dispersal abilities, peripheral communities present a much lower richness and higher beta diversity than central communities. Moreover, metacommunities from real landscapes host more biodiversity than randomized landscapes, a result that is determined by the heterogeneity in the geographic isolation of communities. In a dispersal gradient, mass effects systematically increase the local richness and decrease beta diversity. However, the spatial arrangement of patches only has a large importance in metacommunity processes at intermediate dispersal abilities, which ensures access to central locations but limits dispersal in isolated communities. The ongoing reduction in spatial extent and simplification of the landscape may consequently undermine the metacommunity processes that support biodiversity, something that should be explicitly considered in preserving and restoring strategies.

Introduction

Dispersal among communities is increasingly recognized as a main determinant of biodiversity structure through a range of mechanisms captured in the metacommunity theory (Cadotte, 2006; Leibold and Chase, 2018; Lu, 2021; Cathey and Brown, 2022). In general, metacommunity mechanisms become important when dispersal is strong enough to affect local dynamics but not too high to homogenize the system (Leibold et al., 2004; Loreau, 2010; Leibold and Chase, 2018). At intermediate dispersal levels, both local processes and regional dispersal are determining community assembly, none of them being strong enough to override the effect of the other, a scenario that promotes species coexistence by a balance between patch dynamics, species sorting, and mass effects (Thompson et al., 2017; Suzuki and Economo, 2021). This generates a hump-shaped relationship between the local organization of communities and the dispersal rate (Mouquet and Loreau, 2003; Leibold and Chase, 2018; Suzuki and Economo, 2021). Congruently, local richness, functional diversity, food web complexity, and species co-occurrence have been shown to peak at intermediate dispersal rates among communities (Pillai et al., 2012; Bender et al., 2016; Mougi and Kondoh, 2016; Leibold and Chase, 2018; Suzuki and Economo, 2021). Dispersal is also a force that can homogenize the structure of different communities, causing a decrease in beta diversity (Loreau, 2010; Gianuca et al., 2017; Leibold and Chase, 2018). These trends in alpha and beta diversities were consistently reported in both theoretical and empirical studies (Grainger and Gilbert, 2016; Leibold and Chase, 2018).

The metacommunity theory mostly stands on the analyses of spatially implicit systems, in which all patches are equally connected to all others (but see Economo and Keitt, 2010; Ai et al., 2013; Suzuki and Economo, 2021), and species have equal dispersal rates (Borthagaray et al., 2015a,b; Grainger and Gilbert, 2016). Furthermore, spatially explicit models and experiments usually focus on the idealized spatial arrangement of local communities—random locations, chains, stars, and grids (Economo and Keitt, 2008, 2010; Borthagaray et al., 2014; Arim et al., 2016; Grainger and Gilbert, 2016; Häussler et al., 2021; Suzuki and Economo, 2021). Therefore, it is assumed that dispersal gradients equally affect all species and communities (Grainger and Gilbert, 2016; Suzuki and Economo, 2021). However, in real species pools, there are large variations in species' abilities to disperse (De Bie et al., 2012; Jones et al., 2015), and real landscape communities are not uniformly distant from each other (Economo and Keitt, 2010; Ai et al., 2013; Grilli et al., 2015; McManus et al., 2021; Suzuki and Economo, 2021). In addition, the relative location of patches in the landscape determines a significant heterogeneity in the geographic distance among them that affects source–sink dynamics (Hanski, 1999; Muneepeerakul et al., 2008; Economo and Keitt, 2010; Grilli et al., 2015; McManus et al., 2021; Suzuki and Economo, 2021). This gradient in the geographic isolation of patches was associated with trends in local richness (Urban and Keitt, 2001; Economo and Keitt, 2010; Borthagaray et al., 2015a, 2020; Henriques-Silva et al., 2019), numerical abundance (Cunillera-Montcus et al., 2020b), functional diversity (Arellano-Rivas et al., 2016; Harvey and Altermatt, 2019; Borthagaray et al., 2020), adaptive response of species to environmental change (McManus et al., 2021), spatial and temporal turnover (Carrara et al., 2012; Lu et al., 2019; Jabot et al., 2020), size structure (Castle et al., 2011; Borthagaray et al., 2012), trophic structure (Chase and Shulman, 2009; Arim et al., 2016; Ryser et al., 2019), and ecosystem functioning (Piñeiro-Guerra et al., 2014; Maureaud et al., 2019; Gonzalez et al., 2020).

Consequently, a mass of evidence has accumulated indicating that both the geographic isolation gradient of patches and the gradient in species dispersal abilities are the main determinants of biodiversity patterns in metacommunities (Fortin et al., 2021; Suzuki and Economo, 2021; Cathey and Brown, 2022). Furthermore, their effects are interrelated. At low dispersal rates, local communities are essentially disconnected from each other, no matter if they are in a central or isolated location (e.g., Phillipsen and Lytle, 2013; Cañedo-Argüelles et al., 2015; Castillo-Escriv et al., 2017; Cunillera-Montcus et al., 2021). At high dispersal rates, geographically distant patches are easily reached, and all communities are equally affected by dispersal (e.g., Sarremejane et al., 2017). Thus, it is at intermediate dispersal rates that the geographic isolation gradient of communities may play the most important role in determining biodiversity (e.g., Fortuna et al., 2006; Economo and Keitt, 2010; Cañedo-Argüelles et al., 2015). Notably, we have no clear understanding of how the spatial structure of a real landscape interacting with species dispersal abilities may determine metacommunity biodiversity at intermediate dispersal rates (Hillebrand et al., 2018; Thompson et al., 2020; Bauer et al., 2021; Fortin et al., 2021; Suzuki and Economo, 2021).

In this study, we focus on six metacommunities of ponds for which reliable information about their spatial arrangement is available. Using coalescent metacommunity models in these real landscapes, we analyzed the following: first, the interaction between the geographic isolation gradient and the dispersal ability gradient as the determinant of metacommunity diversity, and second, the importance of the real spatial arrangement of ponds in the landscape on diversity patterns by means of contrasting metacommunity diversity in real vs. randomized landscapes.

Methods

A total of six metacommunities of ponds from Europe (5) and South America (1) were analyzed here as model systems (Figure 1). These landscapes were selected as model systems because they combine reliable information about the spatial location of local communities, the geographic or microclimatic barrier that constrains dispersal within the metacommunity, and previous studies support the role of metacommunity processes determining their diversity patterns (Arim et al., 2010, 2011; Ballón et al., 2016; Tornero et al., 2016, 2018; Cunillera-Montcus et al., 2020a; Rodriguez-Tricot and Arim, 2020). The set of metacommunities considered covers a wide range of spatial arrangements of water bodies, including metacommunities from 61 to up to 542 small waterbodies and total areas from 2 to 560 km² (see Supplementary Table S1). Information about waterbody areas was not available for all communities, and consequently, we focus here only on their spatial arrangement. The range of dispersal abilities that we used in this study for the species inhabiting these systems (i.e., from 10 to 3,000 meters) is within the range of geographic distances observed among ponds (Arim et al., 2011; Borthagaray et al., 2015a,b; Tornero et al., 2018; Cunillera-Montcus et al., 2020a, 2021; Sarremejane et al., 2020). Hereafter, we refer to these systems as pondscapes.

FIGURE 1

Figure 1. Pondscapes analyzed in this study. Five European pondscapes [Albera, Clots de Guils, Empordà Wetlands (PNAE) and Vila Nova de Milfontes] and one Uruguayan (Rocha). Minimum spanning trees (the network that connects all nodes with a minimum of links) are presented for visualizing the spatial arrangement of waterbodies.

To estimate the geographic isolation gradient for each pondscape, we define a fully connected network, in which the links between ponds were weighted by the inverse of the Euclidean distance between them, that is, $w_{i j} = (\frac{1}{d_{i j}})$ . Based on these weighted networks, degree centrality was estimated as the sum of the weights of the links between one pond and all ponds in the network (Newman C. M., 2018). A weighted degree is a metric of centrality in networks that well represents the relative geographic isolation of local communities in the pondscape (Urban et al., 2009). Note that pond isolation is inverse to pond centrality. In this sense, the larger the degree of a pond, the larger the connectedness of this pond to other ponds in the community. Because of this, we opted to use a degree centrality because first, it has a direct interpretation—it is proportional to the average inverse distance between a focal pond and all others—and second, it provides similar information to alternative metrics such as closeness or eigenvector centrality (Newman M., 2018). The effect of geographic isolation on species dispersal rate will depend on the dispersal ability of the species—see in the next section (Keitt et al., 1997; Bunn et al., 2000; Urban and Keitt, 2001).

Metacommunity assembly

Coalescent and lottery models of metacommunities can capture the effect on diversity patterns of the balance between dispersal, drift, and selection by environmental filters and species interactions (Hubbell, 2001; Rosindell et al., 2011; Borthagaray et al., 2014, 2023; Worm and Tittensor, 2018; Cunillera-Montcus et al., 2021). In this sense, these models provide a methodological framework coherent with the community ecology theory (Vellend, 2010, 2016). These models are particularly suitable for detecting the effect of the landscape on dispersal patterns and biodiversity (Hubbell, 2001; Rosindell et al., 2011; Borthagaray et al., 2014, 2023; Worm and Tittensor, 2018; Cunillera-Montcus et al., 2021). While predictions of biodiversity patterns from coalescent and lottery models are more congruent, coalescent models are much faster, ensuring the analysis of a wider range, or scenarios (Munoz et al., 2018; Worm and Tittensor, 2018). Consequently, here, we use a spatially explicit coalescent model in which community dispersal is affected by their spatial location, and species compete for a fixed amount of resources in local communities—represented by the number of individuals that can be hosted in each spatial patch (see Hubbell, 2001; Economo and Keitt, 2008, 2010; Borthagaray et al., 2014, 2023; Worm and Tittensor, 2018; Cunillera-Montcus et al., 2021).

Each metacommunity was assembled with a coalescent dynamic in which one individual was randomly selected from a metacommunity species pool for colonizing each local community (see Figure 2). Then, the community is filled with 100 individuals sequentially sampled as follows: (1) from the metacommunity species pool, (2) from neighboring communities, or (3) from the same local community with rates m.pool, m.neighbor, and 1 − (m.pool + m.neighbor), respectively (Worm and Tittensor, 2018; see also Munoz et al., 2018; Cunillera-Montcus et al., 2021; Borthagaray et al., 2023). m.neighbor represents the migration rate from all other communities in the landscape (P). Migration from other local communities ( $m . n e i g h b o r = \sum_{i = 1}^{P} m_{i j}$ ) is based on a dispersal kernel between two communities. This dispersal kernel was modeled with an exponential decay function as follows: $m_{i j} = m . m a x * e^{- {b^{*} d}_{i j}}$ where d_ij is the geographic distance between communities i and j in meters, m.max is the migration between ponds 0 meter apart (m.max = 1), and b is a dispersal parameter that describes the decay in dispersal with distance d_ij. This parameter b was estimated considering the distance between communities at which migration decays to half of its maximum value, a distance defined as d₅₀, i.e., the parameter b is proportional to the inverse of the average dispersal distance (Hanski, 1999), that is, when m_ij = 0.5, from the exponential decay function, $b = \frac{- l o g (0.5)}{d_{50}}$ . Herein, we use d₅₀ as the parameter that resumes the organism's dispersal ability. This dispersal ability emerges from the interaction between species traits and the resistance to dispersal imposed by the environment between ponds (e.g., predation risk, wind, temperature, and topography). In this sense, d₅₀ is not only determined by geographic distance but also captures the species' biology and environmental conditions. Then, a community-by-community migration matrix M was estimated by applying the exponential decay function to each element of the geographic distance matrix (D): M = m.max ∗ e^−b*D. The element m_ij of this matrix indicates the dispersal rate from community i to community j. Finally, each column of matrix M plus migration from the species pool, m.pool, represents the set of parameters that determine the total source of individuals to recruitment at each community, i.e., the potential incoming migration rate. These parameters are translated to probabilities standardizing in each column by the sum of its elements plus m.pool–i.e., $m_{i j} s t a n d a r d i z e d = \frac{m_{i j}}{(m_{p o o l} + \sum_{i = 1}^{P} m_{i j})}$ . After this standardization, the sum of the non-diagonal elements of each column represents recruitment probability from other communities, m.neighbour, and diagonal elements represent self-recruitments, 1 − (m.pool + ∑m.neighbour). After standardization, even dispersal rates of 1 (e.g., when zero distance between ponds) do not imply that individuals have to arrive from this only source. It should be noted that in this and previous coalescent models, it is implicitly assumed that there are no biases to internal or external recruitment, more than those originated by the dispersal processes. In all simulations, m.pool = 0.01, a value that ensures a connection with a source of species but represents a system in which dispersal among local patches is the main source of recruitments (see also Worm and Tittensor, 2018). The pool of species available to colonize a local community was then estimated as the matrix product of the metacommunity abundance matrix and the standardized migration matrix. For all simulations, the species pool was defined at 200 species with a log-normal abundance distribution.

FIGURE 2

Figure 2. Conceptual model illustration representing the two main steps of the coalescent model. Step 1: the model starts locating one individual randomly selected from the species pool in each community. Step 2: communities are filled with 100 individuals randomly selected from the (i) species pool, (ii) neighboring communities, or (iii) the same community with probabilities m.pool (= 0.01), m.neighbour, and 1 − (m.pool + m.neighbour), respectively. Migration from neighboring communities is determined by geographic distance (d_ij) and species dispersal ability (d₅₀). Step 3: once all communities are filled, alpha and beta diversities are estimated.

Pondscape, dispersal, and biodiversity

A total of two complementary analyses were considered here for each studied pondscape. First, we evaluated the interaction between the geographic isolation gradient and the dispersal ability gradient (represented in the parameter d₅₀; 50 values ranging from 10 to 3,000 meters) as determinants of local alpha and beta diversities. This dispersal ability gradient surpasses the range of dispersals reported for the species inhabiting these waterbodies (Borthagaray et al., 2015a; Cunillera-Montcus et al., 2021). To this aim, alpha and beta diversities were estimated from coalescent models for each community along the dispersal ability gradient of species (d₅₀). Alpha diversity was estimated as the number of species in the community. Beta diversity was estimated as the average of Jaccard dissimilarity between one community and all the other communities in the network. Note that each dispersal ability d₅₀ defines one migration matrix M, and the pool of species available to colonize a local community in a pondscape. Complementary, the ratio between the alpha (beta) diversity value in the most central community (i.e., quantile 0.95 of the degree centrality distribution) and the alpha (beta) diversity value in the most peripheral community (i.e., quantile 0.05 of the degree centrality distribution) was estimated in the dispersal ability gradient (d₅₀) (hereafter, denoted as Ratio_alpha and Ratio_beta).

In the second approach, we built a null model to explore the importance of the spatial arrangement of waterbodies—the real spatial location of ponds along the landscape—for diversity patterns along the gradient of dispersal abilities (d₅₀). This null model preserves the total number of individuals and the number of communities in each metacommunity but randomizes their spatial location—removing the spatial arrangement of waterbodies (see Fortuna et al., 2006). The average of alpha and beta diversities predicted for the real pondscape ( ${\bar{S}}_{o b s}$ ) was estimated along the dispersal gradient (d₅₀). Then, this ${\bar{S}}_{o b s}$ is contrasted with the average diversity expected following the coalescent dynamic running on 200 randomizations of the spatial location of waterbodies ( ${\bar{S}}_{n u l l}$ ). Then, a Z − value for alpha and beta diversities was estimated as $Z - v a l u e = (\frac{{\bar{S}}_{o b s} - {\bar{S}}_{n u l l})}{S D (S_{n u l l})}$ , for each one of the dispersal abilities considered (d₅₀) (Gotelli and Graves, 1996). This Z − value represents the positive or negative effect of the pondscape on the average alpha and beta diversities. Complementary, the difference between the coefficient of variation of S_obs and the coefficient of variation of S_null was estimated along the dispersal gradient (d₅₀) for alpha and beta diversities. This difference reflects the heterogeneity in the local diversity among communities. Finally, we explored the association between pondscape features—i.e., number of ponds, area, and diameter (see Supplementary Table S1)—and the value of the dispersal ability (d₅₀) at which a maximum effect of the spatial arrangement of waterbodies (Z − alpha, Z − beta, Ratio_alpha, and Ratio_beta) was observed. Similarly, the relationship between pondscape features—area, ponds mean distance, and ponds density—and the magnitude of the effect of landscape arrangement on alpha and beta diversities was also explored.

Gamma diversity was not explicitly considered in the previous analysis. We explored the contribution of local communities to global diversity along the dispersal ability gradient (d₅₀). To this aim, we used the following equation: $γ = \bar{α} * (1 + \bar{β})$ . We found that this contribution is determined by the trend in α (see Supplementary Figure S1). Consequently, we focus on alpha and beta diversities, while trends in gamma diversity emerge from trends in previous results. All simulations and analyses were performed using R software (R Development Core Team, 2013), while data and codes to reproduce results are available in the online tutorial and the GitHub repository (https://metacommunity-lab.github.io/FEE_Coal_Rand_Pondscape/ and https://github.com/matiasarim/Coalescent_Lottery_Metacommunity_model).

Results

All results highlighted the importance of the interaction between the spatial arrangement of communities and species dispersal abilities (d₅₀) as a determinant of metacommunity diversities (Figure 3). As dispersal ability increased (d₅₀), local richness also increased, and communities were homogenized (i.e., beta diversity decreased). However, this interaction between biodiversity patterns and dispersal ability was markedly different among peripheral (low degree centrality) and central communities (high degree centrality) (Figures 3A, B). At low dispersal abilities (d₅₀ < 100 meters, left side of the curves), communities were locally assembled from rare colonization events ending with low local richness (ranging from 1 to 15 species depending on the pondscape) and a large beta diversity among them (ranging from 1 to 0.8 Jaccard dissimilarity index depending on the pondscape). At intermediate dispersal levels, differences in alpha and beta diversities between peripheral and central communities became greater. Thus, it was at this level of dispersal that the effect of the spatial arrangement of communities on biodiversity patterns was evident. Central communities (continuous line in Figures 3A, B) presented greater alpha and lower beta diversities, whereas peripheral communities (dashed line in Figures 3A, B) had the opposite pattern (presenting differences of almost 20 species and from 0.3 to 0.5 Jaccard dissimilarity index among central and peripheral communities in most pondscapes). Finally, at high dispersal abilities (d₅₀ > 1,000 meters, right side of the curves), the effect of the pondscape was lost, and communities were homogenized, presenting similar values of alpha and beta diversities in the geographic isolation gradient.

FIGURE 3

Figure 3. Interplay between the spatial arrangement of waterbodies and species dispersal ability as determinants of metacommunity biodiversity. Rows correspond to metacommunities of waterbodies. (A, B) Alpha and beta diversities along a dispersal ability gradient. Dispersal ability in x-axes is the distance at which dispersal falls to half of its maximum rate (d₅₀). The colors of points represent the isolation-centrality gradient (weighted degree centrality) of communities. A jitter was applied to better visualize the range of waterbody isolations. The black curves represent the set of more isolated (dashed) and central (continuous) communities in the pondscape. The inset figure is the ratio for alpha or beta diversity in the most peripheral vs. the most central community in the gradient of dispersal ability—i.e., Ratio_alpha and Ratio_beta. (C, D) Importance of the spatial arrangement of ponds for alpha and beta diversities along the gradient of dispersal abilities. Z − values correspond to standardized deviations of diversity predicted for the real pondscape structure ( ${\bar{S}}_{o b s}$ ) and the average diversity expected by the null model prediction when the location of communities was randomized ( ${\bar{S}}_{n u l l}$ ): $Z - v a l u e = (\frac{{\bar{S}}_{o b s} - {\bar{S}}_{n u l l})}{S D ({\bar{S}}_{n u l l})}$ , for each dispersal ability (d₅₀).

Randomized landscapes in the null model significantly reduced the heterogeneity in the geographic isolation of waterbodies (Supplementary Figure S2). The comparison of alpha and beta diversities between real pondscapes vs. randomized pondscapes (Z − value in Figures 3C, D) also indicated a strong effect of the spatial arrangement of communities at intermediate dispersal abilities. At very low dispersal abilities (d₅₀ < 50 meters), pondscape did not affect biodiversity patterns (Z − value for alpha and beta diversities closer to zero; Figures 3C, D, respectively). At intermediate dispersal levels, large differences with null model expectations were observed (Figures 3C, D). A randomized pondscape supported between 15 to 30 fewer species on average by the community than the real pondscape (see Figures 3A, B). Finally, at high dispersal rates (usually when d₅₀ is close to 1,000 meters), the landscape loses its relevance, and the differences between the real pondscapes and the randomized ones become smaller (Z − value closer to 0). These trends were related to the heterogeneity in local alpha and beta diversities among communities, where the coefficient of variation for local diversity among communities was larger in real than randomized landscapes (see Supplementary Figure S3). Again, this was true for those organisms with relative intermediate dispersal abilities. Notably, null model simulations were reported for dispersal distances up to 2,000 meters for reducing computation time considering that from this distance, observations were not different from expectations.

Overall, the six pondscapes provided different estimations for alpha and beta diversities, but all shared similar relationships between diversity patterns and the interplay between species dispersal ability (d₅₀) and the geographic isolation of communities. Indeed, organisms with intermediate dispersal ability are the ones strongly affected by the spatial arrangement of waterbodies (Figure 3). On the other hand, there were also differential patterns across pondscapes likely associated with the whole size of the systems, i.e., the number of waterbodies or total area (see Figure 4, Supplementary Figures S4, S5). First, the magnitude of the effect of the spatial arrangement of communities on local diversity ranged from 10 times more richness in central than peripheral communities in Albera (560 km²) to a maximum of three times increase in richness in Guils (15 km²) or Rocha (2 km²) (Figure 3A; Ratio_alpha). Similarly, beta diversity was reduced to half in the central communities of Albera and to 70% in other metacommunities (Figure 3B; Ratio_beta). Furthermore, the dispersal ability values for which the spatial arrangement of communities has a maximum effect also differed among pondscapes, according to their total areas (Figure 4 and Supplementary Figures S4, S5). In Albera, the species that are experiencing the greater effect of waterbody location are those with an average dispersal ability of approximately 170 meters, whereas, in other small systems (e.g., Guils or Rocha), the pondscape effect is more important at lower dispersal abilities (50 m). Indeed, despite the low statistical power, a significant association between the dispersal abilities (d₅₀) more affected by the landscape and system size (area) was detected (Figure 4, Supplementary Figures S4, S5). Second, it should also be noted that null model comparisons also detected differences in dispersal values (d₅₀) for which the real spatial arrangement of ponds was an important determinant of community diversity (Figures 3C, D). When comparing diversity estimated on real pondscapes with the null model expectations (randomized pondscapes), a variation across pondscapes ranging from Z − value = 50 to 9 and Z − value = −15 to −30 for alpha and beta diversities, respectively, were detected. Indeed, the effect of the spatial arrangement of patches in metacommunity diversity increased with the size of the study system (Figure 4, Supplementary Figure S4). Finally, all results in alpha and beta diversities presented opposite patterns but were not symmetric in their significance, the magnitude of effects, or dispersal ability values at which the larger effect of the landscape on biodiversity was observed.

FIGURE 4

Figure 4. Effect of pondscape spatial extent (total area) on the interplay between the diversity of local communities, species dispersal ability (d₅₀), and patch geographic isolation among the six metacommunities considered. (A, B) Relationship between the dispersal distance at which the maximum Ratio_alpha (minimum Ratio_beta) occurs and the area of the pondscape. (C, D) Relationship between the dispersal distance at which the maximum Z − alpha (minimum Z − beta) occurs and the area of the pondscape. Relationship between the dispersal distance (d₅₀) at which the maximum effect of the spatial arrangement of waterbodies on Ratio_alpha and Ratio_beta and the area of pondscapes.

Discussion

The metacommunity theory has consolidated a central role in ecology (Leibold and Chase, 2018; Chase et al., 2020; Record et al., 2021; Cathey and Brown, 2022). Landscape, and particularly, the difference in the relative geographic isolation of communities, is progressively suggested as a determinant of community assembly. Despite this, there are main features of the interplay between metacommunity processes and real landscape structure that are not clearly understood (Suzuki and Economo, 2021). Here, it is shown that those metacommunity processes that determine diversity are strongly enhanced in real landscapes by features seldom considered in theoretical or empirical approaches, that is, the existence of large heterogeneity in the geographic isolation of local patches. Congruently, with previous results, in a dispersal gradient, mass effects systematically increase the local richness and decrease beta diversity (Loreau and Mouquet, 1999; Mouquet and Loreau, 2002; Loreau, 2010; Borthagaray et al., 2015b) (Figure 3). However, the spatial arrangement of patches only has a large importance in metacommunity processes at intermediate dispersal levels. Consequently, landscape determines biodiversity patterns through those organisms with dispersal abilities that ensure access to central locations but limits dispersal in isolated communities. This study used a theoretical approach based on the real spatial arrangement of the pondscapes for advancing the role of landscape structure on metacommunity diversity.

Congruently with previous studies, we observed a monotonic increase in alpha diversity and a concomitant decrease in beta diversity from geographically isolated to central communities and from lower to larger dispersal abilities (see also Loreau and Mouquet, 1999; Economo and Keitt, 2010; Borthagaray et al., 2015b). In the gradient of dispersal abilities, low dispersal levels consistently involved low richness among all communities. Communities also presented a high beta diversity because of random variation in species dominance among communities, e.g., great priority effect (Fukami, 2015). The opposite is true at high dispersal abilities, being consistently high among all communities with a concomitant decrease in beta diversity (see also Loreau and Mouquet, 1999; Borthagaray et al., 2015b). As expected, the metacommunity phenomenon emerges as a determinant of biodiversity patterns at intermediate dispersal levels (Loreau, 2010; Leibold and Chase, 2018). However, it is also in this scenario when the spatial location of communities in the landscape really plays a key role in the metacommunity assembly. The geographic isolation gradient of the pondscape becomes an important determinant of diversity patterns for those species with the proper dispersal ability to be affected by the isolation differences among communities (Urban and Keitt, 2001; Economo and Keitt, 2010). Then, central ponds host a high diversity, enhanced by dispersal from neighboring communities, while peripheral ponds present few species, following beta diversity in the opposite trend. Indeed, the null model indicates that metacommunities having the same total area, number of local communities, productivity, or heterogeneity, but lacking a spatial arrangement covering a large isolation gradient may be prone to harboring lower diversity levels (Grilli et al., 2015; Haddad et al., 2017b; Chase et al., 2020). These results agree with empirical observations reported for the metacommunities herein considered and also for other systems (e.g., Vanschoenwinkel et al., 2008; Waterkeyn et al., 2011; Borthagaray et al., 2015a, 2020; O'Neill, 2016; Bellin et al., 2020). In this scenario, central communities provide a large set of species for recruitment elsewhere, while rare species persist in isolated communities (Scheffer et al., 2006; Chase and Shulman, 2009; Brown and Swan, 2010; Chase et al., 2010). Our findings further suggest that this phenomenon is fostered by the total area covered by the metacommunity and the number of local communities involved. In real landscapes, with the increase in the metacommunity size, the environmental heterogeneity and the number of organisms in the species pool might also increase, further promoting the action of the ecological and evolutionary processes to enhance diversity (McManus et al., 2021). These considerations and the detection of significant effects of the metacommunity size, in spite of the low statistical power, point to a large ecological phenomenon that has to be further considered elsewhere.

This study used a theoretical approach based on the real spatial arrangement of the pondscapes for advancing the role of landscape structure on metacommunity diversity. Our results evidenced the importance of real landscape in a strength and features that may not be evident when mean field models or artificial landscapes are considered, e.g., random fields (Borthagaray et al., 2014; Grilli et al., 2015), random graphs (McManus et al., 2021), tree graphs, chains, or grids (Economo and Keitt, 2008, 2010; Suzuki and Economo, 2021). The randomized landscape that hosted lower diversities removed most features of the spatial arrangement of communities as spatial modularity, clustering, the maximum distance between communities, etc. Features that were related elsewhere to metacommunity diversity patterns (Economo, 2011; Mougi and Kondoh, 2016; Suzuki and Economo, 2021). However, here we highlight the role of the spatial isolation gradient. Firstly, because the interaction between the geographic isolation and the species dispersal herein reported (Figures 3A, B). Secondly, a large reduction in the range of relative isolations is observed after landscape randomization (Supplementary Figure S2). Finally, because this isolation gradient is associated with most other features of landscapes, such as clustering and modularity.

Our use of neutral species pools and the spatial location of patches but not local conditions—area, nutrients, or heterogeneity—evidently left significant biological and environmental factors out of the analysis. However, it has been noted that this approach was fruitful in capturing the single effect of the spatial arrangement of the landscape on biodiversity patterns in several systems (Hubbell, 2001; Economo and Keitt, 2008, 2010; Muneepeerakul et al., 2008; Rosindell et al., 2011; Carrara et al., 2012; Borthagaray et al., 2014, 2023; Cunillera-Montcus et al., 2021). Differences in biology and local conditions may affect the dispersal abilities at which this phenomenon is important or its strength. However, the general biological messages reported here are robust to these deviations. However, the geographic isolation gradient is probably a main determinant of biodiversity patterns, and its effect peaks at intermediate dispersal levels, being dismissible in species with relatively low or high dispersal ability. Despite this being an intuitive expectation, it represents a metacommunity phenomenon that was seldom considered in theoretical models—both spatially explicit or implicit (Loreau, 2010; Leibold and Chase, 2018)—and in empirical studies—both in experiments (Grainger and Gilbert, 2016) or field surveys (Borthagaray et al., 2020). Notably, metapopulation models also identified the heterogeneity in the geographic isolation of patches as a determinant of species persistence (Gilarranz and Bascompte, 2012; Grilli et al., 2015).

The expectation of greater importance of metacommunity processes at intermediate dispersal levels is a cornerstone in the metacommunity theory (Hanski, 1999; Loreau, 2010; Leibold and Chase, 2018). The importance of the isolation gradient on biodiversity patterns has also been acknowledged in theoretical (Economo and Keitt, 2010; Thompson et al., 2017) and empirical studies (Chase and Shulman, 2009; Horváth et al., 2019; Borthagaray et al., 2020; Cunillera-Montcus et al., 2020a). Despite this, little attention was devoted to the conditions that foster or constrain the effect of the isolation gradient on community structure. The results, presented herein, also have impacts on management and conservation strategies (Brooker et al., 1999; Barnett and Belote, 2021). On the one hand, the importance of preserving the spatial arrangement of patches and not only an amount of habitat is herein supported (see also Brooker et al., 1999; Haddad et al., 2017a; Suzuki and Economo, 2021). Considering the change in the landscape may help us to better understand the determinants of biodiversity response to habitat degradation (Horváth et al., 2019; Cunillera-Montcus et al., 2021). On the other hand, combining empirical information about landscape structure and species dispersal abilities with theoretical simulations may represent a straightforward approach to guiding management decisions (Brooker et al., 1999; Haddad et al., 2017b; Resasco et al., 2017). Landscapes in general, and particularly pondscapes, are in a threatening scenario due to global change and human activities (Wood et al., 2003; Calhoun et al., 2017; Hill et al., 2021). The present contribution attempts to advance those landscape features that may support biodiversity and should be explicitly considered both for understanding the ecological mechanisms and preserving biodiversity.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

MA, AB, and DC-M lead the conception of the study. MA led the metacommunity model development. AB and DC-M contributed to the metacommunity model development and organized the database. MA, AB, and DC-M wrote the first draft of the manuscript and performed simulations and statistical analysis. AB, DC-M, IT, DB, MA-P, EO, XQ, SG, and MA contributed to the database. AB, DC-M, JB, IT, DB, MA-P, EO, TM XQ, SG, and MA contributed to the manuscript revision, read, and approved the submitted version. All authors contributed to the article and approved the submitted version.

Funding

This study was funded by the H2020 EU-funded project PONDERFUL (869296), CSIC I+D 2020 and CSIC groups (ID 657725) UDELAR, and FCE_1_2021_1_167009. DC-M was supported by the European Union—NextGenerationEU, Ministry of Universities and Recovery, Transformation and Resilience Plan, through a call from Universitat de Girona. Financial support was also provided by a grant PID2020-114440GB-I00 funded by MCIN/AEI/ 10.13039/501100011033.

Acknowledgments

MA and AB are grateful for the support provided by CSIC groups (ID 657725), CSIC I+D 2020, and FCE_1_2021_1_167009.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2023.1125607/full#supplementary-material

Supplementary Table S1. General characteristics of the six metacommunities studied: Number of waterbodies, Area, Density as the number of waterbodies/area, Degree range is the 0.05 and 0.95 quantiles of the weighted degree distribution, Different degree range is the difference between Q_0.95 − Q_0.05, Mean distance is the average distance among all ponds in the network, Diameter is the longest of the geodesic distances among all pairs of ponds, and Species richness is the total number of animals in the metacommunity.^*Corresponds to the species richness of plants.

Supplementary Figure S1. Gamma diversity along a dispersal ability gradient for the six metacommunities studied. Dispersal ability in x-axes is the distance at which dispersal falls to half its maximum rate (d50). Gamma diversity is estimated as γ = α * (1+β). Note that the trends in gamma diversity emerge from trends in previous results.

Supplementary Figure S2. Frequency distribution of the coefficient of variation in degree centrality estimated for 2,000 landscapes with the spatial distribution of ponds randomized. The vertical dashed line indicates the observed coefficient of variation in ponds degree centrality in the real landscape. Real degree centrality was always significantly larger than those observed in randomized landscapes.

Supplementary Figure S3. Difference between the coefficient of variation of the diversity predicted for the real pondscape structure and the coefficient of variation of the diversity expected by the null model prediction when the spatial distribution of ponds was randomized. Dispersal ability in x-axes is the distance at which dispersal falls to half its maximum rate (d₅₀).

Supplementary Figure S4. Effect of pondscape size (total area or the number of waterbodies) on the interplay between the diversity of local communities, species dispersal ability, and patch geographic isolation among the six metacommunities considered. First row, the relationship between the dispersal distance at which the maximum Ratio_alpha (minimum Ratio_beta) occurs and the area or the number of waterbodies of the pondscape. Second row, the relationship between the dispersal ability at which the maximum Z − alpha (minimum Z − beta) occurs and the area or the number of waterbodies of the pondscape. Note that area and the number of water bodies are highly correlated (r = 0.83).

Supplementary Figure S5. Effect of pondscape size (total area or number of waterbodies) on the interplay between the diversity of local communities, species dispersal ability, and patch geographic isolation among the six metacommunities considered. First row, the relationship between the maximum value of Ratio_alpha (Ratio_beta) and the area or the number of waterbodies of the pondscape. Second row, the relationship between the maximum value of Z − alpha (Z − beta) and the area or the number of waterbodies of the pondscape. Note that area and the number of water bodies are highly correlated (r = 0.83).

References

Ai, D., Gravel, D., Chu, C., and Wang, G. (2013). Spatial structures of the environment and of dispersal impact species distribution in competitive metacommunities. PLoS ONE 8, e68927. doi: 10.1371/journal.pone.0068927

PubMed Abstract | CrossRef Full Text | Google Scholar

Arellano-Rivas, A., De-Nova, J. A., and Munguía-Rosas, M. A. (2016). Patch isolation and shape predict plant functional diversity in a naturally fragmented forest. J. Plant Ecol. 11, 136–146. doi: 10.1093/jpe/rtw119