Genetic Diversity and Population Structure of the Antarctic Toothfish, Dissostichus mawsoni, Using Mitochondrial and Microsatellite DNA Markers

The Antarctic toothfish, Dissostichus mawsoni, serves as a valuable fishery resource around the Antarctic Continent since 1997, managed by the Commission for the Conservation of Antarctic Marine Living Resources (CCAMLR). Although delineating genetic or stock structure of populations is crucial for improving fishery management of this species, its number of genetic populations and genetic diversity levels remain ambiguous. In the present study, we assessed the population genetic and phylogeographic structure of the Antarctic toothfish across 20 geographic localities spanning from Subareas 88 (88.1, 88.2, and 88.3) to Subareas 58 (58.4 and 58.5) by using mitochondrial DNA (mtDNA) cytochrome oxidase I (COI) and 16S rRNA (16S) sequences and seven nuclear microsatellite loci. MtDNA revealed a low level of polymorphism (h = 0.571, π = 0.0006) with 40 haplotypes in 392 individuals, connected only by 1–5 mutational steps, which is indicative of shallow evolutionary history. Microsatellites showed a range of allelic richness (AR) from 6.328 (88.3 RB3) to 7.274 (88.3 RB6) within populations. Overall genetic diversity was generally higher in Subareas 58 than in Subareas 88, suggesting that effective population size (NE) is larger in Subareas 58. The results of population analyses using microsatellites suggest that the sampled populations are likely to comprise a well-admixed single gene pool (i.e., one genetic stock), perhaps due to high contemporary gene flow occurring during the prolonged larval phase of this fish. However, given weak, but significant microsatellite differentiation found in six population-pairs, the possibility of existence of multiple genetic populations could not be completely excluded. The mtDNA AMOVA suggests a genetic break between the Subareas 88 and 58 groups (FCT = 0.011, P = 0.004). Moreover, mtDNA genetic distances (FST) between populations were proportionally greater as geographic distances increase. The patterns of isolation by distance (IBD) shown only in mtDNA, but not in microsatellites might suggest that population differentiation or divergence processes underwent faster in mtDNA than microsatellites, due to its NE being only one-quarter of nuclear DNA. Temporal stability in the genetic structure of D. mawsoni is also indicated by the results of no genetic differentiation between juveniles and adults. The findings of this study will help to design effective stock management strategies for this valuable fishery resource. We suggest that a long-term genetic monitoring is needed to understand the population structure and dynamics of toothfish in response to ongoing climate changes.

The Antarctic toothfish, Dissostichus mawsoni, serves as a valuable fishery resource around the Antarctic Continent since 1997, managed by the Commission for the Conservation of Antarctic Marine Living Resources (CCAMLR). Although delineating genetic or stock structure of populations is crucial for improving fishery management of this species, its number of genetic populations and genetic diversity levels remain ambiguous. In the present study, we assessed the population genetic and phylogeographic structure of the Antarctic toothfish across 20 geographic localities spanning from Subareas 88 (88.1, 88.2, and 88.3) to Subareas 58 (58.4 and 58.5) by using mitochondrial DNA (mtDNA) cytochrome oxidase I (COI) and 16S rRNA (16S) sequences and seven nuclear microsatellite loci. MtDNA revealed a low level of polymorphism (h = 0.571, π = 0.0006) with 40 haplotypes in 392 individuals, connected only by 1-5 mutational steps, which is indicative of shallow evolutionary history. Microsatellites showed a range of allelic richness (AR) from 6.328 (88.3 RB3) to 7.274 (88.3 RB6) within populations. Overall genetic diversity was generally higher in Subareas 58 than in Subareas 88, suggesting that effective population size (N E ) is larger in Subareas 58. The results of population analyses using microsatellites suggest that the sampled populations are likely to comprise a well-admixed single gene pool (i.e., one genetic stock), perhaps due to high contemporary gene flow occurring during the prolonged larval phase of this fish. However, given weak, but significant microsatellite differentiation found in six population-pairs, the possibility of existence of multiple genetic populations could not be completely excluded. The mtDNA AMOVA suggests a genetic break between the Subareas 88 and 58 groups (F CT = 0.011, P = 0.004). Moreover, mtDNA genetic distances (F ST ) between populations were proportionally greater as geographic distances increase. The patterns of isolation by distance (IBD) shown only in mtDNA, but not in microsatellites might suggest that population differentiation or

INTRODUCTION
The Antarctic toothfish, Dissostichus mawsoni (Perciformes: Notothenidae) is endemic to the Antarctic Ocean with a circumpolar geographic distribution (Hanchet et al., 2015). This species can live up to more than 30 years and reach a total body length of 175 cm and weigh over 80 kg (Yukhov, 1971;Hanchet et al., 2015). While adults of the Antarctic toothfish typically use the entire water column over deep water, juveniles prefer to exploit benthic habitats in shallow coastal water. This species is known to have a pelagic larval stage of 4-5 months, and following this developmental period, its larvae drift for up to 9-12 months with oceanic currents until they reach a body length of 4-8 cm (Hanchet et al., 2008). Species with such a long pelagic larval stage is predicted to undergo extensive gene flow among geographic populations via passive dispersal of planktonic larvae by oceanic currents, likely ensuing spatial homogeneity in the population structure (Bohonak, 1999;Boulding, 2007, 2009;Faurby and Barber, 2012). However, genetic stock structure and population connectivity of D. mawsoni remain still poorly understood.
Dissostichus mawsoni has provided a valuable fishery resource around the Antarctic area since 1997, managed by the Commission for the Conservation of Antarctic Marine Living Resources (CCAMLR). Although delineating the genetic or stock population structure and genetic diversity levels of this species is critical for its fishery stock enhancement and improved management, the number of populations or stocks and diversity levels remain incongruent. A single genetic stock hypothesis with high gene flow via the Antarctic Circumpolar Current is originally proposed for the Antarctic toothfish, particularly in the Weddell Sea (Subareas 48) and Ross Sea (Subarea 88.1; Figure 1; Orsi et al., 1995). A recent study of single nucleotide polymorphisms (SNPs) at whole genome-levels supported the hypothesis of weak population structuring with high connectivity among geographic populations (Maschette et al., 2019). Two or more stock hypothesis has also been proposed, based on a combined analysis of population genetics [using random amplified polymorphic DNA (RAPD) markers], parasitology, stable isotopes and otolith microchemistry (Parker et al., 2002;Hanchet et al., 2015). Whether D. mawsoni is comprised of a single or multiple genetic populations surrounding this region remains therefore unresolved.
Although there are several studies that have been performed to elucidate the population genetic structure of the toothfish using mitochondrial (mt) and nuclear DNA markers (Parker et al., 2002;Smith and Gaffney, 2005;Kuhn and Gaffney, 2008;Mugue et al., 2014), the observed patterns of the genetic structure were not consistent across these studies, and also inconclusive, mainly due to low genetic polymorphism in the molecular markers used. Parker et al. (2002) initially suggested that populations from the Ross and Bellingshausen Seas (Subareas 88) have split into two genetically unique groups, indicating two separate genetic populations existed. A further investigation (Kuhn and Gaffney, 2008) revealed considerable genetic structure among the four CCAMLR subareas [Australian Antarctic Territory (Subarea 58.4.2), Ross Dependency (Subareas 88.1 and 88.2), and the South Shetland Islands (Subarea 48.1); see Figure 1], providing some evidence for the presence of multiple stocks of D. mawsoni, particular given genetic differentiation detected between Subareas 58 and Subareas 88. However, a lack of population structure was observed among the three geographic subareas, including Subareas 48. 1, 88.1, and 58.4.2 (Smith and Gaffney, 2005;Mugue et al., 2014). Such discrepancies in the results of previous studies highlight the need for further study on the population genetics and phylogeography of the Antarctic toothfish. More importantly, future investigation with more polymorphic markers, such as microsatellites or a great amount of SNPs would provide better insights on discovering hidden population structure (Le Cam et al., 2019). Nevertheless, using 10,303 SNP markers Maschette et al. (2019) observed the absence of the spatial population structure among the sampling localities spanning nearly an entire geographic range of D. mawsoni. Unisexually inherited markers such as mitochondrial DNA (mtDNA) are presumed to diverge more rapidly than nuclear DNA considering genetic drift effects alone, as effective population size (N E ) of mtDNA is only onequarter of nuclear DNA (Hartl and Clark, 2007). By using a combined analysis of mtDNA and nuclear microsatellites, we here compared our results with those of Maschette et al. (2019) and discussed possible explanations on the population genetic structure of the tooth fish. Temporal stability in the genetic structure of D. mawsoni was also suggested, by using temporally replicated geographic samples (Mugue et al., 2014). However, age or body size was not considered in their study as a proxy for different temporal groups of 2011-2013 years. There could be differences in the population structure between juveniles and FIGURE 1 | Sampling localities of Dissostichus mawsoni and the geographic distribution of mitochondrial DNA COI and 16S haplotypes (1,304 bp) that we determined. Each polygon in the map represents sampling areas. Different colors in circles denote 40 different haplotypes that show relative frequencies within populations. The map was extracted from CCAMLR area (Everson, 2017). adults, possibly due to low reproductive success of dispersers and/or strong genetic drift effects over generations (Hedgecock, 1994;Lee and Boulding, 2007).
In the present study, we examined the population genetic and phylogeographic structure of the Antarctic toothfish across 20 geographic sites around the seas of the Antarctica Continent (Subareas 88 and Subareas 58; Figure 1) by using mtDNA cytochrome oxidase I (COI) and 16S rRNA (16S) sequences and also seven nuclear microsatellites. Although the previous work (Maschette et al., 2019) with large amount of genome-wide SNPs would be expected to outperform the present study, this study had five specific objectives, which were: (1) to examine the level of genetic diversity among 20 geographically separated populations (including 10 populations from the Subareas 88 and 10 from the Subareas 58); (2) to assess the genetic and phylogeographic structure of these populations for the stock identification; (3) to determine the phylogenetic relationships among the populations (and haplotypes/genotypes) and conduct phylogeographic analyses in order to infer its evolutionary and demographic histories; (4) to test whether there is a difference in genetic structure between juveniles vs. adults; (5) to discuss differences in genetic signals between mtDNA and microsatellites with respect to the population structure. To our knowledge, this is the first study to examine genetic differentiation between juveniles and adults in order to test for the temporal genetic stability in D. mawsoni, which has implications for the population stability or turnover rates. Moreover, a comparison of genetic signals in the population structure between different marker types (e.g., mtDNA vs. nuclear DNA) has not been made. The results of this study will contribute to make conservation and management efforts for this valuable fishery resource around the Antarctic Continent.

Study Sites and Sample Collection
A total of 392 individuals of the Antarctic toothfish were sampled from 20 sites in two CCAMLR Subareas by a bottom longline: Australian Antarctic Territory (Subareas 58.4; nine localities, and 58.5; one locality), and the Ross Dependency (Subareas 88.1; three localities, 88.2; one locality, and 88.3; six localities) around the seas of the Antarctic Continent from December 2016 to March 2020 (Figure 1 and Table 1). Tissue samples from Subareas 88 were obtained during commercial fishing operations in research blocks, granted by the Distant Water Fisheries Resources Division, National Institute of Fisheries Science (NIFS) in South Korea. The samples collected from Subareas 58.4 and 58.5 were provided by a research group of the Australian Antarctic Division, Australian Government, Australia (Maschette et al., 2019). Body size (total length; cm) and weight (g) were measured for every fish. Adult and juvenile samples from 10 populations for comparison of differences in genetic composition were classified based on a standard length (SL) of 81 cm [juveniles < 81 cm (mtDNA; N = 40, microsatellites; N = 42), adults >120 cm (mtDNA; N = 110, microsatellites; N = 114); Supplementary Table 1] following a previous study (Near et al., 2003). Fish individuals with a body size range of 81∼120 cm were excluded to avoid the possible confounding effects. The samples used for genetic analysis were comprised of muscle or dorsal fin tissues for the Subareas 58.
Polymerase chain reaction (PCR) amplification was carried out in a reaction volume of 15 µl comprised of 1 × PCR buffer, 25 µM of each dNTP (Bio Basic Inc., Canada), 0.6 µM of each of the forward and reverse primers, 0.2 U of Taq polymerase (Thermo Fisher Scientific, Waltham, MA, United States) and 5.64∼561.98 ng/µl of DNA template. The following thermal conditions were applied: initial denaturation at 94 • C for 3 min followed by 35 cycles of denaturation at 94 • C for 1 min, annealing at 56 • C (COI) and 60 • C (16S) for 30 s and extension at 72 • C for 45 s, followed by a final extension at 72 • C for 10 min. PCR products were verified on 2% agarose gels stained with RedSafe (iNtRon Biotechnology, South Korea). The amplified PCR products were purified enzymatically by incubating them at 37 • C for 15 min with Exonuclease I (New England BioLabs, United States) followed by 85 • C for 15 min with Shrimp Alkaline Phosphatase (New England BioLabs). The purified mtDNA fragments were subjected to direct sequencing in the forward and reverse directions using the same forward and reverse primers as in the PCR and the BigDye Terminator 3.1 Cycle Sequencing Ready Reaction Kit (Applied Biosystems, United States) in an ABI 3730xl automated DNA sequencer (Applied Biosystems). The DNA sequences were edited using CHROMAS v2.01 computer software and aligned with Clustal W v2.0 (Larkin et al., 2007) implemented in BioEdit v7.2.5 (Hall, 1999), and finally verified manually.

Microsatellite Genotyping
Seven polymorphic nuclear microsatellite loci (DELIG1, DELIG3, DELIG9, DELIG12, DELIG14, DELIG22, and DELIG33), which were originally developed for the Patagonian toothfish, Dissostichus eleginoides (Garcia et al., 2019), were used for genotyping D. mawsoni in this study. Microsatellite genotyping was only successfully undertaken for dorsal fin or muscle tissue samples (N = 348) because PCR of the stomach specimens (from Subareas 58.4) was failed, due to a low DNA concentration. The forward primers were labeled with fluorescent dye (FAM and HEX). PCR reactions were performed as described for amplifying the mtDNA. PCR cycling conditions comprised an initial denaturation phase at 95 • C for 5 min, 35 cycles of 95 • C for 1 min, 57-62 • C for 40 s, 72 • C for 30 s and a final extension step at 72 • C for 10 min. The PCR products were electrophoresed on an ABI 3730xl automated DNA sequencer (Applied Biosystems). Fragment sizes were determined with the ROX 500 bp size standard (ABI) using GENEMAPPER software v 5.0 (Applied Biosystems).

Genetic Diversity
For mtDNA, we used mtDNA COI and 16S as a genetic marker in this study. We selected these genes, as they had relatively high polymorphism across the five regions evaluated [COI, NADH dehydrogenase 2 (ND2), cytochrome b (cyt b), 16S and control region (CR)] in our previous study (Choi et al., 2019). GenBank accession numbers for COI and 16S, and CR, cyt b, and ND2 are MW684720-MW684761 and MZ170964-MZ170973, respectively. In particular, the 16S marker showed the highest degree of diversity [haplotype diversity (h) = 0.468; nucleotide diversity (π) = 0.0009] (Choi et al., 2019). In the present study, a concatenated sequence (1,304 bp) of COI and 16S was used for formal analysis. The number of polymorphic sites and haplotypes, h and π for the concatenated fragment of mtDNA COI and 16S were estimated for each locality, each Subarea 58 and 88 and the entire pooled population of D. mawsoni using ARLEQUIN v3.5 (Excoffier and Lischer, 2010). These diversity indices were also calculated for each of the juvenile and adult groups. The rarefaction method was applied using CONTRIB v1.02 (Petit et al., 1998) to calculate haplotype richness (HR), which was corrected for unequal sample sizes across the populations. A Student's t-test was performed to examine whether there was a significant difference in the level of HR, h, and π between the Subareas 88 (10 populations) and 58 (10 populations). Haplotype network was constructed using HAPSTAR v0.7 (Teacher and Griffiths, 2011) to determine the relationships between the mtDNA haplotypes identified.
To assess microsatellite diversity in D. mawsoni, the mean number of alleles per locus (N A ), observed (H O ) and expected (H E ) heterozygosity, inbreeding coefficient (F IS ) and allelic richness (AR) corrected for unequal sample sizes were calculated for the seven microsatellite loci using GENEPOP v4.0 (Rousset, 2008) and FSTAT v2.9.3.2 (Goudet, 2001). Multilocus tests for Hardy-Weinberg equilibrium (HWE), pairwise F ST estimates (allele frequency-based genetic differentiation) and linkage disequilibrium (LD) tests for genotypes among pairs of the seven loci were also performed in GENEPOP. A Mann-Whitney U test was also conducted to test whether there was a significant difference in the level of AR, H E , and H O between the Subareas 88 (10 populations) and 58 (5 populations).
To test for evidence of a genetic bottleneck in each population of D. mawsoni, we used BOTTLENECK v1.2.02 (Piry et al., 1999) with the two-phase mutation model (TPM) of microsatellite evolution. A population bottleneck can be identified by the occurrence of a mode-shift (i.e., allele distribution shift) and/or a significant heterozygosity excess (Wilcoxon sign rank test) (Luikart et al., 1998). Contemporary effective population sizes (N E ) were also calculated for each of the Antarctic toothfish samples based on the LD method in NeESTIMATOR v2.01 (Do et al., 2014).

Genetic Structure
Both exact tests for population differentiation (Raymond, 1995) and calculation of pairwise estimates of F ST (Weir and Cockerham, 1984) were performed with ARLEQUIN and GENEPOP to estimate geographic differentiation between populations, and also between juvenile and adult groups. These analyses were done for the identification of genetic stock and temporal structure in our samples. The 95% significance levels for the pairwise population comparisons were adjusted using a Bonferroni correction. A hierarchical analysis of molecular variance (AMOVA) was further performed in ARLEQUIN to assess geographic population structure or a genetic break between the Subareas 88 and 58 populations, based on variation in mtDNA COI and 16S and seven microsatellite loci. The spatial AMOVA was undertaken by grouping populations according to two subareas: the Subareas 58 group ( Figure 1). These AMOVA partitioned the total molecular variances between the Subareas 88 and 58, and the Subareas 58, 88.1, and 88.2/88.3 (F CT = "inter-subarea" genetic variation), among populations within the subareas (F SC = "intrasubarea" genetic variation) and within populations regardless of grouping (F ST ). Furthermore, we tested for isolation by distance (IBD) among populations using the Mantel test in GenAlex v6.502 (Peakall and Smouse, 2006). Geographic distance was calculated as the shortest coastline distance between populations (e.g., a central point of each sampling site) in kilometers (Maier et al., 2005;Fontaine et al., 2007).
The microsatellite population structure of D. mawsoni was further evaluated using an individual-based Bayesian population assignment test in STRUCTURE v2.3.1 (Pritchard et al., 2000) under a model of admixed ancestry among populations and correlated allele frequencies. STRUCTURE calculates a likelihood score when the data are forced into a given number of genetic clusters (K) = 1-15. We tested 10 iterations at each K = 1-15, with 100,000 burn-in steps followed by 1,000,000 Markov Chain Monte Carlo (MCMC) generations. STRUCTURE analyses were also carried out separately for the populations in the Subareas 58 and Subareas 88 at K = 1-5 and K = 1-10, respectively. The analysis was further performed to examine if there is a difference in the population structure between juveniles and adults. The K value was estimated from the Structure Harvester 1 website (Earl and vonHoldt, 2012), based on the rate of change in the probability of data between successive K values. Structure analysis ran three times independently to check for the convergence of similar values and found that the three runs arrived at identical K values. The bar plot from STRUCTURE shows individual assignment probability for a given K-value and each vertical bar represents one particular individual. Moreover, three different factorial correspondence analyses (FCA) for Subareas 58, Subareas 88 and all populations were performed using GENETIX v4.05.2 (Belkhir et al., 2004), based on genetic relationships between individuals with multi-locus genotypes.

Demographic History and Phylogenetic Relationships
Population expansion of D. mawsoni was tested using two different analyses including mismatch distribution and Tajima's D and Fu's Fs neutrality tests in ARLEQUIN (Tajima, 1989;Fu, 1997). The Tajima's D and Fu's Fs analyses were performed for the Subareas 58, Subareas 88 and all populations. The demographic history of this species was further determined using Bayesian Skyline Plot (BSP) analysis in BEAST v2.4.5 (Bouckaert et al., 2014). However, we could not resolve its demographic history (BSP data not shown) because our effective sample size (ESS) estimates could not reach a convergence (ESS = 32.1) due to low mtDNA variation, although ESS >200 is suggested as a rule of thumb for Bayesian phylogenetic analysis (Lanfear et al., 2016).
To investigate phylogenetic relationships between the 40 mtDNA haplotypes detected, mtDNA COI and 16S sequences of D. eleginoides (Accession No.: EU074420) and Notothenia coriiceps (Accession No.: MG729470) were used as an outgroup for the phylogenetic analysis. These two species were chosen as an outgroup based on a reconstructed phylogeny for the Antarctic fishes using mtDNA markers (Bargelloni et al., 2000). Phylogenetic tree was reconstructed with a maximum parsimony (MP) method using as implemented in MEGA v6.0 (Tamura et al., 2013) and statistical support was estimated by 1000 bootstrap replicates. Using POPTREE2 (Takezaki et al., 2010), phylogenetic analysis was also carried out with F ST distances among the 15 populations from the Subareas 58 and 88 based on microsatellite allelic variation calculated.

MtDNA
We found only 40 mtDNA COI and 16S haplotypes among 392 samples with one to maximum five mutational steps between haplotypes (Figures 1, 2). Estimates of h and π values for the entire toothfish population were 0.571 and 0.0006, respectively ( Table 1). The most two common haplotypes (Hap1, Hap2) were shared by 235 and 104 individuals (59.9 and 26.5%), respectively and the second most common haplotype Hap3 had only six individuals (1.5%) and Hap4 three individuals (0.8%) of the total samples (Figure 2). Eight haplotypes were observed in two individuals (0.5%) and the remaining 28 haplotypes occurred only once (i.e., singleton). Haplotypes, Hap4, Hap5, Hap7, Hap11, and Hap12 were detected with 2-3 individuals, which were present only in the Subareas 58, and Hap10 was detected with two individuals, which occurred exclusively in the Subarea 88.1 RBH (Figures 1, 2), suggesting that these particular haplotypes might represent locally unique matrilineal lineages.
Allelic frequency distribution showed a normal L-shaped distribution with a probability of 0.688 (every population ranged from 0.297 to 1.000) (Supplementary Figure 2), assuming that all loci fit the TPM based on mutation-drift equilibrium. Therefore, null hypothesis could not be rejected with a Wilcoxon sign rank test, suggesting the study population is unlikely to have experienced a recent genetic bottleneck.  Table 2). However, no significant genetic differentiation was observed in the pairwise F ST estimates using mtDNA COI and 16S between 20 populations from the Subareas 58 and 88 (Table 3). In addition, no significant mtDNA differentiation was observed in the pairwise F ST estimates between adults and juveniles (F ST = −0.014, P = 0.982), albeit weak but significant results found in microsatellites (F ST = 0.001, P < 0.05) (Supplementary Table 2). The Mantel test of mtDNA, however, revealed a significant correlation between genetic (F ST values) and geographic distances (kilometers) for the pairwise comparisons among the 20 populations (r = 0.318, P = 0.010) ( Figure 4A). The observed patterns of a genetic break between the Subareas 58 and 88 and also IBD suggest non-negligible phylogeograhic or population structure over space at least from a mtDNA perspective.

Microsatellites
The AMOVA for microsatellites revealed differences in the genetic structure among populations within the subareas  Table 2), which is opposite to what we found from mtDNA analysis. We also found weak but statistically significant   Table 3). These results might suggest the possibility of genetically different populations of D. mawsoni present in the Subareas 58 and 88.
The results of STRUCTURE analysis of the 15 populations from all populations (combined Subareas 58 and 88) suggested that they are most likely to form three genetic clusters (K = 3) (Figure 5 and Supplementary Figure 3A). Genetic composition appeared to be largely admixed across the 15 populations with homogeneous distribution of individual genotypes, despite they are located far apart at both ends of the Antarctica (Subareas 58 and 88), suggesting that all these populations are genetically indistinguishable as one group in general, albeit some local structure also present. Approximately similar proportions of the inferred three genetic clusters suggest that all populations are genetically undifferentiated and so could not form clusters most likely due to the very low level of F ST values (overall mean F ST = 0.0053) (Latch et al., 2006). When STRUCTURE analysis was performed separately for each Subarea 58 and 88, the observed most likely number of genetic clusters corresponded to two genetic clusters for the Subareas 58 (Supplementary Figures 3B, 4A), and three clusters for the Subareas 88 (Supplementary Figures 3C, 4B). Moreover, FCA results of seven microsatellites supported the well-admixed panmictic population constituting the Subareas 58 and Subareas 88 (Figure 6). When FCA analysis was performed separately for each Subarea 58 and 88, similarly results of no genetic clustering were found (Supplementary Figure 5). The Mantel test of microsatellites showed a non-significant correlation between genetic (F ST values) and geographic distances for the 15 populations (r = 0.010, P = 0.510) (Figure 4B). The results of STRUCTURE analysis of adult-juvenile comparisons showed that they are most likely to form a single genetic group (K = 1) (Supplementary Figure 6), suggesting no genetic difference between juveniles and adults.

Demographic History and Phylogenetic Relationships
A sudden population expansion of D. mawsoni was supported by the results of the mismatch distribution and neutrality tests. Mismatch distribution of D. mawsoni showed smooth and unimodal shapes albeit statistical significance (sum of squared deviations = 0.1148, P = 0.01, raggedness index = 0.2884, FIGURE 3 | Comparisons of the level of mtDNA COI and 16S (A), seven microsatellite (B) diversity between the Subareas 58 (mtDNA: 10 populations, microsatellites: 5) and 88 (mtDNA: 10, microsatellites: 10). Levels of mtDNA diversity (HR, h, and π) were significantly higher for the Subareas 58 than for the Subareas 88. Levels of microsatellite diversity (AR, H E , and H o ) tended to be higher for the Subareas 58, albeit statistical non-significance. **: P < 0.01, ***: P < 0.001. P < 0.01; Figure 7), suggesting that sudden population expansion had occurred (Rogers and Harpending, 1992 With respect to phylogenetic relationships, forty mtDNA haplotypes linked by only 1-5 mutational steps could not provide any evidence for diverged clades or lineages among haplotypes and also across localities (Supplementary Figure 7A), suggestive of their shallow evolutionary history. Similarly, NJ tree based on microsatellite distances did not show phylogenetic patterns of distinct clades locally evolved in the subareas (Supplementary Figure 7B).

Genetic Diversity
Our study shows that an overall mitochondrial diversity in D. mawsoni is quite lower than the reported levels of mtDNA diversity in temperate or tropical fishes (Manel et al., 2020). The mtDNA COI nucleotide diversity levels (sequences longer than 500 bp) as estimated from the mean number of mutations per base pair (π) are found to be <0.001∼0.002 in temperate fishes, which is substantially lower than tropical fishes (π > 0.002) (Manel et al., 2020). Our estimated π value for D. mawsoni in    COI (685 bp) is 0.00016 ± 0.00029, which is even approximately more than ten-fold lower than those levels. There are two hypotheses that have been proposed to explain the observed low mtDNA gene diversity in D. mawsoni than those of temperate and tropical fishes. Small N E or historical bottleneck events may lead to the low mtDNA diversity in this species (Parker et al., 2002). Alternative explanation would be that extremely cold environments drove low rates of micro-evolutionary change owing to much lower metabolic rates, resulting in the low mtDNA diversity in the Antarctic toothfish relative to temperatezone fish (Parker et al., 2002). The first hypothesis could be rejected by our findings that none of the 20 populations from the Subareas 58 and 88 tested show a sign of population bottlenecks. However, star-like tree topology of the haplotype network, where the most common, two ancestral haplotypes are linked with more recent ones connected with only one or two mutations, typically represent the scenario of population expansion from a small number of founders around the Antarctic Ocean once the glacial maximum had passed and the ice cover sheets retreated (Figure 2). Marine organisms with high dispersal potential in the Antarctic Ocean have often shown the "star-like" pattern of haplotype network (Allcock and Strugnell, 2012). Results of Tajima's D and Fu's Fs statistics and mismatch distribution of D. mawsoni could also reject the hypothesis of the constant population size under mutation-drift equilibrium (Figure 7). A demographic expansion during the last glacial maximum (LGM) was documented for benthic fishes, Trematomus bernacchii and Trematomus pennelli (Janko et al., 2007). Nevertheless, the second "evolutionary speed" hypothesis is more likely to account for low mtDNA diversity, given a strong positive correlation between the rate of metabolism and DNA substitution rate for both nuclear and mtDNA in vertebrates (Avise et al., 1992;Bargelloni et al., 1994). Metabolic rates of Antarctic fish are known to be typically lower than those of temperate fish (Bargelloni et al., 1994). Such an association may be caused by oxidative damage by radicals, which are yielded as by-products of aerobic metabolism (Richter et al., 1988). A recent review of mitochondrial diversity in 3,815 marine fish species at the global scale showed that sea surface temperature is the major determinant of the global patterns of diversity FIGURE 7 | Results of mismatch distribution of nucleotide differences between pairs of mtDNA COI and 16S haplotypes (1,304 bp) in Dissostichus mawsoni. There was a significant difference between the observed and expected frequencies (sum of squared deviation = 0.1148, P = 0.01). (Manel et al., 2020). The causal link between temperature and diversity levels can be explained by the evolutionary speed hypothesis that cold temperature environments like the Antarctic region could lengthen generation times and slow-up mutation rates, thereby likely diminishing genetic diversity (Oppold et al., 2016). This interesting hypothesis of "the low metabolic rate driving slow molecular evolution in mtDNA" needs further a detailed investigation.
At the intraspecific levels, by comparison to a previous study of the Antarctic toothfish (Kuhn and Gaffney, 2008) using different mtDNA markers (ND2: h = 0.182, π = 0.0003; CR/tRNA: h = 0.032, π = 0.0001; cyt b: h = 0.153, π = 0.0003), our study shows that 20 geographic populations inhabiting the two main Subareas 58 and 88 in the Antarctic Ocean have relatively higher levels of mtDNA diversity (h = 0.571, π = 0.0006), although much larger sample sizes used in this study. The results of our study can therefore provide more conclusive information for delineating the population and phylogeographic structure, compared to previous studies using other mtDNA markers (Smith and Gaffney, 2005;Kuhn and Gaffney, 2008). The Patagonian toothfish (D. eleginoides), a sister species of our study organism, also showed a very low mtDNA diversity (h = 0.1193; π = 0.0007) in the 12S rRNA region (249 bp) for 151 individuals (Rogers et al., 2006). Among 40 mtDNA haplotypes detected in the present study, 21 haplotypes occurred exclusively at Subareas 58, whereas 13 haplotypes existed only at Subareas 88, although the remaining six haplotypes (Hap1,  Hap2, Hap3, Hap6, Hap8, and Hap9) were shared between the two Subareas (Figures 1, 2). In particular, Hap4, Hap5, Hap7, Hap11, and Hap12 with 2-3 individuals, have evolved only in the Subareas 58 and Hap10 with two individuals has diverged only in the Subarea 88.1 RBH. The presence of locally and/or regionally unique matrilineal lineages might provide some evidence supporting the previous hypothesis that D. mawsoni populations form two genetically or phylogenetically distinct clusters or lineages (Ross and Bellingshausen Seas) based on RAPD markers (Parker et al., 2002). SNPs in the COI and 16S could distinguish distinct lineages between the Subareas 88 and 58, and also among the Subareas 58, 88.1, and 88.2/88.3, given significant differences detected in the genetic structure between/among these regions suggested by mtDNA AMOVA analyses, albeit no such a pattern shown in microsatellites ( Table 2). Relatively higher mtDNA diversity in the Subareas 58 compared to Subareas 88 might in part contribute to the observed mtDNA population structure. The higher levels of genetic diversities can be interpreted as larger N E for the populations from Subareas 58 than for those from Subareas 88, given a positive correlation expected between population sizes and genetic diversity level according to population genetics theory (Freeland et al., 2011). Mean of a lower bound of N E is higher in the Subareas 58 (N E = 125.1) than the Subareas 88 (N E = 90.9) after excluding 88.3 RB6 with a lower bound of ∞ (Table 1). However, a mtDNA footprint of the phylogeographic structure (e.g., a genetic break from AMOVA, significant IBD) might have been erased by strong contemporary connectivity via ongoing gene flow occurring among the populations during their pelagic larval stage, as shown by weak population differentiation, absence of IBD, and large population admixture in the results of microsatellite analyses. The faster genetic differentiation or divergence processes undergone in mtDNA relative to nuclear DNA due to its one-quarter smaller N E (Avise, 2012) would be the most plausible explanation of how and why some detectable structure could only exist in mtDNA (see below). Although mtDNA genes have frequently been used as genetic markers for inferring evolutionary or demographic histories of organisms (i.e., phylogeography), microsatellite loci are suggested to be more suitable for elucidating "contemporary" processes of populations as they typically evolve faster (Avise, 2012).

Population Structure
In general, the results of our population structure analysis based on microsatellites support the notion that the sampled populations of D. mawsoni are likely to represent a genetically well-admixed panmixia, suggested by multiple lines of analyses such as pairwise F-statistics and STRUCTURE. Previous studies of mtDNA and nuclear intron markers also showed a negligible level of population structure over space and time in the Antarctic toothfish populations from the CCAMLR Subareas (Smith and Gaffney, 2005;Mugue et al., 2014). A recent population genomics study of D. mawsoni from the Subareas 48, 58, 88, and the South Pacific Regional Fisheries Management Organisation (SPRFMO) area north of Subarea 88.1 suggested a complete lack of the genetic structuring across the Southern Ocean (Maschette et al., 2019). The weak stock structure is most likely to result from high population connectivity through extensive larval dispersal by the Southern Ocean circulation systems during the free-swimming larval stage of D. mawsoni. The epipelagic period as planktonic larvae of D. mawsoni has been reported as up to 9-12 months (Hanchet et al., 2008). The findings of genetic homogeneity are therefore due to high ongoing gene flow across geographic localities with oceanic currents (e.g., Antarctic Circumpolar Current, Antarctic Coastal Current) during the planktonic larval stage of this fish. Long-distance dispersal of adults is also plausible, given adults can migrate up to over 4,000 km, based on mark-recapture tagging studies (CCAMLR-Secretariat, 2017).
However, we should be precautious about the interpretation of the single genetic stock (i.e., gene pool) of these populations, considering detectable microsatellite differentiation apparent for six population-pairs, mtDNA spatial structure present between the Subareas 58 and 88 and also among the Subareas 58, 88.1, and 88.2/88.3, and also the significant IBD patterns in mtDNA (see below). Some studies revealed significant geographic variation among populations, suggesting the possible existence of multiple genetic stocks in D. mawsoni in this region (Parker et al., 2002;Kuhn and Gaffney, 2008). Despite geographic isolation among the 15 populations in the Subareas 58 and 88, we found there was weak but statistically significant genetic differentiation only in six population-pairs, such as 58.4. Although there is generally no evidence for mtDNA differentiation among populations between the two Subareas 88 and 58 as indicated by F-statistics, it should be cautious about managing them as one stock (i.e., a panmictic population), given some microsatellite differentiation and also locally specific mtDNA lineages (Subareas 58: Hap4, Hap5, Hap7, Hap11, and Hap12; Subareas 88: Hap10) detected. The relatively higher betweenand among-subarea variation, with reference to within-subarea variation suggested by multiple AMOVAs also indicates weak, but some detectable mtDNA variation exists between the Subareas (Table 2). Furthermore, a Mantel test of mtDNA suggests moderate evidence of positive relationship in genetic distances (F ST ) between populations and geographic distances, which is consistent with a recent genomic study of Maschette et al. (2019). These results suggest that while genetic differences between populations are rather small, they tend to accumulate with increasing geographic distances, supporting the nonnegligible population structure. The observed patterns of IBD only present in mtDNA, but not in microsatellites might suggest the detectable phylogeographic structure shaped by evolutionary history of D. mawsoni, while this evolutionary signature has been wiped out by the strong contemporary processes (connectivity) among the populations via ongoing gene flow during its pelagic larval stage. This pattern was also indicated in the Patagonian toothfish (D. eleginoides) (Shaw et al., 2004). Population differentiation or divergence processes are predicted to undergo more quickly, in theory, four-times faster in mtDNA relative to nuclear DNA due to its N E being one-quarter of the nuclear genome (Moore, 1995). This could explain why we observed the significant IBD only in mtDNA, but not in microsatellites. Another possibility would simply be that the microsatellite loci used in this study have a poor resolution for discriminating geographic populations of the toothfish. No significant genetic differentiation between juveniles and adults suggests the temporal stability in the genetic structure of D. mawsoni, which is consistent with a previous finding of Mugue et al. (2014). No significant temporal genetic variation indicates that their N E values are large enough to prevent genetic drift effects over generations (Lee and Boulding, 2007). Our estimated N E values with a range of 1,514∼∞ may support this notion.
The population/stock structure of Antarctic toothfish provides important information for improving stock management/enhancement programs for at least two particular reasons: (1) it gives detailed information on the number and distribution of the genetic stocks of the species; (2) it offers information on the comparative genetic diversity levels. Based on our findings, we can provide two suggestions for the management and conservation of this fishery resource. The Antarctic toothfish fishery is now managed by CCAMLR at the Subarea and Division levels coordinated by different nations. Its primary aim is to make decision rules to set catch limits or "quota" at the management levels. In this regard, we suggest directions to "reduce the catch (harvest) in Subareas 88" in comparison to Subareas 58, given their significantly lower diversity. In theory, low genetic diversity within populations implies that N E is small and is thereby susceptible to the effects of inbreeding, genetic drift or a combination of the two (Hartl and Clark, 2007). While we observed only weak population structure of D. mawsoni across the Subareas 58 and 88 regions as suggested by Maschette et al. (2019), we cannot argue that the current fishery management policy needs to be changed as we are still uncertain that stock structure between the Subareas 58 and 88 may exist. Genetic analysis alone has, of course, a limitation in determining the stock structure, as information on ecological differentiation would be necessary. Assessments of ecological divergence between populations from Subareas 58 and Subareas 88 could compromise the limitations of genetic analysis alone in identifying the population/stock structure of D. mawsoni. For example, evaluating trophic niche differentiation between the regional populations through the analysis of their food compositions will help to understand the stock structure of this species (see below).
Given our results are inconclusive, adding more microsatellites (particularly D. mawsoni specific primers) to the existing dataset might help to more accurately determine the population or stock structure of D. mawsoni in the seas of the Antarctic Continent. Although the study sites and genetic approaches are different, our findings of a signal of some genetic differentiation or structure are inconsistent with the previous results of a complete lack of structure using 10,303 SNPs that were developed across the whole nuclear genome (less than 0.1% of the genetic variation) (Maschette et al., 2019). Based on the findings of the present study, it should be careful to manage fisheries for the genetically identical stock, although overall patterns suggest they are likely to share a common gene pool. Nevertheless, a mantel test of mtDNA for IBD indicates that the degree of genetic differentiation was positively correlated with geographic distances. Ecological features of the Antarctic toothfish, such as dietary composition and feeding strategy (Yoon et al., 2017), possible locations of its spawning grounds and effects of the Antarctic Ocean currents should be considered more carefully to better understand how the population/stock is shaped and structured. Different dietary composition between Subareas 58 and Subareas 88 represents ecological divergence between the two regional populations, which may facilitate the population structuring (Yoon et al., 2017). Using additional tissue samples from the Subareas 48 and also from the Subareas 88 will allow for defining the population/stock structure of this fish more precisely. A continuous long-term genetic monitoring will be crucial for understanding how well these valuable fishery populations will sustain particularly under current climate changes.

DATA AVAILABILITY STATEMENT
The mtDNA COI and 16S, and CR, Cyt b and ND2 sequences obtained for this study have been deposited in GenBank under the accession numbers MW684720-MW684761 and MZ170964-MZ170973, respectively. The datasets (distribution of mtDNA haplotypes among the 20 sampling localities, microsatellite genotypes) are also included in Supplementary Tables 3, 4, respectively.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because Antarctic toothfish individuals were sampled from around the seas of the Antarctic Continent by a bottom longline for fishery and we used tissues of these animals that were already dead for our genetic analysis.

AUTHOR CONTRIBUTIONS
SC, S-GC, H-WK, and HJL conceived and designed the study. H-kC, JEJ, SYB, YRK, DM, and HJL performed the experiments and analyzed the data. SC, S-GC, and DM supported the field expedition for Antarctic toothfish sampling. H-kC, JEJ, SYB, YRK, and HJL drafted the manuscript. All authors read and approved the final manuscript.

FUNDING
This work was supported by a grant from the National Institute of Fisheries Science, South Korea (NIFS; R2019021, R2020023, and R2021029).

ACKNOWLEDGMENTS
We thank members of the National Institute of Fisheries Science for organizing to collect Antarctic toothfish samples in the field. We also thank Soo Rin Lee and Gun Wook Baeck for helping to obtain tissues or stomach samples for DNA analysis.