Role of trait combinations, habitat matrix, and network topology in metapopulation recovery from regional extinction

We studied the role of oceanographic conditions and life history strategies on recovery after extinction in a metapopulation of marine organisms dispersing as pelagic larvae. We combined an age‐structured model with scenarios defined by realistic oceanographic conditions and species distribution along the Irish Sea coast (North Europe). Species life history strategies were modeled combining the dispersal behaviors with two levels of fecundity. Recovery times were quantified after simulating extinction in four regions. Two alternative strategies (high fecundity or larval tidal transport) led to short recovery times, irrespective of the effects of other drivers. Other strategies and low larval survival exacerbated the effects of oceanographic conditions on recovery times: longer times were associated with for example the presence of frontal zones isolating regions of extinction. Recovery times were well explained by the connectivity of each focal population with those located outside the area of extinction (which was higher in the so‐called small world topologies), but not by subsidies (direct connections with populations located nearby). Our work highlights the complexities involved in population recovery: specific trait combinations may blur the effects of the habitat matrix on recovery times; K‐strategists (i.e., with low fecundities) may achieve quick recovery if they possess the appropriate dispersal traits. High larval mortality can exacerbate the effect of oceanographic conditions and lead to heterogeneity in recovery times. Overall, processes driving whole network topologies rather than conditions surrounding local populations are the key to understand patterns of recovery.

Disturbance is a fundamentally important process in all ecological systems, modifying resource availability and causing disruption to population, community, or ecosystem structure (White and Pickett 1985). A disturbance event by definition occurs over a relatively discrete period of time, but the spatial scale over which it operates can vary from that of an individual through to whole ecosystems. For instance, in the marine environment, regional scale disturbances covering areas of the order of 10-10 3 km may be produced by a range of drivers including summer anoxia at the sea bed (Bishop et al. 2005;Diaz and Rosenberg 2008), storms (Woodley et al. 1981), heat waves, extreme temperatures (Glynn 1993;Coma et al. 2009), and pathogens (Miller and Colodey 1983;Lessios 2016). Many disturbance regimes are currently changing, with profound shifts expected in the coming decades as consequence of climate change and increasing encroachment of human activity over previously pristine habitats (Turner 2010).
Current climate change projections, for example, suggest that extreme weather events, including heat waves and storms, are likely to increase in frequency and magnitude (e.g., Burrows et al. 2014;Di Lorenzo and Mantua 2016) leading to regional scale levels of mass mortality. Current coral bleaching, as consequence of the recent El Nino, is an example of the spread and importance of regional scale events (Tollefson 2016).
For marine species, determining the fate of pelagic larvae is central to understanding the consequences of mass mortalities in terms of population recovery. Benthic or demersal populations are patchily distributed, often in association with specific habitats, but populations are connected through a dispersive larval stage (Cowen and Sponaugle 2009). It is therefore appropriate to use the concept of metapopulation, for such populations. Metapopulation theory was developed with the idea of modeling local extinction and recovery in populations with limited connectivity and has provided a valuable theoretical framework for conservation ecology (Hanski and Gaggiotti 2004), but the concept has been expanded to consider cases where populations are well connected. Hence, one may define marine metapopulations as a set of local populations of adults connected through larval dispersal (Armsworth 2002) with the direction and magnitude of larval dispersal pathways determining patterns of connectivity *Correspondence: Luis. Gimenez@awi.de This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Additional Supporting Information may be found in the online version of this article. and hence the extent to which different locations receive larval subsidies (Levin 2006;Cowen and Sponaugle 2009).
Given the major logistical challenges of directly quantifying larval connectivity, most effort has focused on modeling patterns of larval transport through hydrodynamic models (e.g., Cowen et al. 2006;Paris et al. 2007;Ayata et al. 2010;Robins et al. 2013). In such models, the connectivity and retention coefficients represent the main characteristic of the habitat matrix (i.e., the habitat through which organisms disperse and migrate : Wiens 1997;Joly et al. 2001;Shima and Swearer 2009). These models highlight the importance of local hydrodynamic conditions and species-specific larval behaviors in driving population persistence (e.g., Cowen et al. 2006;North et al. 2008;Botsford et al. 2009); hence, they should also drive the rate of recovery from extinction. A critical output of such models is that in situations where the model domain covers sufficiently large spatial scales (e.g., Cowen et al. 2006;Robins et al. 2013) limitations in dispersal define regions that are weakly connected with each other. Such spatial patterns may result in a reduced capacity to recover from mass mortalities if the scale of disturbance matches the scale of connectivity.
Models coupling dispersal with local processes have helped to understand the conditions of persistence of populations (Armsworth 2002;Hastings and Botsford 2006;Artzy-Randrup and Stone 2010). Such work has recognized pressures on local populations but given the increasing regional scale of disturbance events, there is, in addition, an urgent need to understand the drivers of recovery from regional extinction. Marine metapopulation models should help us to understand how quickly local populations may recover from mass mortalities especially if they are applied to realistic metapopulations. Here, we applied an age-structured metapopulation model to a realistic scenario by modeling the dynamics of a coastal species distributed in fragmented populations in the Irish Sea (Northwest Europe). The Irish Sea, like many other coastal areas worldwide, is impacted by regional and global scale phenomena (see Robins et al. 2016) exposed to a range of anthropogenic stressors leading to local and regional mortality events of benthic species (e.g., Malham et al. 2012). We develop a model for coastal species restricted to sheltered bays and estuaries (habitat patches) in order to examine patterns of recovery from extinction in well-defined regions throughout the Irish Sea. The following questions were addressed: Given an event of regional extinction: (1) Do specific life histories (i.e., combinations of traits such as fecundity and larval behavior) enhance recovery? In particular, is there any optimal strategy enhancing recovery or are strategies contextdependent, that is, do they depend on the region within the habitat matrix? (2) What is the role of the habitat matrix in setting the time scale of recovery? More specifically, would larval survival and variations in oceanographic conditions (as captured by the connectivity matrix) lead to significant spatial or temporal patterns in recovery? (3) What is the importance of the network topology? In particular, is recovery explained by retention or direct subsidy to a focal population or is it driven by whole network connectivity to populations located outside the region of extinction?

General procedure
We constructed a metapopulation model of a species occupying shallow, sheltered habitats in 40 populations throughout the Irish Sea and dispersing during the larval stage. Specifically, we consider a model species living up to 5 yr and starting to reproduce after the first year of life. The direction and magnitude of connection among discrete populations are given by larval transport matrices, accounting for the physical transport of larvae. The larval transport matrices were obtained from particle tracking models incorporating realistic three-dimensional (3D) hydrodynamic conditions. Three main larval strategies were considered : (1) passive transport (no vertical migration); (2) diel vertical migration (upward swimming during the night and downward swimming during the day); and (3) flood tidal migration (upward swimming during the flood phase and downward swimming during the ebb). Vertical swimming speed was set to 3 × 10 −3 m s −1 , representing ciliated larvae (Chia et al. 1984). Robins et al. (2013) showed that dispersal distance, retention, or connectivity did not vary in the range of speeds 1-5 × 10 −3 m s −1 but that speeds < 1 × 10 −3 m s −1 such patterns would resemble those obtained for passive transport.
Regional extinction events were simulated based on total loss of populations in one of four regions (Fig. 1). Overall, we ran 96 simulations differing in the combination of larval strategy (three levels), combinations of fecundity and larval mortality (two levels, controlled through a composite parameter, see below), timing of larval release (two levels-spring and summer), region of extinction (four regions), and the strength of densitydependent mortality in the benthic phase. For each simulation, the model was run for 200 cycles (= years) followed by an event of regional extinction. Recovery time was then quantified after running the model for an additional 400 yr. The model was run in Matlab R (see Section S1 Supporting Information for code).

Metapopulation model and dispersal matrices
The model is a modification of the one developed by Armsworth (2002). Using five age classes and 40 populations, the model is as follows: In Eq. 2, f is the fecundity, η is the larval survival during pelagic dispersal (i.e., due to sources other than overdispersion, e.g., predation, stress, and food limitation), l p p is the fraction of larvae that would return to the original population if mortality were zero, and σ is the probability of settled individuals reaching the first year of age. Notice that we separate physical transport from larval mortality; hence, larval dispersal would be given by η × l p p (or η × l i j see below). In addition, connectivity will be given by σ × η × l p p (σ × η × l i j : see, e.g., Lett et al. 2015). Here for simplicity, we refer to l p p and l i j as larval retention and larval connectivity, respectively. In Eq. 1, a series of Q matrices (see Eq. 3) define the number of larvae originating in each population that settle and survive over the first year in the local population. Each Q matrix is a square matrix of 5 × 5 cells (corresponding to the age classes) given by Eq. 3.
For the upper row, a ij = σ × l i j × η × f represents the product of the larval production and survival, and the larval connectivity coefficient (l j i ).
The larval connectivity and retention coefficients were those given by Robins et al. (2013). Robins et al. (2013) used a coupled 3D hydrodynamic and Lagrangian particle-tracking model to simulate scenarios where larvae are released at one time in the year and allowed to disperse for 28 d under realistic wind, temperature, tidal, and photoperiod conditions. Importantly, the connectivity and retention coefficients were obtained from realistic distribution of populations throughout the region, coastal geography, and sea bottom topography, all contributing to the patterns of dispersal.
The number of larvae competent to settle at the end of the pelagic period is a function of the number produced (fecundity rate) and the number of these which survive (mortality rate). We lack any evidence to realistically vary either production or loss of larvae across geographic locations and age class of adult. Hence, fecundity and larval survival rates were modeled as density-independent and constant over the whole model domain and were combined to form a single parameter, ω = f × η; ω may be interpreted as a component of the maternal fitness, that is, the product of offspring number and survival. Variation in larval survival has been shown to modify patterns of connectivity (Paris et al. 2007); hence, the metapopulation model was run using two values of ω (10 and 10,000 larvae per reproductive adult). The values of ω are chosen in first instance to cover a wide range in the parameter space. Second, such values cover the range of fecundities and known mortality rates observed in marine invertebrates. For instance, instantaneous daily mortality rates for marine larvae range from 0.22 (Rumrill 1990) to 0.14 (White et al. 2014a), resulting in ca. 0.13-1.50% survival over 28 d (used to model the connectivity matrices). Combining these estimates (η = 1.3 × 10 −3 to 1.5 × 10 −2 ) with the values of the term ω used in the model (10, 10 4 ), we obtain a realistic range of fecundities (from f = 10 3 for ω = 10 to f = 10 7 for ω = 10 4 ).
Overall, the model incorporates two sources of mortality, one that was captured in the parameter η, and one that was caused by dispersal away from suitable habitat (subsequently termed overdispersion). Overdispersion is reflected in the coefficients of connectivity and retention since at the end of the simulation any larvae that do not reach the population from which they arose or reach one of the other 39 populations are considered dead. Overdispersion is affected by larval behavior and by temporal changes in oceanographic conditions ). In the present models, we assume that overdispersion and η do not covary, but we recognize potential sources of covariation; for instance, diel vertical migration may lead to specific patterns of overdispersion while minimizing mortality by predation.
The number of larvae arriving to a focal population p, (S t,j = p ) is defined by the contribution of the focal population, accounted for in the first row of the M matrix (Eq. 2), and the subsidy from other populations, accounted for in the first row of the Q matrices (Eq. 3). These contributions are calculated as: The first term in the right hand side of Eq. 4 defines the total number of individuals retained as the product of the number of adults, fecundity, and survival (= ω), and the larval retention coefficient. The second term defines the input of larvae from other populations (subsequently termed subsidy), similar to retention, but based on the larval connectivity coefficients between each particular population, j, and the focal population, p. For simplicity, we assume that larval retention or connectivity coefficients are constant across ages and through each event of larval release (no t or k subscripts in Eq. 4). Note that when the focal population, p, goes extinct, recovery is in the first instance governed by the subsidy and hence the first term of Eq. 4 is zero. Eventually, once individuals reach the reproductive age, recovery will also depend on retention coefficients.
Following larval settlement, and over the first year of benthic life, the survival rate, s, in the M and Q matrices, was modeled as a density-dependent process according to the Beverton-Holt equation. The number of individuals surviving the first year of life (N 1 ) was: where S t is the number of larvae arriving to the nursery habitat (Eq. 5) and α and β are parameters; the subscript 0 indicates that these correspond to young of the year. We defined α as the density-independent survival parameter (since N 1 /S!α when S!0); β is the density-dependent parameter (if β = 0, N 1 is proportional to S t ).
In the following years, all individuals experience mortalities depending on the total number of individuals present in that habitat. The survival rate in the adult habitat (g i in the M matrix) was also modeled by a Beverton-Holt equation. Therefore, the number of individuals of age k surviving to the age k + 1 is: Giménez et al.

Metapopulation recovery after disturbance
where the subscript "a" indicates that the parameters of the function that correspond to the survival of individuals of one or more years of age. The density-dependent coefficient β for both the first year of life, and subsequent years, was varied between 0.1 and 0.0001 to explore the role of benthic survival on recovery. For simplicity, we assume that for all populations, α 0 = α 1 = 1 and β 0 = β 1 = β in Eqs. 4, 5. In such a model, for each individual simulation, the variation in the response variables among populations depended solely on the coefficients of the connectivity matrix. Taken together, the simplifications of α 1 , β, and ω would result in a "stage-structured" model composed of a juvenile stage (with zero fecundity and 1 yr duration), and an adult stage (with nonzero fecundity and 4 yr duration). However, we prefer to present the model in an "agestructured" form shown in Eqs. 1-3 because it is easier to interpret each iteration as equivalent to 1 yr of duration.

Simulations of extinction
Initially the model was run for 200 yr, covering 24 different scenarios, varying in all possible combinations of larval strategy (diel, tidal, passive), season (spring, summer), the term ω (= 10 or 1000 larvae per reproductive adult), and density-dependent coefficient β (= 0.1 or 0.0001). Each model simulation was initialized with 10 individuals per age class at each population. In each cycle, the model computed, at time t, the number of larvae produced and settled at each site, using Eq. 3. The number of settlers was then used, also at the time t, to update S t in Eq. 4. Then, the model computed, at the time t + 1, the number of individuals surviving the first year of life according to Eq. 5. The number of "first years" was then used to update the term N k in Eq. 6, which is used to compute, at the time t + 2, the number of individuals surviving to age = 2 yr. At a given year t + i, the terms in the local matrix (Eq. 2) depending on Eqs. 4-6 were updated simultaneously according to S t and the total number of individuals in the adult habitat.
After the initial 200 yr run of each of the 24 simulations, we then simulated extinction and measured recovery for four regions: (1) Cardigan Bay (populations 4-8 in Fig. 1), (2) Anglesey (populations 9-12), (3) Liverpool Bay (populations 13-22), and (4) central Ireland (populations 33-37). Each extinction and recovery simulation was run separately for each of the four regions, giving a total of 96 simulations (= 24 × 4; i.e., only one region suffered an extinction event at any one time, in a given simulation). We did not simulate extinction in any population located at the border of the model domain since the recovery of such populations should be affected by subsidy populations outside the model domain. Extinction was simulated by setting abundance to zero for all age classes of the populations in the target "extinction" region. The model was then run for a further 400 "yr" and the rate of recovery quantified as T 50 , the time (in years) required by populations to reach 50% of the asymptotic population size. We run a series of preliminary simulations in order to check the behavior of the model: these results are detailed in Section S2 in Supporting Information.

Statistical analysis
We used a statistical approach to quantify the average effect of each driver (region of extinction, month of larval release, β, and ω) on recovery time (T 50 ). We also studied the time needed for a population to double its size when rare (Section 2.2 in Supporting Information) but found that T 50 was more useful as a descriptor of the recovery rates. Statistical analyses (on T 50 ) were run in R (R Core Team 2013). We followed the recommendations of White et al. (2014b) and focused on effect sizes, since applying significance testing to simulation outputs would not be appropriate in a modeling exercise. Effect sizes were quantified using two techniques, boosted regression trees (BRTs) (James et al. 2013) and general least squares (GLS) models (Zuur et al. 2009;Galecki and Burzykowski 2013); such techniques have been used to identify key drivers of metapopulation connectivity (Treml et al. 2015).
BRTs were carried out following Elith et al. (2008). BRT is a technique of nonlinear model fitting based on so called "decision trees." Decision trees partition the predictor space (defined here by our drivers of connectivity) into regions of similar values of the response variable (recovery time) and then fit a constant to each region (Elith et al. 2008). Boosting is a numerical optimization technique minimizing the error through the cumulative fitting of additional trees to the data (each tree is fitted to the residuals obtaining from the fit of a previous tree). We use BRT as way to quantify the importance of each predictor, given as their relative influence (i.e., the proportion of the number of trees where a given predictor is fitted). BRTs were fitted in R, using the package dismo and the function gbm.step; this function enables the use of a crossvalidation method, based on testing the models of the fraction of the data ("bag fraction" = 0.5 in our case) in order to select optimal number of trees for the model. We fitted four models differing in the "learning rate" (range 0.05-0.0005), a parameter controlling the contribution of each tree to the model. All models were fitted with normal residuals and a tree complexity of 5, that is, considering the highest (five-way) interaction. All model fits led to similar patterns in the relative influence of the predictors on the recovery time; we present the results corresponding to a learning rate of 0.01.
In order to interpret the output of the BRT, we present plots of averaged recovery rate in response to all combination of parameters as well as plots of parameter estimates obtained from a GLS model. The GLS was fitted using the package nlme (Pinheiro et al. 2018) considering variance heterogeneity (VarIdent constructor function) and correlations among sites (CorCompSymm function). Although our design was fully replicated, our attempts at fitting the full model for the variance structure led to situations of nonconvergence; when this occurred, we reduced the complexity of the variance structure in the starting model. Models were fit using restricted maximum likelihood method. The GLS technique was applied to the logarithmically transformed values of the recovery time. The full model was a fifth order full factorial for both the variance and the fixed structure (i.e., for the fixed structure: T 50~ω × β × strategy × month × region) and parameter estimates were extracted.

Network topology
In a separate group of models, we evaluated the role of network topology. First, we evaluated how well recovery time was predicted by local subsidy and retention. In this case, we used the coefficient of determination (R 2 ) as a metric of effect sizes because our focus was on the importance of subsidy or retention in explaining variation in recovery time. Subsidy was defined as the sum of the connectivity coefficients indicating input of larvae toward a focal population. Note that subsidy is calculated as the sum of transport coefficients directly connecting the focal population to others many of which that may be inside the area of extinction. For instance, for a focal population located in the center of the area of extinction, the most likely scenario is that subsidy is entirely dependent on adjacent populations located inside the area of extinction. Larval connectivity to the outside source may occur indirectly, that is, through one or more local populations (e.g., in a stepping-stone pattern); it will depend on transport coefficients linking the focal and other populations with those outside the area of extinction and it is calculated as a product. Hence, for each focal population, we used two different indices of connectivity to the sources located outside the region of extinction, the total connectivity (CT) and the one provided by the path giving the maximum connectivity (CM  Fig. 2. Flow diagram summarizing the relationships between factors or natural drivers defining the habitat matrix, the traits, and the terms or coefficients determining the recovery time (T 50 ). The habitat matrix is characterized by biotic and abiotic factors, affecting larval survival (manipulated here variations in the term ω) and by oceanographic conditions driving larval transport and hence the coefficients of the connectivity matrix. The traits are the fecundity (contributing to the term ω) and larval behavior (contributing to the coefficients of the connectivity matrix). The term ω is the product of larval survival and fecundity; l ij denotes the coefficients of the transport matrix. Difference between parameter estimates at predictor levels vs. the reference, estimated from GLS model (from summary output, total of 96 parameters). In (b), the references correspond to ω = 10, β = 0.0001, month = April, dispersal = passive, and region = Cardigan Bay. Each dot corresponds to the difference between the reference level and another level, for a given combination of predictors. For instance, for ω, the additional level is ω = 10,000 and there are 48 dots corresponding to the combinations of levels of all other predictors (region, β, month, and strategy: 48 = 4 × 2 × 2 × 3).
For ω, β, and month, there is a single column of symbols because there is only a single level other than the reference. Notice in (a) that β is the predictor with less relative influence; in (b), this coincides with almost no difference between parameter estimates obtained at the reference (β = 0.1) vs. β = 0.0001 irrespective of remaining predictors (most dots are on the zero line). By contrast, ω has a strong relative influence (a); (b) shows that differences vary depending on other predictor combinations, but they are always negative, indicating that higher ω drives consistently reduced recovery time. Region (as month and strategy) had high relative importance (a) but the magnitude depended on the combinations of other parameters (in b, differences are negative or positive).
located outside the region of extinction either direct or indirect (the latter is calculated as a product of coefficients).
In Eq. 7, each l v ij is a larval transport coefficient between populations forming a path v between the focal population and the source populations located outside the region of extinction. For CT, all paths are considered; one such path, the one used to calculate CM, has the maximum product of the associated transport coefficients: For example, if populations were connected to a source S 0 through a single path (a stepping-stone pattern): S 0 !P 1 ! P 2 !P 3 and larval connectivity were l 1 0 = 10 −1 , l 2 1 = 10 −1 , l 3 2 = 10 −2 , respectively, then CM 1 = 10 −1 , CM 2 = 10 −2 , and CM 3 = 10 −4 . If by contrast S 0 were also connected to P 2 with l 2 0 = 10 −1 , then CM 1 = 10 −1 , CM 2 = 10 −1 , and CM 3 = 10 −3 , because such alternative path leads to higher CM for populations 2 and 3. In many cases, transport between populations i and j is bidirectional because of nonzero coefficients occurring in both directions (l j i > 0 and l i j > 0). Bidirectionality in transport between populations is common; although in most cases, there are strong asymmetries because currents flow in a predominant direction. Where transport between two populations was more symmetrical (i.e., where differences in coefficients were not large), the calculation of CM is based in the first instance on connections with populations from that region, and in the second instance, on the largest coefficient connecting two populations. For example, if we have S 0 !P 1 !P 2 $P 3 and l 3 2 < l 2 3 then, for P 2 , l 3 2 would still be the appropriate coefficient because l 2 3 would only enable subsidy from S 0 to P 2 once P 2 subsidizes P 3 through l 3 2 . However, in a case of two sources, for example, S 0 !P 1 !P 2 $P 3 S 1 , then CM for each population was calculated from the source giving the largest coefficient: if for P 2 , CM-S1 = l 3 S1 Á Á Ál 2 3 > l 1 S0 Á Á Ál 2 1 = CM-S0 , then we chose CM-S1 for population P 2 because CM-S1 > CM-S0 . Section 2.3 in Supporting Information gives detailed information about connectivity coefficients used to calculate CM.
Importantly, CM is influenced by the position of the focal population downstream of the sources and on whether the network is either "stepping-stone" (Carr and Reed 1993) or "small world" type (Watts and Strogatz 1998; that is, networks where populations are highly connected with each other). Here, we also used the coefficient of determination (R 2 ) as a metric of effect sizes based on both the raw data and log-transformed recovery times, that is, log(T 50 ) and log-transformed values of CM.

Life history and habitat matrix
The density-dependent coefficient β did not have any important influence on recovery times (Fig. 3) and it is not considered further. The term ω (the product of fecundity and survival: Fig. 2) had a strong effect on the predicted recovery times (Fig. 3). At high ω (= 10 4 ), the predicted recovery times were much shorter and had a lower degree of variation among larval strategies and time of larval release (April vs. August) and regions (Fig. 4). The term ω also seems to influence variability in recovery time within a region. For instance, for Cardigan Bay (Fig. 5), the simulation resulted in short recovery times (T 50 < 25 yr) for all strategies with ω of 10 4 , but at ω of 10 such times varied considerably among larval behaviors or sites (T 50 varied from < 25 to > 100 yr). Overall, the model predicted that increased fecundity and larval survival mitigates the effects of larval strategies and time of release on recovery times. This conclusion is logical because in Eq. 4, the term ω operates on settlement through a multiplicative effect (the effect of the transport coefficients, l j!p , on settlement is increased ω-fold).
In the model, larval behavior drives recovery times through changes in the coefficients of the transport matrix (Fig. 2). Larval behavior had an important effect on recovery time but this effect depended on other predictors (Figs. 3, 4). In most simulations, tidal migration led to short recovery times (T 50 < 25 yr); under a tidal migration strategy, the effects of ω (month of release) were smaller than under other migration strategies, at both the scale of regions (Fig. 4) and within regions (Figs. 5, 6). Passive and diel migration let to regional scale variation in recovery times (Fig. 6, see also Supporting Information Figs. S4-S7) from short (Irish coast: T 50 < 5 yr) to longer times (e.g., Liverpool Bay: average T 50 = 20-40 yr and Cardigan Bay:~4 0-70 yr both ω = 10). Taken together, tidal migration and high fecundity/larval survival (i.e., high ω) minimized the average and the spatial and temporal variability in recovery time.
There was also important regional scale and temporal variation in recovery times (Figs. 3, 4), driven by oceanography. The Irish coast showed consistently short recovery times (Fig. 6); that is, they were short irrespective of tidal strategy and month of larval release. By contrast, recovery in other regions was largely affected by tidal strategy and month of release. In addition, recovery times were slightly shorter in simulations of release in April as compared with those in August (Figs. 4, 5).

Network topology
Network topology refers to the geometric arrangement of the network (Kininmonth et al. 2010) including its nodes and connections. A key question was whether local subsidies and retention coefficients, constituting direct connections to each population, were able to explain the patterns of recovery time. Subsidy and retention were poor as predictors of recovery time (R 2 < 0.10), although the transport coefficients defining such connections determined recovery times by design (Fig. 2). In general, populations characterized by high subsidy recovered rapidly but recovery time varied considerably when subsidy was low (Fig. 7a); examination of scatterplots showed that subsidy (or retention) had a weak relationship with recovery

Giménez et al.
Metapopulation recovery after disturbance time (exception: combination tidal strategy, August release, ω = 10: Supporting Information Figs. S12, S13). By contrast, total connectivity (CT) or the maximum connectivity to the source populations outside the area of extinction (CM) were a strong predictor of recovery times (R 2 > 0.7; Fig. 7b,c). CM represents the connectivity provided by one of the paths considered in CT. Close examination of populations within regions shows how CM reflects the overall topology of the regional networks as modified by larval behavior and month of larval release. For instance, the subnetwork of Cardigan Bay (populations 4-8; Fig. 5) was characterized by a stepping-stone pattern where populations 6-8 (P 6-8 ) usually received larvae from P 1-3 , outside the network, through P 4-5 . For the spring simulation, in spite of differences associated with larval behavior, all connectivity coefficients were moderate to high (mostly > 10 −7 ), and populations recovered quickly from extinction. However, for the summer simulation, low connectivity coefficients (mostly < 10 −7 ) characterized the passive and diel strategies. In summer, the differences in recovery times among populations were high as extremely low recovery rates were predicted for some locations but not others. For example, for passive dispersal, P 7 has a long predicted recovery Giménez et al. Metapopulation recovery after disturbance rate (T 50 > 100 yr) compared with the adjacent site P 5 (rapid recovery: T 50 < 5 yr) despite similar levels of subsidy (10 −7 vs.~10 −8 ). However, P 5 is strongly and directly connected with a population (P 3 ), located outside the area of extinction (connectivity P 3 -P 5~1 0 −7 ) as compared with P 6 (connectivity P 3 -P 6 < 10 −18 ). Stepping-stone patterns are not reflected in the subsidy but CM incorporates both the subsidy from populations that are connected to the source and the position in the network. Low values of CM are found in populations that are not directly connected to the sources outside the region of extinction and instead are organized stepping-stone patterns.
Stepping-stone patterns were not present in all regions; for instance, the subnetwork of Anglesey (Supporting Information Fig. S5) resembles a "small world" type, especially for the passive strategy (i.e., with populations closely connected with each other and with the source). Whether the network structure resembled a stepping-stone or a small world type also depended on the larval strategy (Supporting Information Figs. S6, S7: Liverpool Bay and West Ireland): in particular, the tidal strategy which led to short recovery times gave rise to several connections between the source and the populations located inside the region of extinction.

Discussion
We have carried out a modeling exercise in order to better understand the role of life history strategies, the habitat matrix, and the network topologies in determining the capacity of local populations to recover from extinction. Recovery in some cases took many decades; this is consistent with empirical observations of metapopulations over large regions (> 20 yr : Lessios 2016); however, in our case, they may reflect the fact that we simulated recovery based only on a single larval release event each year. Recovery strongly depended on the interactive effect of trait combinations and the temporal and spatial variations in the habitat matrix. Real spatial or temporal patterns will look blurred compared with model outputs because many species will produce several batches of larvae each year and hence profit from temporal changes in hydrodynamic conditions. In addition, effects of larval strategy might be blurred because swimming speeds may not be constant over time or may vary intraspecifically according to the physiological state of larvae (e.g., body size, nutritional condition). Recovery will also be driven by Allee effects and the extent to which regional extinction lead to regime shifts in the local habitat (see Lessios 2016 for discussion). Given these points, we use our results only as a guide for understanding processes driving recovery from the standpoint of larval survival and dispersal. In our model, the characteristics of the habitat matrix were driven by oceanographic conditions and larval survival (manipulated through the term ω). The characteristics of the habitat matrix contributed to the properties of the regional subnetwork, which in turn were driving recovery time; the correlation between recovery times and the connectivity to the source populations outside the region of extinction suggested that recovery is driven by whole network properties.

Life history strategies
High fecundity and tidal transport minimized recovery times. High fecundity is known to increase connectivity (Johansson et al. 2012;Treml et al. 2012) but fecundity varies considerably among marine species (Ramirez Llodra 2002). Under the conditions of the model, only the most fecund species (producing~10 6 larvae per female) would be able to produce sufficient survivors to ensure quick recovery (T 50 < 5 yr) irrespective of region and season. Less fecund species may however avoid low recovery rates by producing multiple broods and releasing larvae over a protracted period.
Our model predicts that not only high fecundity but also tidal migration can lead to short recovery times. In the sea, larval vertical migration varies across species (Shanks and Brink 2005;Epifanio and Cohen 2016) even within the same habitat (Lindley et al. 1994;Garrison 1999). While tidal migrations usually promote transport, diel migration can promote retention (e.g., Queiroga et al. 2007;Shanks 2009); hence, one would expect that variations in larval strategies among species would lead to important intraspecific variation in recovery times. Our models however do not include the fact that diel migrations reduce mortality by predation (Hays 2003) and would indirectly contribute to shorter recovery times.
The finding that dispersal strategies can be as important as high fecundity is relevant from the evolutionary standpoint. First, K-strategists, characterized by low fecundity, but possessing the appropriate pattern of larval migration, would be able to quickly recolonize habitats postdisturbance. Second, larval behavior may provide an evolutionary routes to maximize connectivity (and hence fitness), free from constraints associated with increased fecundity. Because of energetic constraints, fecundity is linked through a trade-off with per-offspring investment (Marshall et al. 2007;Kindsvater and Otto 2014). Because maximizing fecundity will come at the costs of reducing larval survival, species may not be able to increase fitness through increments in ω (= fecundity × larval survival). However, if specific behavioral strategies are not linked to energetic constraints, then fitness can be increased through increments in larval transport coefficients. The phenotypic links characterizing the life histories of marine organisms (Marshall and Morgan 2011) predict that selective pressures in the adult habitat may well drive the evolution of larval behavior, or alternatively that selection for specific larval behaviors (e.g., for diel migration) may drive the evolution of fecundity and offspring size. Advance in this field needs information about the covariation between per offspring investment, larval swimming speed, and behavioral strategies.
The effect of life history variation on recovery times modeled here is also relevant to understand the structure of metacommunities. When local communities are structured through priority effects, such trait combinations may determine which species are established first, and which may inhibit or promote the establishment of a second species (Almany 2003;Chase 2007). Depending on trait combinations, conditions of the habitat matrix may lead to trait-based environmental filtering (Lebrija-Trejos et al. 2010). Environmental filtering would occur because disturbance would select species according to traits promoting rapid recovery, for example, high fecundity (Ponge 2013;Seifan et al. 2013) or specific migration strategies. Because the importance of migration strategies depends on properties of the habitat matrix (e.g., oceanographic conditions in our specific case), our models indicate that structure of metacommunities may depend on landscape-dispersal interactions. Landscapedispersal interaction has been highlighted as an overlooked but potentially important driver of metacommunity structure (Ryberg and Fitzgerald 2016).

Role of habitat matrix and network topology
Our model outputs reinforce the finding by others concerning the role of the habitat matrix in driving recovery (Hanski 1999;Joly et al. 2001;Haynes and Cronin 2004;Fisher et al. 2005;Goodsell and Connell 2005). A component of the habitat matrix is given by those factors driving larval survival (Shima and Swearer 2009), incorporated as η, in the term ω. Our findings are consistent with arguments in Paris et al. (2007) on including overall larval survival rates to understand the role of connectivity on population recovery and highlight the necessity to quantify larval survival in the field (Vaughn and Allen 2010, but see White et al. 2014a,b).
Species traits, the quality of the habitat matrix, and the spatial configurations of the local populations contribute to the characteristics of the network topology through effects on the transport coefficients. We found that effects of network topology on recovery were captured in the connectivity to the sources outside the region of extinction, either as the total connectivity or as the path providing the maximum connectivity (CM) because high values of CM would occur more frequently in small world networks than on those characterized by steppingstone patterns. The fact that CM had a much higher predictive power than subsidies and retention coefficients suggest focus should be in understanding the ecological factors driving "emergent" network properties rather than (only) on conditions surrounding local populations. Network topology varied at two scales, defining regions linked by weak connections (see, e.g., Cowen et al. 2006 as a similar example) and groups of weakly connected locations within regions. Reduced larval connectivity among regions and some stepping stone patterns within regions occurred at the time of formation of frontal zones and thermoclines in summer, which act as conduits of larval transport ). On the other hand, strong currents promoted a small world type of network in East Ireland. Hence, it seems that the nature of the habitat matrix is such that it leads to context-dependent recovery.
Overall, we have found high levels of contingency in attempting to determine which biological and physical factors drive recovery from extinction. Key drivers of contingency were the temporal variation and spatial heterogeneity of the habitat matrix driven by variation in hydrology; the effect of habitat heterogeneity on recovery was exacerbated under low fecundity or high larval mortality. Having the right larval behavioral strategy may be as important as high fecundity or low mortality rates in achieving quick recovery times in a heterogeneous habitat. Quick recovery was obtained by a strategy providing sufficient connectivity but limited overdispersion (exemplified by the tidal strategy) which led to small-world type networks. A strategy maximizing retention (here exemplified by diel vertical migration) at the expense of larval connectivity may not ensure quick recovery, unless it is coupled with high fecundity or low larval mortality. Mixed approaches, based on the application of general metapopulation models to situations characterized by realistic seascapes (exemplified by the Irish Sea), might further contribute to understand the mechanisms driving recovery from extinction after disturbance events. This is also relevant for conservation and for understanding invasion dynamics. For instance, the optimization of networks of protected areas, which depend on understanding patterns of retention and connectivity ), would require knowledge of species trait combinations. In addition, conservation should address the quality of the habitat matrix (see also Shima and Swearer 2009), which depends on stressors (e.g., pollutants, toxic algae: Shaber and Sulkin 2007;Vasas et al. 2007;predatory jellyfish: Purcell 2011;Lee et al. 2013). These are the ecological and anthropogenic factors that codetermine whole network topological properties and influence recovery from extinction.