Samples of archived otoliths for large yellow croaker were collected through a series of research surveys by the Institute of Oceanology, Chinese Academy of Science (IOCAS, China) in 1960s and 1980s. Otoliths were preserved dry in paper envelopes at room temperature. Contemporary samples were collected from 2016 to 2020 and fin clips or muscular tissue of each individual were preserved in 95% ethanol. All the samples were collected from different locations in the East China Sea ( Figure 1 ). A total of 130 archived otoliths and 67 contemporary specimens were collected ( Table 1 and Figure 1 ).

Genomic DNA was extracted from the archived otoliths using the EDTA-SDS extraction protocol (Firidin et al., 2017) and subsequently stored at -20°C. For contemporary samples, DNA was extracted from fin or muscular tissue by standard phenol-chloroform extraction and ethanol precipitation (Taggart et al., 1992). Microsatellite discovery and primer design were performed with QDD 3.1.2 (Meglecz et al., 2010) by using the large yellow croaker reference genome (GenBank assembly accession: GCA_000972845.2) on a local Galaxy plarform with default parameters (Afgan et al., 2016). In order to enhance the ability to amplify degraded historical DNA and ensure polymorphism of microsatellite loci, sequences with microsatellites (pure tandem repeats of di- nucleotide motif with at least eleven uninterrupted repeats) and targeted length ranging from 100-180 bases were selected for primers design. To ensure that the analysed microsatellite loci were unlinked, we selected loci on different chromosomes in the large yellow croaker genome ( Table S1 ). Finally, 7 polymophic microsatellite loci were selected for subsequent analysis ( Table S1 ). The amplification of the 7 microsatellite loci was performed by the PCR protocol described in Xue et al. (2017). All the amplified loci were analysed through an automated capillary sequencer ABI (Applied Biosystems) and allele sizes were determined with the GS500-ROX size standard GeneMarker 2.2.0 (SoftGenetics, State College, USA).

Historical DNA is particularly prone to genotyping errors due to its degradation and potential contamination from exogenous DNA, which increases the risks for allele dropout and null alleles. To address these issues, a standardized protocol was used for DNA extraction and genotyping. Firstly, all extractions and preparation of PCRs for historical samples were carried out in a laminar flow clean bench equipped with periodic sterilization. Negative controls were included at all stages. Secondly, DNA extractions and PCRs were performed independently for different samples, in order to avoid potential contamination of DNA between samples. In addition, extraction and genotyping of contemporary samples were conducted after historical samples. Specimens that were failed in genotyping after three trials were discarded.

Microstellites were tested for potential presence of genotyping errors due to large allele dropout, stuttering and null alleles by using Microchecker 2.2.1 (Van Oosterhout et al., 2004). Genetic diversity indices including number of alleles (Na), average observed heterozygosity (Ho), average expected heterozygosity (He) were calculated for each sample using the GENALEX 6.5 (Peakall and Smouse, 2006). Measures of genetic diversity including allelic richness (Ar) based on the smallest population size and inbreeding coefficient (F_IS ) were calculated using FSTAT 2.9.3 (Goudet, 2001). In order to test whether there was temporal changes of genetic diversity, one-way analysis of variance (ANOVA) were performed on Ar, He, F_IS among temporal samples. Departure from Hardy-Weinberg equilibrium (HWE) were tested using Genepop 4.3 (Rousset, 2014). A sequential Bonferroni correction of p-values was applied to account for an increase in type-I error from multiple comparisons (Rice, 1989).

As a precautionary measure against low sample sizes and the limited number of microsatellite loci employed in this study, we used POWSIM 4.1 (Ryman and Palm, 2006) to assess the statistical power of our dataset ( Table S2 ). Different levels of F_ST values (0.001 to 0.05) were generated using an effective population size of 2000 with 200 runs of 1000 iterations. The null hypothesis of genetic homogeneity was tested with Fisher exact test and chi-square test. The proportion of significances (p < 0.05) represented an estimate of the power.

Global and pairwise F_ST values were calculated using Genepop 4.3 (Rousset, 2014) and GENALEX 6.5 (Peakall and Smouse, 2006) respectively. The significance of each pairwise F_ST was assessed by 10,000 bootstrap permutations. The p value of pairwise F_ST was adjusted using a Bonferroni correction (Rice, 1989). Population genetic structure was assessed through Bayesian clustering of microsatellite genotypes as implemented in STRUCTURE 2.3.4 (Pritchard et al., 2000). Ten independent runs were conducted for different numbers of genetic clusters (K=1 to K=7), each with a burn-in period of 50,000 chains followed by 500,000 MCMC iterations, and with no prior knowledge of sampling locations. Admixture models with correlated allele frequencies were applied. To determine the most probable value of K, the Delta K method was used as implemented in Structure Selector (Li and Liu, 2018). Subsequently, the software CLUMPAK (Kopelman et al., 2015) was used to create a graphical representation of STRUCTURE outputs. Principle Components Analysis (PCA) was also conducted to assess genetic structure based on allele frequencies by R package adegenet 1.2-2 (Jombart, 2008).

Analysis of molecular variance (AMOVA) was performed in GENALEX 6.5 (Peakall and Smouse, 2006) to check how the genetic variations was partitioned temporally (among periods) and spatially (among locations). We conducted two AMOVA tests, which partitioned genetic variation (1) among the three sampling periods (1960s, 1980s and contemporary samples), emphasizing temporal variation, and (2) among the two geographic groups (historical Lvsiyang, historical Zhoushan, historical Yangtze River Estuary, contemporary Lvsiyang, contemporary Zhoushan vs. historical Dongtou, contemporary Wenling), focussing on spatial variation.

A redundancy analysis (RDA) was performed to investigate the relationships between genetic structure and spatio-temporal factors (latitude and longitude, year of sample collection) (Lepš, 2003). A set of spatial factors (Spatial1, Spatial 2 and Spatial 3) were generated using the distance-based Moran’s eigenvector maps (dbMEM), based on Cartesian coordinates generated from latitude and longitude coordinates. In addition, temporal factor was confirmed by using sampling time directly. Variation partitioning analysis (VPA) was conducted to address the relative contribution of the spatial and temporal factors and their conbined effect on shaping genetic structure (Peres-Neto et al., 2006). The relative contribution of both components could be explained by pure spatial variables (S|T), pure temporal variables (T|S), spatial variables (S), temporal variables (T) and the combined effects of both spatial and temporal variables (S∩T), respectively. The residual proportion represents the unexplained variance. Significance testing for each variable were conducted with permutation test. All these analyses were performed in the R software (R Core Team, 2019).

Effective population size (Ne) is an important element in efforts to conserve genetic diversity (Poulsen et al., 2006). Ne determines levels of inbreeding and magnitude of genetic diversity erosion for populations due to random genetic drift (Poulsen et al., 2006). Therefore, obtaining reliable estimates of Ne for large yellow croaker populations is of primary importance to evaluate whether populations are able to maintain sufficient genetic variation.

Ne was calculated based on linkage disequilibrium (LD) method implemented in NeEstimator v2 (Waples and Do, 2008; Do et al., 2013), which was considered to be an estimator of inbreeding Ne (Luikart et al., 2010). The method calculates the correlation among alleles at unlinked loci and corrects for downward bias associated with small sample sizes (Waples, 2006). Confidence intervals were determined by jackknifing on loci and the allowed frequency of observed alleles was set at 0.05 because alleles occurring at low frequencies may lead to biased results (Waples and Do, 2008).

Missing genotypes accounted for 0.073% of the overall dataset. There was a lack of allele dropout and stuttering throughout the entire dataset. The occurrence of null alleles was detected in 15 of 49 tests, involving both contemporary (6 instances) and historical samples (9 instances). Since no loci showed systematic presence of null alleles, all loci were used for subsequent analyses. Out of 49 locus-population pairs, 10 pairs were significantly deviated from HWE after Bonferroni correction (n=49, α=0.00102).

The seven microsatellite loci analysed were highly polymorphic across all the samples ( Table 2 ). Average Ar ranged from 8.704 in historical Yangtze River Estuary to 10.725 in historical Dongtou. Average expected heterozygosity (He) across all loci ranged from 0.725 in historical Zhoushan to 0.837 in contemporary Wenling. Average observed heterozygosity (Ho) ranged from 0.875 in historical Zhoushan to 0.707 in contemporary Lvsiyang. The inbreeding coefficient (F_IS ) varied from -0.02 in historical Yangtze River Estuary to 0.169 in contemporary Lvsiyang. No statistically significant difference were detected for the average Ar among temporal samples (p = 0.211 for 1960s vs. contemporary samples; p = 0.106 for historical vs. contemporary samples). Average He was similar among all temporal samples (p = 0.204 for 1960s vs. contemporary samples; p = 0.153 for historical vs. contemporary samples). In addition, we found no significant temporal changes for F_IS (p = 0.452 for 1960s vs. contemporary populations; p = 0.249 for historical vs. contemporary samples). These results indicated that there was no depletion of genetic diversity for large yellow croaker stocks in the East China Sea despite overexploitation.

The POWSIM analyses showed that the combination of the 7 microsatellite loci and sample sizes provided a statistical power sufficient to detect a low level of differentiation. Indeed F _ST = 0.005 or larger could be detected with a probability of at least 96% across all samples (Chi and Fisher exact test; Table S2 ) and the estimated F _ST in our data set ranged from 0.005 to 0.027. Thus the microsatellite loci were powerful enough to differentiate populations.The global F_ST was 0.0144, which reflected a shallow genetic differentiation across samples in the East China Sea. The pairwise F_ST were low, ranging from 0.005 to 0.027 and significant genetic differentiation was detected for only six comparisons ( Table 3 ). The STRUCTURE analysis showed highest support for two (K = 2 with the highest Delta K) clusters among the samples, with three historical samples (LSYhis, ZShis and CJKhis) forming one genetic group and the rest four samples grouping into another one ( Figure 2 ). Most of the significant genetic differentiation were from comparisons between historical and contemporary samples, which was consistent with the estimates of F_ST . No significant genetic differentiation was observed among contemporary populations, suggesting the existence of a single panmictic population across its distribution. However, the PCA analysis showed that all the 7 samples formed a single cluster and overlapped with each other, which was consistent with the low F_ST estimates ( Figure 3 ). The two hierarchical (temporal grouping and geographical grouping) AMOVA demonstrated that over 98% of genetic variation was within populations and only 0.76% and 0.72% of the total genetic variation was attributed to temporal structure and spatial structure, respectively ( Table 4 ). F_CT was non-significant among either geographic groups (p = 0.076) or periods groups (p = 0.134).

The results of RDA suggested that temporal factors played a similar role in explaining genetic structure compared with spatial factors ( Figure 4A ). VPA ( Figure 4B ) revealed that the proportion of pure temporal variation (2.1%) among all samples basically coincided with that of pure spatial factors (1.6%). Remarkably, a large amount of variation (96%) could not be explained by spatial or temporal factors, implying that spatio-temporal factors did not have a significant impact on genetic structure. Although the likelihood of not detecting temporal genetic differentiation could not be excluded, the real degree of temporal genetic differentiation must, nevertheless, had been very small.

Estimates of effective population size using LD method were listed in Table 2 . The Ne estimates for 1960s samples ranged from 120.8 to infinite with lower bounds of confidence interval (CI) ranging from 28.7 to 388.3. Ne estimate for historical Yangtze River Estuary in 1980s was 95.7 with a lower bound of CI of 28.6. For the three contemporary samples, Ne estimates was infinite for Lvsiyang and Wenling, and 176.1 for Zhoushan, with the lower bounds of CI ranging from 34.6 to 82.6. The upper bounds of CI were infinite for all the 7 samples.

Estimates of genetic diversity from both historical and contemporary samples enable us to understand the evolutionary effects of overfishing for overexploited fishery resources. The historical samples in 1960s and 1980s corresponded to large yellow croaker stocks under initial lower fishing stress and population exhaustion after higher fishing stress respectively. Despite sharp decrease in stock size and rapid changes of life-history traits over several decades, no signal of depletion for genetic diversity (Ar, He and F_IS ) was detected for large yellow croaker in the East China Sea. The results suggested that neither genetic drift nor gene flow from allopatric populations had strongly altered the genetic composition of stocks in the East China Sea over this time-span. Consistent with our results, a recent genetic study with three molecular markers also detected relatively high genetic diversity in wild populations of large yellow croaker, which was comparable to those of othermarine fishes (Yuan et al., 2021). Our result was also consistent with previous results of temporal genetic analysis in Atlantic salmon, Brown trout, North Sea sole and Copper redhorse, which also documented temporal stability of genetic diversity over decades of overexploitation (Hansen et al., 2002; Lippe et al., 2006; Cuveliers et al., 2011; Pinsky et al., 2021). For example, genetic erosion were detected neither in the Icelandic nor North Sea stocks of the North Sea sole despite stock crashes in the last century (Hoarau et al., 2005).

Changes of life-history traits together with population declines have been observed in many overexploited fishery stocks (Law, 2000). By contrast, loss of genetic diversity most likely happens when populations are both very small and isolated (Palstra and Ruzzante, 2008). For marine populations, the fishery would likely collapse before the decreased genetic diversity could be detected (Boehlert, 1996). In the classical conservation genetic context, an Ne≥50 is enough to maintain short-term genetic diversity, and an Ne≥500 is sufficient to maintain a longer-term genetic stability (Hoarau et al., 2005). No apparent decrease of Ne was detected for large yellow croaker stocks in the East China Sea over the last 60 years. All of the Ne estimates showed broad confidence intervals with infinite upper boudance of confidance interval. The contemporay effective population sizes were still large enough and genetic drift was probably not strong enough to deplete the genetic diversity. In addition, a mechanism of genetic compensation acting against the negative genetic effects from overexploitation might exist (Palstra and Ruzzante, 2008), which might be interpreted by the reduction of competition for mating partners or spawning sites with lower census population size, leading to more equal reproductive output and a more balanced sex ratio. In general, reproductive success in marine fish of high fecundity strongly depends on larger and older females, which could produce more eggs of good quality (Chambers and Leggett, 1996; Palumbi, 2004; Wright and Trippel, 2009). The removal of those older and larger females by a size-selective fishery might give many other less productive smaller females opportunity of becoming ‘reproductive winners’ (Cuveliers et al., 2011). This might in turn reduce the variance resulting from differences of reproductive success among individuals, so that the population could maintain more genetic diversity. Moreover, the large yellow croaker has high migration ability which facilitates high connectivity among populations (Ma et al., 2016), and gene flow from adjacent populations may be enough to buffer any loss of genetic variation.

Determining the genetic structure is essential for stock improvement and conservation plans (Danish et al., 2021). Results of population genetic analyses indicated that population genetic structure was shallow for large yellow croaker stocks in the East China Sea. The sporadic significant F_ST together with the result of STRUCTURE suggested marginal genetic differentiation among some comparisons between temporal samples. In contrast, the results of PCA indicated that all the temporal populations were clustered together, and no spatio-temporal genetic structure were observed, suggesting a low resolution power of this method for the shallow population structure as revealed by the low F_ST values. AMOVA together with RDA and VPA analysis proved that the proportion of variance explained by temporal factors was small and similar with that of spatial factors. The marginal temporal genetic variation was possibly due to a random genetic drift. It is likely that temporal variations in reproductive success led to random slight allele frequency shifts among generational cohorts and that cumulative genetic variation may be observed after several decades (Han et al., 2010). Especially under the impact of overfishing, the effects of genetic drift might be amplified due to the reduced census population size. However, technical factors could not be excluded. Statistically significant but biologically insignificant differentiation may easily occur at hypervariable markers such as microsatellites (Hansen et al., 2002). Besides, technical errors resulting from mis-scoring or the presence of null alleles remain a concern in historical samples and may lead to biased results. In conclusion, the spatio-temporal genetic structure of the large yellow croaker stocks has essentially been stable. Sporadic genetic differentiation among temporal samples might occur due to random genetic drift or technical errors. To a large extent, these results were in accordance with previous studies on populaiton genetic structure of large yellow croaker in which a single panmictic population were suggested across its distribution (Han et al., 2015; Yuan et al., 2021). Our results indicated that the population genetic structure of large yellow croaker was generally spatio-temporal stable and had not been changed substantially due to the fishing pressure. This had important implications for the potential local adaptation for large yellow croaker and its stock management.

Evolutionary consequences of harvesting must be incorporated into conservation and management strategies (Pandolfi, 2009). Neither significant loss in genetic diversity nor apparent changes in spatio-temporal genetic structure were detected in large yellow croaker stocks in the East China Sea. Genetic diversity plays an important role in phenotypic reversal (Pinsky et al., 2021). The maintanence of genetic diversity in large yellow croaker stocks suggested that potential for stock recovery might still be maintained. So far, no severe threat of losing evolutionary potential due to genetic drift was expected in large yellow croaker populations. However, this optimistic conclusion needed to be treated with caution. Results in the present study based on microsatellite markers mostly reflected neutral variations and the possibility of genetic erosion in adaptive genomic regions could not be excluded. If effective population size of large yellow croaker was indeed large, role of selection could be strong so that the adaptive genomic regions might be more prone to be selected (Cuveliers et al., 2011). Thus, potential selective changes of genetic composition imposed by fisheries might be potent factors making irreversible changes in large yellow croaker populations.

Due to its long generation length (about 15.5 years) and late sexual maturity (about 2-4 years), stock recovery of the large yellow croaker should be difficult after heavy exploitation (Liu et al., 2020). Benefiting from success of artificial mariculture since 1985, stock enhancement of large yellow croaker by release of artificial seedlings was extensively conducted in the East China Sea (Wang et al., 2013). Stock enhancement has been one of the conventional strategies for stock recovery and management of large yellow croaker (Yuan et al., 2021). However, previous studies have shown that the ecological and genetic effects of hatchery releases are of real concern because of the potential in facilitating inbreeding depression and reducing the fitness of wild stocks (Grant et al., 2017), such as Atlantic salmon (Hagen et al., 2020), Salvelinus fontinalis (Marie et al., 2010), and Acanthopagrus schlegelii (Shan et al., 2020). Hatchery-reared large yellow croaker showed low genetic diversity, relaxed purifying selection and increased genetic loads, which might facilitate inbreeding depression and reduce the fitness of wild stocks (Yuan et al., 2021). So the effects of stock enhancement by release of hatchery seedlings should be carefully evaluated in the future. Our results suggested that the genetic diversity of contemporary large yellow croaker wild populations was still in an ideal condition. Considering its potential for self-resilience and the unknown effects of hatchery releases, avoidance of overfishing and habitat restoration should be the best management strategy for the recovery of large yellow croaker. In detail, fishery control measures, including controls of fishing effort and catch, temporal and spatial closures aimed at protecting spawning, nursery and over-wintering grounds and populations, should be regarded as specific actions toward the recovery of large yellow croaker stocks (Liu and Mitcheson, 2008).