The genetic structure of populations is not necessarily reflected in the geographical proximity of individuals, because environmental gradients such as those of vegetation or climate can function as cryptic barriers to gene flow. We examined polymorphisms at nine microsatellite loci to determine and discuss whether a distinctive genetic structure was detectable in a spatially continuous population of the sika deer (Cervus nippon) on the Boso Peninsula of central Japan. Spatially explicit Bayesian analysis revealed that two genetically distinctive clusters exist in the Boso population. The spatial boundary of the two clusters approximately conformed to the border defined previously from a mitochondrial DNA dataset. By combining information on the geomorphic features surrounding the boundary and that on the lineage of 1970s population, we propose a schematic scenario for characterizing the population genetic structure to the present. The current population consists of genetically different lineages, and spatially discontinuous clusters have come into contact in the vicinity of a local road running along a steep-walled ravine that could act as principal barrier to gene flow. Biological factors such as distribution of vegetation and philopatric behavior might also have helped strengthen the cryptic genetic structure of the Boso population.
The use of molecular genetics to detect population structure is a well-established, practical alternative to direct methods such as mark-and-recapture or biotelemetry, and it has the added advantage of giving a temporal perspective to gene flow (Slatkin, 1985; Neigel, 1997). Wright's F statistics are used extensively for estimating the levels of gene flow and genetic dissimilarity among predefined patches or populations. However, there is not always much point in using such predetermined definitions of groups, because the genetic structure of a population is not necessarily reflected in the geographical proximity of individuals. Populations that are distributed in a spatially continuous way can be genetically structured along environmental gradients, such as climate or vegetation gradients, that function as cryptic barriers to gene flow (e.g., Zannèse et al., 2006).
Recent advances in the development of landscape genetic tools have enabled us to detect dissimilar individuals within spatially continuous populations. These tools were developed under the framework of Bayesian statistics, especially for approximations of the posterior probability of assigning individuals to their groups of origin; the updating of plural parameters is estimated in a deterministic order (Pritchard et al., 2000; Falush et al., 2003; Guillot et al., 2005). In addition to being used to find the number of genetically distinctive clusters, these approaches are particularly useful for mapping the borders of clusters and detecting migrants; they could therefore be powerful tools for monitoring the illegal translocation of nonnative individuals (Frantz et al., 2006) and for inferring the level of admixture in a genome (Corander and Marttinen, 2006).
Studies of genetic structure based on mitochondrial D-loop variation have demonstrated that the Boso population of the sika deer (Cervus nippon) consists of two major haplotypes (Nagata et al., 1999; Yoshio et al., 2008). Marked heterogeneity has been detected in the distribution of these haplotypes, and spatial discontinuity in the constitution of the haplotypes has been detected at the edge of a local motorway, Route 81; for details see Yoshio et al. (2008). Because Route 81 on the Boso peninsula passes along a steep-walled ravine, the natural landscape here may have functioned as a cryptic barrier to gene flow (Yoshio et al., 2008). Because the resolution of mtDNA polymorphism analysis is generally lower than that of other polymorphic marker analyses, such as AFLP or microsatellite analyses, more details of the genetic structure can be detected by using such other markers.
Our primary purpose was to dissect the genetic discontinuity by using microsatellites, which are codominant nuclear markers with higher resolution than mtDNA. If potential configurations or other environmental factors suggested by the analysis of mitochondrial D-loop variation were to impede migration and promote the genetic differentiation of underlying groups, then we expected that a similar pattern of genetic discontinuity would be detected by using other molecular markers. We conducted a spatially explicit Bayesian analysis using nine microsatellite markers to delineate the heterogeneity of the genetic structure of the Boso population.
MATERIALS AND METHODS
Specimens and DNA extraction
Whole blood, tail or ear tissue, or hair was obtained from 456 individual sika deer (260 females and 196 males) culled by hunters at Kimitsu, Otaki, Ichihara, Kamogawa, and Katsu-ura in Chiba Prefecture, central Japan, between January 2005 and September 2006. Details, including geographic information on capture points (latitude and longitude), date of collection, and sex of the deer, were individually recorded. Whole-blood samples were stored on FTA Mini Cards (Whatman, Middlesex, UK), and other tissues were stored in a freezer until the extraction of total DNA. DNA was extracted from the FTA Mini Cards in accordance with the manufacturer's specifications and from the tissues by using a DNeasy Tissue Kit (Qiagen, Hilden, Germany) in accordance with the manufacturer's protocol.
Amplification and genotyping of microsatellite DNA
Nine microsatellite loci (BOVIRBP, OarFCB193, BL42, BM203, BM888, BMC1009, Cervid14, ETH225, and BM4107) were used for genotyping (Table 1). These loci, apart from Cervid14, which is derived from the white-tailed deer (Odocoileus virginianus; DeWoody et al., 1995), are derived from cattle or sheep (Bishop et al., 1994; Crawford et al., 1995) and are mutually unlinked on a genetic map of red deer (Cervus elaphus), which is closely related to the sika deer (Slate et al., 2002). All of these loci have been revealed as polymorphic markers in the Kinkazan Island population of the sika deer (Okada and Tamate, 2000; Tamate et al., 2000).PCR was carried out in reaction volumes of 25 µl containing 2.5 µl of 10×NH4 buffer (Bioline, London, UK), 2.5 mM MgCl2, 0.2 mM each dNTP, 0.4 µM fluorescently labeled forward primer, 0.4 µM reverse primer (see Table 1), and 1.25 units of BIOTAQ DNA polymerase (Bioline). PCR cycling conditions consisted of an initial denaturation at 94°C for 1 min, followed by 30 cycles of denaturation at 94°C for 30 s, annealing at an appropriate temperature (see Table 1) for 30 s, and extension at 72°C for 30 s. After PCR, 0.1 µl of the product for each locus was mixed with 10 µl of deionized formamide and run on an ABI 3730 capillary sequencer (Applied Biosystems, Foster, CA, USA) with 0.2 µl of GeneScan-500 LIZ size standard (Applied Biosystems). Genotyping was performed with GeneMapper software, version 3.5 (Applied Biosystems).
The nine microsatellite loci used; annealing temperatures for amplification; allele sizes (in base pairs); numbers of alleles observed; observed and expected heterozygosities (H o and He, respectively); FIS (coefficient of inbreeding) values, and deviation from Hardy-Weinberg equilibrium (P).
a Significance levels for P value after Bonferroni correction: *, α<0.05; **, α<0.01; ***, α<0.001.
Genetic diversity and detection of genetic discontinuity in the Boso population
To document the extent of genetic diversity, the numbers of alleles and observed and expected heterozygosities (Ho and H e) were calculated for each locus. Significant deviations from linkage equilibrium between pairs of loci and from Hardy-Weinberg equilibrium at each locus were tested by the method of Raymond and Rousset (1995). The significance of P values was adjusted by a Bonferroni method (Rice, 1989). These calculations were carried out in GENEPOP, version 4.0 (Raymond and Rousset, 1995).
Before the spatial analysis of genetic structure, the number of genetically identified clusters (K) within the population was determined with a Bayesian approach implemented in the software STRUCTURE, version 2.2 (Pritchard et al., 2000). We chose the admixture model and the option of correlated allele frequencies, as this configuration is considered best in cases of subtle population structure (Falush et al., 2003). Similarly, we let the degree of admixture alpha be inferred from the genotyping data. In a pilot study, we found that a burn-in and Markov chain Monte Carlo (MCMC) length of 104 was sufficient. Longer burn-in or MCMC did not change the results significantly. As different runs could produce different likelihood values, even with much longer chains (e.g., 106; Evanno et al., 2005), five runs were performed to quantify the degree of variation of the likelihood for each K. From the results of the preliminary study, the range of possible K values tested was set at 1 to 4. We fitted a linear regression and second-order polynomial to the likelihood scores and compared Akaike's information criterion (AIC) (Akaike, 1987) scores to choose the best K.
After determination of the best K, the two-dimensional population genetic structure was investigated by using GENELAND software (Guillot et al., 2005); this method uses Bayesian statistics to partition individuals into groups that are as close as possible to Hardy-Weinberg and linkage equilibria (but see Frantz et al., 2006), in consideration of the spatial information explicitly associated with each individual. Partitions of spatially proximal individuals were given higher priority than partitions of individuals farther apart (Guillot et al., 2005). We ran five independent Monte Carlo chains with 106 iterations each (parameter values: maximum rate of Poisson process, 456; no coordinate noise; maximum number of nuclei in the Poisson-Voronoi tessellation, 1368; thinning, 10) (Guillot et al., 2005) and looked for consistency in the results. We adopted both the spatial Dirichlet model (D-model) and the Falush model (F-model) to determine whether the results varied with the priors for allele frequencies; for details see Guillot et al. (2005). Finally, to define the spatial boundary between different clusters, we first applied local polynomial regression fitting (LOESS) to dummy scores (e.g., “0” assigned to individuals of one group and “1” to those of the other when K=2) and then regarded the median of the highest and lowest predicted values from the fitting as the most appropriate boundary between two distinctive clusters. To show the levels of uncertainty of the boundary, we also estimated the standard errors of the predicted values. Data fitting was performed with the R statistical package, version 2.6.2 (R Development Core Team, 2008). The locations of all the samples belonging to each group were plotted, and artificial structures such as golf courses and major roads that could act as barriers to gene flow (e.g., Forman and Alexander, 1998) were overlaid on the same geographic information system (GIS) map using ArcGIS (version 9.1; ESRI Inc., Redlands, CA, USA).
Genetic variation in nine microsatellite markers
In 36 individual tests for linkage disequilibrium, four pairs of loci deviated significantly from equilibrium after Bonferroni adjustment (BOVIRBP and ETH225, P<0.001; OarFCB193 and BMC1009, P<0.005; BL42 and ETH225, P<0.005; Cervid14 and BM4107, P<0.005). Five of the six loci in the first three pairs are unlinked on the genetic map of red deer, which is closely related to the sika deer (Slate et al., 2002), and Cervid14 is known to be unlinked with any of the other loci that we examined here (Tamate et al., 2000). Therefore, all these loci were regarded as independent and were used in subsequent analyses.
The deviation from Hardy-Weinberg equilibrium was not significant at six loci, whereas a significant heterozygote excess was observed at OarFCB193 and significant heterozygote deficits at both BM203 and BM888, after Bonferroni adjustment (Table 1). These three loci were highly polymorphic (Table 1). The largest number of observed alleles was 13, on the OarFCB193 locus, and the smallest was two, on the BOVIRBP locus. Ho ranged from 0.06 at BOVIRBP to 0.71 at OarFCB193 and BMC1009; average observed and expected heterozygosities in the Boso population were 0.52±0.21 (mean±SD) and 0.53±0.20, respectively. In a previous study, H o and the number of alleles observed were 0.08 and 2, respectively, at the BOVIRBP locus and 0.39 and 4 at OarFCB193 (Nagata et al., 1998). The large discrepancy in Ho at the OarFCB193 locus might have resulted from a difference in sample sizes.
Spatial genetic structure of the Boso population
Five replicated runs using STRUCTURE showed that the mean posterior probability [LnP(D)] was highest when K=2 (Fig. 1). We found that the second-order polynomial fitted better (AIC=195 for the linear regression and AIC=169 for the second-order polynomial) and that the posterior probability was ∼1 when K=2. Along these lines of evidence, we concluded that K=2 was the best K. We then estimated the genetic population structure by using GENELAND (Fig. 2). The boundary estimated between identified clusters was situated at 140.11°E±0.0221° (median±SE) longitude; no distinctive latitudinal boundary was detected. The position of the boundary detected between the two clusters was fairly consistent, irrespective of the different priors used in the D-model and the F-model (data not shown). The distribution boundary of the two clusters is surrounded by Route 81 and two dirt forest roads (the Nishimine and Matsumori routes). Six of the nine loci contained private alleles, and most of these were found in individuals in the eastern area (see Fig. 2). Weir and Cockerham's θ, an estimator of genetic differentiation between the western and eastern individuals that was identified from the discontinuity line described above, was 0.017, deviated significantly from zero (P<0.001, 5×104 permutations), implying that a distinguishing genetic differentiation existed in the population.
All nine microsatellite loci tested were polymorphic in this isolated sika deer population on the Boso Peninsula, central Japan. Significant heterozygote excess and deficit were found at a single locus (OarFCB193) and at two loci (BM203 and BM888), respectively (Table 1). The deficits in heterozygosity might have been due to a bottleneck effect and subsequent inbreeding resulting from a severe reduction in population size in the past (Koganezawa et al., 1976; Chiba Prefecture and Deer Research Group on Boso, 1994).
There is serious concern worldwide about the ecological impact of an overabundance of deer in relation to direct economic losses in agriculture and forestry and the transmission of several animal and human diseases (Côté et al., 2004). From an ecological viewpoint, deer overabundance has indirect negative effects on the whole community structure and nutrient cycle (Rooney and Waller, 2003). On the Boso Peninsula, intensive studies of the impacts of overabundance of the sika deer have revealed a marked increase in damage to the richness of plant and arthropod communities with increasing population density of the deer (Miyashita et al., 2004; Suzuki et al., 2008). In the 19th century, sika deer were distributed throughout the Boso Peninsula (Chiba Prefecture and Deer Research Group on Boso, 2004). A severe population reduction was observed around the 1940s, probably as a result of human settlement of the foothills and indiscriminate hunting (Chiba Prefecture and Deer Research Group on Boso, 2004). In 1973–74, the lineage population became very small and was restricted to one area in the southern part of the peninsula (Koganezawa et al., 1976). In addition to the relaxation of hunting pressure when protection of the deer was instigated, extensive deforestation for conifer planting, which provided new feeding sites, was probably the impetus for the population resurgence that began in the late 1980s. This range expansion of the deer population has continued up to the present. The estimated number of deer had increased to around 3000 by 2001, and the distribution area changed from an estimated 40 km2 in 1973–74 to about 440 km2 in 2001 (Chiba Prefecture and Deer Research Group on Boso, 2004). Our principal purpose in this study was to examine whether a genetic population structure was detectable in this spatially continuous sika deer population on the Boso Peninsula, and, if so, to compare the boundary positions of clusters as inferred by both mtDNA and microsatellite analyses.
We found two meaningful clusters, in the western and eastern areas flanking Route 81. In a study of mitochondrial D-loop variation, spatial genetic discontinuity between two main haplotypes was detected in the area surrounding Route 81 (Yoshio et al., 2008); the results of this mtDNA analysis were in agreement with those of our present study, although the mtDNA analysis also suggested that a network of artificial features (two golf courses and two dams) might have contributed to the restriction of deer dispersal. Although further details of the genetic structure were not revealed by using the microsatellites, such incongruence in the results obtained by using the two different types of markers is not extraordinary, because mtDNA is generally inherited from the dam while microsatellite polymorphism arises from both the sire and the dam (e.g., Elmer et al., 2007; Ronka et al., 2008). The finding of an approximate concordance in the spatial patterns of genetic discontinuity, despite the use of disparate genetic markers, is intriguing in terms of dispersal patterns of males and females.
The sika deer might have been fragmentarily distributed in the southern area in the 1970s, because at this time field signs were found only on the margins of this area (Koganezawa et al., 1976). One conceivable explanation is that the present sika deer population on the Boso Peninsula consists of genetically different lineages, and that spatially discontinuous clusters have come into contact in the vicinity of Route 81, where a steep-walled ravine exists. Not Route 81 itself (because of its low volume of traffic), but two rivers and three ridges encompassing Route 81, are also potential candidate geographic barriers to gene flow. Recent studies of red deer have revealed that some landscape features, such as roads and slopes, can impede migration on a local spatial scale, supporting the above view (Pérez-Espona et al., 2008).
Because the mountains on the Boso Peninsula are relatively low and thus probably do not hamper migration, other potential factors, such as the distribution of suitable foods (Takada et al., 2002), are worth considering as facilitators of population structuring. In terms of plant communities, the landscape of the Boso Peninsula is a mixture of fine-scale elements, including broad-leaved forests, Japanese cedar forests, and agricultural fields (Miyashita et al., 2007). In fact, the forest edges distributed sporadically in the study area are generally known to act as large sources of food for the deer (e.g., Sinclair, 1997); furthermore, a wide area of agricultural fields in the southern part of the distribution, to the west of Route 81, supplies the deer with high-quality food (Takada et al., 2002; Miyashita et al., 2007). These spatially biased distributions of food tend to affect migration distance (Kareiva and Odell, 1987; Dwyer and Morris, 2006); this variation in migration distance would, in turn, promote the formation of a cryptic genetic structure in the population.
It is also worth mentioning the effect of social behavior in the sika deer. In general, sika deer are potentially capable of dispersing great distances. According to knowledge of the white-tailed deer and roe deer, for example, dominant females occupy a home range that is at the center of the social group, with younger individuals occupying home ranges that overlap and extend the periphery (Porter et al., 1991). Populations expand as an array of family units composed of females and showing strong philopatry (“the petals of a rose” metaphor; see Porter et al., 1991), while males tend to disperse further than females (Purdue et al., 2000; Wang and Schreiber, 2001). It remains an open question whether such social behavior also facilitates the formation of a global spatial genetic structure; such a characteristic of dispersal could be a factor in the acceleration of genetic structuring on a restricted spatial scale (e.g., Mathews and Porter, 1993; Purdue et al., 2000). In a future study we will need to examine the extent to which the factors mentioned above have contributed to the organization of spatial genetic structure in the Boso population.
We are grateful to Yuji Nagahata (Kamogawa City Hall), Ikuo Yamanaka (University of Tokyo), Takashi Ishihara and Ken Kumagawa (Deer Research Group on Boso), Makoto Kurokawa (Katsu-ura City Hall), and Yukihiro Takahashi (Kamogawa City) for their great help in collecting specimens and to Keiko Ikeda (National Institute for Environmental Studies) for technical assistance. Thanks are also due to all members of the Boso Sika Project and to K. Ryo Takahasi (National Institute of Genetics) and Hirohisa Kishino (University of Tokyo) for discussions. This study was supported by an Environmental Technology Development Grant (No. 60029) from the Ministry of the Environment of Japan.