Genetic structure, phylogeography, and demography of Anadara tuberculosa (Bivalvia) from East Pacific as revealed by mtDNA: Implications to conservation

Abstract Wild populations of the pustulose ark, Anadara tuberculosa (Bivalvia), an emblematic species of the East Pacific mangrove ecosystem declined in South American countries (Colombia, Ecuador, and Peru) mainly due to overharvesting and habitat loss or degradation. Understanding the genetic aspects of geographic variations and population structure of A. tuberculosa, currently unknown, appears as a priority to fishery authorities in order to elaborate integrated and collaborative conservation policies for fishery management, aquaculture, and stock enhancement programs. We used mtDNA sequence data to investigate haplotype diversity, genetic structure, and demography of A. tuberculosa. Results indicate genetic homogeneity of populations distributed north and south of the equator, respectively. However, statistically significant differentiation emerged between northern and southern populations with pairwise ф ST values ranging between 0.036 and 0.092. The oceanic current system acting in the area (Panama Current and Humboldt Current) might play a role in limiting the larval dispersal of the species, still poorly understood. Demography reconstruction supported recent population expansion, possibly started after last glacial maximum. Our results would suggest separate and independent management of populations north and south of the equator.

ecosystem that has been ancestrally collected by coastal populations as a staple food and that is still consumed for traditional meals in several tropical Latin American countries. Extraction is an essential activity for numerous families whose economy relies on A. tuberculosa trade. Most of the natural stocks of A. tuberculosa are over exploited, and some populations are close to collapse in several countries (Lucero et al., 2012;Mora, Moreno, & Jurado, 2011). In Peru, the A. tuberculosa population of Tumbes Region has been reduced by 6.4-fold between 1988 and 2008 in unprotected mangrove areas as well as in the protected National Mangrove Sanctuary of Tumbes (SNLMT) (Ordinola, Montero, Alemán, & Llanos, 2010;Vivar, 1996). The diminution of A. tuberculosa population, like the majority of bivalve natural stocks worldwide, results from several factors such as overexploitation, habitat degradation, and nonidentified mortalities.
To recover depleted aquatic animal stocks, strategies generally rely on quota implementation with, on one hand, permanent collec- Strategy for restocking and enhancing stocks of natural aquatic populations has also been focused on the mass release of hatchery-produced animals (Arnold, 2008;Bell, Rothlisberg, & Munro, 2005). Such a strategy is attractive for aquatic species, in particular mollusks, considering their extremely high fecundity and the subsequent possibility to produce in hatchery millions of spat from a numerically limited wild broodstock. In fact, aquaculture has started in El Salvador and Costarica (FAO, 2017). However, such a strategy could also lead to a quick reduction of population genetic diversity, generating negative effects. Within this context, understanding the genetic aspects of geographic variations and population structure of A. tuberculosa appears as a priority to fishery authorities in order to elaborate integrated and collaborative conservation policies for fishery management, aquaculture, and stock enhancement programs in a concerted way.
At present, genetic studies of A. tuberculosa still lack and no information is available about the genetic structure of the species at a large geographic scale. In particular, gaining insight on the possible role of marine equatorial currents in shaping the genetic structure of the South American populations could prove very useful for the identification of sanctuaries aimed at the preservation of the genetic variation of the species. In this study, we used mtDNA sequence data (partial COI gene) to investigate haplotype diversity, genetic structure, and demography of A. tuberculosa from the southern edge of its geographic distribution. Samples were collected from two sampling sites north of the equator (Colombia and Ecuador) and three south of the equator (Ecuador and Peru).

| DNA extraction and amplification
Total genomic DNA was individually extracted from approximately 100 mg of tissue that could be spats, gill, and mantle, using a standard CTAB protocol (Folmer, Black, Hoeh, & Lutz., & Vrijenhoek, 1994).
Primers used were as proposed by Folmer et al. (1994), to obtain a partial sequence (620 bp) of the COI mitochondrial DNA gene.
Haplotype and nucleotide diversities were estimated by using DnaSP ver. 6.10.01 (Librado & Rozas, 2009). Because the number of haplotypes and private haplotypes depended of sample size (Supporting Information Figure S1), we also recalculated all statistics after randomly extracting 24 individuals from each population, to match the minimum sample size. Among-and within-group F ST , F SC , ф ST , ф SC , and pairwise F ST and ф ST estimates were calculated by using Arlequin ver 3.5.2.2 (Excoffier & Lischer, 2010). The same software was used to perform Mantel tests (Mantel, 1967). The tests were run to investigate possible isolation by distance (IBD) by estimating correlation between the matrix of pairwise geographic distances between sampling locations and two correspondent matrices of estimates of genetic differentiation: F ST , based on haplotype frequency, and ф ST , that takes into account genetic distances among haplotypes (Bird, Karl, Smouse, & Toonen, 2011). Geographic distances were calculated as the shortest pathways along the costal line.
Demographic histories were studied using DnaSP ver. 6.10.01 to estimate three different classes of statistics under the assumption of neutrality. We estimated class I D*, F* (Fu & Li, 1993), and D (Tajima, 1989) test statistics, which use information of the mutation frequency (segregating sites). We also estimated Fs (Fu, 1997), which uses information from the haplotype distribution (class II). Finally, we calculated Harpending's raggedness (r) index (Harpending, 1994), based on the distribution of the observed pairwise nucleotide site differences (mismatch distribution, class III), and the expected values in populations with constant population size and in growing populations (Rogers & Harpending, 1992). Such statistic is expected to show lower values in mismatch distributions of expanded populations but has little power to detect population expansions. For all statistics, we used the coalescent algorithm implemented in DnaSP to estimate the probability of obtaining values that, under the tested demographic model, would be lower than the observed (one-tailed test).
Following a different approach, different demographic models (constant population size, exponential growth, and Bayesian skyline -BSP) were also investigated using Beast 2.4.5 (Bouckaert et al., 2014). Because the BSP makes few a priori assumptions about the historical demographic trend of the population, it can guide to formulate more specific demographic hypotheses (Crandall, Sbrocco, DeBoer, Barber, & Carpenter, 2012). We used the software jModeltest 2.1.7 (Darriba, Taboada, Doallo, & Posada, 2012) and selected the TIM3+I+G substitution model (AC = CG, AT = GT; unequal base frequencies, invariant sites and gamma distribution), based on the Akaike Information Criterion (Akaike, 1973). For both constant and exponential growth models, a lognormal prior was selected as coalescent population size parameter and a strict clock We therefore considered the best grouping methods those in which among-group heterogeneity was maximized (largest and significant F ST or ф ST estimates) and within-group heterogeneity minimized (smallest and least significant F SC or ф SC estimates). In this case, the analysis was restricted to adjacent groupings, under the implicit assumption of possible IBD. Further support was also sought by comparing pairwise F ST and ф ST estimates.
In the absence of prior information on time of divergence, the to a long-term rate of 0.65%/Myr. Whereas the second rate is certainly too low to effectively trace recent events, the first one might not sufficiently account for purifying selection that keeps removing variation even after lethal mutations have been eliminated. However, if divergence is recent, purifying selection is not expected to have played a major role yet; therefore, a high mutation rate is conceivable.

| Genetic variation and differentiation
Overall, 69 different haplotypes were found. The number of haplotypes decreased to 47 after reducing sample size throughout sampling sites (N = 120, see Table 1). Haplotype diversity ( and linearly correlated to sample size (Supporting Information Figure   S1). Haplotype and nucleotide diversities were not strongly affected by the reduction of sample size, which in turn strongly influenced both the number of haplotypes and private haplotypes ( Table 1).
Values of F ST, ф ST, F SC, and ф SC associated to different assemblages of populations in 2, 3, and 4 groups are reported in

| Haplotype network
Genealogical relationships between haplotypes are reported in

Tumaco (Colombia) Esmeralda (Ecuador) Guayas (Ecuador) El Oro (Ecuador) Tumbes (Peru)
Tumaco ( Pie charts indicate the proportion at which each haplotype occurs at each location substitution. Most of them were found at sites south to the equator. Longer branches departed from Haplotypes 1 and 8, some of them leading to groups of haplotypes geographically more localized at north or south to the equator. The unequal sample size biased the graph in favor of Tumbes (Peru), but did not alter the general topology of the network (Supporting Information Figure   S2).

| Demographic histories
Based on ф ST analysis, demographic histories were investigated by considering two groups of populations: Group 1 that included Tumaco and Esmeralda, and Group 2 that comprised Guayas, El Oro, and Tumbes. In all tests, D*, F*, Fs, and D were negative, but  Table 4.
The BSP (Figure 4) indicated a modest increase in N e μ for Group

(a). A clearer evidence of demographic expansion was shown by
Group 2 (b). Bayes factors indicated the BSP as the best representation of the demographic history of Group 1, whereas the exponential model was preferred for Group 2 (Table 5).

| Geographic variation
Our study showed high level of genetic variation in equatorial A. tuberculosa, as estimated using maternally inherited COI sequences. Haplotype diversity (h) and nucleotide diversity (π) were comparable to or higher than values observed in other natural populations of bivalves (Xue, Wang, Zhang, & Liu, 2014). Despite the potentially high larval dispersal, the species appears significantly structured and two groups of populations seem to exist, which comprise populations at north and south of the equator,  Additionally, all statistical class I tests used here increase the power to reject the constant size model with increasing the degree of expansion, and consequently, large samples are needed to detect small population growth events. This is consistent with the sample size for Group 1, smaller than for Group 2, and with results of the Bayesian analyses that indicated population growth for both Groups 1 and 2, but with different characteristics. In fact, Bayes factors supported the skyline model as best fit for Group 1, whereas exponential growth was preferred for Group 2. Interestingly, the slope of the skyline for Group 1 (Figure 4)

| Dating demographic histories
Despite the degree of the demographic events was different in Groups 1 and 2, they seem to have initiated approximately at the same time, around 5E-4 mutations/site. In this case, the application of the two rates provided by Crandall and collaborators would pinpoint the start of the demographic growth between approximately 77Kya (lower rate) and 9.4Kya (higher rate), both estimates suggesting that the event was recent. If the appropriate mutation rate is comprised between 5.3%/Myr and 0.65%/Myr, it is very likely that the demographic process may have started at the end of the Last Glacial Maximum (LGM), which is dated approximately between 19 and 20 Kya in the northern hemisphere and 14-15 Kya in the West Antarctic (Clark et al., 2009) and Central Pacific (Blard, Lave, Pik, Wagnon, & Bourles, 2007). In the tropics, glaciers could have reached their greatest extent 34Kya and were retreating approximately 21 Kya (Smith, Seltzer, Farber, Rodbell, & Finkel, 2005 (Bromley et al., 2009;Shakun et al., 2015).
If so, the proper μ is most likely more proximate to the higher value suggested by Crandall and collaborators and would corre-

| CON CLUS IONS
In conclusion, we showed a high level of genetic variation in equatorial A. tuberculosa. Despite the potential long-range larval dispersal, such variation is geographically structured with populations at north and south of the equator not being completely connected.
As in other bivalves, a combination of life history traits and oceanic prevalent currents may explain the pattern observed. Climatic changes occurred at the end of the LGM may be responsible for the documented demographic expansions. While awaiting for more data becomes available from other locations, we suggest separate management of populations at north and south of the equator.

ACK N OWLED G M ENTS
Results presented in this research were partially funded by the NGO Vicario for his valuable criticism on an early version of this paper. The authors take the opportunity to express their respect to each of the concheros and concheras of Colombia, Ecuador, and Peru, in recognition of the hard work they do to support their families, as well as efforts to continue improving the activity, protecting the resource and the mangrove ecosystem for future generations. We thank Rita Castilho and an anonymous reviewer for their constructive criticism.

CO N FLI C T O F I NTE R E S T
The authors declare no competing interest.

AUTH O R CO NTR I B UTI O N
BD and VC designed the work; BD, KP, RA, and CC collected samples and contributed genetic data; GG and BD contributed data analysis; GG, BD, and VC drafted the manuscript. All authors critically revised and approved the final manuscript.

DATA ACCE SS I B I LIT Y
The raw data underlying the main results of the study are archived in Genbank (accession numbers: MK043084 -MK043325).  Note. Bayes factor for Group 1 (below the diagonal) and Group 2 (above the diagonal). BSP is selected for Group 1 and exponential growth for Group 2. Group 1 includes Tumaco and Esmeralda, whereas Group 2 comprises Guayas, El Oro, and Tumbes