INTRODUCTION
Careproctus is one of the most species-rich genera in family Liparidae (Orr et al., 2015, 2019), and includes more than 100 species, more than 50 of which are found in the North Pacific (Chernova et al., 2004). Of these, species that are morphologically similar to the salmon snailfish Careproctus rastrinus Gilbert & Burke, 1912 collectively constitute the Careproctus rastrinus species complex. Kai et al. (2011) conducted the first genetic study of this species complex, and molecular phylogenetic analyses based on the 16 S ribosomal RNA and cytochrome b (Cytb) mitochondrial genes revealed that this species complex consisted of nine monophyletic lineages, including one from the Pacific Ocean, one from the Arctic Ocean and Gulf of Alaska, three from the Bering Sea, one from the Okhotsk Sea, and three from the Japan Sea.
Orr et al. (2015) revised the species in the Careproctus rastrinus complex based on their morphological characteristics, and concluded that the complex contained eight species. Of these, the pellucid snailfish Careproctus pellucidus Gilbert & Burke, 1912, from the Cytb lineage PAC1 described by Kai et al. (2011), and C. rastrinus, from the lineage OKH1, were found to be the sole and endemic species of the Pacific Ocean and the Okhotsk Sea, respectively. The diagnostic morphological characteristics distinguishing these two species, as described by Orr et al. (2015), are summarized in Table 1. However, no objective standard is available for color characteristics, and the quantitative characteristics overlap between the two species.
Table 1.
Diagnostic morphological characteristics for discriminating C. pellucidus and C. rastrinus.
img-ApY1_01.gifTable 2.
List of sampling sites.
img-z2-32_01.gifFollowing Kai et al. (2011), Hashiguchi (unpublished data) further identified individuals with the C. rastrinus-type mitochondrial DNA (mtDNA), from the Pacific Ocean, off the Sanriku Coast, in northeastern Japan, based on the nucleotide sequences of the Cytb gene. Therefore, the phylogeographic relationships between C. pellucidus and C. rastrinus need to be confirmed. In this study, the genetic relationships between these two species were analyzed, based on the nucleotide sequences of the mitochondrial Cytb gene, the first intron region of the nuclear S7 ribosomal protein gene (S7), and 11 microsatellite loci, from a total of 156 specimens. In addition, their color characteristics were evaluated quantitatively. The results of this study provide important information for revising the Careproctus rastrinus species complex.
MATERIALS AND METHODS
A total of 100 Careproctus specimens were collected from the Pacific Ocean and the Okhotsk Sea (Table 2). The specimens were preserved either by freezing or in ethanol, at the Atmosphere and Ocean Research Institute, The University of Tokyo (AORI), and the Hokkaido University Museum (HUMZ). Total DNA was extracted from the muscle of each preserved specimen, using a DNeasy Blood and Tissue Kit (Qiagen, Valencia, CA, USA), according to the manufacturer's protocol.
A part of the mitochondrial Cytb gene and the first intron region of the nuclear S7 gene were amplified by polymerase chain reaction (PCR), using the primer sets CareCBfwd (5′-GCTTCCTGCACCCTCTAATATCTCA-3′) and CareCBrev (5′-ATGAAGCCAGCTAGGGGGAATATT-3′) for Cytb (Kai et al., 2011), and S7RPEX1F (5′-TGGCCTCTTCCTTGGC-CGTC-3′) and S7RPEX2R (5′-AACTCGTCTGGCTTTTCGCC-3′) for S7 (Chow and Hazama, 1998). The PCR amplification was performed in a final reaction volume of 20 μl containing < 10 ng of the template DNA, 2 μl of 10 ×ばつ Ex Taq buffer, 2.0 μl of dNTPs (2.5 μM), 0.125 μl of Takara Ex Taq Hot Start Version DNA polymerase (5 unit/μl; Takara Bio Inc., Shiga, Japan), 2.5 μl of 10 μM forward and reverse primers, using a Thermal Cycler (Applied Biosystems, Foster City, CA, USA). The PCR conditions were as follows: an initial denaturation at 94°C for 3 minutes (min), followed by 35 cycles of denaturation at 94°C for 1 min, annealing at 58°C for 1 min (Cytb) or 45 s (S7), and extension at 72°C for 2 min (Cytb) or 1 min (S7), with a final extension at 72°C for 10 min. The PCR products were separated electrophoretically in a 2% agarose gel at 100 V for 25 min. Subsequently the gel was stained with 0.5 μg/ml ethidium bromide for about 20 minutes, and the DNA fragments were visualized under ultraviolet light. To degrade the remaining primers and nucleotides, 1.5 μL of each PCR product was mixed with 0.3 μL of ExoSAP-IT (United States Biochemical, Cleveland, OH, USA) and incubated at 37°C for 40 min and then at 80°C for 15 min. Purified PCR products were used for the cycle sequencing reactions, carried out using the BigDye Terminator Cycle Sequencing Kit, version 3.0 (Applied Biosystems) with the same PCR primers. The nucleotide sequences were determined bidirectionally, using an ABI 3130 automated DNA sequencer (Applied Biosystems).
As the nucleotide sequences of some S7 PCR products could not be determined by direct sequencing, TA cloning was performed. The PCR products were subcloned into the pMD20-T vector of the Mighty TA cloning kit (Takara Bio Inc.), following the manufacturer's protocol. Eight clones per sample were used for direct PCR with forward and reverse primers (M13M4 and M13RV, Takara Bio Inc.). The same protocols were used for direct sequencing. No sequence variation was observed among clones obtained from individual specimens.
The obtained sequences were assembled using the MEGA X version 10.0.5 software package (Kumar et al., 2018). The nucleotide sequences obtained in the present study were deposited in the DNA Data Bank of Japan (DDBJ), European Molecular Biology Laboratory (EMBL), and National Center for Biotechnology Information (NCBI) GenBank databases under the accession numbers LC534911–534925 (Cytb) and LC534926–534940 (S7).
To clarify the phylogenetic relationships among the haplotypes, a haplotype network was constructed using the TCS method in PopART software (Leigh and Bryant, 2014). The Cytb sequences of C. pellucidus (PAC1) and C. rastrinus (OKH1), determined by Kai et al. (2011), were retrieved from the DDBJ/Genbank/EMBL databases.
The unbiased fixation index, FST (Weir and Cockerham, 1984), was estimated, and its significance was tested using a nonparametric permutation approach (10,000 permutations) performed in Arlequin version 3.5.1.2. (Excoffier and Lischer, 2010). Differences in the haplotype frequencies between the populations were examined using the exact test of population differentiation (Raymond and Rousset, 1995a) in Arlequin.
Eleven microsatellite loci developed for Acanthopterygii (Orla2-91, Orla3-65, Orla5-131, Orla4-222, Orla8-113, Orla9-204, Orla18-49, Orla20-134, Orla21-231, Orla22-135, and Orla23-61; Gotoh et al., 2013) were used in the microsatellite analyses. Three multiplex PCR panels were designed using Multiplex Manager 1.0 (Holleley and Geerts, 2009); universal probe target sequence (Blacket et al., 2012) was added to the 5′-end of the forward primer and the pigtail sequence (5′-GTTTCTT; Brownstein et al., 1996) was added to the reverse primer, to reduce noise peaks during the fragment analysis. The PCR amplification was performed using a Thermal Cycler (Applied Biosystems). Each 5 μL of the multiplex PCR mixture contained 0.005 μL of each universal probe (100 μM), 0.0004 μL of forward primer (50 μM), 0.0008 μL of reverse primer (50 μM), 2 μL of 2 ×ばつ QIAGEN Multiplex PCR Master Mix (Multiplex PCR Kit, QIAGEN), and 1 μL of DNA. The following conditions were used for PCR: initial denaturation at 95°C for 15 min, followed by 35 cycles of 94°C for 30 s, 53°C for 30 s, and 65°C for 30 s, with a final extension at 60°C for 30 min. Fragments were analyzed by electrophoresis, using an ABI3130 automated DNA sequencer (Applied Biosystems), with GeneScan 500 Liz (Applied Biosystems) as a size marker. Allele size was determined using Peak Scanner Software version 2.0 (Applied Biosystems). Deviation from the Hardy-Weinberg equilibrium (HWE) for each locus and linkage disequilibrium were tested using Genepop version 4.7 (Raymond and Rousset, 1995b; Rousset, 2008), using an exact probability test (Markov chain parameters: 10,000 dememorisations, 100 batches, and 1000 iterations per batch), with the sequential Bonferroni's correction (Rice, 1989).
A Bayesian clustering approach was adopted in Structure version 2.3.4 (Pritchard et al., 2000). The burn-in period was set to 100,000 iterations, following which a Markov chain Monte Carlo (MCMC) program was run 1,000,000 times. The number of clusters (K) was set to 1–7, and 10 independent runs were conducted for each K value. The K value showing the highest ΔK was determined to be the most suitable (Evanno et al., 2005), using the online program Structure Harvester version 0.6.94 (Earl and vonHoldt, 2012). In the run with the most suitable K, the result showing the highest posterior probability (Ln P(D)) was used. The results of the independent runs for each K value were merged using a greedy algorithm in Clumpp 1.1.2 (Jakobsson and Rosenberg, 2007), and a bar plot for each K value was then generated using Distruct (Rosenberg, 2004).
A principal coordinate analysis (PCoA) was conducted based on genetic distances between individuals (DA) (Nei et al., 1983), using Populations version 1.2.32 (Langella, 1999) and GenAlEx version 6.5 (Peakall and Smouse, 2012). The genotypes were clustered by discriminant analysis of principal components (DAPC; Jombart et al., 2010) executed using the program adegenet version 2.1.0 (Jombart, 2008) for R version 3.4.4 software (R Core Team, 2016).
To evaluate the genetic differentiation between the mtDNA lineages, pairwise FST values were calculated based on microsatellites. The significance of the values was tested using a nonparametric permutation approach (10,000 permutations) in Arlequin.
Peritoneum and stomach colors of specimens of sufficiently sized specimens collected in 2019 were evaluated quantitatively based on three indices in the CIE L*a*b* color space, i.e., lightness (L*), red/green value (a*), and yellow/blue value (b*), which were measured using a CR-300 color difference meter (Konica Minolta Inc., Tokyo, Japan). Differences in the average values were tested, using the t-test if the F-test showed no significant difference in variance or using the Mann-Whitney U test if the F-test showed a significant difference in variance. A principal component analysis (PCA) was conducted based on peritoneum and stomach colors. All of these analyses were conducted using BellCurve for Excel version 2.02 (Social Survey Research Information Co., Tokyo, Japan).
RESULTS
From the 100 Careproctus specimens, 692 bp of the Cytb sequences were determined, resulting in 15 haplotypes with 23 nucleotide substitutions and no indels (insertion/deletion) ( Supplementary Table S1 (zs200011.s1.pdf)). Based on the haplotype network (Fig. 1), the specimens could be classified into the two lineages described by Kai et al. (2011), namely, PAC1 (n = 50) and OKH1 (n = 50). Specimens of the OKH1 lineage were obtained from both the Pacific Ocean (n = 22) and the Okhotsk Sea (n = 28), and are hereafter referred to as OKH1pac and OKH1okh, respectively. Only a single haplotype was obtained from the OKH1pac individuals, and it was the most dominant haplotype in the OKH1okh individuals. The test of FST value (FST = 0.01990) and the exact test based on Cytb did not indicate any significant genetic difference between OKH1pac and OKH1okh individuals (P > 0.05). Half (50.0%) of the individuals collected at stations off Miyagi Prefecture (Stations P1–5) were OKH1pac individuals. However, only a single OKH1pac individual (9.1%) was collected at stations off Iwate Prefecture (Stations 8–9), which were to the north of the Miyagi Prefecture stations, and no OKH1pac individuals were collected off Fukushima Prefecture (Stations 6–7) or Ibaraki Prefecture (Station 10), which were to the south of Miyagi Prefecture ( Supplementary Table S1 (zs200011.s1.pdf)).
Fig. 1.
A statistical parsimony haplotype network of Careproctus pellucidus and Careproctus rastrinus based on nucleotide sequences of the cytochrome b (Cytb) mitochondrial genes. Haplotype names are the same as in Supplementary Table S1 (zs200011.s1.pdf). * and ** denote the same sequences as PAC1 and OKH1 individuals in Kai et al. (2011), respectively. The areas of the circles are proportional to the frequency of the occurrence of the haplotypes. Black and white sections denote the relative occurrence frequency of individuals from the Pacific Ocean and the Okhotsk Sea, respectively. Small closed circles are intermediate haplotypes that were not observed in the present study.
Fig. 2.
A statistical parsimony haplotype network of Careproctus pellucidus and Careproctus rastrinus based on nucleotide sequences of the first intron region of the nuclear S7 ribosomal protein gene. Haplotype names are the same as in Supplementary Table S2 (zs200011.s1.pdf). The areas of the open circles are proportional to the frequency of the occurrence of the haplotypes. Black, gray, and white sections denote the relative occurrence frequency of PAC1, OKH1pac, and OKH1okh individuals, respectively.
Table 3.
Pairwise FST values based on the nucleotide sequences of the first intron region of the nuclear S7 ribosomal protein gene. The value indicated in parentheses is based on the individuals collected off Miyagi Prefecture only. Significant differences between groups are indicated by * (P < 0.05) and ** (P < 0.001).
img-z4-6_01.gifFig. 3.
Influence of the number of clusters (K) that best fit the data. ΔK is the standardized second order rate of change of the posterior probability of each K.
Fig. 4.
Bayesian clustering of 84 Careproctus individuals from the Pacific Ocean and the Okhotsk Sea into two genetically distinguishable clusters.
The S7 nucleotide sequences (422 bp) were determined for 49 specimens collected in 2006, 2011, and 2012. They contained 15 haplotypes with nine nucleotide substitutions and no indels ( Supplementary Table S2 (zs200011.s1.pdf)). The phylogenetic relationships between the S7 haplotypes are shown as a haplotype network in Fig. 2. Haplotypes of the two mitochondrial lineages, i.e., PAC1 and OKH1, did not form exclusive clusters in the network, and five haplotypes were shared between the lineages. The tests of pairwise FST values based on the S7 sequences indicated significant genetic differences between OKH1pac and OKHokh individuals (P < 0.01; Table 3). Genetic differences were also detected between Pacific Ocean individuals with PAC1 and OKH1pac mtDNA (P < 0.05), although the FST value was considerably lower (Table 3). The exact test revealed significant genetic differentiation between OKH1pac and OKHokh individuals (P < 0.01), but not between Pacific Ocean individuals with PAC1 and OKH1pac mtDNA (P > 0.05). In the area off Miyagi Prefecture (Stations P1–5), no significant genetic difference was observed between PAC1 and OKH1pac individuals, as evident from the results of FST value (P > 0.05; Table 3) and exact (P > 0.05) tests.
Fig. 5.
Principal coordinate analysis (PCoA) plot based on 11 microsatellites from 84 Careproctus individuals from the Pacific Ocean and the Okhotsk Sea. Open circles, closed circles, and open triangles denote PAC1, OKH1pac, and OKH1okh individuals, respectively.
Fig. 6.
Results of the discriminant analysis of principal components (DAPC) showing genetic similarities between PAC1, OKH1pac, and OKH1okh individuals.
A total of 84 specimens were used for microsatellite analyses. The numbers of alleles per locus and the observed and expected heterozygosities are summarized in Supplementary Table S3 (zs200011.s1.pdf). Significant deviations from HWE after the sequential Bonferroni's correction (P < 0.05) were observed in Orla9-204 for PAC1, Orla18-49 and Orla9-204 for OKH1pac, and Orla3-65 and Orla23-61 for OKH1okh. Significant linkage disequilibrium was not found for any pairs. A Bayesian clustering analysis based on 11 microsatellites was implemented to estimate the genetic population structure of Careproctus individuals from the Pacific Ocean and the Okhotsk Sea. The most probable number of clusters (K) was estimated to be two (Fig. 3). In the Bayesian clustering of individuals, OKH1pac individuals exhibited genetic characteristics intermediate to PAC1 and OKH1okh individuals (Fig. 4). PCoA also revealed a trend similar to that observed in the Bayesian clustering (Fig. 5). The first and second principal coordinates explained 9.14% and 6.20% of the total variation, respectively. In the scatter plot of the first and second scores, the plotted areas of PAC1 and OKH1okh individuals did not overlap with each other; however, plotted areas of OKH1pac individuals overlapped with those of both the aforementioned individuals. The DAPC analysis showed that OKH1pac individuals were more closely related to PAC1 individuals than to OKHokh individuals (Fig. 6). Pairwise FST values based on microsatellites indicated significant genetic differentiation between all pairs among the three groups (P < 0.01); however, the FST value between PAC1 and OKH1pac individuals was slightly lower (Table 4). In the area off Miyagi Prefecture (Sts. P1–5), no significant genetic difference could be detected between the PAC1 and OKH1pac individuals by the tests of FST value (P > 0.05; Table 4) or the exact test (P > 0.05).
Fig. 7.
Peritoneum (A) and stomach (B) colors of Careproctus individuals from the Pacific Ocean and the Okhotsk Sea, expressed as lightness (L*), red/green value (a*), and yellow/blue value (b*). Open circles, open triangles, and closed circles denote PAC1, OKH1pac, and OKH1okh individuals, respectively.
Table 4.
Pairwise FST values based on 11 microsatellites. The value indicated in parentheses is based on individuals collected off Miyagi Prefecture only. Significant differences between groups are indicated by ** (P < 0.01).
img-A0r6_01.gifPeritoneum and stomach colors of 9 PAC1, 1 OKH1pac, and 17 OKHokh individuals were quantitatively evaluated (Fig. 7). Significant differences were observed between the individuals from the Pacific Ocean and the Okhotsk Sea in terms of the stomach L* values (Mann-Whitney U test, P < 0.001), and the peritoneum a* (t-test, P < 0.05) and b* (t-test, P < 0.001) values. As described by Orr et al. (2015), lightness (L*) of the stomach of individuals from the Okhotsk Sea was significantly higher than that of individuals from the Pacific Ocean, although the two ranges overlapped. The OKH1pac individual exhibited greater color similarity to most PAC1 individuals than to most OKHokh individuals (Fig. 7).
DISCUSSION
Our results showed that the OKH1 mitochondrial lineage is distributed not only in the Okhotsk Sea, but also in the Pacific Ocean (Fig. 1), which is different from the findings of Kai et al. (2011). It is worth noting that 22 of the 72 Pacific individuals sampled in the present study were classified as OKH1 individuals. These individuals were collected almost exclusively from the area off Miyagi Prefecture, while few of these individuals were collected from areas north (off Iwate Prefecture) and south (off Fukushima and Ibaraki Prefectures) of Miyagi Prefecture. Kai et al. (2011) only used C. pellucidus specimens collected from Iwate and Fukushima Prefectures, which might explain why they did not detect individuals belonging to the OKH1 mitochondrial lineage from the Pacific Ocean.
In contrast to the OKH1 mitochondrial lineage, individuals with the PAC1-type mtDNA have not yet been obtained from the Okhotsk Sea. According to Orr et al. (2015), C. pellucidus and C. rastrinus have been collected from depths of 145–300 and 114–217 m, respectively. As the Okhotsk Sea and the Pacific Ocean are connected by deep straits, namely, the Kruzenshterna Strait (1900 m deep) and the Bussol Strait (2300 m deep), Careproctus habitats in these areas would not have been isolated from each other, even during the Last Glacial Maximum (LGM). In the Bussol Strait, water flow from the Okhotsk Sea to the Pacific Ocean is in excess of the opposing flow, and most of the water mass flowing into the Okhotsk Sea is estimated to flow back to the Pacific Ocean (Katsumata et al., 2004). Because of the unidirectional migration occurring due to the stronger current flowing from the Okhotsk Sea to the Pacific Ocean, individuals with the OKH1-type Cyt b genes could have migrated to the Pacific Ocean, while those with the PAC1-type Cyt b genes would have found it difficult to migrate to the Okhotsk Sea. Unidirectional migration from the Okhotsk Sea to the Pacific Ocean has also been suggested for the Japan Sea eelpout Bothrocara hollandi (Jordan and Hubbs, 1925), based on the contrasting genetic diversity between the populations of the two areas (Kojima et al., 2011). However, Tohkairin et al. (2016) stated that the population of the barred snailfish Crystallichthys matsushimae (Jordan and Snyder, 1902) in the Okhotsk Sea was most likely formed by immigrants from both the Pacific Ocean and the Japan Sea. As suggested by Tohkairin et al. (2016), this difference in the migratory history may be attributed to a shallower distributional depth of C. matsushimae compared to that of the Careproctus rastrinus complex and B. hollandi.
Although significant genetic differentiation between PAC1 and OKH1pac individuals was detected for the nuclear S7 gene and microsatellites, this differentiation was not significant when the analyses were limited to the individuals off Miyagi Prefecture, where OKH1pac individuals were almost exclusively collected (Tables 3 and 4). The DAPC analysis showed that OKH1pac individuals were more closely related to PAC1 individuals than to OKHokh individuals (Fig. 6). Unfortunately, only a single OKH1pac individual was available for color evaluation, and its peritoneum and stomach were more similar to most PAC1 individuals than to most OKHokh individuals (Fig. 7). If complete reproductive isolation had been established between C. pellucidus and C. rastrinus, OKH1pac individuals would exhibit characteristics that are more similar to OKH1okh than to PAC1 individuals. Thus, the present results suggest that C. pellucidus and C. rastrinus are inbreeding in the Pacific Ocean off Miyagi Prefecture.
Another possibility is that the OKH1-type mtDNAs within the Pacific population might have originated from past hybridization between C. pellucidus and C. rastrinus, resulting in the introgression of mtDNA of the latter into the former. The limited geographical distribution of OKH1pac individuals and the assumed dispersal direction previously discussed suggest that hybridization occurred between migrating C. rastrinus and indigenous C. pellucidus within the Pacific Ocean. As organelle genes are often preferentially introgressed during interspecific hybridization, and as introgression almost exclusively occurs from the local to the invading species, irrespective of the relative densities of the two species (Currat et al., 2008), the possibility of an introgression event from C. rastrinus to C. pellucidus seems to be very low.
In the present study, OKH1pac individuals were collected almost exclusively from the area off Miyagi Prefecture. Such a geographic trend might reflect an interspecific difference in the preference of physical environments. However, it can also be attributed to annual variation, as sampling off Iwate and Ibaraki Prefectures was carried out 8 years after that off Miyagi and Fukushima Prefectures. Systematic sampling of additional specimens, especially from areas around deep straits connecting the Pacific Ocean and the Okhotsk Sea, as well as more comprehensive molecular and morphological analyses are necessary in order to reveal the evolutionary relationships between C. pellucidus and C. rastrinus.