Metacommunity resilience against simulated gradients of wildfire: disturbance intensity and species dispersal ability determine landscape recover capacity

effect of size and intensity. These results not only illustrate the mechanisms shaping the studied metacommunity but also more generally stress the strong role of metacommunity mechanisms and landscape structure in biodiversity resilience. Finally, this study highlights the importance of using theoretical approaches rooted in empirical data to determine metacommunity dynamics and the need to preserve and build connected and heterogeneous landscapes to address future disturbance scenarios.


Introduction
Disturbance has been recognized as a main driver of species turnover, coexistence and successional changes in natural communities (Gleason 1926, Clements 1936, Connell and Slatyer 1977, Connell 1978, White and Pickett 1985, Pickett and McDonnell 1989, Turner and Gardner 2015. The role of disturbances in landscape diversity has typically been linked to disturbance frequency, intensity and extent, which determine post-disturbance succession (Connell 1978, Turner and Gardner 2015, Crabot et al. 2020). More recently, it has been noted that species order of arrival after disturbances is also a main determinant of the differences in the assemblies of communities (Fukami et al. 2010, Vannette andFukami 2017). Wildfires represent a natural disturbance in ecosystems, but their frequency and intensity are expected to increase because of global change across several regions around the globe (Bowman et al. 2009, Moritz et al. 2012, Pausas and Fernández-Muñoz 2012, Young et al. 2017, Turco et al. 2018. This change in wildfire dynamics is expected to impact ecosystem resilience (Wright and Bailey 1982, Christensen et al. 1989, Prepas et al. 2009, Chang et al. 2013, Hegeman et al. 2014, Miller et al. 2019). However, we lack a deep understanding of the determinants of community resilience after wildfires, particularly regarding the role of landscape structure and metacommunity organization (Day et al. 2020).
The role of wildfires in ecosystem dynamics has attracted substantial attention in grasslands, forests and the transition between these ecosystems (Bernardi et al. 2016, Miller et al. 2019. However, less attention has been devoted to other ecosystems that are equally impacted, such as rivers (Gresswell 1999, Kelly et al. 2006) and wetlands, including temporary ponds (Bixby et al. 2015). Temporary ponds create a spatial mosaic of water bodies scattered throughout the landscape, forming a network linked by the dispersal of individuals, which defines a metacommunity (Leibold et al. 2004, Urban 2004, Louette and De Meester 2005. Furthermore, metacommunities are assembled by species belonging to different dispersal groups with contrasting dispersal abilities and dispersal networks (Urban andKeitt 2001, Borthagaray et al. 2015a), which may determine differences in the recovery capacities of species and consequently in community assemblies. Wildfires generate a mixture of both disturbed and undisturbed patches, changing the current metacommunity network (Wang andCumming 2010, Zozaya et al. 2012). Therefore, the metacommunity network and its changes due to wildfires will determine species recruitment after disturbance (Gresswell 1999, Dunham et al. 2003, Altermatt 2013, Verkaik et al. 2013, Whitney et al. 2016, Crabt et al. 2020, Cunillera-Montcusí et al. 2020b.
Mediterranean regions are frequently subject to wildfires, which are expected to intensify in the future because of global change (Pausas 2004, Kovats et al. 2014, Turco et al. 2018. In these regions, networks of water bodies, including temporary ponds, contain diverse faunal communities displaying a wide range of spatial requirements and dispersal abilities, ranging from passive dispersers to active flying dispersers (e.g. Anostraca, Gastropoda and Coleoptera). Field studies have shown that these macrofauna have a high resilience to the impacts of natural wildfires, observing that one hydroperiod after a wildfire, burned pond communities were similar to the unburned ones (Cunillera-Montcusí et al. 2019, 2020a. Wildfire impacts on species were related to the organism's dispersal capacity, which also determined recovery to pre-disturbance levels, with the exception of passive dispersers (Cunillera-Montcusí et al. 2019). Thus, during the postfire hydroperiod, traits fostering dispersal ability were positively selected (Cunillera-Montcusí et al. 2020a). From these empirical results, we hypothesized that metacommunity resilience is tightly linked with its network structure, which is determined by an organism's dispersal ability and the landscape alterations generated by wildfires. The connection between network structure and biodiversity dynamics is a cornerstone of metacommunity theory (Vellend 2016, Leibold andChase 2018), from which we can hypothesise that gradients of wildfire size (i.e. fire area) and intensity (i.e. fraction of burned communities) will produce a change in metacommunity network and consequently in local community resilience. Furthermore, it is expected that this change will interact with the dispersal ability of the considered species. However, the actual relationship between community resilience and disturbance size or intensity is not evident for temporary pond or other systems (Miller et al. 2019, Crabot et al. 2020, Day et al. 2020. For example, it is not evident whether resilience will present gradual or sharp transitions along disturbance gradients (Scheffer 2009, Ratajczak et al. 2018, at which level of disturbance size or intensity this transition could be expected (Moore 2018), what the strength of these disturbance properties could be as determinants of resilience (Seidl et al. 2017), or how the differences in dispersal ability among taxa might determine future patterns (Sorensen and Shade 2020). The empirical evaluation of these relationships is limited to a few ecosystems or taxa (Sorensen and Shade 2020), making theoretical evaluations rooted in real systems particularly valuable to advance our knowledge of disturbance and landscape interactions (Miller et al. 2019). Here, we analyse these general relationships between disturbance and resilience with a focus on a model study system -a water body network located on the NE Iberian Peninsula -for which key information about pondscape arrangement, species biology, community structure and response to disturbance is available (Cunillera-Montcusí et al. 2019, 2020a. To address this aim, we developed a metacommunity lottery model in which the pondscape was captured in a network inferred from species dispersal ability and regional water body locations. Throughout this model, we analysed the whole gradient of disturbance size and intensity, exploring its effect on community resilience, understanding its determinants and drawing its limits against future changes in wildfire disturbance regimes.

Study site
The studied network contains 542 water bodies found in the Albera region on the NE Iberian Peninsula (Fig. 1A). This network was partially affected by a wildfire in July 2012, which burned 13 000 hectares (Serveis territorials de Girona 2012) and affected 97 water bodies (Fig. 1A). It should be noted that within the burned area, not all water bodies were affected, creating a heterogeneous landscape of burned and unburned patches. Cunillera-Montcusí et al. (2019) for more information about wildfire impacts on some water bodies in the region.

Metacommunity model
The simulation model, summarized in Fig. 2, was a lottery model that was presented in several graphs that represent the landscape structure. Each node in the graphs represented a local community (i.e. water body), and each link represented the flow of individuals between the connected communities (Borthagaray et al. 2014). Species dispersal ability was used to estimate the associated graph networks. Resource overlap between dispersal groups was also included in the model as a sensitivity analysis considering a fraction of the resources of local communities (α) as shared among all dispersal groups and a fraction (1 − α) exclusive to each group. In the lottery dynamics, shared individuals could be replaced by any dispersal group, but exclusive individuals were only replaced by individuals from the same dispersal group. The sensitivity of the obtained results to α was evaluated considering a range of resource overlaps (α = 0.1, 0.3, 0.5, 0.7 and 0.9).

Dispersal networks
The lottery model considers the fact that metacommunities are assembled by species belonging to different functional groups with contrasting dispersal networks (Borthagaray et al. 2015a). Each of these dispersal groups has a maximum distance that determines its movement between local communities (following Urban and Keitt 2001, Estrada and Bodin 2008, Economo and Keitt 2010, Gilarranz and Bascompte 2012, Dale 2017). In the model, each dispersal group can only disperse between water bodies that are at a distance below this threshold. Thus, for a determined dispersal group, the water bodies that are below their corresponding threshold distance will be connected, conforming to the group-specific dispersal network. This allowed us to incorporate the different landscape scales that could be observed in the real metacommunity network (Borthagaray et al. 2015b). Based on the literature regarding the dispersal ability of the taxa occurring in the studied ponds and their reported dispersal groups, nine dispersal networks were constructed (Louette andDe Meester 2005, Shurin et al. 2009), which were related to four dispersal ability groups (Heino 2013, Cunillera-Montcusí et al. 2019, 2020a. Dispersal distances of 250 and 500 m corresponded to the dispersal ranges reported for species with passive dispersal ability and aquatic adults, Group 1 (e.g. Oligochaeta, Hirudinea and Gastropoda; Van De Meutter et al. 2007, Kappes andHaase 2012); 1000 and 1500 m corresponded to the dispersal range reported for weak active dispersers with flying adults, Group 2 (e.g. Diptera: Ceratopogonidae and Chironomidae; Delettre and Morvan 2000, Dettinger-Klemm 2003, Vallenduuk and Moller Pillot 2007; 2000 and 2500 m corresponded to moderate active dispersers with flying adults, Group 3 (e.g. Ephemeroptera and Trichoptera; Kovats et al. 1996, Nilsson 1996, Astorga et al. 2012; dispersal distances greater than 3000 m corresponded to species with strong active dispersal and flying adults, Group 4 (e.g. Odonata, Heteroptera and Coleoptera; Nilsson 1996, Conrad et al. 1999, Angelibert and Giani 2003, Chaput-Bardy et al. 2010, Miller and Bergsten 2016. According to this criterion, nine dispersal networks were considered representing the following dispersal distances: 250, 500, 1000, 1500, 2000, 2500, 3000, 3842 and 5000 m. It should be noted that in the study system, a distance of 3842 m corresponded to the percolation distance, defined as the minimum linkage distance at which all metacommunity network nodes become connected to at least one neighbour node (Keitt 1997, Urban and Keitt 2001, Solé 2012. This distance implied a transition in global connectivity from a fragmented to a totally connected landscape in which all water bodies are potentially reachable.

Lottery model
Metacommunities were assembled using a lottery model (Hubbell 2001, Borthagaray et al. 2014, Worm and Tittensor 2018. The model generated the deletion and replacement of individuals in each local community (i.e. water body) 10 000 times following these steps ( Fig. 2): 1) the order of the communities to be updated was randomized; 2) 50 individuals were randomly removed from a local community; 3) individuals were replaced with a probability of m with individuals from communities linked to the focal community and a probability of 1 m for individuals randomly selected from the same community (m = 0.01 for all species and all simulations); and 4) steps b to c were repeated until all communities are updated. Ten thousand repetitions of this cycle were enough for biodiversity stabilization in our simulations (Supporting information). The simulation started with a metacommunity of 200 species, each having 20 individuals in every local community, so a total of 4000 individuals constituted an individual local community (i.e. water body). Figure 2. Conceptual representation of the simulation model. Dispersal networks are defined according to species dispersal abilities. The lottery model consisted of 10 000 repetitions of replacing dead individuals in each community by randomly sampling individuals of the same or neighbouring communities (i.e. neighbours were defined by each species dispersal network). Wildfire scenarios were simulated by setting to zero the community abundance of a randomly selected proportion (fire intensity) of ponds within a burned area (fire size). Recolonization dynamics were based on a coalescent cycle that stopped when burned ponds were recolonized (number of initial individuals) or after 100 iterations of coalescent cycles. This species richness was close to that observed in the studied system (Cunillera-Montcusí et al. 2019).
Once the metacommunities were created, the wildfire simulation started (Fig. 2), eliminating all the individuals in the selected burned local communities. The burned water bodies were randomly selected from the 542 available in the landscape by defining different wildfire scenarios of varying intensities and sizes (Supporting information). After the wildfire, a coalescent recolonization dynamic was implemented in the following steps: 1) communities linked to the burned water bodies were identified according to whether they contained individuals (i.e. unburned or previously recolonized water bodies); 2) one individual was randomly selected for colonization from the linked communities; and 3) the community was then filled via a coalescent dynamic with the same fixed immigration probability (m). These steps were repeated until the community reached the initial number of individuals. After the first arrival, an internal recruitment probability (1-m) was also considered in step 'a'. Every repetition was recorded as an individual iteration. This process ended when all burned water bodies were recolonized or after a defined number of iterations of recolonization dynamics. We considered a total of 100 iterations following the coalescent filling of burned water bodies connected to unburned or recolonized bodies. Lottery and coalescent metacommunity models may provide similar results (Worm and Tittensor 2018). We opted for cycles of coalescent-community filling because it involved an arbitrary time scale that we could use later as a temporal measure of community resilience.

Metacommunity resilience analysis
We used a set of complementary metrics to elucidate the recovery capacity of communities after a disturbance (Table 1). First, we calculated the average proportion of the species pool that was available for the assembly of local water bodies immediately after disturbance (T1) and after 100 iterations of coalescent filling (T100). These proportions represented the available diversity for community assembly. Second, we estimated the number of coalescent iterations required to reach a burned community. This value based on time complemented the first value based on richness. The maximum number of iterations was fixed to 100 because it ensured that all dispersal groups reach burned water bodies if their associated metacommunity networks connected with them. Complementarily, these previous metrics were expressed in terms of resilience as the temporal rate of community recovery and the final extent of return to previous conditions after disturbance. These two metrics are conceptually related to engineering resilience and ecological resilience, which evaluate a system's time or rate to return to previous structural attributes and the amount of change needed to change the system state or equilibrium, respectively (Angeler and Allen 2016). Specifically, we estimated the return rate (hereafter RR) as the inverse of the number of coalescent iterations required to reach burned ponds ( Table 1). The return rate represents the speed at which the system recovers (Pimm 1991, Table 1). An RR closer to one implies a fast recovery, and an RR closer to zero implies a slow recolonization -or any recolonization capacity if RR = 0. Finally, we calculated the diversity recovery (DR), defined as the proportion of richness recovery between burned communities after 1 cycle and after 100 cycles of community refilling, quantifying how much change in diversity the communities experience throughout colonization, and thus, this value is linked to the ecological resilience concept (Holling 1973, Angeler and Allen 2016, Table 1). A DR value of 100 implies a high resilience (i.e. small diversity changes throughout recolonization), where with only one iteration, the community has Table 1. Metrics used to evaluate community recovery and their definitions. Diversity metrics are based on the average richness of burned communities (T1 and T100), time-based values (IT) and system resilience metrics: return rate (RR) and diversity recovery (DR). Note that an iteration corresponds to one cycle of community refilling by a coalescent process.

Metric Definition
T1 and T100 Average number of species in burned water bodies after 1 and 100 recolonization iterations. Define the state of the system right after the disturbance and after the whole recolonization process, respectively. Inverse of the number of iterations, which is the temporal rate of community recovery. RR measures the speed at which burned ponds return to the state before the disturbance, reflecting the engineering resilience of the system.
Diversity recovery (DR) Proportion of species richness immediately after the disturbance (T1) to recovered richness after community recolonization (100 coalescent iterations after the disturbance).
DR T T =( ) 1 100 100 recovered all the available richness. However, lower values quantify the deviation from the former richness level (i.e. high diversity changes throughout recolonization). The resilience metrics RR and DR, together with the mean richness at T1 and T100, provide complementary information about the recovery dynamics of disturbed communities. All calculations and model simulations were carried out using R software (<www.r-project.org>; Supporting information).

Model simulations
Simulations using the defined model were performed across gradients of wildfire size and intensity, including the specific case of the 2012 wildfire. The gradient of wildfire size was designated as every km in a buffer from 2 to 21 km from the starting location of the 2012 wildfire, and historically, these gradients have shown a north-south direction (Castellnou et al. 2010). The gradient of intensity was generated by simulating the percentage of burned water bodies within the burned area, with intensities ranging from 5 to 100% considering 5% increases. Through this process, we simulated 400 fires combining 20 levels of fire size with 20 levels of fire intensity. The response to the simulated wildfires was evaluated by calculating the average number of species present in burned water bodies after one (T1) or 100 cycles of community refilling, i.e. number of coalescent iterations, (T100) and the two metrics of resilience, return rate (RR) and diversity recovery (DR). The original 2012 wildfire information was obtained from the geographic information system of the Firefight Department of the Catalan Government (Serveis territorials de Girona 2012) and analysed using ArcMap 10.0 (ESRI 1999) to obtain the corresponding size and intensity of the fire following the current study definition.

Results
When considering the whole metacommunity, we reported a marked non-linear transition from an almost complete recovery (most species pool recovered) to a metacommunity having almost none of the available species (Fig. 3). This abrupt change implies that a small change along the wildfire gradient, mainly in its intensity, promoted a large decrease in diversity (Fig. 3A). Nevertheless, this transition vanished over time (Fig. 3B) as recolonization advanced through connected communities. Congruently, resilience metrics showed this nonlinear transition across dispersal groups for both the original wildfire scenario (Fig. 4) and the 400 simulated scenarios (Fig. 5).
Within the whole range of simulated scenarios, the original 2012 wildfire corresponded to a 10.5 km wildlife size and a 30% intensity disturbance ( Fig. 1B and black dots in Fig. 3, 5). Under this scenario, only some dispersal distances had low resilience (Fig. 4, Supporting information). In particular, when dispersal distances were greater than 2500 m, species presented return rate (RR) values of one, which meant that only one iteration was needed to reach burned ponds. For lower dispersal distances, the recovery rate decreased, and when the dispersal capacity was below 1000 m, the RR was particularly low, indicating that recolonization was not effective. This result was supported when the whole range of disturbance sizes and intensities were considered.
Within the simulated disturbance regimes, recolonization dynamics qualitatively differed among dispersal groups. Recolonization was only effective when the network was sufficiently connected to allow both direct arrival or arrival through the stepping-stone recolonization processes (Fig. 5, Supporting information). The return rate (RR) was small for functional groups with short dispersal distances across all disturbance scenarios (Fig. 5). Among dispersal groups with long dispersal abilities, a strong interaction between disturbance size and intensity was observed. The recovery rate was close to one when the disturbance area or intensity was low. For a large disturbance with a low intensity, the RR was close to one, and it was also close to one for a disturbance with a high intensity but a reduced size. RR values close to one along the intensity and size gradients were more frequently found with greater dispersal distances (Fig. 5, Supporting information). On the other hand, diversity recovery (DR) also presented an abrupt transition across taxa but was mostly related to the disturbance intensity (Fig. 5). After a threshold in disturbance intensity, DR dropped from almost 100% (richness recovered in only one iteration) to approximately 10% (communities assembled after disturbance had 10% of their final richness value). The wildfire intensity at which this change in DR was observed and the quickness of this transition were also determined by the dispersal ability of the species (Supporting information).
The average richness of the communities after one cycle (T1) and 100 cycles of coalescent filling dynamics (T100) as well as the number of iterations needed for diversity recovery showed a similar non-linear pattern dependent on the dispersal distance (Fig. 5, Supporting information). Once a threshold of wildfire intensity was surpassed, the species richness in the disturbed communities significantly decreased for all dispersal groups. However, short dispersal-distance groups (250 and 500 m) were unable to reach most burned ponds regardless of the number of iterations; therefore, the coalescent dynamics were unable to transport individuals of these organisms to burned water bodies (Fig. 5, Supporting information). Furthermore, for all simulations, the sensitivity analysis based on resource overlap showed the same trends in all different tested levels of α, indicating model robustness (Supporting information).
In summary, our results showed that the increase in the disturbance regime, mostly driven by intensity rather than size, was followed at some point by a transition from high to low resilience. This transition occurred along the intensity axis for diversity recovery (DR), richness immediately after disturbance (T1) or after 100 cycles of coalescent recolonization (T100). However, the rate of species recovery (RR) presented sharp transitions in community resilience related to both disturbance size and intensity.

Discussion
A mechanistic understanding of the relationship between wildfires and system resilience is urgently required at a time where the intensity and area of this type of disturbance are already increasing above historical records around the globe (Pausas and Ribeiro 2013, Kovats et al. 2014, Banks et al. 2017, Young et al. 2017, Turco et al. 2018, Day et al. 2020. Our study identified a gradual decline in resilience capacity with the increase in burned area but a sharp loss in resilience capacity when the fire affected more than half of the waterbodies (at approximately a 50% intensity). This sharp transition was also observed for diversity recovery in response to disturbance size. This transition highlighted that the recovery velocity was also attenuated for large disturbances, and thus, larger disturbances, though having low intensities, would also compromise diversity recovery. Consequently, an increase of larger disturbances may preclude community recovery because the stepping-stone recolonization dynamics of the metacommunity will not have enough time to recolonize impacted communities between disturbance events. This interpretation is further supported by the main differences in how wildfire intensity and size change the metacommunity network through which biodiversity is recovered. An increase in wildfire intensity disconnected burned communities, constraining species arrivals. However, stepping-stone recolonization is possible, even after large fires, if a determined fraction of water bodies remains undisturbed at the landscape level, a process that may involve several generations within community dynamics (Urban andKeitt 2001, Borthagaray et al. 2014). Species dispersal ability strongly determined the recovery patterns that appeared and were as important as the wildfire regime. Thus, the realistic range of considered dispersal distances and the reported loss of some of the dispersal abilities throughout this study represent threatening future scenarios, which may lead to an impoverishment of ecosystem functionality (France and Duffy 2006, Fox and Harpole 2008, Ptacnik et al. 2008, Valiente-Banuet et al. 2015. Our model simulated biodiversity dynamics at the regional level, something unattainable with an exclusively empirical approach. Lottery and coalescent models have been particularly suitable for analysing the effects of landscape structure on biodiversity patterns here and elsewhere (Economo and Keitt 2008, Muneepeerakul et al. 2008, Borthagaray et al. 2014, Sokol et al. 2017, Munoz et al. 2018, Worm and Tittensor 2018. Moreover, our model considered trait-mediated dynamics accounting for species dispersal abilities (Louette andDe Meester 2005, Ruhí et al. 2013), and resources overlapped across them (Shipley 2010, Sokol et al. 2017. In this sense, the current model focuses on reporting the potential biodiversity that is locally available after a wildfire. The subsequent assembly of colonized communities may also involve trait-mediated mechanisms (Weiher and Keddy 1999, Shipley 2010, Vellend 2016, Cadotte and Tucker 2017. In summary, we provided information regarding resilience through both neutral and niche mechanisms (Muneepeerakul et al. 2008, Halley and Iwasa 2011, Rosindell et al. 2011, May et al. 2015. Resilience metrics (RR and DR) and post-disturbance average species richness in burned ponds (in T1 and T100) for groups of species with contrasting dispersal abilities. Dispersal distances are indicated on the right side of the plot: Group 1 (represented by Gastropoda), Group 2 (represented by Chironomidae), Group 3 (represented by Ephemeroptera) and Group 4 (represented by Odonata). Size is the spatial extent covered by the simulated wildfire, and intensity is the proportion of burned water bodies within this extent. All other dispersal groups can be seen in the Supporting information. The number of iterations has not been plotted, as it is equivalent to the RR values, although these values can be seen in the Supporting information. Black points correspond to the location of the studied actual wildfire. All these results correspond to α = 0.5, other α values together with standard deviation values are presented in the Supporting information.
Cunillera-Montcusí et al. 2020a). The theoretical analysis of the present study captured the mechanisms suggested by empirical studies (Cunillera-Montcusí et al. 2019, 2020a, allowing us to explore a range of disturbance scenarios, which is not possible with empirical approaches. Therefore, we advanced the understanding of metacommunity dynamics in response to disturbances with a theoretical-empirical analysis of a model system with a focus on general theory. Spatial extent, intensity and frequency are the three main characteristics of disturbances (White andPickett 1985, Turner andGardner 2015). The studied region experienced a wildfire in 2012 that covered a large area, but it had a low intensity (Serveis territorials de Girona 2012). Unaffected water bodies acted as refuges during the wildfire and then as sources of recolonization, improving postfire recovery (Dunham et al. 2003, Uys et al. 2006. Connectivity has been identified as a key determinant of community resilience in rivers and streams, favouring fish recolonization of disturbed sites after wildfires (Fagan 2012, Whitney et al. 2015. In Mediterranean temporary ponds, connectivity is determined by the dispersal ability of different species and thus strongly varies among taxa (Boix et al. 2016). Consistent with previously observed patterns (Cunillera-Montcusí et al. 2019, 2020a, our simulations also presented distinct recovery patterns in different dispersal groups across wildfire regimes. Classic theories linking disturbance with biodiversity focus on species' abilities to access suitable resources once released (Connell 1978, Sheil andBurslem 2003). However, the idea of an implicit species pool with equal potential to arrive at any local community has been replaced by the recognition that the arrival of a species depends on the network configuration and the organism dispersal potential (Borthagaray et al. 2015a, Cañedo-Argüelles et al. 2015, Leibold and Chase 2018, Tornero et al. 2018. Strong, moderate and weak active dispersers (e.g. Coleoptera, Ephemeroptera, Heteroptera, Odonata and Trichoptera) would be able to rapidly reach burned sites (presenting high resilience values), and consequently, burned or unburned areas would not affect the role of colonization in community assembly. On the other hand, passive dispersers (e.g. Gastropoda, Hirudinea and Oligochaeta) may have difficulties recolonizing affected habitats, as they appear disconnected from the metacommunity species pool. In particular, the slow rate of species recovery indicated a great amount of time was needed for recolonization to be feasible. Therefore, the recovery of the last group after a wildfire will need to rely on other life strategies, such as dormancy in the egg bank or burying in the sediment (Wiggins et al. 1980, Chittapun 2011, Johnston et al. 2014, Strachan et al. 2014, Cunillera-Montcusí et al. 2019). In the end, the studied metacommunity was strongly resilient to the original 2012 wildfire, recovering its composition in only one hydroperiod (Cunillera-Montcusí et al. 2019, 2020a. However, this promising observation is challenged by our current theoretical analysis, which indicates that this resilience would be highly compromised with the expected increase in wildfire regimes (Miller et al. 2019).
Community recovery after a disturbance mainly relies on the ability of the species to move through the landscape. At specific disturbance levels, our analysis predicts an abrupt collapse that would be characterized by a decrease in resilience metrics and a transition from high to low species richness after the disturbance. This pattern is consistent with several results already found in landscape studies from a wide range of perspectives (Urban and Keitt 2001, Scheffer 2009, Gilarranz and Bascompte 2012, Borthagaray et al. 2014, Lenton et al. 2019). Thus, it has been shown that across a gradient of species dispersal distances, there is a network 'percolation distance' that determines an abrupt transition in species potential to move locally or globally (With and Crist 1995, Urban and Keitt 2001, Borthagaray et al. 2014, Bertuzzo et al. 2016). This scenario occurs because when dispersal is high enough, the potential for species to reach habitat patches increases in a stepping-stone manner, allowing them to reach any local community in the landscape (Keitt 1997, Urban andKeitt 2001). In the current study, the transition from low to high resilience was observed at shorter distances than the percolation distance (approximately 3800 m for the study system). This result implies that species that disperse below this distance, which were experiencing a globally fragmented network, were still able to recolonize affected habitats. Consequently, the dispersal distance at which the transition in resilience was observed did not match the percolation distance (Borthagaray et al. 2014). Furthermore, we observed that the arrival of species to disturbed communities depended on stepping-stone dispersal through unburned and recolonized communities (Keitt 1997). Thus, in essence, our results show that landscape disconnection, due to the loss of occupied water bodies in the network, broke the colonization process when there were high fire intensities and/or when species had short dispersal distances, driving metacommunity resilience towards collapse. Consequently, ensuring the persistence of undisturbed communities in the landscape would appear key for improving biodiversity recovery after disturbances, in both the studied region and other ecosystems. Thus, conserving landscape heterogeneity and mitigating the fraction of burned areas may represent a key management strategy to cope with future scenarios (Gongalsky et al. 2012, Chia et al. 2015, Banks et al. 2017.
The metacommunity network that determines a system's resilience is based on the interactions among species with dispersal abilities (De Bie et al. 2012, Jones et al. 2015, Schuler et al. 2017, Tornero et al. 2018, the relative spatial locations of local communities (Urban and Keitt 2001, Economo and Keitt 2010, Borthagaray et al. 2014, 2015b, and landscape resistance to individual movement (Cañedo-Argüelles et al. 2015). These three determinants have seldom been considered jointly in management strategies (Turner and Gardner 2015). However, a system's resilience is determined both by the interaction among these three determinants and by disturbance characteristics. Species dispersal and landscape structure are progressively compromised by human activities (Shanafelt et al. 2018), and the intensity and size of wildfires are increasing with global change (Pausas and Ribeiro 2013, Kovats et al. 2014, Turco et al. 2018. Systems that until now were perceived as resilient might experience an exacerbated loss in their capacity to recover after a disturbance. In this sense, the metacommunity network is key to maintaining a system's resilience against increasing threats, as it maintains greater biodiversity levels (Horváth et al. 2019). Thus, future disturbance regimes will result in metacommunities being in more unstable scenarios, and in line with other globalscale processes (Steffen et al. 2018, Lenton et al. 2019, the reported non-linear change in resilience will compromise the recovery of the entire system. Considering such dynamics with a landscape perspective is key to fully comprehending the determinants of resilience at wider scales (Chase et al. 2020). Coping with future disturbance scenarios is a pressing challenge in a changing world, especially in fire-prone regions, which may become affected by even larger and more catastrophic wildfires (Moritz et al. 2012, Turco et al. 2018. Promoting connected landscapes and reinforcing a system's resilience considering the regional network of the metacommunity constitutes one main path to pursue in the future (Chase et al. 2020), especially for shallow aquatic systems, which are already threatened by human activities (Wood et al. 2003, Pekel et al. 2016, Bastin et al. 2019.