We examined the nucleotide sequences of 500 bp variable portion from the 5′ end of mitochondrial (mt) DNA control region in about 500 individuals from 12 populations that were captured in 11 rivers, six in Hokkaido and five in Honshu, Japan. Comparison of the sequences showed 10 variable sites, defining a total of 12 haplotypes in the examined individuals. All the 12 haplotypes occurred in seven Hokkaido populations, whereas only six haplotypes were found in the five Honshu populations. Among these haplotypes, two were common in all the Hokkaido and Honshu populations. The AMOVA analysis inferred a genetic differentiation among three geographic regions, i.e. Hokkaido, Pacific Ocean coast in Honshu, and Japan Sea coast in Honshu. Haplotype diversity was higher in the populations of Hokkaido than those of Honshu, indicating a greater genetic variation in the Hokkaido than the Honshu populations. The estimates of pairwise population FST suggested that the regional differentiation was mostly ascribed to the divergence between populations in Hokkaido and the Pacific coast in Honshu.
INTRODUCTION
Chum salmon (Oncorhynchus keta) is the most widely distributed species of salmon in the North Pacific Ocean (Salo, 1991). Spawning adults, like other anadromous Pacific salmon species, return to the natal river. Such restricted homing behavior will lead geographically distinct populations to partial genetic isolation.
Stock identification of salmon has been attempted by using tagging, scale characteristics, parasite tagging, thermal otolith marking, and/or protein genetic (allozyme) characters (Ishida et al., 1989; Winans et al., 1994; Willmot et al., 1998; Urawa et al, 1998). Although variability of allozymes is beneficial for genetic stock identification (GSI) of chum salmon, allozyme analysis requires careful collection and handling of tissues (Park et al., 1993). Furthermore, resolution of allozymes remains mostly at the regional- to continental-levels (Wilmot et al., 1998).
Recently developed molecular techniques are a powerful replacement which compensates deficits of the use of allozymes with an increase in accuracy and resolution. Because mitochondrial (mt) DNA has higher sequence variability than single copy nuclear genes (Brown et al., 1979) and clonal haploid inheritance, analysis of mtDNA has become a method of choice in phylogenetic and population genetic studies (Moritz et al., 1987). In fact, mtDNA studies, mostly based on the analysis of restriction fragment length polymorphisms (RFLPs), has been conducted so far in many fish species including salmon (Meyer, 1993). Low levels of mtDNA sequence variation shown by RFLP analysis were reported in chum salmon (Park et al., 1993) and other species of Oncorhynchus (Wilson et al., 1987), providing the estimates of genetic divergence similar to those obtained from allozyme analyses (Seeb and Crane, 1999a).
The control region is considered to be the most variable portion of mtDNA, showing two to five times higher rates of nucleotide substitution than the protein-coding regions (Moritz et al., 1987; Meyer, 1993). However, a low variation in the sequence of the mtDNA control region of chum salmon was reported after RFLP (Cronin et al., 1993) and sequence analyses (Park et al., 1993), probably because the samples in these studies were small or collected from small populations of restricted area.
In the present study, we collected more than 500 chum salmon specimens from 11 rivers in the northern part of Japan to examine the variation of nucleotide sequence of mtDNA control region. Then, genetic variability within or among populations was estimated based on the sequence variation.
MATERIALS AND METHODS
Samples
Liver or blood samples were collected from 537 individuals of 12 populations which returned to their natal river in the latter half of breeding season from 1996 to 1999 (Table 1). Allozymes or other characters of chum salmons homed in these rivers were examined for more than a decade to monitor a possible annual difference in the stock (S. Urawa, personal communication). The homing migration of northern population is earlier than that of southern populations, and the peak runs in Hokkaido occur about one month earlier than those in Honshu in Japan (Salo, 1991). Early run may include descendants of populations which were transferred from rivers of other areas, particularly, in Honshu. Therefore, the samples from fish returned in the latter half of breeding season were used in the present study. Mostly late run fish were obtained from six rivers in Hokkaido. A possible genetic difference, if any, between early and late run fish was examined for those collected in the Tokoro River. Late run fish were obtained from five rivers in Honshu, i.e. two in Iwate and one in Miyagi Prefecture in the Pacific Ocean coast, and one each in Akita and Yamagata Prefectures in the Japan Sea coast (Table 1). The samples from the Otsuchi River were collected at random from juveniles of late run fish in the hatchery, which were obtained by the artificial fertilization of gametes from arbitrarily selected multiple parents and kept in large aquarium. The collected samples were stored at −80°C until DNA extraction.
Table 1
Sampling location, date of collection, and the number of chum salmon samples used for mtDNA analysis
DNA extraction
DNA was isolated from the stored specimens following the routine phenol-chloroform method (Sambrook et al., 1989). About 50 μl of whole blood or liver homogenate was added to 500 μl sodium tris EDTA buffer (0.1 M NaCl, 10 mM Tris-HCl, and 1 mM EDTA, pH8.0) containing 500 μg/ml proteinase K and 0.5% SDS, and incubated at 37°C overnight. DNA was extracted with a mixture of phenol (250 μl) and 24:1 chloroform:isoamylalcohol (250 μl) three times and then twice with 500 μl of 24:1 chloroform:isoamylalcohol alone. DNA in aqueous phase was recovered by ethanol precipitation, dried in air, and dissolved in TE (10 mM Tris-HCl, 1 mM EDTA, pH 7.5).
PCR amplification
The control region of mtDNA was amplified by PCR in a 100 μl of reaction mixture containing 25–100 ng of template DNA, 10 mM TrisHCl (pH 8.3), 50 mM KCl, and 1.5 mM MgCl2, 250 nM each dNTP 250 nM of forward and reverse primers, 0.001% gelatin, and 1.25 U of Taq DNA polymerase (TaKaRa, Tokyo). PCR primers were designed based on the reported sequences of rainbow trout O. mykiss mtDNA (Zardoya et al., 1995) and chum salmon control region (Shedlock et al., 1992). They included forward primers, tRNAglu, 5′-AAC CAC CGT TGT TAT TCA ACT A-3′; tRNAthr-2, 5′-TCT TGT AAT CCG GAA GTC GGA-3′; and reverse primers, tRNAphe, 5′-CA(G/T) CTT CAG TG(T/C) TAT GCT TT-3′; tRNAphe-2, 5′-AAC AGC TTC AGT GTT ATG CT-3′. The tRNAglu is located at 5′ to cytochrome b gene, but the others are in tRNA site flanking the control region. The condition of PCR amplification using a GeneAmp PCR System 9700 was as follows: preheating at 94°C for 5 min, followed by 35 cycles of denaturation at 94°C for 45 sec, annealing at 58°C for 45 sec, and elongation at 72°C for 1 min, and completed with final extensions at 72°C for 7 min. The PCR products were purified by the QIAquick PCR Purification Kit (QIAGEN, Hilden, Germany) after examining the size and quality with 1.5% agarose-gel electrophoresis and ethidium bromide staining.
Nucleotide sequence analysis
The mtDNA fragment amplified with tRNAglu and tRNAphe was cloned into vector pCR2.1-TOPO with the TOPO TA cloningTM system (Invitrogen, Carlsbad, CA). After cloning, the insert was analyzed for the sequence using a Hitachi SQ-5500L DNA sequencer (Hitachi, Tokyo). The nested fragment amplified with tRNAthr-2 and tRNAphe-2 was used as a template for direct sequence analysis. Sequence reaction was preformed using a Thermo Sequenase™ pre-mixed cycle sequencing kit (Vistra Systems, Sunnyvale, CA) according to the manufacturer′s instruction. The forward and reverse sequence primers labeled with Texas Red included M13 reverse and T7 primers for cloned products and four primers for direct sequence analyses, which were designed based on own sequence data of the mtDNA control region. They included forward primer, tRNAthr-3, 5′-GGT TAA AAC CCT CCC TAG TG-3′; Okdl-L1, 5′-AGC TTG CAT ATA TAC AAG TGC A-3′, and reverse primers, Okdl-H1, 5′-ATG GGT TCT CTG GAA TTC AA-3′; Okdl-H2, 5′-TGG GTA ACG AGC AAT AAG AT-3′. The 5′ and 3′ halves of the control region were analyzed using a pair of tRNAthr-3 and Okdl-H1 (or Okdl-H2) and a pair of Okdl-L1 and tRNAphe-2, respectively.
Data analysis
The sequence data were aligned by GENETIX-WIN version 4.0.1 (Software Development Co., Ltd, Japan) to find out nucleotide variations among the sequences of control region, from which the haplotypes of control region were defined. Phylogenetic relationship among the haplotypes was estimated by the neighbor-joining method (Saitou and Nei, 1987), based on genetic distance estimated by Kimura's two-parameter method (Kimura, 1980). A parsimony network which connects the observed haplotypes was drawn after Bandelt (1994).
The heterogeneity of the haplotype frequencies within and between geographic regions was evaluated using the contingency χ2test (Roff and Bentzen, 1989), with 17,000 Monte Carlo simulations generated by the CHIRXC program (Zaykin and Pudovkin, 1993). Some intra- and interpopulation parameters were calculated using Arlequin version 1.1 program package for the mtDNA sequence data (Schneider et al., 1997). Haplotype diversity, being equivalent to expected heterozygosity for diploid data, was calculated for each population according to Nei (1973). In order to assess the extent of genetic differentiation at the different level of geographic hierarchy, the overall molecular variance was partitioned into components corresponding to the population divergence within and among regions using the analysis of molecular variance model (AMOVA; Excoffier et al., 1992). For AMOVA, populations were grouped geographically by the above neighbor-joining algorithm, in which the obtained-topology was tested by a bootstrap analysis with 10,000 pseduoreplicate trees. Significance of the variance components and FST values was tested with a permutation method.
Pairwise FST values were calculated to estimate the genetic distance between populations, according to Slatkin (1995). The amount of gene flow (Nm) between populations was estimated using the approximation FST = 1/(2M+1), where M = Nm for haploid data. The estimation of Nm, however, was not conducted if its validity was not supported by the neutrality test after Tajima (1989) and by the regional equilibrium test after Hutchison and Templeton (1999).
RESULTS
Nucleotide sequences of the mtDNA control regions and flanking regions
About 2300 bp fragment was successfully amplified with a primer pair of tRNAglu and tRNAphe from DNA samples of two individuals. The PCR product was cloned and analyzed for the sequence by fluorochrome-labeled M13 reverse, tRNApro-3, T7, and Okdl-H1 primers. The base homology analysis with reported rainbow trout mtDNA sequence (Zardoya et al., 1995; EMBL/GenBank accession number L29771) revealed that the insert contains a 1002 bp control region sequence. It was flanked by a 3′ part of cytochrome b, tRNAThr, and tRNAPro at the 5′ side and with a 5′ part of tRNAPhe at the 3′ side (DDBJ/EMBL/GenBank accession number AB039956), as shown in Fig. 1. This gene order around control region is in keeping with the previous findings in salmon (Shedlock et al., 1992; Zardoya et al., 1995). The three chum salmon tRNAs could be folded into a cloverleaf secondary structure which was assured with the tRNAscan-SE program (Lowe and Eddy, 1997), verifying that the obtained sequence is truly from mtDNA. The observed sequence of the control region was identical with the previously reported chum salmon sequence (Shedlock et al., 1992), except for a few substitutions and indels probably reflecting an interindividual variation.
Nucleotide sequence variation in chum salmon control region
About 1200 bp fragment was amplified with a primer pair of tRNAthr-2 and tRNAphe-2 on DNA samples from 30 individuals from the Chitose and the Otsuchi populations (Fig. 1). The direct sequence analysis on the fragment containing the entire control region and flanking tRNA sites was performed by fluorochrome-labeled tRNAthr-3, Okdl-L1, Okdl-H1, OkdlH2, and tRNAphe-2. Multiple alignment of the obtained sequences revealed a cluster of variable sites (six of seven sites) within 500 bp from 5′-end of the chum salmon control region. Estimation of the 481 bp sequence in the 5′ variable portion disclosed 10 variable sites in a total of 537 individuals from 12 populations, defining 12 haplotypes designated as OKDL-1 to OKDL-12, as shown in Table 2. The observed variation included base substitutions and indels with an excess of trans-version substitutions (six of eight substitutions). The sequences of entire chum salmon control regions were registered in the DDBJ/EMBL/GenBank with accession numbers AB039890 to AB039901. Nucleotide variation in the 3′ downstream to the variable portion occurred only in the OKDL-4 haplotype (data not shown).
Table 2
Variable nucleotide positions in the 5′ half of mtDNA control region of chum salmon
A parsimony network which connects the 12 control region haplotypes is presented in Fig. 2a-a. The observed haplotypes of chum salmon could be grouped into three clades based on the nucleotide variation shown in Table 2, i.e., OKDL-1 to -5 in clade 1, OKDL-9 to -12 in clade 2, and OKDL-6 to -8 in clade 3. The T to C transition at nucleotide 30 separated the clade 2 from the clade 1, and a deletion at nucleotide 386 and the C to A substitution at nucleotide 395 discriminates clade 3 from clades 1 and 2, respectively. The T to C transition at nucleotide 231 in OKDL-3, -8, and -12 was homo-plasmic. The neighbor-joining consensus tree, which was resulted from the bootstrap analysis of the 481 bp sequence and rooted with the corresponding region of the rainbow trout mtDNA, showed haplotype clustering that is identical to that observed in the network, but with less than 50% of bootstrap support for the nodes (data not shown).
Genetic differentiation among chum salmon populations
The distribution of 12 haplotypes among 12 populations of chum salmon is presented in Table 3 and Fig. 2-b. All the 12 haplotypes occurred in the seven Hokkaido populations, whereas only six haplotypes were found in the five Honshu populations. All the phylogenetic clades were found in the Kawabukuro, Gakko, and Koizumi populations, whereas the clade 3 was not found in the Tsugaruishi and Otsuchi populations (Fig. 2-b). The OKDL-1 and -9 were common in the 12 populations. The OKDL-7 haplotype was seen in all Hokkaido and three Honshu populations. Early and late run populations from the Tokoro river in Hokkaido showed essentially similar haplotype distribution, but differed in haplotype frequencies. The Hokkaido populations showed a higher level of haplo-type diversity than the Honshu populations (Table 3).
Table 3
Distribution of mtDNA control region haplotypes and haplotype diversity within the 12 populations calculated from mtDNA haplotype frequencies
Highly significant heterogeneity in the haplotype frequencies was revealed for the entire set of populations (p<0.001), seven Hokkaido populations (p<0.001), and Honshu populations (p<0.005) by contingency χ2 test. Such significant heterogeneity was also observed for the set of populations from Hokkaido and Pacific coast or Japan Sea coast in Honshu (p<0.001), although no significant heterogeneity was shown for the populations within the two regions of Honshu (p>0.3). These findings suggest a greater genetic difference among populations between regions than within regions.
The extent of genetic differentiation at the different levels of geographic hierarchy was examined by the AMOVA. Using the neighbor-joining method, three population clusters, Pacific coast in Honshu, all Hokkaido and the Kawabukuro, and the Gakko, were obtained on the consensus tree, with more than 90% of bootstrap support for the clusters of Pacific coast in Honshu and the Gakko, respectively (Fig. 3). Early and late runs from the Tokoro River were separated from each other, but with less than 50% of the nodal value, and the heterogeneity in their haplotype frequencies was not significant by the χ2 test (p>0.5). Taking together these observations and the results of the above contingency χ2 test, grouping of three geographic regions, i.e. Pacific coast in Honshu, Hokkaido, and Japan Sea coast (including both the Kawabukuro and Gakko) was considered to be most appropriate for AMOVA. The AMOVA analysis indicated that the magnitude of variance among three regional groups was moderate but significant (11.32%, P<0.001). However, the variance among populations within regions was very low (<1%, P>0.1). This fact suggests that the populations within regions are genetically indistinguishable. Most of the variation occurred within populations (87.70%, P<0.0001). The results of AMOVA imply that the three regional groups of populations are genetically differentiated from each other.
The increased level of genetic isolation among populations between regions was suggested for Hokkaido, except for the Yurappu, and Pacific coast in Honshu by high FST values (0.103 to 0.328) compared with other pairwise estimate, as shown in Table 4. These results imply a low level of gene flow among populations between Hokkaido and Pacific coast in Honshu. Tajima's test (1989) supported the neutrality for the mtDNA variations in all populations (p<0.01, D= −2.092 to 2.840 or p<0.05, D= −1.734 to 2.130), except for the one in the Nishibetsu population (p>0.1, D= 3.7568). In addition, the relationship between FST and geographic distance in the seven Hokkaido populations inferred a lack of regional equilibrium, i.e. absence of isolation by distance (data not shown). From these observations, estimation of Nm was not conducted in the present study. Correlation between genetic and geographic distance was not tested for other regions, because the number of populations was small and pairwise FST values between populations were mostly not significant.
Table 4
Pairwise estimates of Fst based on mtDNA sequence data
DISCUSSION
The present mtDNA analysis demonstrated: 1) 10 variable positions in the nucleotide sequences of about 500 bp from the 5′ end of control region (Table 2), defining a total of 12 haplotypes among more than 500 individuals from 12 Japanese populations of chum salmon (Table 3); 2) probable genetic differentiation among three regional groups of Japanese populations of chum salmon, i.e. Hokkaido, the Pacific coast in Honshu, and Japan Sea coast in Honshu; 3) greater genetic variation in the populations of Hokkaido than those of Honshu; and 4) moderate but significant genetic isolation of the populations between Hokkaido and the Pacific coast in Honshu.
The sequence variation detected in the control region of Japanese chum salmon was similar to or even higher than the diversity in other Oncorhynchus species including O. mykiss, O. kisutch and O. tshawytscha (Nielsen et al., 1994). Thus, the number of observed haplotypes is apparently greater than those obtained from the previous RFLP analyses on chum salmon mtDNA (Cronin et al., 1993; Park et al., 1993; Bickham et al., 1995; Seeb and Crane, 1999b). These findings indicate an increased potential of mtDNA sequence analysis to estimate the genetic structure of chum salmon populations.
The 12 haplotypes were phylogenetically classified into three clades (Fig. 2-a). Because of their prevalence in the examined populations, the haplotype OKDL-1, -9 and -7 are considered to be ancestral one within the clade 1, 2 and 3, respectively, and hence parental to the other haplotypes less frequent in the each clade. However, the frequency of OKDL-7 was low in Honshu populations and this haplotype was not detected in the Tsugaruishi and Otsuchi (Fig. 2-b). The decreased number of haplotypes in these two populations suggest that they might have gone through the decrease of population size, that resulted in extinction of the lineage with the OKDL-7 and other haplotypes.
Results of AMOVA were corroborated by the contingency χ2 test for the heterogeneity in the haplotype frequencies among populations. Genetic differentiation among the three regional groups of Japanese chum salmon, as inferred from the AMOVA, was also suggested by the previous allozyme analysis (Kijima and Fujio, 1979). Little genetic differentiation among populations within regions may partly be associated with a level of straying rate (about 14%) in spawning migration of chum salmon (Salo, 1991), which could cause a gene exchange between rivers. The estimates of haplotype diversity and pairwise population FST strongly suggest that the regional differentiation is mostly ascribed to the divergence between populations in Hokkaido and the Pacific coast in Honshu. Probable low gene flow between Hokkaido and the Pacific coast in Honshu may be contributed by the differences in the route of spawning migration (Kijima and Fujio, 1982; Okazaki, 1986) and run timing (Salo, 1991). In fact, the previous allozyme study suggested that separation by at least 600 km is required for genetic differentiation to occur (Kijima and Fujio, 1982). This estimation well fits the minimum distance between the populations of Hokkaido, expect for the Yurappu, and the Pacific coast in Honshu.
The estimated genetic structure of Japanese chum salmon populations may reflect either the effect of historical factors such as the glacial advancement and resulting shift of the species range to southward or contemporary forces such as human-mediated extensive transplantation of stocks, or both. In fact, the Hokkaido and Honshu populations have undergone extensive hatchery operations and transplantation of eggs and fry from one river population to another in history of raising commercial salmon production (Kijima and Fujio, 1982). Such an extensive transplantation of stocks from Hokkaido might be a cause of the Kawabukuro being positioned within a cluster of Hokkaido populations (Fig. 3). Absence of isolation by geographic distance in Hokkaido populations implies a possible contribution of an artificial transplantation of stocks among rivers. The lack of such regional equilibrium also suggests little gene flow and divergence mostly by drift in the formation of Hokkaido populations (Hutchison and Templeton, 1999). The regional equilibrium test using a larger number of populations in Honshu should be performed to examine the factor(s) involved in the formation of their genetic structure. In any case, mtDNA analysis on foreign populations such as Russia, USA, and Canada would be helpful to understand how population structures in Japanese chum salmon were established.
Acknowledgments
We would like to thank Drs. Dmitri Churikov and Anthony Gharrett, Juneau Center, School of Fisheries and Ocean Sciences, University of Alaska Fairbanks, for valuable advice on statistic analysis. We also would like to thank the Tokoro fisherman's cooperative association and the Otsuchi fisherman's cooperative association for their help in sampling. This study was supported in part by Grants-in-Aid from the Fisheries Agency and the Ministry of Education, Science, Sports and Culture, Japan.