Post-Disturbance Genetic Changes: The Impact of the 2010 Mega-Earthquake and Tsunami on Chilean Sandy Beach Fauna

Earthquake/tsunamis can have profound impacts on species and their genetic patterns. It is expected that the magnitude of this impact might depend on the species and the time since the disturbance occurs, nevertheless these assumptions remain mostly unexplored. Here we studied the genetic responses of the crustacean species Emerita analoga, Excirolana hirsuticauda, and Orchestoidea tuberculata to the 27F mega-earthquake/tsunami that occurred in Chile in February 2010. mtDNA sequence analyses revealed a lower haplotype diversity for E. analoga and E. hirsuticauda in impacted areas one month after the 27F, and the opposite for O. tuberculata. Three years after the 27F we observed a recovery in the genetic diversity of E. analoga and E. hirsuticauda and decrease in the genetic diversity in O. tuberculata in 2/3 of sampled areas. Emerita analoga displayed decrease of genetic differentiation and increase in gene flow explained by long-range population expansion. The other two species revealed slight increase in the number of genetic groups, little change in gene flow and no signal of population expansion associated to adult survival, rapid colonization, and capacity to burrow in the sand. Our results reveal that species response to a same disturbance event could be extremely diverse and depending on life-history traits and the magnitude of the effect.


Results
The sequence alignment length and concatenation of both genes was 875 bp for T1 and T2, respectively, in E. analoga, 1070 bp (T1 and T2) in E. hirsuticauda, and 1078 bp (T1 and T2) in O. tuberculata (Genbank Accession number MK914652 -MK917441). No evidence of stop codons were detected in the COI gene, no substitution saturation was found for any species or sampling time, indicating that data sets were suitable for all of the subsequent genetic analyses.
Genetic structure. The Bayesian approach of spatial clustering based on Geneland analyses showed changes in the genetic structure of all of the studied species through time . The number of genetic clusters decreased throughout time in E. analoga (from K = 5 to K = 2) and slightly increased in E. hirsuticauda (from K = 4 to K = 6) and O. tuberculata (from K = 7 to K = 8; Figs 3a,b-5a,b). Probabilities for each genetic cluster were greater than 0.5 in most cases, except for E. hirsuticauda in T2, where the probabilities ranged between 0.2-0. 4. Emerita analoga populations in the impacted area were well differentiated in T1, while they were part of the largest cluster in T2 (Fig. 3a,b). Excirolana hirsuticauda and O. tuberculata in the impacted area maintained the same genetic differentiation in southern and northern non-impacted areas between T1 and T2 (Figs 4a,b and 5a,b). The AMOVA comparing the zones with different levels of impact showed that E. Genetic variation. Mean haplotype diversity (Hd) in E. analoga populations from the impacted area was lower than Hd from northern and southern non-impacted areas at T1 (Table S2). At T2, Hd from impacted areas became very similar to both non-impacted areas at T2. Hd in E. analoga was significantly different between sites (p < 0.001), time (p < 0.002), and site*time (p < 0.002; Fig. 3c). For E. hirsuticauda, no significant differences Figure 5. Genetic clusters identified with GENELAND (each color represents a genetic cluster) at Time 1 (T1; a) and Time 2 (T2; b), mean haplotype diversity (c) and nucleotide diversity + SD (d) including genetic differences and results of AMOVA analyses among sites, times (T1 and T2), and sites-times, mismatch distribution (e), and gene flow estimated with MIGRATE-n for O. tuberculata at T1 (f) and T2 (g). Arrows' width represents the amount of immigrant numbers, circles' size represents the value of Ne. NNI: northern nonimpacted area, IA: impacted area and SNI: southern non-impacted area. (2019) 9:14239 | https://doi.org/10.1038/s41598-019-50525-1 www.nature.com/scientificreports www.nature.com/scientificreports/ were found in the mean value of Hd for either of the two factors (site and sampling time) nor for their interaction site*time (p < 0.330, p < 0.230, p < 0.480; Fig. 4c; Table S2). Mean Hd in O. tuberculata in the impacted area was higher than Hd in northern and southern non-impacted areas at T1 (p < 0.02; Fig. 5c; Table S2). At T2, the mean Hd in the impacted area was still higher than in both northern and southern non-impacted areas, but northern non-impacted areas showed less Hd than southern non-impacted areas. The lowest mean value of nucleotide diversity (π) for E. analoga was observed in the impacted area at T1, and increased significantly between T1 and T2 (Fig. 3d). At T2, π of this species in impacted and northern non-impacted areas was similar, while it was higher in southern non-impacted areas. Nucleotide diversity (π) in E. analoga showed significant differences for both the factor area and sampling time (p < 0.001 and p = 0.004 respectively), but not for their interaction (Fig. 3d). For E. hirsuticauda, π did not differ for either of the two factors or for the interaction between them (Fig. 4d). Nucleotide diversity of O. tuberculata was lower than in the other two species (Fig. 4d). At T1, π was higher among populations in impacted than in non-impacted areas. At T2, a decrease in π in populations in impacted were observed while values of π were maintained in northern non-impacted areas. Contrastingly, an increase in π among southern non-impacted areas between T1 and T2 was observed. In O. tuberculata, differences in π values were non-significant among sites, times and sites*times (Fig. 5d). Differences between pairs of sequences estimated using the mismatch analyses showed a unimodal distribution for E. analoga at T2, a signal of population expansion (Fig. 3e). Contrastingly, a bimodal distribution for E. hirsuticauda and O. tuberculata in T1 as well as in T2 was observed (Fig. 4e,e).

Gene flow analyses.
In all cases, the migration analysis supported the panmictic model (BF (Panmictic-Full Model) >5.0) and showed temporal changes in effective population size (Ɵ) and migration rate among clusters (M). In E. analoga the cluster composed of populations in the impacted area had the lowest value of effective population size in T1 (Ɵ: 0.43E-3 ± 0.37E-3; Fig. 3f). This magnitude was approximately 8 and 6 times lower than the maximum values of Ɵ at T1, which corresponded to the Cocholgüe (Ɵ: 3.42E-3 ± 3.39E-3) and Boca (Ɵ: 3.42E-3 ± 3.39E-3) clusters, both located in the northern non-impacted area. At T2, the effective population size of the cluster formed by the population of Colcura (Ɵ: 9.49E-3 ± 3.28E-3) was approximately 9.1 times greater than the Ɵ value of the cluster combining the remaining populations (Ɵ: 1.04E-3 ± 1.09E-3; Fig. 3g). Gene flow was moderate to high among all populations at T1, with the highest number of immigrants occurring in non-impacted areas in the southernmost localities. In T2, gene flow was high and mostly symmetric for the two identified clusters. Despite the difference in Ɵ, the migration rate from the largest cluster to Colcura (M: 3901 ± 1539 migrants*generation −1 ) was 12% greater than the gene flow from Colcura to the largest cluster (M: 2731 ± 1421 migrants*generation −1 ; Fig. 4g).
In E. hirsuticauda, the two clusters that contained the populations of the impacted area, the Lebu-Rumena and Lobería-Nigue-Morhuilla clusters, presented the lowest values of effective population size in T1 (Ɵ: 2.94E-3 ± 3.27E-3 and Ɵ: 2.54E-3 ± 4.08E-3 respectively; Fig. 4f). These Ɵvalues were about 75% lower than the effective population size of the cluster formed by the extreme populations from both the northern and southern non-impacted areas (Cocholgüe-Boca-Cheuque-Pichicuyín cluster), and approximately 24% lower than Ɵ of the Chivilingo-Colcura-Lota cluster. At T2, the two northernmost locations (Cocholgüe and Boca) constituted the cluster with the largest population size (Ɵ = 4.92E-3 ± 5.35E-3), followed by three clusters combining localities from the affected area and the southern non-affected area. Of these three clusters, the individuals of the affected populations in Lebu and Rumena were genetically grouped with Nigue, forming the cluster with the second largest population size Ɵ: 2.96E-3 ± 2.50E. At T1, gene flow among all populations was moderate to low and mostly symmetric ( Fig. 4f), while it was moderate and symmetric at T2 (Fig. 4g).

Discussion
Disturbance processes, such as earthquakes associated with coastal uplift and tsunamis, can shape the genetic patterns of coastal species (i.e. population structure, gene flow and diversity) 43 . When populations are suddenly impacted by a mega disturbance, their genetic diversity tends to decline and the rate of that decline depends mostly on the effective population size, the number of mutations occurring after the disruptive event, and the amount of new genetic variants in a population as a result of gene flow [44][45][46] . Nevertheless, not only the genetic attributes of populations can help to understand the genetic responses of species to a disruptive event; ecological traits also play an important role in explaining the direction and magnitude of changes observed in the genetic composition of species under perturbation. As a result of regional and local variations in habitat, magnitude of the disruptive event and ecological attributes characterizing affected communities, each population is expected to react uniquely to perturbations, revealing a unique history of bottleneck and recovery 18 . Our results, comparing the genetic consequences of the 27F earthquake/tsunami in three species with different history traits suggest that species' responses to disruptive events can be extremely diverse, even when exposed to the same disturbance. These responses are highly depending on the combined effects of the main species traits and modalities such as movement type, mobility degree, and intertidal zone habitat occupied (Fig. 1). www.nature.com/scientificreports www.nature.com/scientificreports/ The negative impacts of the combined effects of earthquake and tsunami on impacted areas through time is evident in the genetic parameters of E. analoga. The genetic diversity of E. analoga was lower among populations from impacted areas at T1, showing a significant increase by T2, reaching values similar to those found in non-impacted areas. After the earthquake, E. analoga showed a significant reduction in the number of populations; nonetheless, throughout time, as a consequence of gene flow and population expansion, E. analoga showed a higher genetic diversity at T2 when compared to T1. At T2 the largest cluster, including individuals from 13 of the 14 sampled localities, presented the lowest effective population size (N e ) compared to the populations comprised of individuals from the single locality of Colcura. This lower N e in the largest cluster could be explained by the observed population expansion signal and increase in the gene flow, leading to the homogenization of the 13 populations at T2, a pattern that has been well documented in theoretical and empirical studies 47 . Ecological studies analysing the effects of the 27F earthquake in E. analoga have suggested that this species was one of the most affected by the coastal uplift and tsunami during the 27 earthquake 4,41 . Consequently, the reduction of the genetic diversity of E. analoga through time could be the result of the combined effect of the removal of recruits and adults during the 27F tsunami as well as a source-sink dynamic. In addition, it has been noticed that the higher recruitment of E. analoga throughout the sampled area occurs near February 41 , thus the disturbance events from the 27F might have reduced the available habitat for new recruits. Provided that the massive adult mortality was a consequence of the coastal uplift, which may have wiped out most of the rare genetic variants in adults, a lower genetic diversity and higher population differentiation could be expected after the earthquake. However, our results suggest that three years after the disruptive event of the 27F, populations might have recovered, with new recruits settled, reducing the genetic differentiation, increasing the genetic diversity and revealing signals of population expansion in the mismatch analyses. This last finding could be the result of a recolonization of individuals from distant non-impacted areas, which can be explained by the high dispersal potential of E. analoga given its intermediate larval dispersal stage in areas like Polynesia and Antarctic waters 30,31 .
The genetic patterns observed in E. hirsuticauda are in line with the ones of other intertidal species 23 . This species showed no significant differences in either genetic diversity index at the two sampling times, in impacted and non-impacted areas. However, the analysis of population differentiation throughout time revealed a slight increase in population differentiation, no signal of population expansion, and an increase in gene flow along with a reduction of N e in the southern non-impacted area at T2. A low or nil impact of mega disturbances has also been reported in other marine species. Though this is the first study assessing the effects of earthquakes/tsunamis on crustacean's species, studies in other species can also be a matter of interest for contrasting results. For instance, along the coast of Japan, authors have noticed that the intertidal mud snail Batillaria attramentaria experienced a reduction in its effective population size (N e ) after the 2011 Tohoko tsunami 23 . Nevertheless, the effect of this mass mortality produced no significant reduction in the allelic richness of this mud snail. These authors suggested that the lack of impact on this genetic diversity descriptor may have been the result of a combination of very large populations before the tsunami and a high recruitment of juveniles after the disturbance.
Although the general patterns revealed in our study found that at T2 E. analoga and E. hirsuticauda species decreased their Hd in impacted areas and increased their Hd in non-impacted areas, in O. tuberculata the highest Hd was recorded in the impacted area at T1. The same pattern was observed at T1 for π in O. tuberculata, but the increase in π in northern non-impacted areas was so high that it turned out to be the area with the highest Hd at T2. Similar to that observed in E. hirsuticauda, an increase in the number of populations was detected in O. tuberculata at T2, with the differentiated populations located within the southern non-impacted area. Although there are only a few of studies available linking ecology and genetics in sandy beach organisms exposed to disturbances, our results are on the line of previous results with talitrids [48][49][50] . Ecological studies have found a similar trend, with the species E. hirsuticauda and O. tuberculata showing similar abundances before and after the 27F earthquake 4 . Other authors have found that both species quickly recovered their abundance after the tsunami 51 , suggesting that the particular life history traits of these species may help to buffer the impacts of the mega-perturbation on their local populations. For instance, a high recovery of Orchestoidea tuberculata has been suggested to be related with its feeding strategy as scavengers and detritivorous (Fig. 1) 51 . For E. hirsuticauda, the lack of significant reductions in its density and genetic diversity could be explained by the high mobility of adults, and high dispersal potential of planktonic larvae. The swimming capacity of adults may allow them to resist the tsunami wave and then, relocate to a new site. In addition, the free larval stage has the capacity to rapidly colonize impacted areas as was suggested previously 42 . Contrary to E. hirsuticauda, O. tuberculata does not present an intermediate dispersal stage given that females incubate their offspring until they reach the juvenile stage (Fig. 1). However, adults have the capacity to jump and burrow up to 30 cm in the sand for shelter (Fig. 1) 34 , which together with their high resistance to desiccation could help this species to survive after disrupting events. Thus, it appear that O. tuberculata shows a high genetic and ecological resilience, possibly providing this species with advantageous adaptations to evade the direct effects of the tsunami and survive a coastal uplift and tsunami such as the 27F mega-earthquake.
Our study clearly reveal that species' response to a major disturbance, in this case an earthquake and tsunami, depends on the ecology of each species; and that these distinct responses can influence the genetic diversity of species. The use of highly variable DNA sequences allowed to observe these changes in genetic diversity over a short period of time, only a few generations within each species. Whether DNA sequences can be used for recent processes has been the debate of multiple studies throughout the last decade (e.g. [52][53][54], with no decisive answer. Here we have shown that recent massive disruptive events can actually be detected with DNA sequences, most likely due to the magnitude of the event. Nevertheless, it will be necessary to compare these results with other markers characterized by faster mutation rates, such as microsatellite or SNPs 53 , in order to confirm our conclusion. In a broader perspective, no hypotheses has been made yet on the ecologically meaningful units at which genetic patterns are revealed for sandy beach resident organisms. Our results suggest that we can detect changes in the genetic patterns of sandy beach organisms at population level after a disruptive event, even between relatively close populations. Although more observations are needed, this study can help future researchers to evaluate what major drivers explain the genetic diversity of species (with different life history traits) inhabiting areas frequently exposed to disruptive events. www.nature.com/scientificreports www.nature.com/scientificreports/ Although areas suffering constant disturbances present important challenges for human populations, they also present interesting opportunities, since they act as natural laboratories to study changes in the ecological and genetic patterns of populations during recovery process after strong perturbations. It could be crucial to submit these areas to continuous sampling and monitoring in order to understand how species respond to changes in a landscape under recurrent perturbation regime. This information could then be used to predict what may happen to species facing natural or human-induced habitat modifications and thus provide better information for the conservation of species undergoing environmental disturbances.

Methodology
Individuals' collection and laboratory protocols. Given that no samples were available before the mega-earthquake that occurred in February 2010, we used an impact-control approach. For this, we sampled 14 localities (sandy beach localities between 36°S and 39°S) with different levels of impacts, according to the uplift level and tsunami strength 4,14 . The magnitude of the impact was categorized in 3 areas: northern non-impacted area (NNI), impacted area (IA) and southern non-impacted area (SNI) (Fig. 2; Table S2). Samplings were carried out in March 2010 (T1), one month after the mega-earthquake, and repeated three years later in 2013 (T2), between February and June. All sampled localities were inside non-protected areas. Whole DNA was extracted from collected individuals using an already published protocol 55 and its quality was evaluated in a ND 1000 spectrophotometer (NanoDrop Technologies). DNA templates were amplified for the mitochondrial genes COI and 16S using the PCR protocols described by other authors [56][57][58] . Amplified products were then sequenced in Macrogen Inc. (Seoul, Korea). Sequences of both genes were concatenated with Mesquite 3.4 software 59 and translated to determine whether stop codons were present in the COI gene.
Characterization of the sequence matrices. The potential existence of saturation 60 and adjustment to a neutral model of mutation-drift equilibrium 61 were assessed with the software DAMBE 6.4.100 62 and DnaSP 6.0 63 , respectively. DnaSP was also used to calculate both haplotype and nucleotide diversity (Hd and π). For each species, mean Hd and π values were compared using a two-way Analysis of Variance (ANOVA), with two fixed factors: area with three levels (northern non-impacted area, impacted area and southern non-impacted area), and time with two levels (T1 and T2). When significant effects were detected, Tukey a posteriori tests were performed. All two-way ANOVAs were run in the STATISTICA 7.0 program 64 .
Analysis of spatial genetic structure. GENELAND 4.08 was used to detect genetic discontinuities under a spatially explicit genotypic context 65 for each species and sampling time (T1 and T2). Short analyses were carried out with K = 12 for all species, with 5,000,000 MCMC iterations, thinning each 1,000 iterations and burn-in of 10%, using the Akaike Information Criteria in MCMC context (AICM) 66 . Final simulations were carried out by fixing the number of most probable clusters (K) obtained in the first simulations for each species and sampling time, using 30,000,000 iterations of MCMC, thinning each 1,000 iterations and burn-in of 10%. In each analysis, we computed the Poisson-Voronoi tessellation maximizing the value for the rate of the Poisson process (λ) 67 . The number of genetic clusters was assessed through probability density distribution graphs, for each species in T1 and T2.
Genetic comparison and population expansion. From the variations in the observed haplotype frequency, the genetic differentiation for each species and sampling time was compared applying an Analysis of Molecular Variance (AMOVA). The design of these analyses was carried out by grouping localities depending on the level of impact of the tsunami/earthquake (i.e., impacted area vs. northern and southern non-impacted areas). The genetic differentiation was assessed through AMOVAs based on matrices of estimated genetic distances according to the better-fitted sequence evolution model 68 using the software Arlequin 3.5.2.2 69 . Evolution models of the sequences (estimated for each gene independently and for the concatenated data set), and the underlying parameters, were estimated with the program jModelTetst 2.1.10 70,71 . To investigate signals of population expansion in the three species at each sampling time (T1 and T2) we computed the distribution of pairwise differences from the segregating sites of the concatenated data set by calculating the mismatch distribution analyses implemented in DnaSP.
Gene flow. It has been recently noted that Migrate-n can provide accurate results of recent migration rates, despite the fact that it has most frequently been used to explain long-term historical patterns 72 . Therefore, we used this software to calculate the genetic diversity and migration patterns of populations with the parameters effective population size (Θ) and migration rate among clusters (M) scaled by the mutation rate (μ) in Migrate-n 3.6.11 [73][74][75][76] . The analyses were performed using 10 replicates of short chains with 50,000 iterations and one long chain with 1,000,000 iterations and 10% burn-in. The obtained values were used to parameterize the priors of Bayesian runs using six independent replicates for each species and sampling time. Each Bayesian run was performed with the transition/transversion ratio value of each gene obtained from the previously inferred evolution models. In each replicate we computed 5,000,000 iterations with samples taken every 1,000 steps and 10% burning. We compared the parameters obtained based on: (1) better adjustment considering constant or variable mutation rates for both genes, and (2) better adjustment to an exponential or uniform model of a posteriori distribution. For each type of mutation rate and distribution model, the six replicates were merged, and the associated marginal likelihood was compared under the Bayes Factor criterion for the selection of models 77 implemented in Tracer 1.6.0 software 78 with 10,000 Bootstrap replicates. The mean and standard deviation of the obtained Θ and M distributions were represented in gene flow diagrams using the Excel add-ins PopTools 3.2 79 . Additionally, the inferred population genetic structure was contrasted with a panmictic model, where it was assumed that all of the clusters form part of the same genetically interconnected unit. The adjustment was determined using Bayes Factor criterion estimated from the marginal likelihoods of both models using Tracer 1.6.0 software with 10,000 Bootstrap replicates.