The white-spotted charr (Salvelinus leucomaenis) is a coldwater-adapted fish distributed in far-eastern Asia. To assess phylogeographic patterns of this species over most of its range in the Japanese archipelago and Sakhalin Island, Russia, we examined nucleotide sequences of the mitochondrial DNA (mtDNA) cytochrome b region (557 bp) in 141 individuals from 50 populations. A total of 33 (5.5%) nucleotide positions were polymorphic and defined 29 haplotypes. Phylogenetic analysis assigned the observed haplotypes to four main clades, which were characterized by the idiosyncrasies and discontinuity of geographic distributions. The nested clade analyses revealed that the geographical distribution patterns of some haplotypes and clades were explained by historical event such as past fragmentation. Although substantial genetic differentiation was found among the four main clades, their geographic distributions overlapped extensively in several regions. Since white-spotted charr can potentially use both freshwater and marine environments, coexistence among different lineages can be attributed to secondary contact through range expansion by migratory individuals during multiple glacial periods after interglacial isolation. Finally, our data demonstrate that the current subspecies designation does not reflect the phylogeography of this species based on mtDNA analysis. Hierarchical analysis (AMOVA) also showed that genetic variation was far more pronounced within subspecies than among subspecies (i.e., among discrete regions). These results suggest that each population, rather than each subspecies, must be treated as an evolutionarily significant unit.
INTRODUCTION
Recent theoretical and empirical advances in the field of molecular genetics have provided many insights into the assessment of evolutionary genetic relationships among populations for a number of species (Avise, 2000). Glacial-interglacial cycles during the Pleistocene have had profound effects on the evolutionary history of northern temperate salmonid fishes, especially those distributed in North America and Eurasia (e.g., Bernatchez and Dodson, 1991; Turgeon and Bernatchez, 2001; Brunner et al., 2001; Bernatchez, 2001). However, the phylogenetic or population structure of salmonids in Asia, where effects of glaciation were relatively weak, remain largely unknown (Sato et al., 2001).
The white-spotted charr (Salvelinus leucomaenis) is a widespread species in far-eastern Asia (Kawanabe, 1989). Like other salmonids, the white-spotted charr is adapted to coldwater habitats, which restricts its distribution to higher altitudes in lower-latitude regions (Nakano et al., 1996). At higher latitudes, on the other hand, this species uses a range of different habitats within a river, from the headwaters to the mouth, and often to coastal sea water (Fausch et al., 1994). Across the species' range, the white-spotted charr has recently been designated into four subspecies based on zoogeographic patterns and morphological characteristics (Nakabo, 2000): S. l. leucomaenis (Japanese Name: amemasu), S. l. japonicus (yamato-iwana), S. l. pluvius (nikko-iwana), S. l. imbrius (gogi). Roughly, populations north of northern Honshu Island, including Hokkaido Island, Japan, and the Sakhalin Island and Kamchatka Peninsula, Russia, are classified as S. l. leucomaenis. They are characterized by large white spots on their sides, and generally have a migratory (anadromous) life history (Yamamoto et al., 1999). However, this subspecies is easily landlocked in rivers following the construction of barriers that prevent upstream migration (Morita et al., 2000; Shimoda et al., 2002). The other three subspecies are endemic to Honshu Island, and are composed mostly of non-migratory (fluvial) individuals. These subspecies have rather small white spots on their sides and are flecked with red, reddish-orange, or yellow spots. Salvelinus l. imbrius, which is distributed in only a few rivers of the westernmost Honshu Island, is characterized by large spots on or around the head. The distributions of S. l. imbrius, S. l. japonicus, and southern populations of S. l. pluvius are restricted to the high mountainous regions of Honshu Island (Nakano et al., 1996). The southernmost populations of S. l. japonicus as well as S. l. imbrius are currently listed among the “Threatened Local Populations” (Ministry of the Environment, Japan, 2003) due to their locally declining population sizes. The other regions, mainly central Honshu Island, are inhabited by S. l. pluvius.
Although many biologists have attempted to assign aspects of intraspecific phylogeny of the S. leucomaenis complex using several morphological characters (reviewed in Kawanabe, 1989), its taxonomic status remains controversial due to highly variable morphological patterns among and even within populations (Nagasawa, 1989; Yamamoto et al., 2000; Nakamura, 2003). In this study, we examined the geographic distribution pattern of mitochondrial DNA (mtDNA) sequence polymorphisms in 50 populations of white-spotted charr over most of its range in the Japanese archipelago and Sakhalin Island, Russia, to assess how its spatial population-genetic structure has been moulded. Genetic data were also used to elucidate the congruence of phylogeographic patterns with current taxonomy of the four subspecies based on morphologic and zoogeographic patterns.
MATERIALS AND METHODS
Sample collection
Fifty populations of white-spotted charr (141 individuals) were sampled from their geographic range in the Japanese archipelago (49 populations) and Sakhalin Island, Russia (one population; Fig. 1; Appendix). We found no records that charr have ever been artificially released at the selected sampling sites, so we considered all caught individuals to be native. Dolly Varden (S. malma) from Hokkaido Island and bull charr (S. confluentus) from North America were used as the outgroup, because these two species are considered to form the sister group of white-spotted charr (Phillips et al., 1994; Phillips et al., 1995). Each white-spotted charr population was further classified into four subspecies in accordance with zoogeo-graphic distinctions made by Kawanabe (1989) and Nakabo (2000): LE (abbreviation of the subspecies name leucomaenis), JA (japonicus), PL (pluvius), and IM (imbrius). Exact sampling locations and taxonomic designations are presented in the Appendix.
Fig. 1
Sample localities of white-spotted charr (Salvelinus leucomaenis). Details of the sample sites with a list of abbreviations and sample sizes are provided in the Appendix.
DNA and data analyses
Total genomic DNA was extracted from the adipose fin by proteinase K digestion, phenol/chloroform extraction, and ethanol precipitation. The cytochrome b region of mtDNA was partially amplified using the primers H15915 (5′-ACCTCCGATCTYCGGATTACAAGAC-3′; Aoyama et al., 2000) and L15285 (5′-CCCTAACCCGVTTCTTYGC-3′; Inoue et al., 2000). Thirty cycles of amplification with 30 sec at 94°C, 30 sec at 55°C, and 60 sec at 72°C were preceded by a predenaturation step (11-min at 94°C) and followed by an additional 7-min extension step at 72°C. Amplifications were performed in 50-μl volumes containing 5 μl of 10[.dotmath]PCR Buffer, 4 mM MgCl2, 0.2 μM of each dNTP, 8 pmol of each primer, and one unit of AmpliTaq DNA polymerase (Applied Biosystems Inc.). A nucleotide sequence data reported here is available from DDBJ, EMBL and Genebank accession number AB111031.
DNA sequences were aligned with the multiple sequence editor CLUSTAL W (Thompson et al., 1994). Inter- and intraspecific relationships of haplotypes were inferred by the maximum parsimony (MP) analyses using PAUP version 4.0b10 (Swofford, 2000). Maximum parsimony analyses were performed using heuristic searches with tree bisection-reconnection branch swapping algorithm. Bootstrap analysis (Felsenstein, 1985) with 1000 peudoreplicates was used to measure support of the resulting topologies. To examine geographic and genetic structuring of white-spotted charr populations, minimum spanning network for nested clade analysis was calculated using the program TCS version 1.13 (Clement et al., 2000). Haplotypes were parsimoniously connected with the 95% probability limits, and then converted into a nested design following the procedures of Templeton et al. (1987), up to the final level of nesting comprising the entire network. Two statistical parameters, the clade distance, Dc, which measures the geographical range of a particular clade, and the nested clade distance, Dn, which measures how a particular clade is geographically distributed relative to its closest evolutionary sister clades (Templeton, 1998) were calculated at all hierarchical levels using the program GEODIS (Posada et al., 2000). These statistics were recalculated after each of 1000 random permutations of clades or haplotypes against sampling locality to test the null hypothesis of no geographic association. Geographical distances among populations were measured as the minimum sea-shoreline distances between the river mouths. The observed clade and nested clade distances were then contrasted to the null distribution to infer which distances are statistically significantly large and which are significantly small. When the two distance statistics and interior-tip (I-T) contrasts for each clade were significantly small or large, then the results were interpreted with the inference key provided on the GEODIS web page ( http://inbio.byu.edu/Faculty/kac/crandall_lab/geodis.htm), modified from Templeton (1998).
We also assessed the phylogenetic relationships among the four subspecies using the constraint-tree option in PAUP 4.0b10 (Swofford, 2000). Differences in tree topologies were compared between maximum parsimony (MP) trees unconstrained and constrained by current taxonomy, where tree-length differences were statistically evaluated using the Templeton test (Templeton, 1983). Finally, a hierarchical analysis of molecular variance (AMOVA) was performed using Arlequin version 2.0 (Schneider et al., 2000) to compare the component of genetic diversity for the variance among the four subspecies to that observed within each subspecies.
RESULTS
Geographical distribution of haplotypes and nested clade analysis
Of the 557 bp of the mitochondrial cytochrome b gene from 141 sequenced white-spotted charr, a total of 33 (5.5%) nucleotide positions were polymorphic and defined 29 haplotypes (Table 1). Of the 33 nucleotide position variations, 27 (82% in variable sites) were in third-codon positions, five (15%) were in first positions, and one (3%) was in second position. The estimated ratio of transitions to trans-versions was 4.7 (28/6). The most divergent haplotypes (Hap-2, Hap-29) differed by 12 substitutions (2.1%). Of the 29 haplotypes identified, 20 were found to be specific to a river, the other 9 were distributed across rivers (see Appendix).
Table 1
Mitochondrial sequence variation from 557 bp of the cytochrome b region of 141 samples from white-spotted charr and two species of outgroups, Dolly Varden (Salvelinus malma) and bull charr (S. confluentus). A nucleotide sequence data of Hap-1 is available from DDBJ, EMBL and Gene Bank accession number AB 111031.
A maximum-parsimony tree, which connects the 29 haplotypes with two outgroup haplotypes of Dolly Varden and bull charr, was supported by a high consistency index (CI=0.861) and bootstrap probabilities (54–100% at intern-odes; Fig. 2). The root of the haplotypes, estimated by the outgroup sequences, was located between Hap-1 and Hap-15. All 29 haplotypes could be parsimoniously connected in a single network with 95% probability (Fig. 3). The whole network was included in four 2-step clades (clade 2-1, 2-2, 2-3, 2-4), each clade being composed of a different number of 1-step clades: two clades (clade 2-1), six (2-2), two (2-3), and two (2-4). Of these, clade 2-2 was characterized by forming a star-like topology radiating from a haplotype (Hap-3). Fig. 4 shows the geographic distribution of each clade. Haplotypes belonging the clades 2-1 and 2-2 were widely distributed throughout the investigated range, including Honshu, Hokkaido, and Sakhalin Island, where encircled the geographical areas of all four subspecies. Of these, Hap-1 (clade 2-1) and Hap-3 (clade 2-2) were found across an immense geographic area, from central Honshu Island north to northern Hokkaido Island and Sakhalin Island (see Appendix). Clades 2-3 and 2-4, on the other hand, had relatively narrow distributions, which were restricted to central and southwestern regions of Honshu Island. Clade 2-3 (Hap-28 and Hap-29) and 2-4 (Hap-17, Hap-18, Hap-19) were found only in two rivers [Kiso (JA-4) and Kumano (JA-6)] and in three rivers [Ane (an inlet stream of Lake Biwa; JA-5), Kuzuryu (PU-10), and Tenjin (PU-12)], respectively. Coexistences of each higher-level clade were observed in seven of the 50 studied populations. The rivers Pilenga (Sakhalin Island; LE-1), Ishikari (Hokkaido Island; LE-12), Kurobe (Honshu Island; PU-8), and Tone (Honshu Island; PU-3) had haplotypes belonging to clades 2-1 and 2-2; Kumano (Honshu Island; JA-6) had clades 2-1 and 2-3; Tenjin (Honshu Island; PU-12) had clades 2-1 and 2-4; and Kuzuryu (Honshu Island; PU-10) had clades 2-2 and 2-4.
Fig. 2
Maximum parsimony tree of haplotypes detected from 50 populations of white-spotted charr and two species of outgroups, Dolly Varden (Salvelinus malma) and bull charr (S. confluentus). Haplotype abbreviations correspond to those listed in the Appendix. Bootstrap probabilities of 1000 resamplings are shown on the internodes.
Fig. 3
The 95% parsimoniously set of cladogram for mtDNA cytochrome b haplotypes detected from 50 populations of white-spotted charr. Each connection represents one mutation step. 0 indicates an interior node in the network that was not in the sample. Boxes represent clades of increasing number of steps. Haplotype abbreviations correspond to those listed in the Appendix.
Fig. 4
Geographic distributions of each clade defined by the parsimony network shown in Fig. 3. The broken lines indicate the general boundary among the distribution areas of the four subspecies.
Of the twelve 1-step clades identified in the network, four clades (clade 1-8, 1-9, 1-10, 1-12) were not informative due to the limited amount polymorphism and the restricted distribution of haplotypes in clades. Three 1-step clades (1-2, 1-3, 1-4) and two 2-step clades (2-1, 2-2) and a total cladogram showed significant geographical associations in the nested contingency tests (Table 2). According the inference chain the geographical arrangements of these haplotypes and clades could be best interpreted as past fragmentation (Table 3). The restricted gene flow/dispersal with some long-distance colonization was only inferred for total cladogram.
Table 2
Results of the nested contingency tests of geographical associations for clades. Nested design and clade designation are given in Fig. 3.
Table 3
Results of the nested geographical analysis of the S. leucomaenis mtDNA haplotypes, following the inference key modified from Templeton (1988). Nested design and clade designation are given in Fig. 3. Following the name of haplotypes or clade number are the clade (Dc) and nested clade (Dn) distances. A bold S means that the distance measure was significantly small at the 5% level, and a bold L means that the distance measure was significantly large at the 5% level. The average difference between interior versus tip clades for both distance measures is given in the row labelled I-T.
Relationships among subspecies
Overall haplotype (h) and nucleotide (π) diversities in the white-spotted charr complex were 0.890±0.018 (SD) and 0.0049±0.0029 (SD), respectively. The indices of diversity among the four subspecies are shown in Table 4. The small sample size of S. l. imbrius (three populations) might reflect the small number of haplotypes and low genetic diversity. For the other three subspecies, S. l. japonicus had the highest diversity indices, while S. l. leucomaenis had the lowest. Analysis of molecular variance (AMOVA) revealed significant spatial patterns of genetic structure both among sub-species (P<0.05) and within subspecies (P<0.01; Table 5). Most of the genetic variance was distributed within subspecies (84.4%), whereas among subspecies, variance was only 15.6%. The Templeton test showed a statistically significant difference between the unconstraint MP tree and the MP tree from the constraint analysis of four subspecies (Wilcoxon rank test; Z=−3.21, P=0.0013).
Table 4
Haplotype diversity and nucleotide diversity of white-spotted charr populations grouped by four subspecies.
Table 5
Hierarchical analysis based on the genetic distance among subspecies.
DISCUSSION
General Phylogeographic Patterns
The phylogeographic survey of white-spotted charr throughout most of its distribution range in the Japanese archipelago and Sakhalin Island revealed the existence of the four main clades. They were characterized by the idiosyncrasies and discontinuity of geographic distributions, and showed significant geographical associations in nested contingency tests, suggesting that there has been a strong historical component to the organization of white-spotted charr mtDNA.
Of the four main clades identified, clade 2-1 was considered to be an ancestral (older) mtDNA type of the white-spotted charr, owing to its closest position to the sequences of Dolly Varden and bull charr. In addition, this group had the widest geographic distribution throughout the Japanese archipelago and north to Sakhalin Island. The current distribution pattern of white-spotted charr is thought to reflect species-specific preferred temperature regimes (Takami et al., 1997), and is determined by high-temperature thermal barriers in rivers (Nakano et al., 1996). As a result, individuals are now restricted to higher altitudes in southern regions (mainly central Honshu Island), whereas those of Hokkaido and Sakhalin islands occur along entire river courses, often to coastal sea water (Fausch et al., 1994). Lower among-population mtDNA diversity was observed in S. l. leucomaenis, which was distributed in the more northern portion of the range compared to the other three sub-species (Honshu Island). This may partly explain that high-latitude populations have a high degree of diadromy, and combined with continued gene flow among populations, results in the constraint factor of genetic divergence (McDowall, 1999). In glacial periods, on the other hand, we hypothesize that access to and use of lower-altitude habitats could be possible even in southern populations as thermal barriers moved downstream (Nakano et al., 1996), and eventually allowed colonization by dispersal via seaward migrations. The existence of a widespread haplotype (Hap-1, see Fig. 4) in this group suggests that the extensive dispersal occurred through seaward migration, and that fish were subsequently isolated in some rivers in the southern range of the distribution where water temperatures had risen during the interglacial period. Nested clade analyses revealed that past fragmentation was inferred for the geographical distributions of haplotypes within clade 1-2 (Hap-13, 15, 21). This might result from their long-term isolation, because they are now patchily distributed at the southwestern edge of the range [Kumano (JA-6), Sufu (IM-2), and Takatsu (IM-3)], or at higher altitudes [Tone (PU-3) and Shinano (PU-7)].
The most dominant group identified was clade 2-2. This group was characterized by a star-like appearance in the parsimony network, with many relatively short branches radiating from one major haplotype (Hap-3). In addition, the phylogeographic patterns of some haplotypes and lower-level clades within clade 2-2 in the nested clade analysis were described by past fragmentation. The rapid range expansion of this group, as well as clade 2-1, is indicative by the widespread haplotype Hap-3, which would have been separated and isolated long enough to fix diagnostic haplo-types in each river.
In contrast to clade 2-2, clade 2-3 was geographically very restricted, and was distributed only in the Kiso (JA-4) and Kumano (JA-6) Rivers, on the Pacific side of central Honshu Island. This area corresponds to the southernmost distribution of the genus Salvelinus. Although clade 2-3 implicitly originated from clade 2-2, the former differed by five or more substitutions from the closest haplotype of the latter. This may indicate that gene flow from adjacent rivers has been prohibited extensively in that region, possibly due to being at the edge of the distribution range, and would result in isolation long enough to accumulate the substitutions.
Similar to clade 2-3, clade 2-4 also had a restricted distribution, and was found only in the Kuzuryu (PU-10) and Tenjin (PU-12) Rivers, in west-central Honshu, and the Ane River (JA-5), an inlet river of Lake Biwa, indicating that clade 2-4 had also been isolated for a long time in that restricted region. Despite geographic proximity, the Ane River, which flows into the Pacific Ocean through Lake Biwa, belongs to a different watershed than the Kuzuryu and Tenjin Rivers, which flow into the Sea of Japan. The Ane River population, therefore, seems unlikely to have been established through the dispersal process via diadromous migration; rather, they might have colonized as a result of paleogeographic processes around the Lake Biwa basin, such as stream-capture between adjacent rivers (Matsuura, 1999). Another example of gene flow between the Sea of Japan side and the Lake Biwa basin (Pacific Ocean side) across a watershed is the loach Cobitis taenia (Kimuzuka and Kobayashi, 1983).
Although substantial genetic differentiations were found among the four main clades, their geographic distributions overlapped extensively in several regions. For instance, Hap-21 (Clade 2-1) and Hap-29 (clade 2-3), which differ at nine sites (1.6%) of all positions examined, coexist in the Kumano River in southernmost Honshu Island. A previous calibration derived from phylogeographic patterns for salmonids and successions of Pleistocene glaciation events led to an average estimate of mtDNA substitution rates of 1–2% per million years per nucleotide site (see Bernatchez, 2001). The estimated divergence time between the Hap-21 and Hap-29, hence, corresponds to approximately 0.8–1.6 million years. Thus, it seems plausible to suppose that the coexistence of different clades has been attributed to secondary contact through range expansion during multiple glacial periods after interglacial isolation.
The glacial-interglacial cycles of the Pleistocene have had enormous impacts on the phylogeographic structure of several coldwater-dwelling freshwater fishes, especially in European and North American regions where the impacts of glaciation on hydrology were particularly severe. For instance, it has been suggested for lake whitefish Coregonus clupeaformis (Bernatchez and Dodson, 1991), Coregonus artedi (Turgeon and Bernatchez, 2001), Arctic charr Salvelinus alpinus (Brunner et al., 2001), and brown trout Salmo trutta (Bernatchez, 2001) that the development of ice-cover during the glacial periods forced these species into habitat refugia that resulted in significant genetic differentiation. Although the extent of glaciation in far-eastern Asia was much more limited (Yonekura et al., 2001), climate oscillations might also affect the dispersal and vicariance processes of freshwater fishes through changes in the hydrological network and/or thermal responses of species (Takahashi et al., 2001; Yokoyama and Goto, 2002). Such an inference seems especially applicable to stenothermal, coldwater stream fishes such as the white-spotted charr. This study suggests that a feasible range expansion of white-spotted charr resulted mainly from an invasion by migratory individuals during glacial periods, followed by their persistence in rivers during the interglacial periods. This process would have occurred repeatedly over time and space in relation to climate fluctuations.
Relationships among the four subspecies
Several genetic studies of Japanese freshwater fishes have demonstrated the relatively clear spatial patterns in population genetic structure (e.g., freshwater sculpin Cottus nozawae, Yokoyama and Goto, 2002; threespine stickleback Gasterosteus aculeatus, Higuchi and Goto, 1996, Yamada et al., 2001, Watanabe et al., 2003; medaka Oryzias latipes, Sakaizumi et al., 1983). Population genetic analysis of Gasterosteus aculeatus (Higuchi and Goto, 1996), for instance, revealed two salient genetically divergent groups that comprised populations of both the Pacific Ocean and the Sea of Japan. Our data, however, showed that the geographic distribution of each clade overlapped broadly; in addition, we found only limited congruence between the mtDNA phylogeographic patterns and the distributions of the four currently recognized subspecies. Moreover, the Templeton test revealed that mtDNA lineages of these four subspecies were not monophyletic. These results, in turn, may suggest that the current subspecies designations are not compatible with differentiations at the mtDNA level.
As mentioned above, the lack of clear geographic patterns in mtDNA among subspecies could be caused by both the large-scale dispersal events and by multiple colonizations of different lineages. Historical intergradations of mtDNA, compared to nuclear genes, might easily be achieved by secondary contact because the effective population size of mtDNA for salmonids is much smaller than that of nuclear genes (Hansen and Loeschecke, 1996; Laikre et al., 1998). Another possible explanation for the incongruence between mtDNA and morphological differentiation could be a rapid adaptive response to local environments. Given that the morphological differences among extant sub-species have originated recently (i.e., after the last glacial period), differences detectable with mtDNA sequences would not have accumulated sufficiently over a limited time span. Further detailed examination, e.g., using high resolution nuclear gene markers and large sample sizes at fine geographic scales, would be helpful for understanding the morphological divergence, genetic differentiation, and post-glacial distribution of the species.
Hierarchical analysis (AMOVA) showed that genetic variation was far more pronounced within subspecies than among subspecies (i.e., among discrete regions). Most white-spotted charr populations, particularly in Honshu Island, are substantially reproductively isolated from other conspecific populations and represented by having diagnostic haplotypes. This suggests that each population, rather than each subspecies, must be treated as an evolutionarily significant unit (Waples, 1995). In addition to genetic divergence, some ecological differences among populations, such as egg size (Morita, 2003), growth rate (Yamamoto et al., 1999), and size and age at the time of seaward migration (Yamamoto and Morita, 2002), have been reported for white-spotted charr. Many white-spotted charr populations, as well as other salmonids, have currently suffered from genetic degradation, excessive competition, overfishing, and introductions (Crisp, 2000). The importance of their characteristic population structure, therefore, must be adequately considered when establishing population units for conservation and management.
Acknowledgments
The authors appreciate following individuals for providing samples and technical advices: M. Fukushima, K. Goka, K. Ishikawa, N. Kamata, S. Kubota, H. Kubota, H. Kurachi, Y. Narimatsu, T. Nakamura, K. Oomura, A. Oomori, N. Okada, H. Sato, S. Suzuki, K. Takata, T. Takami, T. Tamate, K. Tamura, J. Tsuboi, G. Yamamoto, M. Yamamoto. This research was supported in part by grants from the Japan Society for the Promotion of Science (No. 13874104), Global Environmental Research Fund, Ministry of the Environment, Japan (F-3), the River Environment Fund in charge of the Foundation of River and Watershed Environment Management, Japan (FOREM), and a 21st Century COE Program on “Neo-Science of Natural History” at Hokkaido University financed by the Ministry of Education, Culture, Sports, Science and Technology, Japan.
REFERENCES
Appendices
Appendix.
Population abbreviation, current taxonomy, locality, Haplotype abbreviation (see Table 1), and sample size.
Notes
[1] Corresponding author Corresponding author: Tel. +81-268-22-0594; FAX. +81-268-22-0544. E-mail: ysho@fra.affrc.go.jp
