Population genomics reveals a mismatch between management and biological units in green abalone (Haliotis fulgens)

Sample collection and study area

Collections were performed from July to August 2017 and May 2018 from 10 locations along the BCP, which belonged to three of four fishing administrative zones (Fig. 1). Abalones above the minimum legal size (n = 185) were collected under SAGARPA permits PRMN/DGOPA-011/2017 and PRMN/DGOPA/016-2018 by hookah diving during fishery evaluation surveys performed by the cooperative fishers. Epipodial tentacle fragments were taken without sacrificing the animals and preserved in 95% ethanol until DNA extraction (Table 1).

DNA extraction and library preparation

Genomic DNA was extracted using a salt-extraction protocol (Sambrook, Fritsch & Maniatis, 1989) with RNase treatment. Genomic DNA integrity was verified using 1% agarose gels and quantified with a Qubit fluorometer 2.0 (Life Technologies, Carlsbad, CA, USA). All samples were normalized to a final concentration of 50 ng µl⁻¹.

Double-digest restriction-associated DNA (ddRAD) libraries were prepared following Peterson et al. (2012) with some modifications; 500 ng of genomic DNA per sample was digested using the restriction enzymes EcoRI-HF (NEB^®) and MspI (NEB^®). A Pippin Prep device (Saga Science) was used to select fragment sizes from 376 bp–408 bp for all the libraries. Two ddRAD libraries were completed, each one with two pools; each pool comprised 48 individuals with a different barcode, and a unique index. A 4-index combinatorial and 48 barcodes described by Peterson et al. (2012) were used. To avoid batch effect, individuals from the same location were pooled in different pools and libraries in unequal proportions. Genotyping accuracy was validated using three replicates with different index-barcodes; then, the concordance genotype percentage between replicates of the same subject genotyped was assessed. The final libraries were sequenced on two lanes of an Illumina Hi-seq 4,000 platform at Novogene (Sacramento, CA, USA), using 150-bp pair-end reads.

Raw data filtering and SNP calling

De-multiplexing and quality filtering were undertaken using the process_radtags pipeline in the STACKS software v1.30 (Catchen et al., 2011). Raw reads were trimmed to 140 base pairs to reduce sequencing errors present at the tail of the sequences. Read-pairs with the expected restriction sites and full barcodes at both ends were identified, allowing up to one error in each barcode. To homogenize the number of reads of each individual sample, they were standardized to 1 million reads each.

Reads were assembled de novo using the denovo_map.pl pipeline in STACKS. After several combinations of the main parameter optimization, the parameter combination was evaluated assessing the parameter effect and yield loci number as recommended in the literature (Paris, Stevens & Catchen, 2017; Rochette & Catchen, 2017) (Supplemental S1). The final parameter combination to build loci was a minimum stack depth of five (−m = 5), and a maximum of three nucleotide mismatches among stacks within individuals (−M = 3). To build the catalog, a maximum of four mismatches between sample loci between individuals (−n = 4) were allowed.

Further data filtering was performed using the population pipeline in STACKS, retaining those loci present in at least 80% of the individuals at each sampling site (−r = 0.8) in all of the localities 100% (−p = 10) and with a minor allele frequency MAF > 0.05. To prevent physically linked loci for subsequent analyses, the write_single_SNP option was used to retain only the first SNP from each RAD tag. Individual samples with >10% missing data across loci estimated with vcftools v.3.0 (Danecek et al., 2011), were excluded from subsequent analyses. The filtered vcf file was converted into the necessary file formats for downstream analyses using PGDSPIDER v.2.1.1.5 (Lischer & Excoffier, 2012). Linkage disequilibrium between all marker pairs was estimated using the function pair.ia of the package poppr v2.8.1 (Kamvar, Tabima & Grünwald, 2014) in R v.3.4.3 (R core Team, 2014) with 199 permutations; p-values were converted to false discovery rate (FDR) q-values (Benjamini & Hochberg, 1995) using the p.adjust function of the package stats v3.5.1 in R v.3.4.3 (R core Team, 2014) with a q-value of 0.05. Linked loci were discarded for further analysis.

Outlier detection and Hardy Weinberg Equilibrium

To identify and eliminate potential candidate loci under selection, two independent methods were used. The first one employed a Bayesian approach to estimate the population specific F_ST coefficients using the BayeScan v3 software (Foll & Gaggiotti, 2008), which integrated population-specific (neutral genetic divergence) and locus-specific (selection) effects (Foll, 2012). The decision criterion to determine whether a locus was likely to be under a strong selection was the q-value (Foll, 2012) analogous to a FDR p-value that must be under 0.05. The BayeScan software was run with its default parameters set at a minimum number of 10,000 iterations, the length of 20 pilot runs to 5,000 iterations, and burn-in length to 50,000 iterations. For the second method, coalescent simulations were used to obtain a null distribution and confidence intervals around the observed values to see if the observed locus specific F_ST values could be considered as an outlier F_ST, conditioned on the globallly observed F_ST value. For this purpose, Arlequin v3.5.2.2 software (Excofier & Lischer, 2010) was used with 1,000 simulated demes and 100,000 coalescent simulations were run to detect outlier loci based on their F_ST and p-value (p < 0.01).

Hardy–Weinberg equilibrium (HWE) was estimated with the software Arlequin v3.5.2.2 (Excofier & Lischer, 2010) by exact tests with 100,000 dememorization steps to remove those loci that were in disequilibrium in five (50%) or more locations (p < 2.3 × 10⁻⁵ after Bonferroni correction).

Summary statistics

Genetic diversity estimates, such as allele frequencies, rarefied allelic richness (Ar), observed heterozygosity (H_o), expected heterozygosity (H_e) and Inbreeding coefficient (F_IS) were estimated with the package PopGenReport v.3.0.4 (Adamack & Gruber, 2014) in R v.3.4.3 (R core Team, 2014). Kruskal–Wallis test and a posteriori Gao test were conducted to examine the differences between genetic diversity values of localities with the packages stats and nparcomp v.3.0 (Konietschke, Noguchi & Rubarth, 2019) in R v.3.4.3 (R core Team, 2014).

Population structure analysis

The F_ST was calculated using the basic.stats function in the R package hierfstat v.0.04-22 (Goudet & Jombart, 2015). The differentiation degree between locality pairs was quantified using θ, an unbiased estimator of F_ST (Weir & Cockerham, 1984) in Arlequin v.3.5.2.2 software (Excofier & Lischer, 2010). The significance level was assessed using 10,000 permutations and corrected for multiple comparisons using the Bonferroni correction.

A discriminant analysis of principal components (DAPC) was performed in the R Package adegent v.2.1.1 (Jombart, Devillard & Balloux, 2010). This method maximize the variance among groups while minimizing the variation within groups and does not make any assumptions about the underlying population genetic model. To assess the optimal group number the function find.clusters was used with the Bayesian information criterion (BIC) method. To estimate the number of principal components (PCs) to be retained, the function optim.a.score was used.

A Bayesian clustering method with the software STRUCTURE v.2.3 (Pitchard, Stephens & Donnelly, 2000) was used to infer population structure and identify the most likely value of K (groups or putative populations), based on the admixture ancestry model and correlated allele frequencies. Runs consisted of an initial burn-in of 10,000 Markov Chain Monte Carlo (MCMC) iterations followed by 100,000 iterations for each inferred cluster (K 1 to 10 with five replicates) and the locprior model. The delta K method was used to determine an optimal K (Evanno, Regnaut & Goudet, 2005).

Analyses of molecular variance (AMOVA) using Arlequin v3.5.2.2 software (Excofier & Lischer, 2010) were performed with 10,000 permutations to evaluate the genetic variation within and between the different groups testing for several hierarchical hypothesized structures: (A) administrative zones as genetic groups; (B) two groups based on STRUCTURE results; (C) three groups based on F_ST results; (D) three groups based on DAPC and the unweighted pair group method with arithmetic mean (UPGMA) results; and (E) four groups based on visual division on the DAPC results.

An UPGMA tree was constructed using the program Phylip (Felsenstein, 1989) based on genetic distance estimates among populations (Cavalli-Sforza and Edwards chord distance, Dc) with 3,000 bootstrap replicates. The consensus tree was visualized using the Treeview software v1.6.5 (Page, 1996).

To further investigate the directional relative migration rates between localities, migration networks were generated using divMigrate function (Sundqvist et al., 2016) of the R package diveRsity (Keenan et al., 2013) with G_st (Nei, 1973) as a measure of genetic distance with 1,000 bootstrap repetitions and an arbitrary filter threshold of 0.60. Gene flow patterns were visualized using network graphics produced with the R package qgraph (Epskamp et al., 2012). Putative source and sink populations were identified by divMigrate as it calculates relative levels of migration by assessing the genetic differentiation between two populations and a hypothetical pool of migrants.

Finally, isolation by distance (IBD) was tested in the R Package adegent v2.1.1 (Jombart, Devillard & Balloux, 2010) using a Mantel test (Mantel, 1967) by correlating Slatkin’s linearized [F_ST/(1 – F_ST)] (Rousset, 1997) vs. geographical distance (shortest waterway distance in kilometers between sampling location pairs) with the function mantel.randtest using a Monte Carlo simulation with 999 permutations. These tests were performed at three hierarchical levels (1) including all locations; (2) using only coastal locations (excluding Guadalupe Island); and (3) within the central-south locations (excluding Guadalupe and San Jeronimo islands and Faro San José). To test whether the correlations between genetic and geographic distances were the result of a continuous or distant patchy cline of genetic differentiation, a 2-dimensional Kernel density estimator was applied with the function kde2d using the MASS package (Venables & Ripley, 2002) in R v.3.6.2 (R core Team, 2014).

Individuals and single nucleotide polymorphism filtering

A total of 1,420,020,148 pair-end raw reads were obtained from two sequencing lanes. After applying the quality filters with the process_radtags program of STACKS (barcode or restriction enzyme site absence, low-quality reads), 96.8% of the reads were retained, resulting in an average of 7,458,559 reads per sample (range: 0.33 M–23.3 M). Individual samples with less than 1 million reads were discarded (n = 8). The remaining samples were standardized to 1 million reads and initial catalog of 116,443 RAD loci with 353,332 SNPs was built with the denovo_map.pl module in STACKS, from which the populations module extracted a dataset of 2,216 loci with an average read depth of 32X. Individuals with >10% of missing genotypes also were excluded (n = 2). After the FDR correction, no evidence was observed for significant linkage disequilibrium between loci (p > 0.05). which resulted in a global database of 175 individuals and 2,216 SNPs.

Outlier detection and Hardy Weinberg equilibrium

Forty-five loci outliers were identified by Arlequin and five by Bayescan (five were common in both methods), leaving a dataset of 2,171 SNPs (see Supplemental S2 for details on the number of SNPs removed after applying quality filters). One locus showed significant deviations from HWE (p < 0.05) in more than half of the locations and it was removed from the database, leaving a neutral data set of 2,170 SNPs. Genotyping accuracy was evaluated using individual replicates, and an average genotyping error rate of 0.04% was recorded.

Summary statistics

The estimates of global neutral genetic diversity were H_o = 0.284; H_e = 0.300; Ar = 1.877. Samples from the northern region (Guadalupe and San Jerónimo islands, and Faro San José) were the least diverse, while Punta Eugenia, Cedros island, and Tortugas South were the sites with the highest diversity (Table 1). Significant differences in genetic diversity values were found among localities, of which Guadalupe and San Jerónimo islands those with more pairwise significant differences (Supplemental S3). The F_IS was 0.004 (p = 0.38) and did not show significant values across locations (F_IS = −0.019 to 0.020; Table 1).

Population structure analysis

The population substructure was significant (F_ST = 0.005, p = 0.006). Pairwise comparisons showed significant differences between the oceanic Guadalupe island (~250 km away from the coast) and the rest of the localities (p < 0.001). Aditionally, the northern coastal San Jerónimo island had significant F_ST estimates with all locations (p < 0.001) with the exception of Faro San José. Sites within the southern region showed no significant genetic divergence values among localities except for Faro San José with Tortugas South (Supplemental S4; p < 0.001).

The DAPC analysis retained 28 principal components and identified three genetic clusters, one composed by the oceanic Guadalupe island, the second by the coastal northern locations (SJi and FSJ), and the third by coastal central-southern localities (CI, PE, TN,TS, AN, AS, and Bo) (Fig. 2). The Bayesian cluster analysis in STRUCTURE revealed two genetic clusters (K = 2) based on the delta K statistic of Evanno, grouping all individuals from Guadalupe Island separately from the rest of the coastal locations (Fig. 3A). However, a visual inspection of STRUCTURE admixture plots for K = 3 (Fig. 3B) showed an additional population subdivision resulting in a grouping pattern similar to the DAPC (Fig. 2) and the F_ST results (Supplemental S4).

Hierarchical AMOVA supported a regional population structure where only the grouping C (GI, SJI, all other localities), D (GI, SJI + FSJ, all other localities), and E (GI, SJI, FSJ, all other localities) resulted in a significant variation among groups (F_CT = 0.0104−0.0148; p < 0.05). However, only the last grouping (E) that consider Guadalupe Island, San José Island, Faro San José and the rest of the localities as separated groups, the differentiation among populations within groups (F_SC = 0.0005, p = 0.22) was no longer significant (Supplemental S5). Neither did any of these regional groupings support the current administrative zones (A) which resulted in no significant variation among groups (F_CT = 0.0003; p > 0.05) (Supplemental S5).

The UPGMA tree identified two to three major groups with a high bootstrap value (100%). The most robust topology separated the northern sample sites (Guadalupe and San Jerónimo island and Faro San José) from the rest (Supplemental S6).

Gene flow analysis

The directional relative migration using divMigrate (Fig. 4A) showed a predominant gene flow from northern locations (Guadalupe and San Jerónimo Island and Faro San José) to the rest of the southern sites. Despite the geographic proximity of San Jerónimo Island from Faro San José, migration between these locations was lower compared to the rest of the localities. A close examination of the connectivity among localities in the southern region showed no apparent predominant directional flow (Fig. 4B).

Isolation by distance

The Mantel test showed a strong relationship between linearized F_ST and geographical distances on the full data set (Mantel r = 0.76, p = 0.001) (Fig. 5). The regression analysis showed this relationship to be positive and linear (R² = 0.57, p < 0.001). This pattern of IBD remains significant without the most isolated location, Guadalupe Island (Mantel r = 0.54, p < 0.017; R² = 0.44, p = 0.005), but it was no longer significant when only the central-south locations were analyzed (Mantel r = 0.34, p = 0.129; R² = 0.23, p = 0.996). This pattern was further explored using 2-dimensional kernel density plots. The dataset—including all locations—showed a discontinuity indicating a patched pattern of genetic differentiation among populations (Supplemental S7A). However, the other scatterplots showed a single patch (Supplemental S7B and S7C). These results support the regional green abalone population structure along the BCP.

Management implications

An essential prerequisite of a sustainable fishery management is the match of biologically meaningful entities and management units. Assessing population status or developing a management strategy might not be possible without first knowing where a population begins and ends (Funk et al., 2012). Unfortunately, fish stocks do not often correspond to true biological populations (Reiss et al., 2009), which may lead to reduced productivity, genetic diversity loss, local population reduction, and in extreme cases, to a local population collapse. The risk for such loss depends both on the structure itself and the size of the component populations (Laikre, Palm & Ryman, 2005). Hence, understanding biocomplexity within stocks is important to maintain their resilience to future environmental change (Hilborn et al., 2003).

The abalone fishery in the BCP has been regulated mainly by minimum legal size, closed seasons and zones, annual quotas, and internal controls of the cooperatives (Diario Oficial de la Federación (DOF), 1993, 1994). A contemporary stock assessment uses reefs as separate stocks, and since 2000, management has mainly used catch quotas, which are assigned by reefs considering an annual growth rate for each reef (Sierra-Rodríguez et al., 2006). Nevertheless, the separate management units that include only a portion of a larger population may lead to population stock dynamics bias (Mullins et al., 2018) and then of the reference points, as well as stock recovery times (Hilborn & Walters, 1992; Caddy & Mahon, 1995; Sierra-Rodríguez et al., 2006). This study addressed a need asked by the Mexican National Institute of Fisheries in assessing and characterizing the genetic composition of H. fulgens’ reefs within and between zones (Sierra-Rodríguez et al., 2006). The results suggest that H. fulgens biological units do not support the administrative areas (Fig. 1) of the abalone fishery, so stock assessment and management areas should be revised. For instance, Guadalupe Island should be considered as an independent stock, separate from any of the coastal locations (putative Zone I). Locations in the northern part of the BPC (north of Guerrero Negro, covering San Jerónimo Island, and Faro San José) would constitute a second subpopulation (putative Zone II). All locations from Cedros Island to La Bocana (and Punta Abreojos) would be a third group (putative Zone III). Locations south of San Ignacio Lagoon, which were not sampled in this study, would keep its present classification as Zone IV, until further sampling and analysis are performed.

Nevertheless, since differences are evident in some biological and morphological characteristics, as well as productivity of H. fulgens banks along the four administrative zones (Ortíz-Quintanilla & León-Carballo, 1996; Vélez-Arellano, 2016; Guzmán-del Próo & Del Monte-Luna, 2017), the genetic results should be complemented with other techniques (morphological, reproductive and recruitment rates, etc.) to improve management decisions about the resource and help to better delineate the stock units.