Mammal Study 33: 59–69 (2008)© the Mammalogical Society of Japan

MSJMammalogical
Mammalthe1343-4152
Original
papersStudy
Society of Japan distribution of mtDNA haplotypes in a sika
Spatially
heterogeneousdeer (Cervus nippon) population on the Boso Peninsula, central Japan
Population genetic structure in sika deer

Masanobu Yoshio1, Masahiko Asada2, Keiji Ochiai2, Koichi Goka3, Kaori Murase4,
Tadashi Miyashita1 and Haruki Tatsuta3,*
1

Laboratory of Biodiversity Science, Graduate School of Agricultural and Life Sciences, The University of Tokyo, Japan
Natural History Museum and Institute, Chiba, Japan
3
Research Centre for Environmental Risk, National Institute for Environmental Studies, Tsukuba, Japan
4
Laboratory of Biometry and Bioinformatics, Graduate School of Agricultural and Life Sciences, The University of Tokyo,
Japan
2

69008
2520
33
Received;??
Copyright
2008
accepted13.7.2007
20.12.2007© 2008 MSJ

Abstract. We used variation in the mitochondrial DNA (mtDNA) control region (D-loop) to
examine the genetic structure of the sika deer (Cervus nippon) population on the Boso Peninsula,central Japan. A total of four haplotypes was found. In order to examine whether or not artificial
barriers such as roads, dams, and golf courses affect the spatial heterogeneity of mtDNA haplotypes,
we implemented two exclusive spatial analyses (SAMOVA and network analysis based on
Monmonier’s algorithm) for searching genetic discontinuities between artificial barriers. Prior to theanalyses, the whole distribution area was divided into meaningful eight blocks. Analysis of molecularvariance (AMOVA) detected significant spatial heterogeneity in the constitution of the haplotypes
among the blocks. The subsequent spatial analyses detected some significant spatial discontinuities
on borders of the blocks. In particular, the largest discontinuity was observed in the area including
motorway Line 81, but the traffic density of Line 81 is generally not very heavy compared to othermajor roads. These findings suggest that roads could be one of major barriers to hamper migration ofsika deer to some extent, but other potential factors such as the location of food resources and/or thehistory of bottleneck event are also likely to more or less contribute to configure the present patternsof haplotype distribution.
Key words: artificial structure, Cervus nippon, gene flow, mitochondrial DNA, population.

Estimation of the patterns of gene flow is one of the cenIntroduction
tral problems in the field of population genetics, because
gene flow determines the extent to which local populations of a species form independent evolutionary units(Slatkin 1985). Two following approaches are possible
to estimate the amount of gene flow. One of them
directly estimates dispersal distances of individuals, forinstance, with the aid of the Global Positioning System(GPS) or by means of traditional mark-and-recapture
methods. The other is founded on a mathematical modelof the interaction of gene flow and other forces (e.g. genetic drift) to predict how much gene flow must havebeen occurring in order for the patterns observed in thedata to be present (Slatkin 1994). Although these two
approaches have different drawbacks for estimating gene
*To whom correspondence should be addressed. E-mail: htatsuta@nies.go.jp

flow, the latter is greatly advantageous because it incorporates the effects of various kinds of dispersal andaverages variations in dispersal over time. Our principal
aims are to research genetic population structure in awild animal using the latter approach.
Movement of individuals is often restricted by physical barriers in wild animals (Avise et al. 1987). Geographic gaps such as rivers, valleys, and mountains havebeen known to contribute to form genetic structure
within or between populations (Avise et al. 1987).
Besides these natural barriers, potential artificial barriers
to gene flow are roads (Forman and Alexander 1998;
Aars and Ims 1999; Mech and Hallet 2001; Riley et al.
2006), dams (Hansen and Loeschcke 1996; Yamamoto etal. 2004), and urbanized areas (Wang and Schreiber

602001). Thus, environmental modification mediated byhuman activities is likely to result in habitat fragmentation in wildlife (Ceballos and Ehrlich 2002) with a consequent reduction of genetic diversity and of the levels ofgene flow (Young et al. 1996). Understanding the practical effects of such artificial structures on gene flow of
wildlife is important in point of detecting of spatial
genetic structure in landscape genetics (Manel et al.
2003).
The population density of sika deer (Cervus nippon)on the Boso Peninsula in Chiba Prefecture, central Japan,
is known to have been relatively high in the latter half ofthe 19th century, but the prompt population shrinkage
occurred in about the 1940s probably because of humansettlement of the foothills and hunting (Chiba Prefecture
and Deer Research Group on Boso 2004). To restore
population size, culling of deer was suspended in 1961.
By 1973–74, when the density was at its lowest, a relict
deer population was distributed in a local area of thesouthern peninsula (Koganezawa et al. 1976; Chiba
Prefecture and Deer Research Group on Boso 2004)(Fig. 1a). For reasons not entirely known, a rapid resurgence of numbers began in the late 1980s. The present
distribution is consequently much wider than ever previously recorded, and its density is also conspicuously highespecially in the southern part of its distribution (Chiba
Prefecture and Deer Research Group on Boso 2004,
Fig. 1). The southwestern limit of the current distribution of sika deer on the Boso Peninsula corresponds wellto the motorway Line 34 (see Fig. 1), implying that thisroad acts as an effective barrier against distribution
expansion to the southward. Geographical barriers are ofimportance as well. A very high mountainous region
and deep gorge are absent, but some roads pass along
the brink, which could affect dispersion of sika deer.
Furthermore, the amount of food resources could cause
non-random dispersal in wild animals (Dwyer andMorris 2006), and thus could be another potential factor
for weakening random dispersal. However, the precise
estimate of the area covered with suitable food plants isgenerally awkward and the survey of plant species
over the area is beyond our scope in the present study.
Therefore, here we focused on the relationship between
artificial structures and population genetic structure,
together with geographical features where artificial structures exist.
In the present study, we assessed the genetic variationand spatial genetic structure in the Boso population using
mtDNA D-loop sequences which have been frequently

Mammal Study

33 (2008)

used to assess such purposes (e.g. Nagata et al. 1998,
1999; Nabata et al. 2004, 2007; Yamada et al. 2006;
Yuasa et al. 2007). Based on the result, we compared theamount of variation in D-loop region in other populations (Yuasa et al. 2007) and discussed effects of staple
artificial structures and other potential factors on geneflow among subdivided populations.

Materials and methods
Specimens and DNA extraction
To date, local governments have started programs tocontrol the population density of sika deer. The present
study was conducted under the operation of the programs. Whole blood, tail or ear tissue, or hair was
obtained from 259 individuals of sika deer culled byhunters at Kimitsu, Ohtaki, Ichihara, Kamogawa, andKatsu-ura, in Chiba Prefecture, central Japan, between
January and October, 2005 (Fig. 1 and Table 1). Details
including geographic information (latitude and longitude) and date of collection of the deer were individually
recorded for each site. Whole-blood samples werestored on FTA Mini Cards (Whatman, Middlesex, UK),and other tissues were stored in a freezer until extractionof 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 sequencing
To amplify the mitochondrial D-loop region by polymerase chain reaction (PCR), we used the primer set
L15926 and H597 (Nagata et al. 1998). Amplification of
the DNA extracted from the FTA Mini Card was carriedout in reaction solution (50 µL in a tube) containing 5 µL
of 10× PCR Buffer II (Applied Biosystems, CA, USA),2.5 mM MgCl2, 0.2 mM each dNTP, 0.4 µM eachprimer, and 1.5 units of AmpliTaq DNA polymerase(Applied Biosystems). Amplification from each tissue
extract was done in 25 µL of the above solution. PCR
was performed under the following conditions: 2 min at
94°C; 40 cycles at 94°C for 30 sec, 55°C for 30 sec, and72°C for 1 min; and a final extension at 72°C for 6 min.
We also used previously developed primers HD2, HD6(Nagata et al. 1998), and original primers, HD2R,
HD6R, and HDN for sequencing. HD2R (5'-CATCTGGTTCTTTTTTCAGG-3') and HD6R (5'-GCAGTCAATGGTCACAGGAC-3') are complementary to HD2 andYoshio et al., Population genetic structure in sika deer61

Fig. 1. (a) Changes in the distribution of sika deer, Cervus nippon, on the Boso Peninsula, central Japan. Dots show locations of observed fieldsigns of sika deer (e.g. footprints, food remains, droppings, wallows) observed in 1973–1974 (Koganezawa et al. 1976). The distributions in
1973–1974 (dashed line) and in 2001 (shaded area) followed the study of the Chiba Prefecture and Deer Research Group on Boso (2004).
(b) Map of the study area on the Boso Peninsula and schematics of boundaries of municipals (bold dashed line) and management units (dashedline), artificial structures such as main roads (solid line), dams or artificial lakes (shaded area), and golf courses (letter “G” in circle). Numbers
in reverse-triangle and in hexagon indicate route (National motorway) and line (Prefecture mortorway), respectively. Eight areas enclosed bythese structures (Block I–VIII) were determined prior to the analysis of spatial genetic structure. Management units were overlaid in this figure(see also Tables 1 and 3).

HD6, respectively. The nucleotide sequence of HDN is
5'-CATAATGTATCTTACGCACCGG-3'. Before the sequencing reaction, the PCR product was purified with a
QIAquick PCR Purification Kit (Qiagen). The purified
PCR product was directly sequenced with a BigDyeTerminator Cycle Sequencing Kit (ver. 3.1, Applied
Biosystems), and analyzed on an automated DNA
sequencer, model 3730 (Applied Biosystems).

Nucleotide and haplotype diversity
All sequences were edited and aligned by CodonCode
Aligner (CodonCode Corp). Definition and alignment oftandem repeat units conformed to the methods of Nagata
et al. (1999). Haplotype diversity (h) and nucleotide
diversity (π) were calculated using DnaSP ver. 4.10.4(Rozas et al. 2003), excluding the gap region and seventh
tandem repeat unit which cause interference in the calculations of the indices (Table 2).

62Mammal Study

Table 1.

33 (2008)

Collecting sites, number of sika deers collected from January to October, 2005 (N), and haplotype frequency at each collecting site
Location

Municipal

Management Unit*

Kamogoawa

A2
A3
A4
A5
G1
G2
G3
G4
G7
T1
T2
T3
T4
T5
T6
T8

Haplotype**

N

1a

1b

2a

2b

51
9
13
5
13
24
18
1
4
5
5
19
4
4
6
6

38
8
12
4
2
17
6
0
3
4
4
14
3
3
6
3

0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

11
1
1
1
10
7
12
1
1
1
1
5
1
1
0
3

2
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0

T9

5

5

0

0

0

Otaki

O2
O3
O4
O5
O6
O7

13
11
7
2
4
12

13
8
6
2
4
12

0
0
1
0
0
0

0
3
0
0
0
0

0
0
0
0
0
0

Ichihara

I1

15

14

0

1

0

Katsu-ura

U3

3

2

0

1

0

Kimitsu* Names followed Chiba Prefecture and Deer Research Group on Boso (2004).
** Definition of each haplotype is shown in Table 2.

Table 2. Comparison of nucleotide sequences (variable positions are shown) and numbers of individuals sharing each haplotype (N) in the Boso
population of Cervus nippon
Nucleotide position
Haplotype

N

1
9
7

2
6
3

2
9
2

3
0
2

3
0
7

1a
1b
2a
2b

193
1
62
3

A
G
·
·

C
·
T
T

C
·
T
T

T
·
C
C

C
·
T
T

3
5
9

–
–
–
*
*

3
9
7

4
8
1

7
5
9

8
1
0

8
2
4

8
5
7

9
5
8

9
7
1

A
·
G
G

A
·
G
G

T
C
·
·

T
·
C
C

T
C
·
·

C
·
·
–

–
·
T
T

Dots indicate nucleotide identity with haplotype 1a, and bars indicate lack of the tandem repeat unit (359–397 bp) or single nucleotide.
Asterisks between 359 and 397 bp show the tandem repeat unit named d8 by Nagata et al. (1999).

Spatial genetic structure of the Boso population
For the purpose of assessing the spatial patterns of
mtDNA haplotypes, we first conducted local interpolation using inverse distance weighting (IDW: Fortin andDale 2005) in order to estimate the intensity of aggregation of haplotypes. A total of four haplotypes was foundin consequence of DNA sequencing, but they consistedof two major haplotypes that differ in the number oftandem repeat unit (see details in Results). Therefore,
the interpolation was performed for the major two haplotypes. The common form of the inverse distance weightisYoshio et al., Population genetic structure in sika deerm

63tion subdivision (Table 1, Fig. 1b). To assess whether or
not staple artificial structures play an important role torestrict gene flow, we divided the population into eightblocks (Block I–VIII) enclosed by main roads, dams orlakes, and golf courses, covered with all sampling points(Fig. 1b, Table 3). Two management units, G7 and T8,
which contained only four and six samples, respectively,
were excluded from subsequent analysis because they
were not included into any of the eight blocks. Based onthese blocks, we examined the relationship between theexistence of artificial structures enumerated above andspatial genetic structure of the sika population. For comparison of the amount of traffic on major roads in study
area, mean number of cars that have passed through eachroad for several years was calculated (Chiba Prefectural
Government 1989, 2000).
Wright’s F statistics, in which the inbreeding coefficient in a subdivided population (FT) can be decomposed
into components due to nonrandom mating within local
populations (FIS) and due to population subdivision(FST), is frequently used (Wright 1921) for assessing thedegree of population substructure and the level of geneflow among (sub)populations. We calculated θ (Weirand Cockerham 1984) as FST estimator, in which thenumbers of samples and locations were entered explicitly, between different blocks described above. θ withunequal sample size was calculated by Arlequin (ver.
3.01, Excoffier et al. 2005).
Two exclusive approaches were used to assess thespatial population structure. First, we conducted Spatial
Analysis of Molecular Variance (SAMOVA) whichallows us to detect groups of (sub)populations that are
genetically and geographically homogeneous using simulated annealing procedure (Dupanloup et al. 2002). In–pz ( x j )d ij∑
j= 1ẑ ( x 0 ) = ---------------------------m
–p
d
∑ ij
j= 1

where ẑ (x0) is an interpolated value of variable z for agiven point x0, z(xj) is the value at the sampled location
j, m is the number of neighboring sampling locations
based on some definition such as being within a searchradius, p is a real value, and dij is the distance between
unsampled location i and sampled location j. In thepresent study, the IDW interpolation was conducted with
p = 2 and the number of sampled location for interpolation is 12.
To find an area in which the constitution of the major
two haplotypes was changed large, local regression
model (Simonoff 1996) was used for fitting data given
by dummy variables (i.e., one major haplotype = 0,
another haplotype = 1). We also estimated Ripley’s Kfunction (Ripley 1981) that was estimated using a
Spatstat package (Baddeley and Turner 2005) in statistical program R, version 2.5. The statistical significanceof K(r) values for given radii, r, was estimated by 95%
confidence intervals generated by 1000 times Monte
Carlo simulations of complete spatial randomness (CSR)
which is implemented as a homogeneous Poisson process (Diggle 1983). Circle radii for the K(r) analysis
increased by 50-m increments up to the maximum distance of 40 km.
Chiba Prefecture Government divided 13 municipalities that were located in and around the distribution ofdeer into 71 units for the deer management (Chiba Prefecture and Deer Research Group on Boso 2004). We
adopted those management units for determining popula-

Table 3. Eight blocks defined by artificial structures, management units containing in each block, the numbers of samples (N),
haplotype frequency, and geographic coordinates of the center of the block
Block
I
II
III
IV
V
VI
VII
VIII

Block composition(Management unit)*
T5, T6
G3, G6
T2, T3, T4, G5
G1, G2, G4
T9, I1
T1, O1, O3, O5,
A2, A3, A4, A5
O2, O4, U2
O6, O7, U3

Haplotype frequency
1a

1b

2a

2b

Geographic coordinates(longitude (E), latitude (N))

10
18
28
38
20
96

9
6
21
19
19
76

0
0
0
0
0
0

1
12
7
18
1
18

0
0
0
1
0
2

(140°01'57'', 35°11'36'')(140°02'47'', 35°08'21'')(140°06'12'', 35°12'04'')(140°06'23'', 35°08'42'')(140°06'59'', 35°15'53'')(140°09'50'', 35°10'50'')

20
19

19
18

1
0

0
1

0
0

(140°11'26'', 35°13'22'')(140°14'14'', 35°12'14'')N* Names of management unit follow Chiba Prefecture and Deer Research Group on Boso (2004).

64

hierarchical AMOVA, FCT (proportion of genetic variation among groups of populations to the total variation)shows the extent of differentiation between groups(Excoffier et al. 1992). SAMOVA calculates the largest
FCT value corresponding to arbitrary number of groups(K) (Dupanloup et al. 2002). We ran SAMOVA program
developed by Dupanloup et al. (2002) with 100-times
simulated annealing process for various number ofgroups (K = 2 to 5) and compared simulated FCT values
for determining the best grouping of blocks. Anotherapproach for assessing population structure was based on
Monmonier’s algorithm which is used to identify the discontinuity. Delaunay triangulation networks between
blocks were constructed using geographic coordinates,and pairwise genetic distance was calculated between
adjacent blocks. Monmonier’s algorithm firstly selected
the edge with the largest distance on the triangulation as
the starting edge to find the discontinuity and thenextended it across the adjacent edge associated with thelargest distance. This procedure was continued until thebarrier reached the limit of the triangulation or formed aloop by itself or met another barrier (Manni et al. 2004).
Pairwise θ between blocks was used as genetic distance,and genetic discontinuities were detected using program
BARRIER 2.2 (Manni et al. 2004). Both methods
require geographic data, namely latitude and longitude,of blocks and sometimes result in somewhat similar
results, while algorithms for finding discontinuities aredifferent (Dupanloup et al. 2002). To determine geographic coordinates of each block, both latitude andlongitude of north, south, east, and west extremity of theblock were found on a 1:25000 scale map, and line connecting north and south extremities and line connecting
east and west ones were each drawn. The point where
north-south and east-west line intersect was defined as
the center of the block in this study (Table 3).

Results
The mitochondrial D-loop region was successfully
amplified from all samples, and nucleotide sequenceswere determined. The length of the compared region
varied from 1051 to 1091 bp. Comparison of all thesequences obtained revealed a total of four haplotypes in
the Boso population (Table 2). The sequences of these
haplotypes have been deposited in EMBL/GenBank/DDBJ (Accession Nos. AB247654–AB247657). Transitions were found at 10 nucleotide sites and gaps at twosites. According to Nagata et al. (1999), the D-loop

Mammal Study

33 (2008)

region of sika deer in Japan possesses tandem repeat
units of 37 to 40 bp, and samples from Chiba containeither six or seven units. Our haplotypes 1a and 1b contained six units, whereas haplotypes 2a and 2b contained
seven (Table 2). Comparison of our results with thoseof Nagata et al. (1999), in which three haplotypes
from Chiba Prefecture, Ama1, Kmo1, and Kmo2, were
described, revealed that haplotype Ama1 was identical tohaplotype 1a, but we did not find haplotype Kmo1. It
was unclear whether haplotype Kmo2 corresponded toeither haplotype 2a or 2b, because the partial sequences
consisting of 169 bp in the left domain of the control
region near tRNA-Phe were not recorded by Nagata et al.
(1999). Haplotype 1b was newly observed in our study.
Haplotype 1a was mostly dominant throughout the Boso
population (74.5% of the total individuals) (Table 1).
The percentage nucleotide difference was 0.29%
between haplotypes 1a and 1b (3/1050 bp), 0.67%
between 1a and 2a/2b (7/1050 bp), and 0.95% between
1b and 2a/2b (10/1050 bp). The haplotype diversity (h)and nucleotide diversity (π) were 0.38 and 0.0025,
respectively.
The distribution of each haplotype demonstrated thatthe constitution of the two dominant haplotypes, 1aand 2a, varied among subpopulations: haplotype 1a was
widely observed throughout the localities, whereas the
majority of haplotype 2a was localized in the south andsouthwest (Fig. 2). Haplotype 1b was observed only in a
management unit “O4”, where the sika deer population
had not been confirmed to be distributed in the 1970s,
and haplotype 2b was found only in “A2” and “G1”
(Table 1). The null hypothesis of complete spatial randomness (CSR) was evaluated for the two major groupsof haplotypes consisting of haplotype 1 (haplotype 1a +1b, see Fig. 3) and haplotype 2 (haplotype 2a + 2b) using
the Ripley’s K-function. CSR was statistically rejectedin the both groups (P < 0.05). Interpolation by local
regression method revealed that major two haplotype
constitution (haplotype 1 and 2 described above)changed largely in an area surrounding motorway Line81 (indicated by bold arrow in Fig. 3) on which a boundary of block IV and VI exists (Figs. 1b and 3).
Analysis of molecular variance (AMOVA, Excoffieret al. 1992) demonstrated that the global FST value overall blocks (θ = 0.184) was significantly greater than zero(P < 0.0001, based on 10 000 permutations), implyingthat distinctive heterogeneity of the haplotype distributions (i.e. population structure) exists among blocks.
SAMOVA demonstrated that FCT value was larger forYoshio et al., Population genetic structure in sika deer

65

Fig. 2. Distributions of haplotypes on the Boso Peninsula. (a) All pooled haplotypes; (b) haplotype 1a; (c) haplotype 1b; (d) haplotype 2a;(e) haplotype 2b. Dots indicate sampling points where individuals were collected. Some dots overlap because many deer were collected at thesame geographic points. Solid lines show the network of main roads at study area (see text for details).

smaller number of groups (K) and that it was the highest
when all samples were divided into two groups, block II
and the others (FCT = 0.320, Table 4). A significant FCT
was also calculated when blocks II and IV were separeted from the others (K = 3, FCT = 0.299, Table 4). Thiscan be interpreted that based on a locus on mtDNA, subdivided populations existing in southwestern part of thewhole distribution were genetically distinct to the othersin the Boso population. Genetic discontinuities among
blocks connected by Delaunay triangulation were delineated using Monmonier’s algorithm (Fig. 4), a result
showed similar tendency to the consequence obtainedby SAMOVA: a main discontinuity was recognized
between a pair of block II and IV and the others, andthird discontinuity was drawn between block II and IV,
while the model of three groups was marginally significant (Table 4). Note that results of these methods were
not completely congruous: the second discontinuity
obtained by latter analytical method did not correspondwith the result of SAMOVA. Such an incongruity would
be attributed to the difference in searching algorithm forgenetically heterogeneous groups (see Dupanloup et al.
2002 for detailed explanation). In conclusion, theseresults support the view that some of artificial structuresconsidered in the present study have tended to constrict
movement of sika deer and also suggest that primaryfactors that restrain migration exist between block II, IV
and the others.

Discussion
Nucleotide and haplotype diversity
One of the main purposes we intended to analyze wasto show the level of genetic variation in the mitochondrial DNA in the Boso population. The percentage difference in nucleotide sequences between haplotypes inthe D-loop region varied from 0.29% to 0.95%. This
gives 0.13–0.43 Myr have passed since haplotypes havediverged, based on the knowledge that a substitution ratein the mitochondrial D-loop is 1.11–1.13% per Myr in
Cervinae (Randi et al. 2001). Two diversity indices,
haplotype diversity (h = 0.38) and nucleotide diversity

66

Mammal Study

33 (2008)

Fig. 3. Distribution of pooled haplotypes (1a + 1b: green circles; 2a + 2b, red triangles). Boundaries between cities (solid black lines), roads thatare larger than 5 m width (bold black lines), golf course (yellow squares) are overlaid, respectively. Large bold arrow indicates Line 81 (seeMaterials and methods and Discussion for details). Local interpolation based on inverse distance weighting (IDW, see Materials and Methods fordetails) is shown in the area of study (probability of existence of haplotype 2a + 2b is shown by black contours). Bottom and right figures indicate
scatterplots of the distribution of pooled haplotypes represented by dummy scores (the score of haplotype 1a and 1b is set as “1” and that of
haplotype 2a and 2b as “0”). The unit of latitude and longitude represents projected coordinate (distance unit: meters) from center (i.e. projecteddistance from the location of N 36°0'0'' E 139°50'0''). These plots were interpolated by local regression model (see Materials and methods). The
regression lines are overlaid on the scatterplots and the area where a large increment of regression coefficients was observed is shown by blue
square.

(π = 0.0025) of the Boso population were included in therange of both diversities in other neighboring populationof southern Kanto-Izu district, Honshu (0.063–0.533 in
haplotype diversity, 0.0001–0.0063 in nucleotide diversity, see Yuasa et al. 2007). Grant and Bowen (1998)
proposed to infer the population history using thesediversity indices in marine fishes. According to Grant
and Bowen (1998) criteria, populations with both low
haplotype diversity (h < 0.5) and nucleotide diversity(π < 0.5) can be interpreted that they recently experienced severe bottlenecks. Since many Japanese populations of sika deer including the Boso one haveexperienced recovery from severe bottleneck, our estimations could be interpreted as a consequence of such ahistory. This inference is not discrepant to previous
reports on the Boso population (Koganezawa et al. 1976;
Chiba Prefecture and Deer Research Group on Boso2004).

Possible effects of artificial structures and other factors
on population genetic structure
While only four haplotypes were detected and the distribution area of sika deer is not very wide, marked heterogeneity of haplotype distribution was detected in the
Boso population. In particular, frequency of haplotype
2 (2a + 2b) was relatively high in the southern and thesouthwestern areas of the whole distribution, whereas
haplotype 1 (1a + 1b) was found over the entire area
(Figs. 2, 3). From the knowledge of substitution rate inmitochondrial D-loop region (Randi et al. 2001), thepossibility that haplotype 2 occurred in the last two or
three decades and has not still spread over the area is
not promising. Therefore, based on the current data, it
is most likely that both major haplotypes existed in thelast bottleneck event and the distribution patterns ofthese haplotypes have gradually varied since the population expansion has started. The mechanism inducing

Yoshio et al., Population genetic structure in sika deer

67

Table 4. Fixation indices corrsponding to the different number of groups (K) of blocks inferred by SAMOVA for the sika deer
population on the Boso Peninsula
Number of groups (K)Group composition

Two groups

1. Blocks I + III + IV +V +VI +VII + VIII
2. Block II
1. Blocks I + III +V +VI + VII + VIII
2. Block II
3. Block IV

Three groups

FSC

FST

FCT

0.119***

0.401***

0.320***

0.021***

0.314***

0.299**P<0.05, **P<0.01, ***P<0.001 (based on 1000 permutations).

Fig. 4. Delaunay triangulation (solid lines) based on geographic
coordinates of the center of the blocks (plots) and three main barriers(bold lines) deduced by pairwise θ between blocks based on
Monmonier’s algorithm.

Fig. 5. Mean number of cars that have passed through each road per
half a day (12 h: 7:00–19:00) with ±1 standard errors (SE) (data from
reports of Chiba Prefectural Government 1989, 2000). Data are basedon the survey for 35 years (20 years for Load 178 and 465).

heterogeneity of the frequency constitution of haplotypes
among subpopulations is a very intriguing issue.
Since artificial structures including roads and dams are
known to hamper animal movement across habitats andalter landscape and spatial patterns, the network of thesestructures is considered to be important to know thepatterns of spatial genetic structure of wildlife. In the
Boso Peninsula, the motorway Line 34 seems to act asan effective barrier against distribution expansion to thesouthward (see Fig. 1). Indeed, the average number ofcars that have passed through the Line 34 was extremely
large among main roads in the study area (Fig. 5). In thisstudy, eight blocks surrounded by roads, dams including
artificial lakes, and golf courses, were treated units to beanalyzed to assess the effect of artificial structures ongene flow (Table 3, Fig. 1b). AMOVA and spatial
analyses detected distinctive genetic structures of sika
deer population on the Boso Peninsula. In particular,
spatial analyses detected the greatest heterogeneity inthe constitution of haplotype frequencies between a pairof block II and IV and the others (Table 4, Fig. 4).
Moreover, interpolation and local regression approachrevealed the area in which haplotype constitution wasconspicuously changed existed in an area surroundingLine 81 (Fig. 3). These concordant results suggest that
the Line 81 is a potential candidate that has left sika
deer in the block II and IV apart from other area. Road
width and vehicle traffic levels could be major determinants of the barrier effect (Forman and Alexander 1998),
but width of the Line 81 is relatively narrow (i.e., theminimum width of road is shortest among comparedroads) and traffic is generally not so heavy comparedwith other roads (Fig. 5). A reasonable explanation forthe discontinuity is that the Line 81 is running througha steep-walled ravine that may practically function as ageographic barrier. It is worth studying whether or nota similar genetic discontinuity is also found in the area
in other animal species.
Our results also suggest that existence of golf courses
on management units G2 and T5 has promoted heteroge-

68

neity of spatial distribution of mtDNA. The constructionof these two golf courses has been started in early 1970s,
approximately corresponding to the beginning of range
expansion of deer population on the Boso Peninsula (seeFig. 1). Another important factor that could affect dispersal patterns in animals is the abundance of resources.
When resources are abundant, consumers tend to disperse shorter distances, whereas the dispersal distancewould be longer when local resource availability is low(Kareiva and Odell 1987; Dwyer and Morris 2006).
These forecast that deer that have migrated from an area
in which food resources are relatively poor would staylonger in areas where the food resources are relatively
abundant, leading to heterogeneous dispersal. These
resource-dependent movements could consequently
affect estimations of migration rate between local habitats. Indeed, in the southern parts of the blocks II and IV,
there exists a wide area of agricultural fields, which is
known to provide high-quality food to deer in this region(Takada et al. 2002; Miyashita et al. 2007). Anotherdetailed study on the distribution of food resource is ofcourse necessary, but the evidence of the agricultural
fields supports the theoretical prediction that individuals
tend to stay in area where food resource is constantly
supplied and as a result, heterogeneous spatial structureof haplotypes arises due to the promotion of local dispersal. Inversely, genetic differentiation will not be
facilitated in such area as the amount of food resource ispoor, which would induce a high frequency of migration.
Not only artificial structures but food resource in block II
and IV may have been enhanced the difference in genetic
constitution in these from other area.
History of population shrinkage in the studied areacannot also be ignored. It is very likely that present
genetic structure is partly derived from the pattern ofrefugia during 1970s. In fact, previous study on a relict
sika population did not take account of the number ofrefugia nor genetic structure (Koganezawa et al. 1976).
The mitochondria are maternally inherited except someevidence (e.g. Ujvari et al. 2007) and therefore, thepresent distributions of haplotypes can be viewed as apattern of past dispersal in females. For the detailed
study of inferring population structure during bottleneck
event, coalescent simulation studies using information of
nuclear markers such as microsatellites are inevitable. It
is important to understand what the extent of various
factors listed above have affected the formation of population genetic structure of sika deer in the Boso peninsulain the subsequent studies.

Mammal Study

33 (2008)

Acknowledgements: This study was supported by theEnvironmental Technology Development Fund (No.
60029) supplied by the Ministry of the Environment inJapan. We are grateful to Yuji Nagahata (Kamogawa
City Hall), Ikuo Yamanaka (University of Tokyo),Takashi Ishihara and Ken Kumagawa (Deer ResearchGroup on Boso), and Yukihiro Takahashi (Kamogawa
City) for their great help in collecting specimens, andto Keiko Ikeda (National Institute for Environmental
Studies) for technical assistance. Thanks are also due
to all members of the Boso sika project for their discussions and encouragement.

References
Aars, J. and Ims, A. R. 1999. The effect of habitat corridors on ratesof transfer and interbreeding between vole demes. Ecology 80:
1648–1655.
Avise, J. C., Arnold, J., Ball, R. M., Bermingham, E., Lamb, T.,
Neigel, J. E., Reeb, C. A. and Saunders, C. 1987. Intraspecificphylogeography: The mitochondrial DNA bridge between population genetics and systematics. Annual Review of Ecology and
Systematics 18: 489–522.
Baddeley, A. and Turner, R. 2005. Spatstat: an R package for analyzing spatial point patterns. Journal of Statistical Software 12:
1–42.
Ceballos, G. and Ehrlich, P. R. 2002. Mammal population losses andthe extinction crisis. Science 296: 904–907.
Chiba Prefecture and Deer Research Group on Boso 2004. Integrated
Science Reports on the Management of Sika Deer on Boso
Peninsula, Chiba Prefecture. Chiba Prefecture, Chiba, 134 pp.
(in Japanese).
Chiba Prefectural Government 1989. Report of Road Traffic Census.
Chiba Prefecture, Chiba (in Japanese).
Chiba Prefectural Government 2000. A New Report of Road Traffic
Census. Chiba Prefecture, Chiba, 1692 pp. (in Japanese).
Diggle, P. J. 1983. Statistical Analysis of Spatial Point Patterns.
Academic Press, New York, 288 pp.
Dupanloup, I., Schneider, S. and Excoffier, L. 2002. A simulated
annealing approach to define the genetic structure of populations.
Molecular Ecology 11: 2571–2581.
Dwyer, G. and Morris, W. F. 2006. Resource-dependent dispersal andthe speed of biological invasions. American Naturalist 167: 165–
176.
Excoffier, L., Smouse, P. and Quattro, J. 1992. Analysis of molecularvariance inferred from metric distances among DNA haplotypes:
application to human mitochondrial DNA restriction data.
Genetics 131: 479–491.
Excoffier, L., Laval, G. and Schneider, S. 2005. Arlequin ver. 3.0: an
integrated software package for population genetics data analysis.
Evolutionary Bioinformatics Online 1: 47–50.
Forman, R. T. and Alexander, L. E. 1998. Roads and their major
ecological effects. Annual Review of Ecology and Systematics
29: 207–231.
Fortin, M-J. and Dale, M. 2005. Spatial Analysis: A Guide for
Ecologist. Cambridge University Press, Cambridge, 365 pp.
Grant, W. S. and Bowen, B. W. 1998. Shallow population histories indeep evolutionary lineages of marine fishes: insighs from sardines
and anchoviesand lessons for concervation. The Journal of

Yoshio et al., Population genetic structure in sika deer
Heredity 89: 415–426.
Hansen, M. M. and Loeschcke, V. 1996. Genetic differentiation
among Danish brown trout populations, as detected by RFLP
analysis of PCR amplified mitochondrial DNA segments. Journalof Fish Biology 48: 422–436.
Kareiva, P. M. and Odell, G. 1987. Swarms of predators exhibit“prey taxis” if individual predators use area-restricted search.
American Naturalist 130: 233–270.
Koganezawa, M., Katai, N. and Maruyama, N. 1976. Distribution ofsika (Cervus nippon centralis Kishida) in the eastern part of the
Boso Mountains, Chiba Prefecture. The Nihonzaru 2: 115–121(in Japanese).
Manel, S., Schwartz, M. K., Luikart, G. and Taberlet, P. 2003. Landscape genetics: combining landscape ecology and population
genetics. Trends in Ecology and Evolution 18: 189–197.
Manni, F., Guérard, E. and Heyer, E. 2004. Geographic patterns of(genetic, morphologic, linguistic) variation: how barriers can bedetected by using Monmonier’s algorithm. Human Biology 76:
173–190.
Mech, S. G. and Hallet, J. G. 2001. Evaluating the effectiveness ofcorridors: a genetic approach. Conservation Biology 15: 467–
474.
Miyashita, T., Suzuki, M., Takada, M., Fujita, G., Ochiai, K. andAsada, M. 2007. Landscape structure affects food quality of sika
deer (Cervus nippon) evidenced by fecal nitrogen levels. Population Ecology 49: 185–190.
Nabata, D., Masuda, R., Takahashi, O. and Nagata, J. 2004. Bottleneck effects on sika deer Cervus nippon population in Hokkaido,
revealed by ancient DNA analysis. Zoological Science 21: 473–
481.
Nabata, D., Kaji, K., Nagata, J. and Masuda, R. 2007. Genetic structure changes of expanding sika deer (Cervus nippon) populations
in central and western Hokkaido, revealed by mitochondrial DNA
analysis. Mammal Study 32: 17–22.
Nagata, J., Masuda, R., Kaji, K., Kaneko, M. and Yoshida, M. C.
1998. Genetic variation and population structure of the Japanesesika deer (Cervus nippon) in Hokkaido Island, based on mitochondrial D-loop sequences. Molecular Ecology 7: 871–877.
Nagata, J., Masuda, R., Tamate, H. B., Hamasaki, S., Ochiai, K.,
Asada, M., Tatsuzawa, S., Suda, K., Tado, H. and Yoshida, M. C.
1999. Two genetically distinct lineages of sika deer, Cervus
nippon, in Japanese Islands: comparison of mitochondrial D-loop
region sequences. Molecular Phylogenetics and Evolution 13:
511–519.
Randi, E., Mucci, N., Claro-Hergueta, F., Bonnet, A. and Douzery, E.
J. P. 2001. A mitochondrial DNA control region phylogeny of
the Cervinae: speciation in Cervus and implications for conservation. Animal Conservation 4: 1–11.

69
Riley, S. P. D., Pollinger, J. P., Sauvajot, R. M., York, E. C., Bromley,C., Fuller, T. K. and Wayne, R. K. 2006. A southern California
freeway is a physical and social barrier to gene flow in carnivores.
Molecular Ecology 15: 1733–1741.
Ripley, B. D. 1981. Spatial Statistics. John Wiley and Sons Inc, 272
pp.
Rozas, J., Sánchez-DelBarrio, J. C., Messeguer, X. and Rozas, R.
2003. DnaSP, DNA polymorphism analyses by the coalescentand other methods. Bioinformatics 19: 2496–2497.
Simonoff, J. S. 1996. Smoothing Methods in Statistics. Springer, New
York, 356 pp.
Slatkin, M. 1985. Gene flow in natural populations. Annual Review
of Ecology and Systematics 16: 393–430.
Slatkin, M. 1994. Gene flow and population structure. In (L. A. Real,ed.) Ecological Genetics, pp. 3–17. Princeton University Press,
Princeton.
Takada, M., Asada, M. and Miyashita, T. 2002. Cross-habitat
foraging by sika deer influences plant community structure in aforest-grassland landscape. Oecologia 133: 389–394.
Ujvari, B., Dowton, M. and Madsen, T. 2007. Mitochondrial DNA
recombination in a free-ranging Australian lizard. Biology Letters
3: 189–192.
Wang, M. and Schreiber, A. 2001. The impact of habitat fragmentation and social structure on the population genetics of roe deer
(Capreolus capreolus L.) in Central Europe. Heredity 86: 703–
715.
Weir, B. S. and Cockerham, C. C. 1984. Estimating F-statistics for the
analysis of population structure. Evolution 38: 1358–1370.
Wright, S. 1921. Systems of mating. Genetics 6: 111–178.
Yamada, M., Hosoi, E., Tamate, H. B., Nagata, L., Tatsuzawa, S.,
Tado, H. and Ozawa, S. 2006. Distribution of two distinct
lineages of sika deer (Cervus nippon) on Shikoku Island revealedby mitochondrial DNA analysis. Mammal Study 31: 23–28.
Yamamoto, S., Morita, K., Koizumi, I. and Maekawa, K. 2004.
Genetic differentiation of white-spotted charr (Salvelinus
leucomaenis) populations after habitat fragmentation: spatialtemporal changes in gene frequencies. Conservation Genetics 5:
529–538.
Young, A., Boyle, T. and Brown, T. 1996. The population genetic
consequences of habitat fragmentation for plants. Trends in
Ecology and Evolution 11: 413–418.
Yuasa, T., Nagata, J., Hamasaki, S., Tsuruga, H. and Furubayashi, K.
2007. The impact of habitat fragmentation on genetic structure of
the Japanese sika deer (Cervus nippon) in southern Kantoh,
revealed by mitochondrial D-loop sequences. Ecological Research
22: 97–106.
Received 13 July 2007. Accepted 20 December 2007.