The sample collection consisted of 295 individuals of Merluccius merluccius from the Southern hake stock, which spans the Cantabrian Sea in northern Iberia ( Supplementary Table S1 ). 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.

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 Table S2 ). 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.

Annual samples were grouped per lustrum as well as per decade to allow for a sufficient sample size to achieve meaningful results ( Supplementary T able 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 Table S2 ). 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 ( Z ¯ _ I C E S ) 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 SSB L t ) was estimated upon the SSB biomass estimated in the first year of the lustrum ( N SSB L 0 ) and the SSB biomass estimated in the first year of the subsequent lustrum ( N 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 e 0 as the N _e from the previous lustrum ( N e L − 1 ) and N e t as the N _e from the proper lustrum N 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 ( Z ¯ _ I C E S ). Mortality rates (Z_LDN _e, Z ¯ _ I C E S , and Z _ N 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.

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.

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 Table S2 ). 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 T able 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.

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 Figure S1 ). 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).

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 Z ¯ _ I C E S was maximal ( Table 3 ). The relative LDN _e size and the ratio LDN _e/ N ¯ 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., Z ¯ _ I C E S 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 SSB t as well as the population effective genetic mortality rate Z_LDN _e were negative from 1996 to 2010, while cohort-based mortality Z ¯ _ I C E S 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 ).

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 _{e t} /N _{e t 0} ( 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 _{e t} = N _{e t 0} implies the conservation of N _e between moments, i.e., zero genetic mortality.

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 ).

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).

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 F igure 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 ).

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 N ¯ SSB ≈ 21.5 mill. (1975–1980) to a minimum of N ¯ SSB ≈ 3.5 mill. (1996–2000) and slowly rebounded thereafter up to N ¯ 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 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 ( Z ¯ _ I C E S ) 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 Z ¯ _ I C E S to infer fishing mortality on age classes provides biased estimates of the real impact exerted by fishing pressure on the whole population. Provided that Z ¯ _ I C E S 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 _{e t} decays with respect to N _{e t 0} and that the inflection point at which N _{e t} = N _{e t 0} 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 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/ N ¯ 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.