Phylogeography of the Eurasian red squirrel (Sciurus vulgaris orientis) on the boreal island of Hokkaido, Japan

Abstract To test postglacial population expansion in small arboreal mammals dependent on boreal and subboreal forests, we used complete mitochondrial cytochrome b gene sequences (1140 bases) to investigate the phylogeography of the Eurasian red squirrel (Sciurus vulgaris orientis) on the Hokkaido Island of Japan. This subspecies is common in the boreal and subboreal forests of Hokkaido. We examined 61 specimens from 27 localities in Hokkaido. Phylogenetic relationships among 29 haplotypes found in the Hokkaido populations were not associated with geographic distribution of sampling localities. There were four mitochondrial DNA phylogroups. Phylogeographic analyses support sudden expansion of S. vulgaris orientis from restricted refugium in the southern part of Hokkaido during the last glacial period. The phylogeographic structure of this subspecies directly reflects the boreal and subboreal forest dynamics occurring in Hokkaido.


Introduction
The refugia theory explains that animal populations were reduced when available habitats decreased during the cyclic glaciers of the alternating cold and warm episodes during the Pleistocene (Haffer 1969;Cracraft and Prum 1988;Cox and Moore 2005). Boreal and subboreal forest dynamics caused by climatic changes have greatly affected continental scale distributions of arboreal small mammals. Resolution of the phylogeographic structure of arboreal small mammals permits indirect tracing of historical boreal and subboreal forest dynamics during the Pleistocene, making arboreal small mammals useful phylogeographic indicators. In fact, phylogeographic structure of the northern flying squirrel (Glaucomys sabrinus) is associated with forest dynamics in North America (Arbogast 1999). In the Eurasian Continent, phylogeography of the Siberian flying squirrel (Pteromys volans) is also influenced by forest dynamics (Oshida et al. 2005).
Boreal and subboreal forest dynamics could have affected the evolutionary history of arboreal small mammals on islands as well as continents. Islands are effective natural laboratories for biogeographic studies (e.g. Whittaker 1998). Evolutionary factors such as geographic isolation, population expansion and restricted gene flow may be more clearly correlated with phylogeographic structure. Populations of arboreal mammals distributed on islands do demonstrate clear phylogeographic traits caused by forest dynamics. Both the Japanese macaque (Macaca fuscata yakui) (Kawamoto et al. 2007) and the Japanese giant flying squirrel (Petaurista leucogenys) (Oshida et al. 2009b) experienced population contraction and expansion on Honshu Island, Japan. Pallas's squirrel (Callosciurus erythraeus) experienced population contraction, expansion and isolation in the mountain ranges of Taiwan (Oshida et al. 2006).
To investigate phylogeographic consequences of boreal and subboreal forest dynamics on arboreal small mammals, we examined the Eurasian red squirrel (Sciurus vulgaris orientis) population on Hokkaido Island, Japan. Hokkaido, whose area encompasses 78,073 km 2 , is the northernmost island of Japan's four main islands. At present, it is mainly covered with boreal coniferous forests characterized by Picea and Abies. There are Tsuga in some montane forests at high elevations. Temperate mixed deciduous forests are present in the southern part of the island (Yamanaka 1979;Morley et al. 1986;Nakanishi et al. 1996; Figure 1A). During the last glaciation, however, montane forests characterized by Pinus and Larix dominated northern Hokkaido, with boreal coniferous forests shifting to southern Hokkaido (Ono and Igarashi 1991;Dobson 1994; Millien-Parra and Jaeger 1999; Figure 1B). This simple change in forest dynamics during the Pleistocene should have affected the phylogeography of arboreal small mammals. Therefore, Hokkaido is an ideal model island to study the association between boreal and subboreal forest dynamics and arboreal small mammals.  Sciurus vulgaris is widely distributed in the Palearctic region (from Iberia and Great Britain east to Kamchatka Peninsula, Korea and the Sakhalin and Hokkaido Islands; south to the Mediterranean and Black Seas, northern Mongolia and western and northeastern China) (Gurnell 1987;Thorington and Hoffmann 2005). The population in Hokkaido is the endemic subspecies S. vulgaris orientis (e.g. Thorington et al. 2012). It occurs commonly from lowlands to mountainous areas (Tamura 2009). In Hokkaido, this squirrel mainly inhabits coniferous and mixed coniferous forests, but not montane forests (Tamura 2009). Therefore, when the boreal and subboreal coniferous forests shifted and contracted with glaciation (Ono and Igarashi 1991;Dobson 1994;Millien-Parra and Jaeger 1999), S. vulgaris orientis would have retreated to southern Hokkaido. After glaciation, this subspecies would have again expanded to northern Hokkaido. During the last glacial maximum, Hokkaido was connected to the Eurasian Continent (Ohshima 1990(Ohshima , 1991, but S. vulgaris could not have migrated from the continent to Hokkaido Island, because of the reduction of boreal and subboreal coniferous forests ( Figure 1B).
To test these hypotheses, we examined the complete mitochondrial cytochrome b gene sequence of S. vulgaris orientis in Hokkaido. This genetic marker has been successfully used in studies of the phylogeographic structure of other squirrel populations (Demboski et al. 1998;Arbogast 1999;Good and Sullivan 2001;Oshida et al. 2005Oshida et al. , 2009bGündüz et al. 2007;Fernández 2012;Moncrief et al. 2012). The phylogeographic structure of regional populations of S. vulgaris in China (Liu et al. 2014) and Europe (Grill et al. 2009) has also been resolved with this marker. In the present study, we test the following three hypotheses.
Hypothesis 1: We hypothesize that the subspecies expanded rapidly from a southern refugium to the whole of Hokkaido Island after glaciation. If true, we expect a single unstructured population with low genetic diversity.
Hypothesis 2: We hypothesize that the subspecies was isolated or fragmented into several refugia during the last glacial maximum in Hokkaido Island. If true, we expect that clear divergences among regional populations will be found.
Hypothesis 3: We hypothesize that the subspecies did not experience radical contraction or expansion and had a fairly stable evolutionary history. If true, we expect that moderate genetic variations, but no clear divergences will be found.
Considering the dynamics of Hokkaido's boreal and subboreal forests, however, we expect the data to support hypotheses 1 or 2. Evidence of population expansion of S. vulgaris orientis would also support forest shifts in Hokkaido Island. Here, we discuss the evolutionary history of S. vulgaris orientis in Hokkaido both during and after glaciation.

Study regions and specimens
Specimens of Sciurus vulgaris orientis examined in this study are listed in Appendix 1. A total of 61 specimens from 27 localities in two regions (southern and northern) of Hokkaido (Table 1 and Figure 2) were used for phylogeographic analyses. During the glaciation, the southern region was covered with boreal coniferous forests (refugia for S. vulgaris orientis), but the northern region was covered by montane forest and tundra (Dobson 1994; Millien-Parra and Jaeger 1999; Figure 1B). Since the southern and northern regions differ in historical forest dynamics, we separated Hokkaido into these two regions for phylogeographic examinations.
Seventeen specimens were provided from the Shiretoko Museum, Shari. Two were provided by the Higashi Taisetsu Museum of Natural History, Kamishihoro. Two specimens were from Dr. M. Asakawa (Rakuno Gakuen University, Ebetsu). Seven specimens deposited in the Museum of Botanic Garden, Hokkaido University, Sapporo were used. Remaining specimens were road-kill collected by many collaborators and deposited in our laboratory. Moreover, to approximate diversity between the Eurasian continent and Hokkaido Island populations, two Sciurus vulgaris specimens (Sc1 and Sc2) from Zaigraevo, Transbaikalia, Russia deposited in the Zoological Institute of the Russian Academy Sciences, St. Petersburg were added to the molecular phylogenetic analyses, although this was not directly related to our purposes.

DNA extraction, amplification and sequencing
Total genomic DNA was extracted with the QIAQuick Kit (QIAGEN K.K., Tokyo, Japan) from muscle tissue, preserved in 99% ethanol and suspended in TE buffer. Whole Table 1. Sampling localities in Hokkaido, Japan, and mtDNA haplotypes of Sciurus vulgaris orientis examined in this study.

Sequence and phylogenetic analyses
Sequence alignment to detect unique haplotypes was with the software program DNASIS (Hitachi, Tokyo, Japan). Modeltest version 3.06 (Posada and Crandall 1998) was used to find the nucleotide substitution model that best fitted the haplotype data set. The hierarchical likelihood ratio tests implemented in Modeltest selected the Hasegawa-Kishino-Yano model (Hasegawa et al. 1985) with a gamma distribution for variable sites (HKY + G). Bayesian inference (BI) was carried out using MrBayes 3.2 (Huelsenbeck and Ronquist 2001). Bayesian analysis was conducted using the HKY + G model selected by Modeltest for our data set. The analysis involved four runs for one million iterations with four Markov chain Monte Carlo chains sampled every 1000 generations and a burn-in of 20%. Then, 50% majority rule consensus trees based on the remaining trees were generated. Posterior probabilities assessed nodal support of the BI tree. Maximum-likelihood (ML) analysis was conducted with the HKY + G model in PAUP* 4.0b10 (Swofford 2001). The ML tree was constructed with the heuristic search option of tree-bisection-reconnection. To assess nodal supports, bootstrapping (Felsenstein 1985) was performed with 300 replicates. To root each phylogenetic tree, we used cytochrome b gene sequences of the Japanese squirrel Sciurus lis as an outgroup: accession number AB192913 in the DNA Data Bank of Japan. Sciurus lis is most closely related to S. vulgaris (Oshida and Yoshida 1997;Oshida and Masuda 2000;Grill et al. 2009;Oshida et al. 2009a;Pečnerová and Martínková 2012), making it a suitable evolutionary root. To grasp phylogenetic relationships among haplotypes, a haplotype network was constructed by statistical parsimony using a 95% confidence limit between haplotypes (Templeton et al. 1992) in TCS version 1.21 (Clement et al. 2005).

Phylogeographic analyses
For the southern and northern regions, including all 27 sampling sites (Figure 2), partitioning of total genetic variation was hierarchically examined by analysis of molecular variance (AMOVA: Excoffier et al. 1992) in ARLEQUIN 2.001 (Schneider et al. 2000). The significance of the fixation indices was tested using permutation procedures described by Excoffier et al. (1992) with 1000 permutations with ΦST as the correlation of random haplotypes between populations. Haplotype and nucleotide diversities (Nei 1987) within groups were calculated with DnaSP 5.0 (Librado and Rozas 2009). To assess the past demographic expansion of this subspecies with three data sets (southern, northern and over all of Hokkaido), we calculated Tajima's D (Tajima 1989) and Fu's F S (Fu 1997) neutralities in DnaSP 5.0 (Librado and Rozas 2009). When a population has recently begun to expand, Tajima's D and Fu's F S have large negative values. In each of the southern, northern and over all of Hokkaido populations, mismatch distributions of haplotypes (Slatkin and Hudson 1991;Schneider and Excoffier 1999) were analysed according to the sudden expansion model (Rogers and Harpending 1992) in ARLEQUIN 2.001. The confidence interval for ι was obtained from 1000 bootstrap replicates. Under this model, recent population expansion was expected to generate a unimodal distribution of pairwise sequence differences (Slatkin and Hudson 1991;Rogers 1995;Schneider and Excoffier 1999).

Sequence composition and variation
Complete sequences (1140 bp) of cytochrome b gene were successfully determined from 61 specimens of S. vulgaris orientis, resulting in 29 unique haplotypes (Table 1). An additional haplotype was found from the two S. vulgaris specimens from Russia (Appendix 1). All sequences (Appendix 1) were deposited in the DNA Data Bank of Japan (DDBJ). Among the haplotypes from Hokkaido, there were 19 segregating sites, of which 11 were parsimony informative.

Phylogenetic analysis
There were one major (A) and one minor (B) phylogroup in the BI tree with high posterior probabilities (Figure 3). Within phylogroup A, there were two minor phylogroups (A-1 and A-2). Phylogroup A was distributed throughout Hokkaido. Also, phylogroup A-1 was widely found in Hokkaido, but phylogroup A-2 was constrained to southern Hokkaido (see Table 2 and Figure 2). Phylogroup B was found in both southern and northern Hokkaido. The ML phylogenetic tree, however, had low bootstrap values, and showed no distinct mtDNA phylogroups between S. vulgaris orientis in Hokkaido and S. vulgaris from Russia (tree not shown). A haplotype (RS1) from two Russian specimens was distantly related to other haplotypes from Hokkaido. The network tree (Figure 4) indicated haplotype HS4 as an ancestral form of S. vulgaris orientis; eight haplotypes (HS3, HS13, HS15, HS16, HS17, HS22, HS24 and HS28) seemed to be derived from haplotype HS4.

Phylogeographic analysis
Haplotype and nucleotide diversities of southern and northern populations are shown in Table 2. Tajima's D and Fu's Fs neutralities are also shown in Table 2. All values were negative. Although Fu's Fs test significantly supported neutrality for southern, northern and over all of Hokkaido populations (P < 0.005, 0.0001 and 0.0001, respectively), Tajima's D test for three data sets did not show significant neutrality (P > 0.1, 0.05 and 0.1, respectively). The test for partitioning of the total genetic variation (Table 3) showed no significant variation between the two regional populations (P > 0.1). Mismatch analyses in northern and over all of Hokkaido populations showed a clear unimodal distribution that fitted the sudden expansion model; in particular, the mismatch distribution of the northern population was almost similar to that of the simulated expansion model ( Figure 5). Also, the mismatch distribution of the southern population seemed to be unimodal, but its shape was not clear with one irregular peak ( Figure 5). The observed mismatch distributions in southern, northern and over all of Hokkaido populations were not significantly different from the simulated sudden expansion model (P = 0.85, 0.29 and 0.66, respectively).

Discussion
In the BI tree (Figure 3), we found four mtDNA phylogroups of S. vulgaris orientis in Hokkaido, but not in the ML tree. In general, it is reported that posterior probabilities of BI trees are often higher than bootstrap values of the ML tree (e.g. Suzuki et al. 2002;Lewis et al. 2005;Randle et al. 2005). Therefore, although our result suggests that the Hokkaido population experienced geographical isolation or fragmentation during glaciation (hypothesis 2), it is difficult to completely discard hypotheses 1 and 3. The AMOVA showed no partitioning among regional populations. Tajima's D and Fu's Fs were negative for all three data sets (southern, northern and over all of Hokkaido), although we did not find significant neutrality in Tajima's D test. In addition, the mismatch distributions for all three data sets showed unimodal shape and insignificant differences from the sudden expansion model. Therefore, S. vulgaris orientis experienced a sudden expansion from a restricted refugium in the glacial periods. The mismatch distribution pattern of the northern population was almost similar to that of the expansion model ( Figure 5). This may be important proof to support the northward  expansion of S. vulgaris orientis. All specimens of phylogroup A-2 were restricted to southern Hokkaido. Therefore, the refugium could have been in the southern part of Hokkaido, where boreal coniferous forests retreated during glaciation (Figure 1), suggesting that hypothesis 1 may be reasonable. Haplotype HS4 seemed to be the ancestor of eight haplotypes (Figure 4). This may suggest that many haplotypes found in Hokkaido were derived from some founders in the refugium. Forest expansion northwards from a southern refugium may have aided population expansion of S. vulgaris orientis. During the Pleistocene, interglacial retreat of the ice sheet resulted in migration of conifer species from glacial refugia into open spaces in North America (Jaramillo-Correa et al. 2004) and Europe (Taberlet et al. 1998). Except for some small parts of mountainous regions, northeast Asia including Hokkaido was virtually free of ice during the Pleistocene (e.g. Hewitt 2000). The lack of ice allowed moderate southward or northward shifting of plant species (Maekawa 1974). The arboreal S. vulgaris orientis probably shifted range with the coniferous forest.
There are a few reports on the phylogeographic structure of forest mammals in Hokkaido. Interestingly, these reports show contradictory phylogeographic patterns. Based on mtDNA control region sequences, Matsuhashi et al. (1999) found three clear lineages of the brown bear (Ursus arctos), which originated from different continental populations. The brown bear has the ability to travel great distances and to overcome geographic barriers. Nevertheless, in Hokkaido, phylogroups of brown bears are geographically separated into three regions: western, central and eastern. Unlike brown bears, the sika deer (Cervus nippon) experienced a population bottleneck in Hokkaido and has extremely low genetic variability (Nagata et al. 1998). The red fox (Vulpes vulpes) has two lineages, but these lineages are sympatrically distributed throughout Hokkaido (Inoue et al. 2007). The long-clawed shrew Sorex unguiculatus in Hokkaido did not show a divergence clear enough for geographic phylogroups (Ohdachi et al. 2001). Thus, Hokkaido's forest mammals do not show congruent phylogeographic characteristics, suggesting different evolutionary histories. During glaciation, it would have been difficult for arboreal forest mammals to move across wide treeless areas, such as grasslands. Therefore, the phylogeographic history of S. vulgaris orientis in Hokkaido may be more tightly correlated with subboreal and boreal forest dynamics, showing evidences of population expansion. To precisely explain its phylogeographic structure, we would need a more complicated hypothesis on the basis of a combination of hypotheses 1 and 2.
Since we analysed mtDNA, our findings may be affected by sex-biased dispersal or mating system. Further study of the population dynamics of S. vulgaris orientis with nuclear genes or microsatellite markers is necessary for a more complete understanding of the phylogeography of this subspecies. Table 3. Analysis of molecular variance (AMOVA) for 61 specimens of Sciurus vulgaris orientis in two populations as defined in Figure 2, in Hokkaido, Japan.