2.1. D-SPOM: coupling spatio-temporal dynamics wetland hydrology and metapopulation patch occupancy
Dynamic spatio-temporal patterns of metapopulation occupancy in patchy habitats
R Soc Open Sci, Jan 13, 2021; DOI: 10.1098/rsos.201309

We developed a dynamic version of the stochastic patch occupancy model, SPOM [6,39,40], to simulate the presence/absence of a given species in a time-varying wetlandscape. SPOM computes the distribution of occupied patches, in discrete time steps, by considering focal species traits (extinction, colonization and dispersal distance) and patch spatial organization (patch gap distances). A binary state variable, pi(t), is set to 1 if site i is occupied at time t, and 0 otherwise. From an initial distribution of occupied patches, dynamics are modelled as a discrete time Markov chain. At each time step, the model allows unoccupied patches to be colonized by surrounding occupied patches with probability

Similarly, species in occupied patches can go extinct with a probability

A discrete time SPOM is a homogeneous first-order Markov chain in which the state of the metapopulation at time t + 1 depends only on the state (occupancy pattern) of the metapopulation at time t [39]. For each patch and each time step, the probabilities of colonization and extinction events depend on colonization and extinction rates with exponential survival probability [41]


where Δt is the simulation time step and Ci(t) and Ei(t) are the colonization and extinction rates for the i-patch at time t. Note that both Ci(t) and Ei(t) are time dependent based on the current conditions of the landscape.

The key difference between the static and dynamic SPOM is that while SPOM assumes static conditions for landscape patches, D-SPOM captures the time variability in habitat structure, embedded by the amount of suitable habitat, Si(t), and patch spatial configurations, dij(t). The colonization rate Ci(t) is specified by the following equation:

where dij(t) is the pairwise distance between the patches i and j, D is the species dispersal distance, c is the species colonization rate and Sj(t) defines patch suitability. The inter-patch distances dij(t) were estimated based on perimeter-to-perimeter distances, which vary in time with expansion and contraction of patch areas and perimeters. The local extinction rate on the i-patch is inversely proportional to the variable Si(t), i.e. Ei(t) = e/Si(t).

In classical metapopulation theory [9,42], the metapopulation capacity, λmax, is derived for a static landscape (i.e. Si(t) = Si) and captures the impact of landscape structure—the amount of habitat and its spatial configuration—on metapopulation persistence. Species can persist if the equilibrium state p0 corresponding to global extinction (i.e. characterized by pi = 0 for any i) is unstable. This condition is met when the leading eigenvalue of the Jacobian matrix J of the system linearized around p0 is positive. By defining a landscape matrix M consisting of elements mij = exp( − dij/D)SiSj for ij and mii = 0, the Jacobian reads J = cM − eI, where I is the identity matrix. The leading eigenvalue of J can be expressed as maxe, where λmax is the leading eigenvalue of matrix M that contains information about the landscape and the quality of the patches. Therefore, the persistence condition reads λmax > e/c [6,9].

The theoretical prediction of the classical metapopulation theory [6,9] is that a species can persist whenever λmax > e/c. However, when simulating the stochastic discrete process using SPOM, a species with λmax slightly above the threshold e/c could go extinct due to demographic stochasticity. The novelty we introduce here with the D-SPOM consists in calculating the metapopulation capacity λmax at each discrete time step because of the different hydrological habitat conditions, dij(t), Si(t) and Sj(t) that are manifested at time t. In a dynamic landscape, λmax(t) represents the survival threshold for the hydrologic conditions at time t. Metapopulation dynamics are estimated in terms of wetlandscape occupancy, Ω(t), which identifies the temporal dynamics of the fraction of patches that are occupied by a focal species.

The advantage of the proposed D-SPOM approach is to explicitly account for patch dynamics influencing habitat availability and accessibility for a focal species. The model requires two fundamental variables of metapopulation dynamics: habitat suitability (e.g. patch areas) and connectivity (e.g. gap distances). The key novelty introduced here is to consider the two quantities as temporally variable, driven by external stochastic forcing. Following Bertassello et al. [26], we estimated temporal fluctuations in attributes (e.g. surface area) of each wetland resulting from the net of precipitation falling over each wetland contributing area, evapotranspiration losses and water exchanged with shallow groundwater (see electronic supplementary material for details). The choice of an appropriate proxy for habitat suitability, Sj(t), depends on species characteristics. For example, amphibians may be associated more with wetland perimeters since most anurans and many caudates lay their eggs in shallow water near shorelines [43], while fish are more linked to the stage inside the water body [44].

In this work, Sj(t) is first set equal to wetland area Aj(t), (see §§3.1, 3.2. and 3.3), and in subsequent analyses we consider only partial wetland area, where each wetland is divided into concentric annular micro-habitats for species that prefer shallow (less than 30 cm), intermediate (30–80 cm) and deep (greater than 80 cm) water depth (see §3.4). In simulating three metapopulations of a hypothetical metacommunity, our key assumptions are that (i) species have particular preferences for each zone, and (ii) they inhabit these zones independent of each other (i.e. no species interactions). By coupling the temporal variability in wetland habitat suitability, Sj(t), and gap distances between wetland habitat patches, dij(t), with fixed species traits (extinction rate e, colonization c and dispersal distance D), we used the static SPOM and D-SPOM to generate time series of metapopulation capacity, λmax(t), and patch occupancy, Ω(t).