Differentiation of the frog sculpin Myoxocephalus stelleri Tilesius, 1811 (Actinopterygii, Cottidae) based on mtDNA and karyotype analyses

Abstract A molecular genetic and karyological study of the frog sculpin Myoxocephalus stelleri Tilesius, 1811 was carried out on an extensive sample from a large area of the species’ range. A total of 42 specimens was sampled from the Sea of Japan, Sea of Okhotsk, and coastal waters off the southern Kuril Islands, which makes this sampling scheme the most comprehensive to date. The level of mtDNA polymorphism was found to be low. The haplotypes of the species formed three phylogenetic groups. The unique M. stelleri haplotype from the coast of Shikotan Island linked all the studied groups, indicating that it is likely ancestral. Robertsonian polymorphism was identified in the species. In all five cytotypes (I – 2n = 44, II – 2n = 43, III – 2n = 42, IV – 2n = 41, V – 2n = 40; NF = 44+2) were identified, all of which were present in the Sea of Japan. Only one (cytotype I) was found in the Sea of Okhotsk, which is probably the closest to the ancestral karyotype. The significant chromosomal polymorphism and the presence of common haplotypes in the studied samples indicate their recent origin from a common ancestor and/or relatively recent contacts within the range. The discrepancies between mtDNA and karyotypes in assigning the ancestral M. stelleri to the coastal waters off Shikotan Island (southern Kuril Islands) and the Sea of Okhotsk, respectively, can be explained by the different inheritance mechanisms and the rates of evolution of molecular genetic and karyological traits.


Introduction
The genus Myoxocephalus Tilesius, 1811 is a large, taxonomically complex group of sculpins of the subfamily Myoxocephalinae (family Cottidae) (Neelov 1979). The modern world catalog (Fricke et al. 2020) lists 14 species and subspecies of this genus, which inhabit coastal and shelf zones of the Atlantic, Arctic, and Pacific Oceans (Neelov 1979;Parin et al. 2014;Mecklenburg et al. 2016). The northern Pacific Ocean is considered the center of diversification of this genus (Nazarkin 2000). The frog sculpin Myoxocephalus stelleri Tilesius, 1811 is one of the most common species of the genus, found across a wide geographic range from the Sea of Japan to the Gulf of Alaska (Mecklenburg et al. 2002). M. stelleri is characterized by significant morphological and ecological plasticity, as it tolerates substantial fluctuations in water temperature and salinity, being found even in river estuaries (Neelov 1979;Sokolovsky et al. 2011). Similar species, including M. decastrensis (Kner, 1865) and M. raninus Jordan et Starks, 1904(Parin et al. 2014Fricke et al. 2020), which are now considered to be junior synonyms of M. stelleri, have been described from various parts of the M. stelleri range.
Cytogenetic and genetic studies of M. stelleri have been described in two publications (Miller 2000;Podlesnykh and Moreva 2014). The only molecular genetic study thus far was based on a small number of specimens from Amur Bay, Sea of Japan, and Odyan Bay, Sea of Okhotsk (Podlesnykh and Moreva 2014). The study highlighted the limitations of using the short fragment of the CO1 mitochondrial gene (525 bp) when constructing an adequate system of the Myoxocephalus species due to the lack of informative characters. The frog sculpin karyotype was first described from Amur Bay, Sea of Japan (Miller 2000). Unlike other Myoxocephalus species from the Sea of Japan and Sea of Okhotsk (Radchenko et al. 2020), M. stelleri exhibits chromosomal variation. Karyological analysis of morphologically indistinguishable individuals from Odyan Bay and Amur Bay revealed their differences not only in the diploid numbers (2n) (2n = 44 vs. 2n = 40; fundamental numbers (NF) 44 + 2), but also in the number and localization of active nucleolar organizers (NORs) (Podlesnykh and Moreva 2014). This led the authors to suggest revising the species affiliation of M. stelleri from the Sea of Japan (Podlesnykh and Moreva 2014). However, the same study showed that the M. stelleri clade is a combination of haplotypes of specimens not only from different geographic localities, but also with different karyotypes. Hence, the discrepancy between the genetic and cytogenetic results, along with the small geographic area of sampling and the limited resolution of the short CO1 fragment, warrant further studies.
We conducted a study to determine the level of genetic and karyological differentiation within and between M. stelleri from the Sea of Japan, Sea of Okhotsk, and the southern Kuril Islands. Using this karyological (N = 42) and molecular genetic data (N = 34), we also aimed to find the centers of species diversification. Our extensive sample included individuals of M. stelleri captured from waters near the site of the original species description -the estuary of the Bolshaya Zapadnaya Kamchatka River (Tilesius 1811) -and also from Chikhachev Bay (De Kastri), which is the type locality of M. decastrensis (Fricke et al. 2020). The results obtained allow assessment of the geographical variation of M. stelleri and identification of its probable causes.  Table 1 show the sampling sites and number of specimens examined. Fish were collected from 2016 to 2019, at 20 localities in coastal waters of the Sea of Japan (7 localities), Sea of Okhotsk (10 localities), and off the southern Kuril Islands (3 localities). Species were identified by morphological traits (Neelov 1979). The Sea of Japan is divided into the northern and western parts according to Appendix to "Nekton of the Northwestern part of Japan (East) Sea" (Shuntov et al. 2004). The fish specimens are stored at the Ichthyological laboratory collection of the Institute of Biological Problems of the North, Far Eastern Branch, Russian Academy of Sciences, Russia (voucher numbers are listed in Table 1).

Ethics statements
This study utilizes samples collected with all applicable international, national and/ or institutional guidelines for sampling, care and experimental use of organisms. The fishes studied here are not included in the IUCN Red List of Threatened Species, nor are they are listed as endangered, vulnerable, rare, or protected species in the Russian Federation. The sampling points are located beyond any protected areas.

DNA analysis
We obtained sequences for three mtDNA markers: COI, cytochrome b, and 16S rRNA. Total DNA was extracted from muscle tissues by standard phenol extraction (Maniatis et al. 1982) following tissue lysing with 1% SDS using Proteinase K (0.2 mg/mL). The following oligonucleotide primers were used to amplify and sequence the DNA markers: for COI, F-33 TCACAAAGACATTGGCACCCTA and R-1421 TTCACGTT-TAGCAGCGAATGCTT ; for cytochrome b, L14795 CAATG-GCAAGCCTACGAAAA and H15844 AGCTACTAGTGCATGACCATC (Radchenko 2005); for 16S rRNA, L2510 CGCCTGTTTATCAAAAACAT and H3080 CCGGTCTGAACTCAGATCACGT (Meyer 1993).  Table 1. DNA sequences were aligned in Clustal W (MEGA version X: Kumar et al. 2018) with default settings, and manually edited after visual inspection. To identify haplotypes and their relationships, a median network was built in SplitsTree4 v4.12.3 (Huson and Bryant 2006) using the median-joining method (Bandelt et al. 1999). Genetic distances (d) between haplotypes were calculated using p-distances in MEGA X. DNA sequences were first analyzed for each gene independently, then a concatenated matrix was created in which different sets of partition scenarios were investigated. For each gene, we used MEGA X under Bayesian Information Criterion (BIC) and Akaike Information Criterion (AIC) to select the optimal models of nucleotide substitutions. Bayesian Inference trees were constructed using MrBayes v3.2.1 (Ronquist et al. 2012), with the prior set to fit the evolutionary models suggested by MEGA X but allowing the parameters to be recalculated during the run. The Markov Chain Monte Carlo process was set for four chains to be run simultaneously for 1,000,000 generations, with trees sampled every 100 generations. Out of the total 10,001 obtained trees, the first 1,001 with unstable parameters of the models of nucleotide substitutions were rejected. The Bayesian analysis dynamics were controlled in Tracer v1.5 (http://tree.bio.ed.ac.uk/software/ tracer/). Nodes with posterior probabilities ≥0.95 were accepted as statistically significant (Leache and Reeder 2002). We used the following four congeneric species of the subfamily Myoxocephalinae as outgroups: Myoxocephalus jaok (Cuvier, 1829), M. brandtii (Steindachner, 1867), M. polyacanthocephalus (Pallas, 1814), and M. ochotensis Schmidt, 1929.

Karyological analysis
Chromosomes were prepared by the air-drying technique (Kligerman and Bloom 1977). Slides were stained with a 4% azure eosin solution (Giemsa, Merck, Germany) in 0,067 M Sorenson's buffer contains 4.53 g KH 2 PO 4 and 4.72 g Na 2 HPO 4 per liter of distilled water. Metaphase plates were examined under a Leica microscope. The best metaphase plates were photographed using an AxioCam HR CCD camera with the AXIOVISION software (Carl Zeiss MicroImaging GmbH, Germany). Microscopy and imaging were carried out at the Far Eastern Center of Electron and Light Microscopy, A.V. Zhirmunsky National Scientific Center of Marine Biology (NSCMB), FEB RAS, Vladivostok.
Chromosomes were classified according to the nomenclature of Levan et al. (1964). Metacentric chromosomes with equal arms and submetacentric chromosomes with unequal arms were referred to as bi-armed chromosomes; subtelocentric chromosomes with a very short second arm and acrocentric chromosomes with an invisible second arm were referred to as uni-armed chromosomes. Submeta-subtelocentrics were also distinguished in the cytotypes. The chromosomes of this pair are morphologically highly variable. In different cells, they may look like submetacentrics or subtelocentrics. For this reason, the number of chromosome arms (NF) is indicated with both possible variants: 44 + 2.

DNA analysis
For M. stelleri, the length of the partial COI was 1,009 base pairs (bp), including 16 variable sites, 8 parsimony informative sites, and 5 non-synonymous substitutions. The length of the partial cytochrome b was 747 bp, including 16 variable sites, 10 parsimony informative sites, and 1 non-synonymous substitution. The length of the partial 16S rRNA was 600 bp, including 1 variable site and 1 parsimony informative site. All sequence data are deposited in GenBank/NCBI (www.ncbi.nlm.nih.gov) (for accession numbers, see Table 1). Based on the parsimony informative sites, 13 haplotypes were found (Table 2).
The haplotype network for M. stelleri is a star-shaped structure with the central haplotype (1b) from Shikotan Island (Fig. 2). There were three haplogroups, each formed largely by the haplotypes from the same geographic area. Group NSJ (the Northern Sea of Japan group) includes three haplotypes, of which 2a is the most common. Haplotypes of this group are found in the northern Sea of Japan and in the coastal waters off the southern Kuril Islands. Group SO (the Sea of Okhotsk group) has six haplotypes, of which 3b is the most common. This group includes all specimens from the Sea of Okhotsk, as well as individuals from the northern Sea of Japan and the southern Kuril Islands. Group WSJ (the Western Sea of Japan group) has four haplotypes and includes all specimens from the western Sea of Japan and the northern part of the Sea of Japan.  The Bayesian tree (Fig. 3) shows three major clades that are congruent with the three haplogroups described above (NSJ, SO, and WSJ). Group WSJ is basal, while groups NSJ and SO are sister clades. The main nodes are supported by significant posterior probability values (≥0.95).

Discussion
Our DNA data show that M. stelleri is genetically diverse at various levels: within the same locality, between geographically close localities, and between geographically distant areas. Each sample had 2 to 4 haplotypes, some of which were unique and others widespread; however, there were no haplotypes common to all studied samples. Several samples had specific haplotypes, e.g. haplotypes 4b (frequency 14.7%) and 6b were found only in the Sea of Okhotsk, while haplotype 4c (frequency 5.9%) was found only in Peter the Great Bay. Similarly, the unique haplotypes 1a and 1b were found in M. stelleri from the coastal waters off the southern Kuril Islands, and haplotypes 1c and 5b were found only in the Zolotaya Bay sample.
The most common haplotype (2a; frequency of 20.6%) was found in M. stelleri from the northern Sea of Japan (Chikhachev Bay, Zolotaya Bay, Dzhigit Bay, and Olga Bay) and from the coastal waters off the southern Kuril Islands. Haplotype 3b (17.7%) was distributed wider across the species range, from the southern Kuril Islands and the northern Sea of Japan to the Sea of Okhotsk. The less common haplotypes (3a, 2b, 2c, and 3c) were also found in more than one sample. The common haplotypes from the Sea of Okhotsk, Sea of Japan, and the southern Kuril Islands can be explained either by their origin from a common ancestor followed by dispersal from the same center, or by recent contact in various parts of the species range. A pattern of haplotype distribution with a few haplotypes being very common and others being rare or unique is frequently found in marine fish (Avise 2000).
In general, M. stelleri from the northern Sea of Japan exhibits higher genetic variation: eight haplotypes, many of which are shared with other geographic areas, and a low frequency of unique haplotypes. The variation between mtDNA sequences of M. stelleri from different localities is low, with the genetic distance between samples being approximately at the same level (Table 4). An exception is the southernmost sample from Peter the Great Bay, which is the most divergent (d = 0.34-0.52%).
The mtDNA haplotypes form three haplogroups, congruent with the three clades formed in the phylogenetic tree. These haplotypes belong to different phylogenetic groups: NSJ, SO, and WSJ. There were more differences within the Sea of Japan (groups NSJ and WSJ: d = 0.47%) than between the Sea of Okhotsk and the Sea of Japan (SO and NSJ: d = 0.36%; SO and WSJ: d = 0.42%). The geographic distribution of haplotypes is not uniform (Fig. 1). M. stelleri from the Sea of Okhotsk have only SO haplotypes, while haplotypes from all phylogenetic groups were found in the Sea of Japan. For example, in the northern part of the sea, NSJ haplotypes were found at a frequency of 50%, SO haplotypes occurred at a frequency of 31%, and WSJ haplotypes at a frequency of 19%. However, in the western part, only WSJ haplotypes were found. In the Pacific coastal waters off the southern Kuril Islands, NSJ and SO haplotypes were found with equal frequencies.
Two unique haplotypes, 1a and 1b, were found in the coastal waters off Shikotan Island. In the Bayesian tree (Fig. 3), these haplotypes have a basal position in clades NSJ and SO; in the haplotype network (Fig. 2), haplotype 1b has a central position connecting all haplogroups, indicating that it is an ancestral haplotype. Such position of the Shikotan (Pacific) haplotypes probably shows the ancestral role that this population played in the formation of the M. stelleri species pattern, as evidenced by the starshaped network with the central "founder" haplotype. This structure is typical when a population historically experienced an exponential increase in numbers after an earlier reduction in effective population size (bottleneck event) (Avise 2000).
The variation between haplotypes from different parts of the geographic range is most clear between the most distant localities: Shikotan Island vs. the western Sea of Japan, or the Sea of Okhotsk vs. the western Sea of Japan. Similar DNA differentiation has been reported for Cottidae species found in the Sea of Okhotsk and the Sea of Japan, documented both at the subspecies level, e.g. Megalocottus platycephalus platycephalus (Pallas, 1814) vs. M. platycephalus taeniopterus (Kner, 1868) (Radchenko and Petrovskaya 2019), and at the species level, e.g. species of the genera Enophrys Swainson, 1839 and Porocottus Gill, 1859 (Moreva et al. , 2019. This trend has also been observed in other fish families: Lycodes matsubarai Toyoshima, 1985 (Zoarcidae; Sakuma et al. 2014), Bothrocara hollandi (Jordan et Hubbs, 1925) (Zoarcidae; Kodama et al. 2008), and the southern and northern forms of Tribolodon hakonensis (Günther, 1877) (Cyprinidae; Ryazanova and Polyakova 2012). In our case, the mtDNA differentiation of M. stelleri from the Sea of Okhotsk and the Sea of Japan is not higher than the intraspecific level (d = 0.21-0.47%).
The results of the karyological analysis are consistent with the conclusion drawn from the genetic data: M. stelleri is a heterogeneous species (Figs 1, 3). We have shown that M. stelleri has a higher chromosomal polymorphism than was previously assumed (Podlesnykh and Moreva 2014). The highest variation was found among individuals sampled from the northern Sea of Japan: 2n = 44 (cytotype I), 2n = 43 (cytotype II), 2n = 42 (cytotype III), and 2n = 41 (cytotype IV). In the western Sea of Japan, the coastal waters off the southern Kuril Islands (cytotype V; 2n = 40), and in the Sea of Okhotsk (cytotype I), individuals are stable in terms of diploid number (Table 3; Figs 1, 4).
All cytotypes of M. stelleri differed in the 2n (Fig. 4). Unlike other cytotypes, cytotype I lacks metacentrics and large submetacentrics. Cytotypes II-V contain different numbers of bi-armed chromosomes. All cytotypes differ in the number of subtelocentrics. Despite these differences, there are also traits common to all chromosome sets (Fig. 4). These similarities, as well as NF (44 + 2), suggest that the metacentrics and submetacentrics in cytotypes II-V were formed through Robertsonian fusions. Karyological studies show that Robertsonian translocations are the main mechanism of structural changes in the chromosomes of marine sculpins of the family Cottidae (Vasilyev 1985;Terashima and Ida 1991;Moreva et al. , 2019Radchenko et al. 2020). Here we report the highest intraspecific variability in 2n (40 to 44) documented to date for this fish group (Arai 2011;Moreva et al. , 2019Moreva 2020).
The differences in the number of subtelocentrics between the chromosome sets ( Fig. 4) indicate that these chromosomes have been involved in the formation of metacentrics and submetacentrics in cytotypes II-V. Cytotypes II-IV lack a pair of large subtelocentrics, which are the additional markers of cytotypes I and V. This fact, along with the size of the non-homologous chromosomes of cytotype II (Fig. 4b, boxed), suggests that the submetacentrics of cytotypes II-IV were formed through the Robertsonian fusion of large subtelocentrics (additional markers of cytotype I) with small subtelocentrics. Other subtelocentrics were involved in the formation of large submetacentrics of cytotype V; these were most likely the medium-sized subtelocentrics of cytotype I. Our data show that the large submetacentrics in cytotypes II-IV and cytotype V formed because of various Robertsonian fusions.
The low level of genetic differentiation in mtDNA between the studied M. stelleri (Fig. 3) confirms their close relationship. The karyological analysis suggests that the cytotype I (2n = 44) is close to the ancestral karyotype of the genus Myoxocephalus in general . Within Group SO (Fig. 3), the frequency of cytotype I was 76%. Among the individuals of this group that had 2n = 44, 77% were from the Sea of Okhotsk and 23% were from the northern Sea of Japan (Aleksandrovsky Bay, Zolotaya Bay). Cytotype I was not found among the individuals from the western Sea of Japan (Fig. 3, Group WSJ), which can be explained by their significant karyological divergence from the individuals from the Sea of Okhotsk and the northern Sea of Japan.
The significant chromosomal variability of individuals from the northern Sea of Japan may indicate their later divergence as compared to individuals from the western part of the sea. The assumption about different divergence times from the Okhotsk Sea is confirmed by the fact that the formation of submetacentrics of similar size in their cytotypes (II-IV and V) occurred because of different Robertsonian translocations. We assume that the polymorphism observed in individuals from the northern Sea of Japan could arise in the relatively recent past in the following way: the chromosomal rearrangement that took place in one or more individuals (Fig. 4b, cytotype II) then became fixed (Fig. 4c, cytotype III) and distributed among M. stelleri from the northern Sea of Japan. Cytotype IV (Fig. 4d) from the northern Sea of Japan could be formed because of subsequent hybridization of karyotypically different individuals from the northern and western Sea of Japan. Despite the lack of isolation and their geographic proximity to individuals from the northern Sea of Japan (cytotypes II-IV; Fig. 1), the modern M. stelleri from the western part of the sea have a uniform number of chromo-somes (cytotype V) (Miller 2000;Podlesnykh and Moreva 2014). This may indicate the lack of secondary contact between M. stelleri from the western and northern Sea of Japan at present.
The frog sculpins from the coastal waters off the southern Kuril Islands deserve special attention. In the geological history of the basins of the Far Eastern seas, several long-lasting barriers existed during the regressions of the level of the World Ocean, separating the unified area of the species. One of them, extending along Sakhalin/ Hokkaido/Honshu, kept the Japanese Sea and South Kuril populations separated for a long time. During this period of geographic isolation, the chromosome sets of individuals from the western Sea of Japan and the southern Kuril Islands could have formed independently of each other, despite their visual identity (cytotype V). One of the southern Kuril specimens has haplotype 1b, which connects all other haplotypes found in M. stelleri, and may thus be the closest to the ancestral haplotype. Contrary to the results of mtDNA analysis, the karyological data point to a significant divergence of M. stelleri (2n = 40) from the coastal waters off the southern Kuril Island. This discrepancy may be caused by the high rate of change in karyological traits compared to that of DNA markers (Lukhtanov and Kuznetsova 2009).
Our results do not support the assumption that the individuals of M. stelleri from the Sea of Japan and the Sea of Okhotsk belong to different species (Podlesnykh and Moreva 2014). Judging by the values of genetic and chromosomal polymorphism, the low level of genetic differentiation, and the frequency and spatial distribution of haplotypes, this species is relatively young. It can be suggested that the fragmentation of the single range of M. stelleri, which took place during the last Quaternary glaciation 20-25 thousand years ago, resulted in the geographic isolation of populations from the Sea of Japan and Sea of Okhotsk. The accumulation of Robertsonian translocations in M. stelleri from the western Sea of Japan presumably occurred during this period. We suggest that, as the glaciation ended, the ocean transgression and restoration of connections between the seas of the Northwest Pacific (Nishimura 1964;Lindberg 1972;Yokoyama et al. 2007;Matsuzaki et al. 2018) allowed the secondary colonization of M. stelleri, which likely resulted in a significant chromosomal polymorphism in the northern Sea of Japan. The paleoclimate events also influenced the distribution of genetic variation in M. stelleri. The formation of the three phylogenetic groups is likely associated with the geographic isolation of ancestral forms in different parts of the species range, and with the limited gene exchange and secondary contact because of migrations after the last ice age (Briggs 2003;Altukhov 2005).

Conclusions
Our data have revealed a similarity in chromosome sets, as well as low levels of differentiation in mtDNA between M. stelleri from the Sea of Okhotsk, the Sea of Japan, and the coast of Shikotan Island (southern Kuril Islands), thus, confirming that these represent a single, yet variable, species across its geographic range. The significant chromosomal polymorphism and the presence of common haplotypes in the studied samples indicate their recent origin from a common ancestor and/or relatively recent contacts within the range. The contribution of different Robertsonian translocations to the formation of cytotypes (II-IV and V) of individuals from the northern and western Sea of Japan allows us to conclude that they diverged from the Sea of Okhotsk M. stelleri independently of each other. The star-shaped structure of the haplotype network with a central ancestral haplotype indicates a connection between all constituent haplotypes and the ancestral position of the southern Kuril Islands haplotype (1b). The discrepancy in assessments of the divergence of individuals from the Sea of Okhotsk and waters off the southern Kuril Islands can be attributed to the different mechanisms of inheritance and rates of evolution of karyological traits and mtDNA markers.
The results of our study demonstrate the necessity of further detailed analysis of the widely sampled M. stelleri populations from the Pacific part of the species range and from the southern Sea of Okhotsk. Such studies should include differential chromosome staining and SNP markers, as well as comparative morphological and osteological analyses using up-to-date methods.