Natural mega disturbances drive spatial and temporal changes in diversity and genetic structure on the toadfish Aphos porosus

Natural disturbances can modify extinction-colonization dynamics, driving changes in the genetic diversity and structure of marine populations. Along Chilean coast (36°S, 73°W), a strong hypoxic-upwelling event in 2008, and a mega earthquake-tsunami in 2010 caused mass mortality within the Aphos porosus population, which is a vulnerable species with low dispersal potential. We evaluated the effects of these two major disturbances on the diversity and spatial-temporal genetic structure of Aphos porosus in two neighboring areas that were impacted on different levels (High level: Coliumo Bay; Low level: Itata Shelf). Thirteen microsatellites (from 2008 to 2015) amplified in individuals collected from both locations were used to evaluate the effects of the two disturbances. Results showed that after the strong hypoxic-upwelling event and the mega earthquake-tsunami, Aphos porosus populations exhibited lower genetic diversity and less effective population sizes (Ne < 20), as well as asymmetries in migration and spatial-temporal genetic structure. These findings suggest a rise in extinction-recolonization dynamics in local Aphos porosus populations after the disturbances, which led to a loss of local genetic diversity (mainly in Coliumo Bay area impacted the most), and to greater spatial-temporal genetic structure caused by drift and gene flow. Our results suggest that continuous genetic monitoring is needed in order to assess potential risks for Aphos porosus in light of new natural and anthropogenic disturbances.

intensity, and extent), and the life-history traits of the species affected (i.e.dispersal potential, population size, and generational time) 11 .
Along the coast of central-southern Chile (34°-39°S), two strong natural disturbances impacted an extension of approximately 500 km of coastline.The first disturbance was a hypoxic-upwelling event (< 0.5 ml O 2 l −1 in the water column), which occurred on January 3, 2008, and caused massive mortality among pelagic and benthic organisms, mainly fish species in the Coliumo Bay 12,13 .The second disturbance was a mega earthquake-tsunami (8.8 Mw, the 6th strongest earthquake since 1900), which occurred on February 27, 2010, impacting the coast along the Central-South of Chile.This event uplifted part of the coastline by more than 3 m, generating a tsunami with waves up to 14 m that devastated extensive areas of the coastal seabed and modified benthic and coastal environments, causing a massive mortality among different species in intertidal and subtidal habitats, especially in the Coliumo Bay [14][15][16] .Both disturbances, presented different degrees of magnitudes throughout the affected area, creating a complex patchy landscape, where the Coliumo Bay was impacted the most, and the surrounding northern area of the Itata shelf to a lesser extent.In the Itata Shelf area, the hypoxic-upwelling event decreased the dissolved oxygen on the surface, but the values remained at normoxia levels (> 2 ml O 2 l −1 ) on the sea floor.In this area, the tsunami during the mega-earthquake was less destructive, with wave heights almost half of those in Coliumo Bay 17 (Fig. 1).
The Aphos porosus toadfish (Valenciennes, 1837) was one of the species severely affected by the 2008 hypoxicupwelling event, as well as by the 2010 mega earthquake-tsunami 18 .Aphos is a monotypic genus, and Aphos porosus is the only species of the Batrachoididae family present in Chile 19 .This benthic-demersal fish is distributed in the Southwest Pacific (3°-53°S) from Puerto Pizarro in Perú to the Strait of Magellan in Chile, and from the coast to the demersal shallow environment (0-100 m) [20][21][22] .Aphos porosus has a generation time frame

Materials and methods
A total of 565 Aphos porosus individuals were collected from the sea floor between the years 2008 and 2015 in Coliumo Bay and Itata Shelf areas mentioned above.Benthic sampling was carried out every 3 months using a modified Agassiz trawl (1 m wide × 1 m long × 30 cm high, lined with 5 mm "knot to knot" netting).All samplings were carried out using R.V. Kay-Kay I and R.V. Kay-Kay II (Department of Oceanography, University of Concepción).The population density (individuals m −2 ) was standardized according to the number of fish collected in each tow, in relation to the swept area 12,16,18,24 .All individuals collected were genotyped using 18 microsatellite markers, previously developed by Silva et al. (2015)  25 .All experimental protocols were approved by bioethical committee of the Universidad Católica de la Santísima Concepción, Chile.Project Fondecyt Regular N° 1130868.All methods used in this study were carried out in accordance with Chilean bioethical guidelines and regulations.A physical euthanasian protocol for moribund collected fish was developed according to the "Manual de Normas de Bioseguridad y Riesgos Asociados de FONDECYT-CONICYT " ("Manual of Biosecurity Norms and Associated Risks by FONDECYT-CONICYT) and its references.No chemical product was used during physical euthanasia, and no laboratory experiments were carried out with collected fishes.Where applicable, all methods were reported in accordance with ARRIVE guidelines 2.0.
Microsatellite genotyping and detection of amplification errors were performed locus by locus using the automatic procedure with manual correction in Geneious v7.1 26 .Presence of null alleles and other genotyping errors were evaluated with MICRO-CHECKER v2.2.3 27 , by comparing observed genotypes to the distribution of randomized genotypes.Deviations from Hardy-Weinberg equilibrium (HWE) and Linkage disequilibrium (LD) were evaluated for each locus by area and year using a Markov chain algorithm 28 with dememorization = 1000, batches = 100, Iterations per batch = 1,000 in GENEPOP v4.7 + 29,30 .Due to multiple comparisons, we applied a Bonferroni correction to HWE and LD p-values, with a new p-value = 0.004.
Deviations from neutral expectations were evaluated with the Bayesian method implemented in the BayeScan v2.1 31 , using Number of pilot runs = 20, Burn-in length = 50,000, Number of outputted iterations = 100,000, and Thinning interval size = 10.To correct for multiple testing, the program computes q-values based on the posterior probability for each locus.q-values < 0.05 and α-values significantly > 0 suggest selection diversification, whereas q-values < 0.05 and α-values significantly < 0 suggest balancing or purifying selection.
Intrapopulation genetic diversity estimated over all loci by area and year were evaluated as follows: no.Alleles (Na) = number of different alleles; No. Alleles Private (Nap) = number of alleles unique to a single population; No. Alleles Common (Nac) = number of locally common alleles (frequency > = 5%) found in 50% or fewer populations, observed heterozygosity (Ho), and unbiased expected heterozygosity (uHe).All indices were estimated in GenAlex v6.5 32 according to Hartl and Clark (1989), as well as Peakall and Smouse (2012) 32 .Comparisons between areas, as well as before and after disturbances years (when possible) were made with the two-sample t-test, assuming unequal variances in Excel, and applying a Bonferroni correction to solve issues resulting from multiple comparisons.Effective population sizes (N e ) for both areas were estimated with NeEstimator v2.1 33 , using Model II of the temporal method proposed by Nei and Tajima (1981) 34 , with a Jackknife confidence interval of 95% and P Crit = 0.020, assuming a generation timeframe of tree years according to a previous study 18 .The individuals collected in 2008 from Coliumo Bay were eliminated from this analysis, so as to allow for comparisons between both areas, given that the number of individuals collected in 2008 on the Itata Shelf were insufficient for genetic analysis.Thus, individuals from 2009 were considered generation 0, and whose from 2015 were deemed generation 2.
To evaluate spatial and temporal genetic structures, a Principal Coordinate Analysis (PCoA) was carried out with GenAlex v6.5 32 , which provided a view of genetic differentiation between the areas under study (i.e., Coliumo Bay and the Itata shelf) as well as over time (2008-2015).Then, a clustering analysis was carried out with STRU CTU RE v2.3 35 , so as to infer differences in population structure between the two areas and years.A mixed www.nature.com/scientificreports/model was used, with k = 1 to k = 15, each with 10 replicates, burning = 50,000 steps, repeat = 1,000,000 steps, and LocPrior = 0. k was selected based on the highest Delta K 36 value in the Structure Harvester 37 .Replicas of the cluster analysis were aligned by using CLUMPP v1.1.2 38, and plotted with DISTRUCT 39 .The number of genetic clusters or populations (k) inferred with STRU CTU RE was then contrasted with an Analysis of Molecular Variance (AMOVA) by using No. of different alleles as a distance method (F ST ), No. of permutations = 999, considering 14 populations within 3 groups, where 1 = Coliumo Bay (with the years 2008-2011 nested), 2 = Coliumo Bay (with the years 2012-2015 nested), and 3 = Itata Shelf (with the years 2009-2015 nested) with Arlequin v3.5 40 .In addition, paired F ST tests were estimated among 14 populations using No. of different alleles as a distance method (F ST ), No. of permutations = 999, and applying Bonferroni corrections to solve issues resulting from multiple comparisons.F IS was estimated for each area and year, following the method proposed by Weir and Cockerham  (1984) with GENEPOP v4.7.5 29,30,41 .
Gene flow between Coliumo Bay and the Itata Shelf, as well as recent migration rates per generation (m), and posterior probabilities for individual ancestral distributions of immigrants were estimated with a Bayesian method implemented on BayesASS v3.0 42 , setting the following parameters: Generator seeds s = 100, iterations I = 10,000,000-100,000,000, burning = 1,000,000, sampling intervals n = 100, allelic frequencies a = 0.30, inbreeding coefficient f = 0.50 and migration rate m = 0.10.Convergence of the CMCMs were reviewed in Tracer v1.7 43 .

Results
Among the 18 loci, 5 presented null alleles and were removed from the analysis, including 118 individuals presenting amplification errors.Therefore, a total of 13 loci and 447 individuals collected between 2008 and 2015 were included in the analyses (305 individuals from Coliumo Bay, and 142 from the Itata Shelf; see Table 1).HWE deviations were observed in Coliumo Bay and the Itata Shelf for all years, except in 2010 on the Itata Shelf.Coliumo Bay showed larger deviations from HWE than the Itata Shelf, which increased between 2014 and 2015 (see Supplementary Table 1).LDs were also observed in Coliumo Bay for all years, except 2011, while on the Itata Shelf they were observed only for 2009, 2014, and 2015 (see Supplementary Table 2).Furthermore, 8 of the 13 loci were shown to be under balancing or purifying selection (see Supplementary Table 3).
Gene flow between Coliumo Bay and the Itata Shelf during the years after 2008 ("the hypoxia year"), was asymmetric from Coliumo Bay to the Itata Shelf.In Coliumo Bay, the rate of immigrants per generation from the Itata Shelf was nearly zero, except in 2013, where it increased to 0.3.On the Itata Shelf, the rate of immigrants per generation from Coliumo Bay was equal to 0.3, except for 2013 when gene flow decreased to 0 (Table 1).

Discussion
Natural strong disturbances, including hypoxic-upwelling events, earthquakes, and tsunamis, have resulted in mass mortality among marine populations, increasing drift strength and affecting gene flow, causing the loss of genetic diversity as well as changes in population structure [2][3][4]10,44 . Undertanding the relationship between demographic and genetic changes over time is of critical importance to predict resistance in at-risk populations, particularly in spatially structured populations (i.e.metapopulations) where spatial arrangements of local populations can modulate both demographic and genetic changes 45 .Considering this, we evaluated the effects of the strong 2008 hypoxic-upwelling event and the 2010 mega earthquake-tsunami on the genetic diversity and genetic structure of the Aphos porosus toadfish in two nearby areas, namely Coliumo Bay and the Itata Shelf, which were affected by the two disturbances to different magnitudes of impact.Results showed that genetic diversity and effective population size of Aphos porosus were lower after the two disturbances, mostly in the Coliumo Bay area, which was impacted the most, and exhibited greater spatial-temporal genetic structure and asymmetric migration during the years analyzed.Biological and ecological traits of Aphos porosus, and particular oceanographic characteristics of the inhabiting areas can increase the probability of genetic isolation and structure, limiting the geographical distribution  of this species 23,46 .The latter in combination with the impact from the disturbances, could explain the genetic diversity pattern observed.Population differentiation occurred despite the short distance between Coliumo Bay and the Itata Shelf, as suggested by the different approaches used in this study.Thus, a combined effect between ecological life history traits of the species, oceanographic conditions, and environmental stressors could explain the spatial structure pattern observed.Firstly, Aphos porosus larvae have limited dispersal potential, and adult mobility is also restricted, similar to what has been observed for Opsanus tau, which is another Batrachoididae 47,48 .Additionally, the literature suggests temporal reproductive isolation for Aphos porosus populations on the latitudinal gradient, which is aligned with what has been observed in other regions of Chile 18,23 . Fo instance, previous studies in Chile have shown gaps of 1-2 months in the reproductive season of this species, specifically in populations in Valparaiso (32°S, 71°W), and Coliumo Bay (36°S, 72°W), which are areas separated by approximately 500 km 18,23 .
Secondly, from an oceanographic perspective, studies have suggested that upwelling 49 , which is a physical process that occurs regularly in the geographical area of south Pacific, can act as a hydrodynamic force, affecting connectivity among local populations, acting as a physical barrier to dispersal, and promoting population genetic structure, as observed in Merluccius capensis 9 , Epinephelus andersoni 7 , and Pomacanthus maculosus 8 fishes in other upwelling areas, for instance.Thirdly, the 2008 hypoxic-upwelling, and the 2010 mega earthquake-tsunami caused a drastic decline in the densities of Aphos porosus populations (> 60%) 18 both in the Coliumo Bay and Itata Shelf areas, and in a short time period (i.e., 2 years).This led to lower genetic diversity and a greater spatialtemporal genetic structure, similar to the effects of a bottleneck 50 .Low dispersal from the Itata Shelf to Coliumo Bay, as inferred from our estimates of gene flow, could generate temporal reproductive isolations, and genetic differentiations between both areas.Consequently, compounding factors including a generation timeframe of 3 years, spatial reproductive gaps, and oceanographic barriers could explain the population differentiation observed in these nearby areas, as well as their dissimilar responses to the same major disturbances 11,51 .
Previous studies showed that Aphos porosus population density in Coliumo Bay was sixty times higher than in the Itata Shelf area before the 2008 hypoxic-upwelling, which suggests that Coliumo Bay could be used as a nursey area for different fish species, including A. porosus 12,18 .Population density decreased in Coliumo Bay and the Itata Shelf after both disturbances.Differences between these areas before and after the disturbances could explain the asymmetries observed in migration and spatial genetic structure.Because genetic drift tends to be higher in smaller populations, the Itata Shelf could be more strongly impacted by drift, generating a differentiated population from Coliumo Bay due to rapid fixation of alleles, and the rapid loss of private alleles.Genetic structure has also been observed in the Poricthys notatus toadfish, the sister genus of Aphos, along the Pacific coast the US, as well as in Opsanus tau and O. beta, other Batrachoididae species along the Western Atlantic coast of the US 46,47,51 .A temporal population substructure was only observed in Coliumo Bay, which may have been related to the stronger intensity of both disturbances in this area.Coliumo Bay exhibited greater genetic differentiation between 2008 (hypoxic-upwelling year) and the years 2013 and 2015 (after the mega earthquake-tsunami).Inter-annual substructures observed in Coliumo Bay may reflect the increase in drift and gene flow after disturbances, probably due to the seasonal (Summer-Autumn) recruitment of new cohorts and immigrants from neighboring localities with new alleles.Unfortunately, this hypothesis could not be corroborated, given that other surrounding areas were not genetically sampled, and the method used in this study did not include migration from unobserved populations 34 .Studies on species affected by mega earthquakes and tsunamis reveal that populations with high densities and high larval dispersal are key to maintain genetic diversity, as occurred with the snail Batillaria attramentaria and the sea urchin Mesocentrotus nudus after the 2011 Japan's earthquake 52,53 .Other genetic studies suggest that disturbances like earthquakes and tsunamis can enable the colonization of new lineages 54,55 from less affected areas.For example, another toadfish species from the genera Porichthys and Opsanus have been observed to raft on logs for 25 miles during the summer period, and on sponges during stormy times 56,57 .Moreover, all samples and years revealed evidence of negative F IS values, which signals outbreeding that could be explained on the Itata Shelf by the high number of immigrants observed from Coliumo Bay.Heterozygosity increased by 10% in Coliumo Bay after 2011, which together with negative F IS values could indicate immigration from an unsampled population, but might also signal the effects of genetic drift resulting from the disturbances.The latter is supported by the high population differentiation observed after 2011.
Drift, immigration and substructure could have also caused another pattern observed in Aphos porosus populations, for instance, HWE deviations and evidence of LD, which in general were higher in Coliumo Bay compared to the Itata Shelf.Deviations from HWE and LD can appear to be consequences of immigrants strongly differing genetically from local populations.This could be due to mixing and finite population abundances with higher frequencies of heterozygotes than expected, but it might also be due to substructure [58][59][60][61][62] .
The effects of the 2010 mega earthquake-tsunami on genetic diversity indices have been studied in other organisms, suggesting that not all species were affected similarly.For instance, the red seaweed Agarophyton chilense showed a population decline and a loss of 10-40% in allelic diversity, recovering 2 years post disturbance due to migration 63 .The intertidal crustaceans Emerita analoga and Excirolana hirsuticauda showed low haplotype diversity after the 2010 mega earthquake-tsunami, followed by a recovery after 3 years.As a consequence of this disturbance, E. analoga revealed a lower genetic differentiation and increased gene flow post disturbance, probably due to its high larval dispersal potential.Unlike Agarophyton chilense or E. analoga, Aphos porosus did not show recovery in genetic diversity nor density population after 5 years of the 2010 mega earthquake-tsunami.The low dispersal potential, together with the generational timeframe of 3 years, could limit their ability to recover fast in impacted areas 11,51 .E. hirsuticauda showed a slight increase in population differentiation, and few changes in gene flow.In this case recovery was associated with adult survival and rapid colonization, due to the strong swimming capacity of juveniles and adults 11 .Similar to what was observed for E. hirsuticauda, Aphos porosus showed an increase in genetic structure in Coliumo Bay after the 2010 mega earthquake-tsunami, which could probably be explained by the combined effect of increasing drift, recruitment of new cohorts from survivors, and asymmetric immigration from neighboring areas introducing new alleles.
The effects of disturbances on genetic diversity depend on both the disturbance regime and life history traits of each species, which as a whole determine the resilience capacity of the populations to resist and subsequently recovery 64 .Hence, the response of each species to the same disturbance could be extremely diverse, as shown previously 11,64 .Studies suggest that species with high population abundances, and high dispersal or mobility have higher resilience capabilities 11,51 .For Aphos porosus, which is a species with low dispersal potential, there are no previous studies regarding its genetic diversity and structure in its entire geographic distribution range, and there is a lack of information about the mobility capacity of adults and juveniles.Our results suggest that after the two disturbances, local populations inhabiting Coliumo Bay and the Itata Shelf decreased in genetic diversity and increased in spatial-temporal genetic structure caused by drift and gene flow.

Conclusions
Several studies have examined how marine species can be affected by hypoxia and earthquakes separately, but not by the combined effect of both disturbances.Here, we focused on both types of disturbances, and provided evidence on how the Aphos porosus population responded sequentially after the two major natural disturbances considered in this study, deferring along short spatial distances with different intensities of impact.These results gain importance under a global perspective, given that natural disturbance regimes such as hypoxic-upwelling events are expected to increase due to climate change 65,66 .Results suggest that monitoring for the conservation and maintenance of genetic diversity should be an essential management priority for ensuring the future resilience and adaptive potential of marine populations worldwide 67,68 .

Figure 1 .
Figure 1.Map of the study area; (A).SW Pacific coast, central-southern Chile.2008 hypoxic-upwelling event is represented with a gradient color, where high dissolved oxygen surface concentration (mL/L) is in blue, and low is in red (data obtained from Hernández-Miranda et al., 2010).The 2010 mega-earthquake uplift (cm) is represented with proportional sized black circled, and the tsunami height (m) is represented with blue proportional sized circles (data obtained from Vargas et al., 2011).(B).Sample stations on the Itata Shelf, and C. Sample stations in the Coliumo Bay, Bío Bío Region (36°S, 72°W).

Figure 2 .
Figure 2. PCoA with genetic data from Aphos porosus, by year and area; Coliumo Bay (in red) and Itata Shelf (in green).The axis 1 and 2 explained 23.8% and 20.2% of the data variation, respectively.

Figure 4 .
Figure 4. F ST population pair-wise comparisons for areas and years.Distance method = No. of different alleles.All pair-wise F ST comparisons were significant after Bonferroni correction (p < 0.004).

Table 2 .
Endogamy coefficients (Fis), for each locus, year and area.Hardy-Weinberg deviations are indicated with an asterisk.