Temporal uncoupling between demographic and genetic metrics in fisheries assessment: the European hake case study

Introduction

An authentic ecological transition for fisheries sustainability requires the integration of all relevant biological criteria in assessment models. Irrespective of what the population census size might be, maintaining the maximum amount of genetic variation needs to become a central goal of long-term management plans for fisheries (McNeely et al., 1990). In this regard, the genetic background of a fishery should be assessed at the same level of relevance as biomass, acting as a guarantee of its biological sustainability. Theoretically, a low census size can bear the same amount of genetic variation as a large census size, thus allowing population rebound from a reduced population size, and vice versa, a large census size can bear an extremely low genetic variation after overharvest or natural bottlenecks (Pinsky and Palumbi, 2014; Gurgel et al., 2020). The erosion of genetic diversity in exploited fisheries can be addressed with population metrics drawn from genotyping approaches, such as the genetically effective population size (N_e). The N_e was defined by Wright (1931) to approach the amount of genetic drift in populations with a binomial distribution of progeny numbers, segregated sexes, or fluctuations in population size, and was later extended to populations with overlapping generations. The genetic assessment of a fishery in terms of N_e should comprise measuring the erosion of genetic diversity caused by overfishing and random drift over several generations. In addition, N_e becomes a critical parameter in fishery management for its relationship to fitness risk. While predictions based on N_e equations may prove inaccurate for a number of reasons, e.g., estimates with a large variance and error and approximate formulation for neutral loci (Lande and Barrowclough, 1987), an insufficient number of SNP markers (Barbato et al., 2015), or small samples not covering the whole range of large populations (Santos-del-Blanco et al., 2022), in-silico trials have shown that a loss trend in temporal estimates of N_e is a warning flag of fishing-induced drift (Waples et al., 2008). That drift can bring about collateral effects such as the lack of both adaptive responses to stressful conditions (Bijlsma and Loeschcke, 2012) and low resilience to human pressure (Lopes et al., 2009; Munsch et al., 2022).

Coastal fisheries of the European hake (Merluccius merluccius, Linnaeus, 1758) extend from Mauritania to Iceland and Norway and from the western Mediterranean Sea to the Southern coast of the Black Sea (Lloris et al., 2005). This hake has been exploited since the 18th century and has been heavily harvested over the last three decades (Casey and Pereiro, 1995; ICES, 2009a; ICES, 2009b), which makes one suspect an erosion of its genetic background as well as a rarefaction of its population genetic structure. Namely, the Southern hake stock spanning the Cantabrian Sea and the Atlantic Iberian Peninsula (ICES management Divisions VIIIc and IXa, respectively) is devoid of genetic structuring, as shown by several genetic and geochemical marker types (e.g., Pita et al., 2011; Tanner et al., 2014). Nowadays, the population genetic challenge in this fishery is to manage the post-overharvesting genetic diversity as a precious asset for the sustainability of this fishery. Previous genetic studies on that fishery using microsatellites have shown postharvest erosion in both the genetically effective population size (N_e) (43-fold reduction) and the spawning stock biomass (N_SSB, 80% loss) (Pita et al., 2017). Such genetic evaluation pointed out that this stock passed through a genetic bottleneck of N_e ≈ 300 in the early 2000s and rebounded to N_e ≈ 3,200 in 2007. Letting alone the hypothesis of environmental facilitation to explain population growth (Díez et al., 2012), the EU regulatory measures adopted on this fishery seem to have succeeded (EC Reg. No. 2166/2005). However, for N_e to rebound from its historical minimum, both the demographic and genetic enrichment from the neighboring Northern stock, as well as the maintenance of the genetic status above a minimum evolutionary N_e threshold, have presumably played a synergic role in protecting the overfishing resilience of this species (Pita et al., 2017).

A general concern exists about the need for a deeper analytical horizon in fisheries assessment, comprising multiple data sources, long-term data sets, and the incorporation of genetic metrics (Cuéllar-Pinzón et al., 2016). Such an approach should allow us to distinguish meaningful genetic changes from sampling errors and environmental variability (Shapiro and Hofreiter, 2012). For instance, the short-term unpredictability of N_e following the post-overharvesting SSB plateau (Pita et al., 2017) might be reduced by increasing time-point estimates of N_e before and after the postharvesting SSB plateau. The adequate dissection of marker sets such as bulks of microsatellites or SNPs into distinct functional categories, prevents the enforcement of marker types not reflecting the extent of population dynamics processes, i.e., Fossil Adaptive Polymorphic Systems (FAPS) sensu Pita et al. (2022). Provided that the number of alleles per locus markedly outperformed other indicators across multiple in-silico situations, increasing the number of neutral microsatellites and sample number and size seems to dramatically increase power (Hoban et al., 2014). This study aimed to provide a less biased estimation of N_e from the Southern hake stock over four decades (1975–2014) as compared to a previous 7-year sample frame. This timely collection was screened with microsatellite types expected to provide a differential erosion signal caused by overharvesting (Pita et al., 2022). Also, historical records of molecular variation in the mtDNA cytochrome b gene (cyt b) were used as a conservative control under purifying selection (Lundy et al., 1999; Pita et al., 2014). The null hypothesis states that the temporal point where the minimum records of demographic and genetic metrics meet should be roughly the same, e.g., the mortality rate upon the global reproductive census (Z_N_SSB), the mortality rate of exploited cohort strata ($\overline{Z}\_ICES$), and the population genetic effective mortality (Z_LDN_e). Because the three metrics are definitively influenced by overfishing, alternate hypotheses suggest that their temporal uncoupling would be due either to methodological constraints or to their sequential interdependence in the time period comprised.

Materials and methods

Sampling and DNA extraction

The sample collection consisted of 295 individuals of Merluccius merluccius from the Southern hake stock, which spans the Cantabrian Sea in northern Iberia (Supplementary TableS1). The sampling range comprised from ICES Subarea VIIIc (San Sebastián, 43°27′00″N–2°00′00″W) to ICES Subarea IXa (La Guardia, 41°55′36.72″N–09°13′53.52″W). Samples were collected during research surveys carried out by the Spanish Institute of Oceanography (IEO) as well as by the authors onboard commercial fishing vessels. Samples from the period 1975 to 1999 consisted of dry-preserved otoliths and samples from the period 2000 to 2014 consisted of gill or muscle tissues preserved in 96% ethanol. DNA extraction from tissue was carried out following the method FENOSALT (Pérez and Presa, 2011), and DNA extraction from otoliths followed protocols previously described (Pita et al., 2016). DNA was stored upon extraction at −20°C in 50 µL TE 1× buffer.

DNA amplification and genotyping

Two microsatellite pentaplexes were PCR-amplified on DNA samples. One pentaplex consisted of five anonymous hake microsatellites, i.e., Mmer-hk3b, Mmer-hk9b, Mmer-hk20, Mmer-hk29, and Mmer-hk34b (Morán et al., 1999), which were coamplified at 55°C and 1.5 mM MgCl₂ using 5 pmol of each primer (Supplementary TableS2). The second pentaplex was formed by five EST microsatellites (Mmer-EST_6.3, Mmer-EST_8.1, Mmer-EST_9.4, Mmer-EST_13.1, and Mmer-EST_13.3) as characterized from an Expressed Sequence Tag (ESTs) library of the European hake (Milano et al., 2011; Pita et al., 2022) and coamplified at 59°C and 1.5 mM MgCl₂ using 3 pmol per primer. The first pentaplex was mostly (three out of five) composed of Adaptive Polymorphic Systems (APS) microsatellites while the second pentaplex comprised EST microsatellites mostly (four out of five) belonging to the Non-Adaptive Polymorphic Systems (NAPS) type (Pita et al., 2022).

The PCR reaction of 15 µL consisted of 0.2 mM dNTPS, 1× reaction buffer, 0.75 U Taq DNA polymerase (Biotaq, Bioline, Meridian Bioscience, Germany), and 50 ng DNA. The PCR routine set in a Mastercycler Gradient Thermocycler [Eppendorf SE (Germany)] consisted of 5 min at 95°C followed by 40 cycles of 70 s at 95°C, 70 s at the specific annealing temperature, and 80 s at 72°C. A final extension step was held for 20 min at 72°C. Genotyping was performed in an automatic sequencer (Applied Biosystems (USA), ABI PRISM 3130) from the Centre for Scientific and Technological Support to Research [CACTI, Universidade de Vigo (Vigo, Spain)] using 9.5 µL deionized formamide, 0.25 µL DNA from PCR amplicons, and 0.25 µL GeneScan 500 ROX Size Standard. Amplicons were sized with GeneMarker v2.6.0 [SoftGenetics LLC (USA)] and cross-checked in two independent readings.

Statistical analyses of microsatellite markers

Annual samples were grouped per lustrum as well as per decade to allow for a sufficient sample size to achieve meaningful results (Supplementary Table S1). The presence of null alleles was estimated with FreeNA (Chapuis and Estoup, 2007) and with Micro-Checker v2.2.3 (Van Oosterhout et al., 2004) (Supplementary TableS2). The existence of selection signatures on marker polymorphism was tested with BayeScan v2.1 (Foll and Gaggiotti, 2008). Allele frequencies were calculated with GenAlEx v6.5 (Peakall and Smouse, 2006). The Hardy–Weinberg equilibrium (HWE) per sample was tested with Genepop v4.3 (Rousset, 2008) through a Markov chain (MC) algorithm of 20 batches of 5,000 iterations each. The heterozygosity observed (H_O) and the one expected (H_E) were estimated with Genepop. The program FSTAT v2.9.3 (Goudet, 1995) was applied to estimate basic genetic diversity parameters such as the allelic richness (R_S) and the fixation index within samples (F_IS) as well as to compare the scores of those parameters among lustrums and among decades using 1,500 permutations. The temporal integrity of the genetic pool under study was assessed with BAPS 6 (Corander et al., 2004), assuming either a mixture model or an admixture model following the recommendations of the manual. The analysis of molecular variance (AMOVA) among hierarchical levels was carried out with Arlequin v3.5 (Excoffier and Lischer, 2010). The F_ST distance between samples was estimated with FSTAT v2.9.3. and the differentiation parameter D_EST (Jost, 2008) was calculated with the R-package DEMEtics v0.8-7 (Gerlach et al., 2010) upon 10,000 iterations. The false discovery rate approach (FDR, Benjamini and Hochberg, 1995) was applied to correct for type I error in multiple test comparisons.

Demographic data from the Southern hake stock was taken from ICES annual reports starting from 1982 (ICES Advice, 2019). The spawning stock biomass (SSB) for lustrum 1976–1980 was replicated from the SSB of 1982, and the SSB for lustrum 1981–1985 was averaged from the period 1982–1985. The population effective genetic size (N_e) was calculated using the linkage disequilibrium method (LDN_e) implemented in NeESTIMATOR v2.1 (Do et al., 2014) as a highly sensitive indicator of recent bottlenecks (England et al., 2010). An Approximate Bayesian Computation was used to assess the potential effects of overharvesting pressure on the genetic diversity of the Southern hake stock. Two scenarios were implemented in DIYABC 2.1 (Cornuet et al., 2014), the first one settled a population with equal N_e intervals at times t1 to t8 (10–10,000; lower and upper limits, respectively) and led the algorithm to decide its fluctuation upon observed diversity values; the second one implemented a reduction of N_e from time t3 to time t5 (i.e., between lustrums 1986–1990 and 1996–2000) with an algorithm decision on the maintenance or decrease of N_e along the four decades. The scenarios were evaluated using 100,000 simulated data sets on the prior distribution of parameters. Significant deviations from mutation-drift equilibrium were explored with program BOTTLENECK v1.2.02 (Piry et al., 1999). A linear interpolation of LDN_e estimates was applied to the open intervals associated with lustrums 1996–2000 and 2001–2005. The calculation of the census size of the SSB (N_SSB) was afforded on both the length-age analytical assessment (GADGET) (ICES, 2014) and the average weight of hakes landed in the sampled years (≈ 2 kg). This latter was calculated as the average of the fresh body weight range 372–4,587 g inferred after the size-weight relationship W = 0.00733L^2.981 and the observed range of landing lengths (37.9–88.0 cm) (Piñeiro and Saínza, 2003). The N_SSB declined per lustrum (t = 1 time unit) with a variable mortality rate ($\overline{Z}\_ICES$) comprising a constant natural mortality rate (M) and a variable per year fishing mortality rate (F). Thereafter, the total demographic mortality per lustrum ($Z\_N_{{\text{SSB}}_{Lt}}$) was estimated upon the SSB biomass estimated in the first year of the lustrum ($N_{{\text{SSB}}_{L0}}$) and the SSB biomass estimated in the first year of the subsequent lustrum ($N_{{\text{SSB}}_{L+1}}$), i.e.,

The instantaneous rate of change of the effective population size due to effective mortality is negative and proportional to N_e at a given time point:

and

and integrating between moments t₁ and t₂ where N_e1 and N_e2 were estimated, respectively,

we have

or

In our case study, a straightforward way to compute the population effective genetic mortality rate (Z_LDN_e) from current data employed $N_{{\text{e}}_0}$ as the N_e from the previous lustrum $(N_{{\text{e}}_{L-1}})$ and $N_{{\text{e}}_t}$ as the N_e from the proper lustrum $N_{{\text{e}}_L}$), i.e.,

Z_LDN_e is defined as the instantaneous rate of effective mortality, i.e., the instantaneous reduction of the dispersion of allele frequencies under random genetic drift as observed in the population under consideration. Intuitively, it measures the reduction of the population effective genetic size LDN_e in period t, i.e., between two infinitesimal LDN_e moments.

Official annual fishing mortality estimates (F) and the natural mortality rate (M) of 0.4 were taken from ICES (ICES Advice, 2019) as based on the inter-annual analyses of landed cohorts in order to calculate the per lustrum mortality average ($\overline{Z}\_ICES$). Mortality rates (Z_LDN_e, $\overline{Z}\_ICES$, and $Z\_N_{{\text{SSB}}_t}$) could not be calculated for the first lustrum (1976-1980) as being dependent on unavailable SSB and N_e from the previous lustrum.

Cytochrome b analyses

The PCR reaction mix for the cyt b gene was prepared in a volume of 25 µL using final concentrations of 2.5 mM MgCl₂, 0.3 mM dNTPs, 0.5 mM of each primer, 1.5 U Taq polymerase (Bioline), and 5% BSA. The volume of the DNA template varied on the sample tissue, i.e., 1 µL DNA from gill or muscle and 3 µL DNA from otoliths. The PCR amplification was performed in a 2720 Thermal Cycler (Applied Biosystems) as described in Pérez and Presa (2008) with subtle modifications, i.e., 3 min at 95°C, 40 cycles of 40 s at 94°C, 50 s at 61°C, and 80 s at 72°C, and a final extension step of 10 min at 72°C. The PCR products were purified with Exonuclease I [Thermo Fisher Scientific (USA)] before sequencing in an ABI PRISM 3130 from CACTI facilities (Universidade de Vigo). The sequences of cyt b were first visualized in CHROMAS Version 2.5.1 (http://www.technelysium.com.au/chromas.html) to discard potential errors and subsequently aligned with BioEdit v7.2.5 (Hall, 1999) using the 3′-end of the tRNA-Glu gen. Nucleotide frequencies were estimated in MEGA v6.0 (Tamura et al., 2013). The number of haplotypes, haplotypes, and nucleotide diversities (Nd and Hd), and the number of variable sites were calculated with DnaSP v5.10.01 (Librado and Rozas, 2009).

The temporal homogeneity of the haplotype distribution across lustrum was tested using a Spearman correlation test implemented in IBM SPSS version 23.0 (IBM Corporation). A hierarchical AMOVA was performed on cyt b variation using Arlequin. The molecular relationships between haplotypes were inferred after NETWORK 5.0.0.0 (Bandelt et al., 1999). The demographic history upon cyt b variation was explored using the mismatch distribution on pooled samples per lustrum, using observed versus simulated values under demographic expansion in Arlequin, e.g., unimodal distribution is expected in a population with rapid demographic growth in the recent past (Rogers and Harpending, 1992). In addition, Harpending’s (Harpending, 1994) raggedness index (Hri) and the sum of squared deviations (SSD) as implemented in Arlequin, were used to evaluate Rogers’ (Rogers, 1995) sudden expansion model (Rogers and Harpending, 1992). Further testing for population expansion was assessed through the Fu’s FS test (Fu, 1997) implemented in DnaSP and its statistical significance was calculated after 10,000 replicates of coalescent simulations.

Results

Genetic diversity of microsatellites

No genetic structuring was detected with BAPS within the Southern hake stock between 1975 and 2014 (data not shown). The polymorphism of the anonymous pentaplexed microsatellite set consisted of three APS (Mmer-hk9b, Mmer-hk20, and Mmer-hk34b), one FAPS (Mmer-hk3b), and one NAPS (Mmer-hk29b), therefore mostly composed of highly polymorphic adaptive polymorphisms. Markers Mmer-hk29 and Mmer-EST_9.4 showed signals of null alleles at frequencies higher than 10%, and further analyses were assayed with or without them (Supplementary TableS2). The polymorphism of the EST-microsatellite pentaplex consisted of four NAPS (Mmer-EST_8.1, Mmer-EST_9.4, Mmer-EST_13.1, and Mmer-EST_13.3) and one APS (Mmer-EST_6.3), therefore mostly composed of nonadaptive or neutral polymorphisms.

A total of 186 alleles were observed in the 295 individuals genotyped. The average number of alleles per locus varied between 6.500 ± 0.756 of locus Mmer-EST_8.1 and 28.500 ± 5.014 of locus Mmer-hk9b (Supplementary Table S3). The number of alleles per lustrum was higher (t-student = −7.41, one-tailed p = 2.02 × 10⁻⁵) in the anonymous maker set (77.00 ± 12.35) than in the EST marker set (42.75 ± 4.27). Likewise, the average number of alleles per locus was higher (t-student = −2.18, one-tailed p = 0.047) in the anonymous maker set (17.25 ± 7.64, excluding Mmer-hk3b, which showed signals of non-neutrality) than in the EST marker set (8.55 ± 2.51).

No correlation was observed between the allelic series of both marker sets across lustrums (R² = 0.0095, F = 0.058, p = 0.818). The highest and lowest average allelic richness (R_S) per locus were observed in lustrums 1986–1990 (R_S = 7.18 ± 2.40) and 1991–1995 (R_S = 6.38 ± 1.74), respectively. The main inflection points were highly influenced by the anonymous marker set whose trend did not correlate with that of EST microsatellites (R² = 0.017, F = 0.104, p = 0.758) (Figure 1). A significant deficit of heterozygotes was observed in 70% and 60% of loci in lustrums 2001–2005 and 2006–2010, respectively. The two loci showing signals of null alleles (Mmer-hk29 and Mmer-EST_9.4) also showed a deficit of heterozygotes across lustrums.

Figure 1 Pairwise comparison of allele richness (R_s) among lustrums for anonymous microsatellites (top) and EST microsatellites (bottom). Significant differences within marker type, between lustrums 1991 and 1995, and the rest are indicated with asterisks (^*p< 0.05; ^**p< 0.01).

Genetic differentiation of microsatellites

The molecular variation of microsatellites was low but significantly distributed among lustrums within decades (Table 1). The genetic distance F_ST was low (range, 0.000–0.016) and nonsignificant between lustrums (Table 2). The genetic differentiation D_EST differed significantly in all pairwise comparisons involving lustrum 1991–1995 as well as in three pairs involving lustrum 2011–2014. The accumulated estimate of the fixation index (F_ST) differed by one order of magnitude from the genetic differentiation per lustrum D_EST (Supplementary FigureS1). Both estimates showed quasi-parallel distributions and a maximum divergence among lustrums in 1991–1995. The correlation between both statistics was significant as per lustrum basis (Y = 5.002X + 0.032; R² = 0.706; F = 62.582; p< 0.001) but not significant as per decade basis (Y = 2.7961X + 0.0192; R² = 0.799; F = 15.907; p = 0.016).

Table 1 Hierarchical locus by locus analysis of molecular variance (AMOVA) on 295 individuals from eight lustrums of the Southern hake stock.

Table 2 Pairwise values of genetic differentiation among lustrums.

Effective population genetic size and demography

Estimates of demographic parameters per lustrum upon official SSB data (ICES Advice, 2019) showed a minimum N_SSB in lustrum 1996–2000, when the official mortality rate per lustrum $\overline{Z}\_ICES$ was maximal (Table 3). The relative LDN_e size and the ratio LDN_e/${\overline{N}}_{\text{SSB}}$ were minimal in lustrum 1991–1995 when the estimated population effective genetic mortality rate Z_LDN_e was maximal. The highest mortality rate calculated after the survival equation occurred in lustrum 1986–1990. The SSB curve adjusted significantly to a polynomic regression (R² = 0.9027, F = 23.1870, p = 0.003), and the temporal trends of N_SSB and LDN_e did not correlate to each other, e.g., the minimum LDN_e was in lustrum 1991–1995 and predated the minimum N_SSB in 1996–2000 (Figure 2). The three mortality rates were temporarily uncoupled, e.g., $\overline{Z}\_ICES$ was maximal in lustrum 1996–2000 when Z_N_SSB was null and Z_LDN_e reached its minimum (Figure 3). The global mortality rate of the population $Z\_N_{{\text{SSB}}_t}$ as well as the population effective genetic mortality rate Z_LDN_e were negative from 1996 to 2010, while cohort-based mortality $\overline{Z}\_ICES$ was always above 0.82. The approximate Bayesian computation observed in both scenarios, but mainly in the second one (i.e., using the inferior limit of observed N_e intervals) suggested a population contraction between lustrum 3 (1986–1990, N_e = 5,270) and lustrum 4 (1991–1995, N_e = 994) and a population expansion flanking such an interval (Figure 4).

Table 3 Estimates of the population effective genetic size LDN_e per lustrum using the linkage disequilibrium method (Waples and Do, 2010) and several mortality rates of the Southern hake stock.

Figure 2 Temporal distribution of the spawning stock biomass (SSB, open circles) and the genetic effective size (LDN_e, closed squares) of the Atlantic European hake Merluccius merluccius in the period 1974–2014. The error bars of the $\overline{\text{SSB}}$ trend indicates its standard deviation within lustrum.

Figure 3 Temporal trends of mortality rates exerted on the Southern Atlantic European hake stock. $\overline{Z}\_ICES$, total mortality estimated by ICES (closed squares); Z_LDN_e, global population effective genetic mortality (open circles); $Z\_{\overline{N}}_{{\text{SSB}}_t}$, global demographic mortality rate calculated after SSB estimates (shaded area); estimated SSB (tons; open triangles).

Figure 4 Main temporal milestones of demographic and genetic parameters estimated on the Southern hake stock from 1975 to 2014. Red cells and green cells indicate the maximum erosion and the minimum erosion intervals of each parameter on the left, respectively (see Table 1). The superimposed trend line shows the summarized posterior median values obtained from DIYABC as implementing the observed N_e population contraction in lustrum 1991–1995 (10–10,000; lower and upper limits, respectively) and the temporal flanking expansions down to lustrum 1975–1980 (10–10,000) and up to lustrum 2011–2014 (10–10,000) with maintenance or decrease of N_e during the last four decades.

The population effective genetic mortality Z_LDN_e comprises the natural effective genetic mortality (M_e) and the fishing effective genetic mortality (F_e). The simulation plot shows the expected behavior of Z_LDN_e as calculated after formula (2) against the ratio between two N_e moments N_et/N_et0 (Figure 5), i.e., positive genetic mortalities occur when N_e at moment t decays with respect to N_e at the initial moment t₀ and vice versa. The inflection point where N_et = N_et0 implies the conservation of N_e between moments, i.e., zero genetic mortality.

Figure 5 Relationship between the effective population genetic mortality Z_LDN_e (the natural effective genetic mortality (M_e) plus the fishing effective genetic mortality (F_e)) as integrated in moments t = 1, 2, 10, and the ratio N_et/N_et0 between N_e moments (N_et0, initial moment; N_et, final moment).

Genetic diversity of cytochrome b

The 465 bp fragment of the cyt b gene from 173 specimens presented 15 variable sites, which consisted of six unique substitutions and nine parsimonious substitutions, defining 22 haplotypes (Supplementary Table S4). The haplotypes Mmcytb.1, Mmcytb.2, and Mmcytb.3 were the most frequent ones across studies and years. The haplotypic diversity (Hd ± SD) was 0.569 ± 0.040, and the nucleotide diversity (Pi ± SD) was 0.0023 ± 0.00026. The mismatch distribution used to explore the demographic history of the Southern hake stock showed a global adjustment of observed values to simulated values under a recent past demographic expansion in the period 1975–2014 (Supplementary Figure S2A). The sum of squared deviations (SSD) and the Harpending’s raggedness index (Hri) used to evaluate Rogers’ sudden expansion model were low (< 0.10) and nonsignificant across lustrums. Global Fu’s F_S parameter was significant under the hypothesis of nonexpansion (F_S = −20.817, p< 0.00001). The mismatch distribution for lustrum 1996–2000 deviated from its expected values under demographic expansion, i.e., SSD = 0.40121, p = 0.002 (Supplementary Figure S2B). The global mutational network of cyt b showed a large number of low-frequency haplotypes in the period 1975–2014 (Figure 6).

Figure 6 Molecular network for cyt b haplotypes as built with NETWORK v10.2.0.0 (http://www.fluxusengineering.com). MM1-27 is the code for observed haplotypes, while code mv1–8 is for haplotypes inferred under sampling drift.

Discussion

Current analyses aim to dissect the evolution of demographic and genetic metrics in the Southern hake fishery from the last 40 years. Demographic data were drawn from official records of ICES, while genetic data were generated in the laboratory on a lustrum-based sampling series between 1975 and 2014. Despite the uncertainty on the upper limit of some confidence intervals, the LDN_e method requires a single population sample and incorporates a correction for the downward N_e bias caused by low sample sizes (Waples and Gaggiotti, 2006; Waples and Do, 2010). We show that type APS microsatellites (most microsatellites of the anonymous pentaplex) were superior to those of type NAPS (most microsatellites of the EST-pentaplex) to detect genetic erosion at the maximum population decline in lustrum 1996–2000. That contraction point was also made patent with the molecular variation of the cyt b. We also show that the maximum demographic mortality Z_N_SSB (1986–1990) predated both the maximum genetic effective mortality Z_LDN_e (1991–1995) and the maximum cohort-based mortality $\overline{Z}\_ICES$ (1996–2000). That lack of temporal congruence between demographic and genetic parameters suggests that official methods based on cohort analyses are insufficient to assess both the real biological status of fish stocks and the genetic consequences of sustained fishing pressure.

Genetic diversity

The single gene pool (k = 1) inferred after Bayesian methods in the period 1975–2014 is congruent with independent data on dispersal patterns of hake juveniles along the Cantabrian shelf (Sánchez and Gil, 2000). Those data indicate that high connectivity exists over the main ground of the Southern hake stock, which spans 1,000 km coast from La Guardia to San Sebastian (northern Iberian coast) (e.g., Pita et al., 2011). Such a structural property is very basic before allowing us to assess shifts in the background gene pool of this stock. A previous genetic study on a 7-year sampling body reported a nonapparent erosion of polymorphism in anonymous microsatellite loci (Pita et al., 2017). The stability of elemental genetic metrics (N_a, H_E, and R_S) was counterintuitive, facing evidence of the historical overharvesting of this fishery (ICES, 2014). The apparent lack of erosion can result from the insensitivity of highly variable microsatellites to biomass reduction due to sampling drift (Pita et al., 2017). Moreover, a change in population size following a moderate bottleneck can have half the power to detect such a bottleneck (Hoban et al., 2013). Here, we show that the capacity to detect erosion depends not only on those two factors but also on both: a saturated sample design along the 40-year period comprising time point estimates of N_e before and after the postharvesting SSB plateau and an adequate selection of the genetic marker set.

The minimum allelic richness (R_S) per locus in lustrum 1991–1995 was due to erosion of polymorphism in the anonymous microsatellite pentaplex (APS marker type), whose amount of variation across lustrums did not correlate with that of the EST-microsatellite pentaplex (NAPS marker type). Moreover, despite genetic parameters showing range values consistent with previous studies in Atlantic hake, e.g., F_ST = 0.008 (Tanner et al., 2014), the systemic low polymorphism of EST microsatellites as compared to APS markers suggests that only these latter are adequate to score useful variation to assess population dynamic processes within fisheries (Pita et al., 2022). The analytical trials showed little sensibility to the inclusion or exclusion of markers Mmer-hk29 and Mmer-EST_9.4, which exhibited a high heterozygote deficit across lustrums. That deficit was putatively caused by null alleles, as suggested in previous studies using those loci, e.g., Mmer-hk29 (Lundy et al., 2000). However, the major deficit of heterozygotes in the temporal series observed in lustrums 2001–2005 and 2006–2010 not only ascribes to the null-allele hypothesis but also to a global deficit of heterozygotes following a postoverharvest population expansion, i.e., an opposite effect to the general excess of heterozygotes expected after a recent bottleneck (Maruyama and Fuerst, 1985; Cornuet and Luikart, 1996; Luikart and Cornuet, 1998). Regarding the temporal landmarks of the genetic markers used to measure overexploitation (see Figure 4), we observed that (1) the molecular erosion of neutral microsatellites predated that of cyt b and that (2) a deficit of heterozygotes could only be detected after the big population contraction. Those results suggest that neither the conservative molecular variation of cyt b nor heterozygosity are useful as immediacy indexes of overharvesting. However, allele diversity and genetic effective size from APS microsatellites seem to be adequate to detect postharvesting genetic erosion in real time (e.g., Yorisue et al., 2020).

Temporal genetic differentiation

Occasionally, the coefficient F_ST on microsatellite data has revealed a significant differentiation between samples over decades, e.g., snappers (Hauser et al., 2002). However, an ocean-wide genetic homogeneity is the scenario currently reported for many marine fishes (Poulsen et al., 2006). The tiny but significant distance (F_SC = 0.007*) among lustrums within this hake stock can be explained by sudden shifts in its genotypic background as a by-product of randomness and overharvesting. Such lustrum-based occasional divergence disappeared upon biomass rebound in successive years, making the temporal genetic divergence unnoticed over decades (F_CT = 0.002). That result highlights the need for long-term fine-grain sampling to detect population sinks that are buffered on a wide temporal scale (Cuéllar-Pinzón et al., 2016). The significant regression between F_ST and D_EST indicates the suitability of both indexes to identify the higher divergence in pairwise comparisons involving lustrums 1991–1995 and 2011–2014. Nevertheless, the statistical significance of D_EST makes it more adequate than F_ST to depict subtle genetic divergence among temporal samples (Jost, 2008) as has also been observed in other marine taxa (e.g., Faria et al., 2018). The genetic divergence between lustrum 1991–1995 and the rest of lustrums was due to a significant allele reduction in that period, likely caused by the maximum overharvesting pressure on this stock (ICES, 2014). Opposite, the cumulative divergence between lustrum 2011–2014 and period 1990–2005 was likely caused by the rebound of the stock due to technical and political EU measures adopted to protect this fishery (Council Reg. No 1162/2001, 2602/2001, and 2166/2005) and putatively helped by migration from the neighboring northern hake stock (Pita et al., 2011).

Despite the observed temporal stability of molecular variation in wide temporal frameworks of exploited marine fishes, e.g., Gadus morhua (Ruzzante et al., 2001; Poulsen et al., 2006; Therkildsen et al., 2010), Raja clavata (Chevolot et al., 2008), Clupea harengus (Larsson et al., 2010), and Solea solea (Cuveliers, 2011), current fishery metrics reflect the demographic erosion as well as the genetic erosion that overharvesting brought about on this hake stock. From the genetic viewpoint, the erosion was patent from nuclear DNA markers as much as from the mitochondrial DNA markers. For instance, the mismatch distribution in the cyt b variation indicates a steady demographic expansion in the period 1975–2014 with a strong population contraction in lustrum 1996–2000 (see Supplementary Figure S2). Such population expansion flanking the contraction is patent in the mutational network of cyt b, where a large bulk of low-frequency haplotypes appeared in that period (Figure 6).

The effective reproductive size

SSB is a demographic metrics used to assess the population status and its reproductive potential. Estimates of SSB are dependent on the algorithm adopted, which is permanently under revision (http://www.ices.dk). Therefore, current SSB data represent a rough approach to the status of the Southern hake biomass. SSB data ranged from ${\overline{N}}_{\text{SSB}}$ ≈ 21.5 mill. (1975–1980) to a minimum of ${\overline{N}}_{\text{SSB}}$ ≈ 3.5 mill. (1996–2000) and slowly rebounded thereafter up to ${\overline{N}}_{\text{SSB}}$ ≈ 7.5 mill. fish (2015–2019). A biomass rebound is good news for this stock and indicates it still holds reproductive plasticity (Pita et al., 2017). However, in absolute biomass terms, the recovery pace was slow under continued exploitation, suggesting the erosion of its reproductive potential by overharvesting. Noteworthy, the SSB rebound up to 15,102 t or ~7.5 mill. fish by 2014 (ICES, 2015) was about one-third of the historical size of this stock as scored in the 1970s. That situation agrees with the outcome of demographic simulations on Atlantic cod, in which one should not be too optimistic about reaching preharvest biomass levels after episodes of sustained overharvest (Kuparinen et al., 2016).

The biomass rebounding potential is not only a property of the census size but also of the genetic background variation aimed at being mirrored by LDN_e. This latter was relatively high until 1990, was down-eroded to its minimum in lustrum 1991–1995, when the effective genetic mortality rate Z_LDN_e approached its historical maximum and rebounded under relaxed fishing intensity from 2000 on. As fishing intensity decreased, LDN_e also increased, likely due to a lower reproductive variance and a higher probability of survival of immigrants from the Northern hake stock, making the ratio N_e/N_SSB rebound very rapidly. Theoretically, the low LDN_e observed at the overharvest maximum in 1991–1995 could be limiting for subsequent generations (Vucetich et al., 1997). Nevertheless, the genetic erosion caused by fishing drift can be partially compensated by gene flow (Fraser et al., 2004; Palstra and Ruzzante, 2008), as previously shown in this species (Pita et al., 2016). Indeed, all Bayesian simulations support each other in rebuilding the evolutionary history of the effective size of this stock and highlight the maximum N_e contraction in the decade 1990–2000 followed by a rapid expansion (see Figure 4) that is fully congruent with other genetic metrics discussed above.

The uncorrelated measures of mortality

The maximum demographic mortality rate (Z_N_SSB) in 1986–1990 predated the minimum of both, the genetic richness R_S and the effective size N_e in lustrum 1991–1995. However, the cohort-based mortality $(\overline{Z}\_ICES$) was maximal only from 1996–2000 on, when the other two mortality rates were negative or low and other lower-pace genetic effects were scored. Namely, the minimum molecular diversity of cyt b and the general deficit of heterozygotes of microsatellites were the latest genetic effects observed after the overharvest maximum (see Figure 4).

The temporal uncoupling of demographic and genetic metrics suggests that they depend causally on each other, i.e., the maximum demographic mortality (1986–1990) predated the overharvesting-mediated genetic erosion in the consecutive lustrum (1991–1995), while this latter predated both the official cohort-based SSB minimum (1996–2000) as well as the population inability to reach preharvest biomass levels (2011–2014). Alternatively, the temporal uncoupling between those population metrics could simply be due to their calculation on distinct population fractions, i.e., a random genetic screening vs. a cohort-based SSB demographic assessment. The global scenario suggests that (1) the large fishing pressure exerted on this stock in 1986–1990 was detected in real time using the classic survival equation approach on Z_N_SSB, (2) the most drastic genetic consequences of fishing pressure were detected in the following lustrum 1991–1995, and (3) the official method applied in $\overline{Z}\_ICES$ to infer fishing mortality on age classes provides biased estimates of the real impact exerted by fishing pressure on the whole population. Provided that $\overline{Z}\_ICES$ only detected deep overharvesting one decade later than other methods (Z_N_SSB), its methodology seems to largely underestimate overfishing consequences for age classes other than the fishing target cohorts. Therefore, a more comprehensive incorporation of meaningful biological data in fisheries stock assessment (e.g., Scott et al., 2016) seems badly needed.

The population effective genetic mortality Z_LDN_e is the theoretical addition of the natural effective genetic mortality (M_e) and the fishing effective genetic mortality (F_e). The simulation plot (Figure 5) shows that positive genetic mortalities occur when N_et decays with respect to N_et0 and that the inflection point at which N_et = N_et0 implies the conservation of N_e between moments, i.e., zero genetic mortality. However, Z_LDN_e = 0 does not preclude that a positive Z_LDN_e is compensated by other factors, such as the effective migration m_e. In turn, negative genetic mortality rates are allowed, implying either a significant effective migration exceeding the genetic erosion or a reduction of the reproductive variance by natural o sampling causes. Noteworthy, a null Z_LDN_e score tells nothing about the quality of the genetic variation, i.e., the substitution of polymorphisms between temporal estimates.

The ratio N_e/N_SSB

The metrics N_SSB is a proxy to the effective reproductive size, and N_e/N_SSB becomes a critical parameter to analyze trends in the genetic effective size of the breeding fraction under population fluctuations. Because of the large uncertainty on population size, overlapping generations, sex ratios, nonrandom mating, variation in family size, and habitat-related productivity (Nunney, 1999; Hare et al., 2011), empirical N_e/N_SSB is much less than expected from demography (Frankham, 1995; Vucetich et al., 1997). That ratio in hake is orders of magnitude less than current estimates of N_SSB, i.e., N_e/N_SSB< (10⁻⁰⁵–10⁻⁰⁴), as also reported in diverse teleosts (10⁻⁰²–10⁻⁰⁶, Hauser and Carvalho, 2008), Cuplea harengus (10⁻⁰⁷–10⁻⁰⁸, Larsson et al., 2010), Gadus morhua (10⁻⁰⁵, Hutchinson et al., 2003), Pagrus auratus (10⁻⁰⁵, Hauser et al., 2002), Raja clavata (10⁻⁰⁵–10⁻⁰⁶, Chevolot et al., 2008), Sciaenops ocellatus (10⁻⁰⁴, Turner et al., 2002), and Solea solea (10⁻⁰⁵, Cuveliers, 2011). We assume the population equilibrium of this hake stock at its pre-overharvest status (1975), i.e., with N_e maximal, N_SSB constant, and N_e/N_SSB at a relative maximum. Massive fishing mortality in 1986–1990 rapidly eroded LDN_e, and N_e/${\overline{N}}_{\text{SSB}}$ reached its minimum in 1991–1995 well before the official SSB minimum in the following lustrum. These data confirm that a stock does not need to be demographically exhausted to exhibit very low levels of genetic diversity (Pinsky and Palumbi, 2014). The LDN_e minimum in 1991–1995 lets us believe that this hake stock traversed a genetic bottleneck in that period due to a drastic increment of the reproductive variance in both sexes. Therefore, fading transient population structures are likely structural sweepstake artifacts (Hedgecock, 1994) that become more pronounced at overharvest (Pita et al., 2017). The absence of correlation between N_e and SSB in the period 1975–2014 denotes the genetic insensitiveness of SSB as affected by compensatory recruitment (Anderson, 1988), as well as the N_e fluctuation with reproductive variance and unequal sex ratios (Waples, 2002). Despite such a lack of correlation, current data suggest that management warnings could have been triggered if the minimum N_e/N_SSB ratio scored on lustrum 1991–1995 (N_e/N_SSB = 3.0 × 10⁻⁶) were taken into consideration.

Conclusions

The population genetic analysis of the Southern hake stock in the period 1975–2014 indicates that its gene pool has been eroded at the overharvest maximum in 1991–1995. Such molecular erosion could be early detected with the large polymorphism of APS microsatellites as opposed to the less within-stock polymorphic NAPS markers. Just-on-time EU fishing policies, the maintenance of a minimum N_e threshold that prevented collapse by genetic factors, and a compensatory effect of gene flow from the Northern hake stock (Pita et al., 2016) likely played synergically on this firm-stock recovery after 2005. However, such recovery is not expected to regain preharvesting levels after episodes of sustained overharvest (Kuparinen et al., 2016). Stock assessment based solely on demographic models or their subrogations (e.g., $\overline{Z}\_ICES$) underestimates the genetic erosion that overfishing exerts on the reproductive potential of a fish stock. The incorporation of different data sources and temporal genetic controls would bring about a more realistic and integral biological assessment (see Basset et al., 2001) and allow for the assessment of the risks of extinction by genetic factors. Long-term, spatially comprehensive, and systematic sampling designs are required to assess the predictability and usefulness of Z_e, N_e, and N_e/N_SSB. Those metrics are promising candidates to assess genetic sustainability under a blue transition agenda for fisheries.

Data availability statement

The data presented in the study are deposited in the Figshare repository, https:/doi.org/10.6084/m9.figshare.22761986.v1.

Ethics statement

Ethical approval was not required for the study involving animals in accordance with the local legislation and institutional requirements because samples were taken from fishing campaigns.

Author contributions

PP designed and conceived the study. AG organized and executed field sampling and tissue preservation of fish samples during CSIC-IEO fishing campaigns. AP conducted the development of EST microsatellites starting from 454 and GAII sequence data. MF-M performed the preliminary in-silico calibration of multiplex assemblages. AP and MF-M genotyped the sample set. PP drafted the manuscript with contributions from the coauthors. All authors contributed to the article and approved the submitted version.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was financed with contract TED2021–132258BI00 from MCIN/AEI/10.13039/501100011033 and The European Union NextGenerationEU/PRTR. The Regional Government of Xunta de Galicia supported the Ph.D. grant (ED481A-2017/362) of MF-M.

Acknowledgments

The authors warmly thank Angie Maribeth Santafé-Muñoz and José Andrés Cuéllar-Pinzón for their technical help in the laboratory.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2023.1214469/full#supplementary-material

References