Mapping the resilience of chemosynthetic communities in hydrothermal vent fields

Hydrothermal vent fields are vulnerable to natural disturbances, such as volcanic activity, and are currently being considered as targets for mineral mining. Local vent communities are linked by pelagic larval dispersal and form regional metacommunities, nested within a number of biogeographic provinces. Larval supply depends on the connectivity of the dispersal networks, and affects recoverability of communities from disturbances. However, it is unclear how the dispersal networks contribute to recoverability of local communities. Here, we integrated a population dynamics model and estimation of large scale dispersal networks. By simulating disturbances to vent fields, we mapped recoverability of communities in 131 hydrothermal vent fields in the western Pacific Ocean. Our analysis showed substantial variation in recovery time due to variation in regional connectivity between known vent fields, and was not qualitatively affected by potential larval recruitment from unknown vent fields. In certain cases, simultaneous disturbance of a series of vent fields either delayed or wholly prevented recovery. Our approach is applicable to a dispersal network estimated from genetic diversity. Our method not only reveals distribution of recoverability of chemosynthetic communities in hydrothermal vent fields, but is also a practical tool for planning conservation strategies.

Chemosynthetic communities in hydrothermal vent fields (HVFs) demonstrate adaptations to extreme environments [1][2][3][4][5] and provide various ecological services [6][7][8][9] . Local vent communities are linked by pelagic larval dispersal and form regional metacommunities, nested within a number of biogeographic provinces [10][11][12] . The incidence of disruptive natural disturbances to vent communities can range from several decades to several hundred years 13 . Faunal adaptations for colonising new vent fields are thus important aspects of the sustainability of these communities, especially since neighbouring vent fields are often separated by 10 s or 100 s of kilometres.
Recoverability, or resilience, refers to persistence of ecosystems in the face of natural or anthropogenic disturbances 14 . It can be quantified as the recovery time to the original state after disturbance 15 . For chemosynthetic communities in HVFs, observations on recovery from disruption caused by volcanic activities suggest that most of the diversity and biomass recovered within five years after the disturbance [16][17][18][19] , reviewed by Gollner et al. 20 . For example, total mega-and macrofaunal species richness at the vents in Juan de Fuca Ridge reached 75% of the pre-disturbance values three years after the 1998 eruption 18 , and 90% two years after the 1993 eruption 16 , representing about 30-60% of species from the larger regional species pool. At the East Pacific Rise (EPR), total megaand macrofaunal species richness reached 69% of pre-disturbance values 4.6 years after the 1991 eruption 17 . After the EPR 2006 eruption, the recovery reached 55% for macrofaunal and 48% for meiofaunal species after 4 years 19 , with 39% of the macro and 42% of meiofaunal species returned. Recoverability varies significantly among communities in active vents, inactive vents, or within the vent periphery 20 . However, differences in recovery time among active vent fields are not well documented. Limited accessibility to HVFs (most of them are remote and in depth of more than 1,000 m) and rare opportunities to observe natural disturbances in less active volcanic areas (e.g., slow-spreading ridge systems and arc-backarc basin systems) hinders the ability of researchers to assess recoverability of chemosynthetic communities in the environment. To address this question, we integrated estimation of dispersal networks of HVFs with a differential equation model that describes the recovery of disturbed populations. We used estimated larval dispersal between 131 HVFs in the western Pacific Ocean 12 . The estimates were based on a physical model of deep-ocean circulation that was validated through a deep-ocean profiling float experiment and considered temperature dependency of larval development that can control duration of pelagic larval stages. For example, larvae that use shallower depths can disperse further because of fast ocean currents and duration of the pelagic larval stage would be shortened when water temperature is higher 21 . Thus, dispersal distance will depend both on the speed of ocean currents and expected duration of larvae at a given depth. We assumed a dispersal depth of 1,000 m in our analysis based on the published data on the hydrographic structure of water columns obtained from observation and simulation models [22][23][24] . The biogeographic studies of larval dispersal at hydrothermal vents have a consensus on the transport mechanism by deep-sea advection and effluent layers 25,26 , where physicochemical parameters, physiographic features of a region and seafloor topography are recognized as possible barriers to dispersal 27 . The water mass below 1,000 m depth generally has stable physicochemical parameters, such as temperature and salinity, and may offer suitable conditions for larval survival in advection above the vent area and lateral transport in the effluent layer of the deep-sea. For example, there is a discontinuity in water temperature and salinity at 500-700 m in the Okinawa Trough, and at 1,000 m in the western and southern Pacific Ocean 28 . Moreover, two chemosynthetic communities found at a difference of more than 1,000 m depth had significant differences in their community structure 29 , suggesting the effect of environmental barriers (but see 30,31 for recent findings).
The vent fields were separated into seven regions by grouping them in terms of their connectivity, i.e., each group had no interconnections ( Table 1). The dispersal networks were implemented as a dispersal matrix A in the equation (1) below (see Methods).
Here, K i is the carrying capacity, (i.e., equilibrium population size of vent field i), r is the reproduction rate defined as the number of larvae that one individual produces per year and A ij is the dispersal rate that a larva produced at vent field j will migrate into i per year, where A ii corresponds to self-recruitment. Our model does not include the duration of pelagic larval phase as a parameter because it is included in the calculation of A 12,21 . We assumed that x i is the population abundance of a species or group of species that share the same niche in vent field i. The species or group of species is assumed to distribute across all vent fields with sufficient abundance. This assumption is realistic because vent communities frequently have a dominant taxon that constitutes most of the biomass at a regional scale, e.g., Bathymodiolus mussels, Shinkaia squat lobsters and Alviniconcha gastropods 29,32,33 , although the dominance-diversity relationships may depend on the environmental conditions such as fluid-flux intensity and sediment types 34 . Amongst different regions however, it would be reasonable to interpret x i as the abundance of different species or group of species that accounts for a similar proportion of biomass and having the same growth and dispersal characteristics. We do not directly consider the effect of disturbances on biological diversity because our model accounts for only one species. However, we expect that recoverability of these representative species would be a proxy for the recoverability of other infrequent community members, thus we regarded it as the recoverability of the community as a whole. For example, this is supported by observations of recovery after eruptions [16][17][18][19][20] which showed a concordance of recovery in total organismal density and species richness. While we used the simplest model for our analysis, availability of more detailed data sets that include, for example, biomass, age structure or trophic interactions, would make more extensive models like those used for fishery stock assessments 35 applicable and may provide more detailed insights. By considering previous observations [16][17][18][19][20] , we defined τ i as the mean recovery time of a vent field i for the ensemble of various spatial distributions of K i , where recovery time is the time required for x i to recover 75% of its original abundance (equal to K i by definition) after it was temporarily reduced to zero. To calculate τ i , we assumed  Table S5. Here, "Mean in-degree" is the mean number of incoming links per vent fields within a region, "Mean self-recruitment" is the mean of self-recruitment (A ii ) within a region and "Mean between vents recruitment" is the mean of A ij (i ≠ j) within a region. "Mean selfrecruitment" and "Mean between vents recruitment" is shown as multiplied by 10 6 to facilitate understanding. that K i is assigned from a uniform distribution P K = 2 q K 0 where q is a uniform distribution between −1 to 1 and K 0 = 10000, and all vent fields except for i are in equilibrium. However, the following results are independent of K i if K i values are identical across vent fields (equation (4)). Furthermore, they can be referred to as the mean recovery time in the ensemble of potential spatial heterogeneity in K i as long as K i of all vent fields follows the same probability distribution (see Methods). For simplicity of notation, in what follows, we regard the set of τ i s for all vents fields as {τ i }. Here, we do not consider the effect of time on K i . However, in nature, K i will not be fixed in time. Our result may underestimate recovery time if frequency of a species is highly constrained by the amount of suitable habitat that increases slowly along with the re-establishment of vent fields after disturbances. For example, the loss and recovery of mussel beds may restrict abundance of associated small invertebrate species 36 .

Results and Discussion
We obtained the spatial distribution of τ i as in Fig. 1. To be consistent with previous observations [16][17][18][19][20] , we set r = 17.4 so that the median of {τ i } was standardised to five years ( Supplementary Fig. S2). Hence, we used five years as a reference representing a standard time scale of recovery of HVFs from large disturbances. Mean τ i in Izu-Bonin was 131 years (Table 1) and the shortest τ i (Mokuyo Seamount) was 56 years ( Fig. 1; see Supplementary  Table S1 for further information). Mean τ i in Manus-Woodlark, Solomon and Mariana was also longer than 20 years. These regions contained multiple HVFs with τ i longer than 40 years: W Syoyo and SW Syoyo in Mariana, New World Seamount, Edison Seamount, Woodlark Basin Segment 5B, Woodlark Basin Segment 3B and Franklin Seamount in Manus-Woodlark, Starfish Seamount and Tikopia Area in Solomon. P values of the τ i were low, suggesting region specific factors were determinants of recovery time. There was a large contrast between regions with high recoverability (Okinawa and New Hebrides-North Fiji-Lau Tonga) and regions with low recoverability (Izu-Bonin, Manus-Woodlark, Mariana and Solomon; Table 1; see Methods). The inequality in recovery time was dependent on the interplay between the ocean circulation and the spatial distribution of HVFs in each region. Within some regions (Mariana, Manus-Woodlark and NewHebrides-NorthFiji-LauTonga), we found a positive correlation between τ i and the distance of HVFs from the mean (see Supporting Table S2). Hence, the recovery time of a HVF would be biased depending on the region to which it belongs and occasionally also on its locality within the region. We did not include Eva, Kaiwo Barat and Teahitia Vents in our results because these HVFs had no known incoming links and thus were unable to recover based on our assumptions. It may be unrealistic to judge that communities in these HVFs cannot completely recover from disturbances because there would be overlap of species between communities in active vent fields and other environments, such as inactive vents or along the vent periphery 37,38 . Some opportunistic/non-endemic species could colonise a recently disturbed vent field through adult migration rather than being limited by larval dispersal 6 . This will contribute to the recovery of communities in these HVFs (and of course, also affect recovery of communities in HVFs having incoming links). However, we did not consider this because the interrelationship between larval dispersal and adult migration for population recoverability is an issue of complexity within the local community. 95% CI of recovery time (thick and thin lines) is shown in Fig. 2a with τ i (shown in white points). While τ i was mostly explained by the total A ij (ΣA ij ) to HVFs (Supporting Fig. S3), there were two factors that affected 95% CI of recovery time. One was K of disturbed HVF (K i ) and the other was that of other HVFs (K j values). Figure 2a shows the 95% CI when K i was fixed at the mean of P K (i.e., K i = 10000; thick lines) and when it was also assigned from P K (thin lines). Thus, Fig. 2a illustrates relative contribution of K i and K j values to the 95% CI as the length of thick and thin lines, respectively. Because the 95% CI was the result of the ensemble of virtual distribution of K values, these lines represent uncertainty in recovery time when we do not have information on both K i and K j values (thin lines) and when we can specify K i (thick lines). Supplementary Fig. S4 shows how much uncertainty can be reduced by specifying K i . The uncertainty in recovery time decreased substantially when ΣA ij was large because HVFs also had large in-degree ( Supplementary Fig. S5), and the law of large numbers reduced the effect of variation of K j values (i.e., the variation of recovery time was explained by the variation of K i itself). It is also worth noting that HVFs in Okinawa and several HVFs in Manus-Woodlark (e.g., DESMOS Cauldron, PACMANUS field and Solwara 13, etc.) had larger uncertainty for K j s than HVFs in New Hebrides-North Fiji-Lau Tonga having similar ΣA ij (e.g., Mata Fa, Mata Fitu and Mata Ono, etc.) (Fig. 2a), which was also explained by the difference of in-degree among these HVFs (Supplementary Fig. S5).  (2) with u = 0). Here, thick and thin lines indicate the 95% CI when K i was fixed at the mean of P K (i.e., K i = 10,000) and when it was also assigned from P K , respectively. (b-d) The same result for u = 1, 2, 3, respectively. In (a-d), points with light colors indicate τ i for other u values. For u = 0, 1, 2, 3, we set r = 17.4, 16.7, 16, 15.6, respectively, to keep median of τ i s as five years. (e) The reduction of mean recovery time within a region, defined as ∆τ + = τ + − τ (u 1) (u 1) (u) i r eg i r eg i r eg is shown. Here, 〈·〉 reg is used to emphasize that it is the mean of recovery time within a region. Several unknown HVFs may exist in each region and provide larvae to known HVFs. The effect of unknown HVFs can be introduced as in equation (2). Here, we introduced a parameter u as the number of unknown HVFs in a region, and set r so that median of {τ i } was kept identical to five years for each u (see legend of Fig. 2).
Addition of unknown HVFs generally reduced both mean and variation of recovery time (Figs 2, S4 and Supplementary Table S1). In particular, unknown fields had a large effect in the regions with large τ i , such as Izu-Bonin, Solomon and Manus-Woodlark (Fig. 2e). However, in Mariana, the effect of unknown HVFs was similar to that of Okinawa or Kermadec despite a higher mean τ i in these region (Table 1). This is due to the fact that Mariana has many HVFs with high inter-connectivity (in-degree was 6.56, the second highest among all regions), but a low dispersal rate per link (mean between vent recruitment was 232.91 10 −6 larvae/adults/year, the second lowest among all regions). Hence, addition of several unknown HVFs had only a small effect on τ i in Mariana. Mean τ i of Izu-Bonin, Manus-Woodlark and Solomon would be largely overestimated considering the potential effect from unknown HVFs. However, their P values still suggested inequality in recoverability among regions (Supplementary Table S3). Hence, presence of unknown HVFs would not qualitatively change the imbalance of recovery time.
On-going development of deep-sea resource mining technologies [38][39][40] raises the possibility that natural and anthropogenic factors will simultaneously disturb multiple HVFs in a region. We evaluated recovery time for disturbance on multiple HVFs as τ C , where C = {i, j, …} is a possible combination of disturbed HVFs in a region. We tested all combinations if the total number of combinations was less than 2,000, and, if the total number exceeded 2,000, we randomly selected 2,000 combinations without duplication. The procedure to calculate recovery time is the same as previous analyses except the abundances of populations in HVFs included in C were simultaneously reduced to zero. τ C is the largest recovery time of HVFs included in C, i.e., τ c = max({τ i } i∈c ). For simplicity, we did not consider variation of K in this analysis and assumed K = 10,000 for all HVFs, i.e., we focused on the effect of disturbances on the mean recovery time.
In Fig. 3, τ c s are shown as a function of the number of simultaneously disturbed HVFs with the 95% CI. The mean recovery time monotonously increases with the number of disturbed HVFs, while the pattern of the increase is different depending on differences in connectivity of each region. We found that some combinations extremely delay recovery compared to others. For example, when more than three HVFs were disturbed in Okinawa, τ C was larger than 20 years if all of Hatoma Knoll, Dai-Yon Yonaguni Knoll and Irabu Knoll were included in the disturbance (this is shown in panel (a) of Supplementary Fig. S6 where a point indicating recovery time longer than 20 years first appears when three HVFs are simultaneously disturbed), while that of other combinations were smaller than 10 years. This is clear from the dispersal matrix in this region (Supplementary Table S4) showing that if all three HVFs are disturbed, recovery will depend on dispersal from Iheya Ridge to Hatoma Knoll and Dai-Yon Yonaguni Knoll, and dispersal via these links is more than 10 times smaller than the mean between vent dispersal in this region (Table 1). In addition, Mariana, Manus-Woodlark and Solomon had combinations that prevent recovery (Supporting Fig. S7). In these cases, no larval supply from the undisturbed HVFs to the disturbed HVFs remained after the disturbance. For example, in Manus-Woodlark, no larval supply remained if Woodlark Basin Segment 5B, Woodlark Basin Segment 3B and Franklin Seamount were simultaneously disturbed. Thus, when disturbance to multiple HVFs is considered, inclusion of some HVFs can either slow down or completely prevent recovery. Further investigation on this observation would benefit by incorporating network analysis, e.g. use of centrality measure 41 to distinguish source and sink HVFs, and would provide relevant insight for the application of our approach to management strategy planning. This application would be further strengthened if knowledge on the frequency of natural disturbances (e.g., cycles of volcanic and hydrothermal activities) in each region could be integrated.
A method using genetic diversity to estimate the dispersal matrix for a species distributed across multiple HVFs has been proposed 42,43 . Using the dispersal matrix of a galatheid crab, Shinkaia crosnieri, in Okinawa Trough (Watanabe et al. unpublished data), we found some consistency between the dispersal matrix estimated from the bio-physical model 12 and the matrix estimated from genetic diversity in terms of recovery time (supplementary material). It is worth noting that the latter does not explicitly include information about ocean circulation. To date, 700 HVFs have been found worldwide 44 . 57% (399) of them are located on the mid-ocean ridge system and the rest are on arc-backarc systems. 88% (249 among 283) of arc-backarc systems are distributed in the Pacific Ocean, of which this study included 52%. In future studies, researchers will reveal dispersal networks in a broader region of the world ocean and range of taxonomy. The coverage of this analysis can be extended provided data are available. Our approach will allow us to map the resilience of different species, enable comparative studies on the relationship between resilience and species traits and, when the model is modified to multispecies systems, it will further help us to understand community dynamics, such as succession. Thus, we think that this study is also an important step towards understanding the processes that form and sustain communities in HVFs. Conservation strategies for deep-sea ecosystems may include creation of new habitats as well as design of ecologically protected areas 38,[45][46][47] . Our methodology would reasonably direct these efforts by predicting the outcome of various strategies in any region worldwide. However, it should be mentioned that there is still a relative lack of empirical data from disturbance-recovery studies that can support our results. All previous studies have been carried out on either the Juan De Fuca Ridge 16,18 or the East Pacific Rise 17,19 . Because both are in a fast spreading ridge where biological communities are frequently disturbed (e.g., ~15 years 48 ), it would be controversial whether these results represent recovery of HVFs in slow spreading ridges or arc-backarc basins where disturbance to communities are relatively infrequent 20 . However, our results still suggest substantial difference of recovery time among HVFs, which can span two orders of magnitude, highlighting the importance of understanding connectivity among HVFs to assess their recoverability.

Methods
Population dynamics. We assumed that populations in a vent field are supported by recruitment from the same vent field (self-recruitment) as well as recruitment from other vent fields. For equation (1), we assumed that K i is determined by the total amount of resource supply in vent field i, and limits maximum population size. rA ii x i represents the larval supply via self-recruitment per year and rΣ i≠j A ij K j is the total larval supply from other vent fields per year, assuming that other vent fields are in equilibrium (x j = K j ).

Estimation of the effect of unknown vent fields.
To estimate the effect from unknown vent fields, we assumed that each unknown vent field supplies larvae to all known vent fields in the region to which it belongs and the amount of larval supply per each link is identical to the mean between patch recruitment in the region ( 0  ; Table 1). By adding the effect from unknown vent fields s as u 0  to equation (1), we obtain, where,  0 is the mean between patch recruitment in the region that i belongs to, and u is the number of unknown vent fields.
Analytical calculation of recovery time. Recovery time can be simply obtained by numerically solving equations with appropriate numerical integration scheme such as the Runge-Kutta method. However, it is more convenient to use their analytical solution if only one vent field is disturbed. For example, for equation (1), we obtain, Moreover, if K i is identical across vent fields, i.e. K i = K, by substituting x i = αK, we obtain, Here, α is the abundance criteria (ratio) to judge recovery, which we set as 0.75 in this paper. Equation (4) shows that recovery time is independent of K if it is identical across vent fields.

Calculation of P values for mean recovery time. P values for the mean recovery time of each region
were obtained by using {τ i } as the empirical distribution. Here, for each region, we simply repeated random sampling from the empirical distribution up to the number of vent fields in the region, then calculated the mean recovery time. This procedure was repeated 10,000 times providing 10,000 bootstrap samples for each region from which we evaluated the P value of the mean recovery time calculated from actual τ i values.
Spatial heterogeneity in population size. If carrying capacity K i of all vent fields follows the same probability distribution Φ, equation (4) holds for the mean of Φ, μ K , i.e., for μ K , recovery time is independent of Φ. of each vent field. As a base line expectation, the best approach to implement spatial heterogeneity in K i is to consider the ensemble. Hence, in this paper, we identify {〈T i 〉} i=1,…,N with τ = … { } i i 1, ,N (we simply denoted τ i as τ i in other parts of this paper).