Phylogeny and genetic structure of the Yellow ground squirrel, Spermophilus fulvus (Lichtenstein, 1823), in Iran

Mammalian Biology 98 (2019) 137–145

Contents lists available at ScienceDirect

Mammalian Biology journal homepage: www.elsevier.com/locate/mambio

Original investigation

Phylogeny and genetic structure of the Yellow ground squirrel, Spermophilus fulvus (Lichtenstein, 1823), in Iran Afsaneh Asgharzadeh a , Mohammad Kaboli b,∗ , Hassan Rajabi-Maham c , Morteza Naderi d a Department of Environment Engineering, Faculty of Natural Resources and Environment, Science and Research Branch, Islamic Azad University, Tehran, Iran b Department of Environmental Science, Faculty of Natural Resources, University of Tehran, Karaj, Iran c Department of Animal Sciences and Biotechnology, Faculty of Life Sciences and Biotechnology, Shahid Beheshti University, G.C., Evin, 1983963113, Tehran, Iran d Department of Environmental Sciences, Faculty of Agriculture and Natural Resources, Arak University, 38156-8-8349, Arak, Iran

a r t i c l e

i n f o

Article history: Received 27 December 2018 Accepted 12 September 2019 Available online 13 September 2019 Handled by Paul Grobler Keywords: Spermophilus fulvus Phylogeny Phylogeography Divergence time Continental interglacial refugia Quaternary climate oscillations

a b s t r a c t Old world ground squirrels (genus Spermophilus) are distributed throughout the Holarctic and Palearctic regions, of which two Iranian species, the Yellow ground squirrel S. fulvus and the Asia Minor ground squirrel S. xanthoprymnus, comprise the southernmost distribution of the genus in the Palearctic. The two species are found in fragmented populations from northeastern to northwestern Iran, with S. fulvus being more common and widespread in the country. The enormous geographical distance (more than 1000 km) between Yellow ground squirrel’s patchy populations in eastern and western Iran has led to ambiguous evolutionary relationships and consequently uncertain conservation planning for the species. We studied the phylogenetic relationships and spatial genetic structure of the isolated populations using the mitochondrial Cytochrome b (cyt b) among 79 individuals. Our phylogenetic analysis found that S. fulvus was divided into three main mtDNA clades in Iran. Molecular dating suggested an ancestral expansion from high latitudes towards Iran during the cold periods (before the Holsteinian temperate period), followed by contraction of populations into interglacial refugia around the Holsteinian temperate period (430-350 Kya), leading to their isolation from their ancestral pool about 427 Kya. Furthermore, we found a non-deep phylogenetic divergence between patchy isolated populations of the species in Iran. A major genetic break was detected between populations of eastern Iran (steppes in Hezar Masjed and Binalud mountains, near Iran-Turkmenistan-Afghanistan borderline) and the others, probably associated with their confinement to interglacial refugia during the Domnitz temperate period (340-300 Kya). However, the recent divergence between northeastern and western populations of Iran probably resulted from their contraction into interglacial refugia in the Saale/Drenthe temperate period (250-190 Kya). These results may suggest that populations of the Yellow ground squirrel in Iran fell under the influence of glacial and interglacial cycles along a longitudinal axis. Finally, we propose the three isolated clades of the Yellow ground squirrel identified in Iran as three ESUs. ¨ Saugetierkunde. ¨ © 2019 Deutsche Gesellschaft fur Published by Elsevier GmbH. All rights reserved.

Introduction Old world ground squirrels (genus Spermophilus Cuvier, 1825) are known as one of the most widespread representatives of the subfamily Xerinae (Thorington et al., 2012). This genus comprises 14 species, mostly including the Eurasian members of ground squirrels, predominantly common to open plains, steppes, and tundra habitats of the Palearctic, although some species are found in semi-

∗ Corresponding author. E-mail address: [email protected] (M. Kaboli).

desert ecosystems across the southern portion of the distribution range (Helgen et al., 2009). The genus Spermophilus is represented by two species in Iran, the Asia Minor ground squirrel, S. xanthoprymnus (Bennett, 1835) and the Yellow ground squirrel, S. fulvus (Lichtenstein, 1823). The latter, S. fulvus, inhabits cold steppes with sparse vegetation and productive clay soils (Karami et al., 2016). It occurs in two fragmented areas in northwest and northeast of Iran (Kryˇstufek and Vohralík, 2012; Karami et al., 2016) at the southern fringe of the distribution area of the species in the Palearctic. Although the phylogeny and phylogeography of ground squirrels have been investigated by a number of studies (e.g., Harrison

https://doi.org/10.1016/j.mambio.2019.09.007 ¨ Saugetierkunde. ¨ Published by Elsevier GmbH. All rights reserved. 1616-5047/© 2019 Deutsche Gesellschaft fur

138

A. Asgharzadeh et al. / Mammalian Biology 98 (2019) 137–145

et al., 2003; Herron et al., 2004; Kapustina et al., 2015; Faerman et al., 2017), the evolutionary history and genetic structure of the Iranian Yellow ground squirrel is still poorly known. The species comprises some completely isolated populations with more than 1000 km distance among them. All field surveys conducted by the experts and guards of the Department of Environment of Iran as well as our extensive field sampling during this study failed to detect any new populations between the northeastern and western populations. Hence, several questions have remained unsolved regarding the evolutionary history of the species, including the genetic structure, phylogenetic relationships among populations, gene flow and factors affecting gene flow of the species in Iran. Provision of such data is essential for conservation of genetic diversity, which enables species to adapt to changes in their surrounding environment (Sgro et al., 2011). Herein, we aimed to clarify the phylogenetic relationships and spatial genetic structure of the patchy populations of S. fulvus in Iran based on a mitochondrial gene (cyt b). As this study is the first mitochondrial analysis of the phylogeny and phylogeography of a member of the genus Spermophilus in its southernmost distribution in the Palearctic, we aimed to address the following questions: (i) Are isolated populations of S. fulvus in northeastern and northwestern Iran genetically and significantly distinct from one another, or do they belong to a single phylogroup? (ii) When did the Iranian ground squirrel disperse into Iran and when did it diverge throughout Iran? (iii) Could the divergence of the Iranian populations be attributed to population expansions and contractions during glacial and interglacial periods, followed by their consequent confinement to glacial refugia?

Material and methods Sample collection During April to July 2014–2017, 79 ear tissue (from livecaptured specimens) or muscle tissue samples (from road-killed carcasses) of the Yellow ground squirrel were collected from 50 localities in Iran (Fig. 1; Table S1). The individuals were captured by net and small ear tissue samples were obtained using a sterile scalpel. The captured animals were immediately released back into the wild after sampling. Our sampling was licensed by the Iranian Department of Environment under permits of 94/6150 and 94/57879.

Sequence analysis All sequences were edited by SeqScape version 2.6 (Applied Biosystems). We included the only available (two) samples from Russia and Kazakhstan downloaded from GenBank (Harrison et al., 2003: AF157913 and AF157908; Table S1). For outgroups, we added two cyt b sequences of S. taurensis (KY938058) and S. citellus (AF157859). Sequences were aligned using Clustal W in MEGA 6 (Thompson et al., 1994; Tamura et al., 2013). The alignments were tested for substitution saturation using DAMBE 6.0.4 program (Xia and Xie, 2001; Xia and Lemey, 2009). Then, nucleotide and haplotype diversities, number of haplotypes and polymorphic sites were estimated by DnaSP 5.0 (Librado and Rozas, 2009). Furthermore, nucleotide composition and genetic distances were calculated using an uncorrected genetic distance in MEGA 6 based on 10,000 bootstrap replications. Phylogenetic analyses Phylogenetic reconstruction was performed using Bayesian Inference (BI) and Maximum Likelihood (ML) methods. The program PartitionFinder (Lanfear et al., 2012) was used to determine the best partitioning schemes and substitution models for the analysis based on the first, second, and third positions within codons. We found three subset partitions including: (i) cyt b pos1, (ii) cyt b pos2, and (iii) cyt b pos3 with the optimal models K80, F81, and GTR, respectively. This partitioning strategy was used for both the Bayesian Inference and Maximum Likelihood analyses. A Bayesian phylogenetic analysis was conducted by implementing the selected scheme in MrBayes v3.2.6 (Ronquist and Huelsenbeck, 2003). Four Markov Chains Monte Carlo (MCMC) were simultaneously run with two replicates for 80 million generations. As the trees were sampled every 4000 generations, the first 25% of the steps were discarded as burn-in. The remaining trees were combined in a 50% majority rule consensus tree. The results of each run were visualized using Tracer v1.5 (Rambaut and Drummond, 2007) to ensure stationarity and convergence were achieved and effective sample sizes (ESS) were higher than 500. Bayesian posterior probabilities were used to assess branch support of the BI tree. IQ-Tree version 1.6.8 (Nguyen et al., 2014) was used to construct a maximum likelihood tree and topology tested by 1000 non-parametric bootstrap replicates. Finally, we performed the Shimodaira-Hasegawa (SH) (Shimodaira and Hasegawa, 1999) likelihood-based test, using PAUP 4.0a165 (Swofford, 2002) in 1000 repeats to confirm both ML and BI trees have the same topology. Spatial and non-spatial analyses of genetic variation

DNA isolation, amplification, and sequencing First, DNA was extracted from tissue samples (muscle and ear) fixed in 96% ethanol using proteinase K and phenol–chloroform standard protocol (Sambrook et al., 1989). Then, DNA pellets were dissolved in standard Tris-EDTA buffer (TE) pH 8.0. We amplified a fragment of the mitochondrial genome including 1140 base pairs (bp) of cyt b for 79 samples. We used specific primers L14725spermo (5 -TGAAAAAYCATCGTTGT-3 ) and H15915-spermo (5 TCTTCATTTYWGGTTTACAAGAC -3 ) from Harrison et al.’s study (2003). In addition, polymerase chain reaction was conducted by the Ampliqon Taq DNA Polymerase 2x Master Mix Red (1.5 Mm MgCl2 final concentration). Each PCR reaction contained 12.5 ␮l Master Mix, 1 ␮l of 10 pmol of each primer, 1 ␮l template DNA, and 9.5 ␮l distilled H2 O in a final volume of 25 ␮l. PCR reaction conditions are shown in Table S2. PCR products were sent to Macrogen and Bioneer (South Korea), where sequences were determined based on the ABI 3730 automated DNA sequencer (Applied Biosystems).

Bayesian analysis of population structure was conducted to identify distinct genetic clusters using BAPS V. 6.0 (Cheng et al., 2013). Values of K (number of clusters) were fixed from one to seven for the analysis. In addition, F-statistics (pairwise FST ) were used with 10,000 permutations to estimate the proportion of genetic variability among populations. Two haplotype networks were constructed on the mtDNA matrix (81 individuals with 1140 bp), to represent genetic relationships amongst S. fulvus haplotypes: (i) TCS (Clement et al., 2000) implemented in PopART software (Leigh and Bryant, 2015; http:// popart.otago.ac.nz), and (ii) SplitsTree v4.6 (Huson and Bryant, 2006) to draw a neighbor-net network based on uncorrected patristic distances and 1000 bootstrap replicates. Molecular dating estimates Molecular divergence dating was estimated using BEAST 1.8.0 (Drummond et al., 2013). Molecular dating is largely affected by the

A. Asgharzadeh et al. / Mammalian Biology 98 (2019) 137–145

139

Fig. 1. Distribution of S. fulvus in Asia (cross-hatched area) and sampling locations in Iran: yellow, green, and purple circles correspond to phylogenetic clades (clades IR1, IR2, and IR3, respectively) of the species in Iran (see Fig. 2). The two blue circles represent the two sample locations from Russia and Kazakhstan. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

quality of calibration points at internal nodes, which are normally based on fossil evidence (e.g., Lieberman, 2002). In this regard, we used three fossil calibration points (nodes 102, 103, and 107: Goodwin, 1995; Harrison et al., 2003) : (i) divergence of the Cynomys subgenera, 2 Mya was considered using a lognormal prior model, (ii) divergence of genus Cynomys, 2.7 Mya was considered using a lognormal distribution model, and (iii) divergence of the genus Spermophilus, 7.3 Mya was implemented using a lognormal prior model. Based on the time-calibrated species tree of Harrison et al. (2003), the present molecular dating was estimated and the cyt b taxa of Cynomys and Spermophilus’s nodes were added to determine the divergence time of Spermophilus in Iran. A birth-death process was implemented as a more appropriate model for sequences involving different species. In addition, three molecular clock models (strict, exponential, and lognormal relaxed) were evaluated, and accordingly, the best clock model was specified based on a Bayes factor analysis (value of 2LnBF; Brandley et al., 2005) as applied in Tracer v1.5. Then, two independent runs of 100 million generations were performed and sampled every 5000 generations, with the first 25% discarded as burn-in. In the next step, the convergence statistics of BEAST runs were evaluated by examining the effective sample size in Tracer version v1.5. Finally, TreeAnnotator v1.8.2 and FigTree v1.4.1 were used to generate and visualize a molecular time tree, respectively. Historical demographic analyses An analysis of past demographic history including Tajima’s D (Tajima, 1989) and Fu’s Fs (Fu, 1997) was estimated by 10,000

replicates under the pairwise difference model in Arlequin v.3.1 (Excoffier et al., 2005). In addition, based on the approach suggested by Rogers and Harpending (1992) and Slatkin and Hudson (1991), a mismatch distribution analysis was inferred by pairwise differences between individuals in order to determine whether a particular clade could represent recent population expansion. In the next step, the distribution was compared with the expected distribution under the sudden expansion model of Rogers (1995) and evaluated with goodness-of-fit tests based on Harpending’s raggedness index (RI) and the sum of the squared deviations (SSD) using Arlequin version 3.1 (Excoffier et al., 2005). Results Sequence analysis The alignments of the complete nucleotide sequence of cyt b (1140 bp) in 79 S. fulvus specimens contained 11 singletons, 24 parsimony informative sites, and 26 unique haplotypes. Nucleotide diversity (pi) and haplotype diversity (Hd) values were estimated 0.005 and 0.916, respectively, for all sample sets (Table 1). No insertions, deletions, or stop codons were observed in our dataset. Finally, transition/transversion bias ratio (R) was estimated to be approximately 6.67. Nucleotide frequencies were T = 34.0%, C = 25.2%, A = 28.3%, and G = 12.6%. Table 1 represents nucleotide and haplotype diversities for Iranian S. fulvus clades obtained from phylogenetic analyses (see Fig. 2). Saturation analyses showed that our cyt b matrix has little substitution saturation (Fig. S1). The observed index of substitution saturation (ISS) value was significantly lower than the critical index

140

A. Asgharzadeh et al. / Mammalian Biology 98 (2019) 137–145

Fig. 2. Bayesian tree reconstructed from cyt b (in congruence with the ML tree in terms of the branching pattern and the positions of the clades), using S. citellus and S. taurensis as outgroups. Nodal supports presented at nodes indicate Bayesian posterior probability (top) by MrBayes, and bootstrap (bottom) by IQ-TREE. Yellow: Clade IR1 (populations from eastern and northeastern Iran along with Iran-Turkmenistan-Afghanistan border), green: Clade IR2 (northeastern and western populations), purple: Clade IR3 (western populations), and blue: two samples from Russia (AF157913: Saratovskaya obl., Krasnokutskii rai., vicinity Diakova village) and Kazakhstan (AF157908: Uralskaya obl., Urdinskii rai., vicinity Tarlykova village) from Harrison et al. (2003). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1 Results of genetic statistics based on cyt b for Iranian lineages of S. fulvus. n: number of sequences, h: number of haplotypes, pi: nucleotide diversity, and Hd: haplotype diversity.

All Clade IR1 Clade IR2 Clade IR3

n

h

pi

Hd

79 18 28 33

26 7 9 10

0.00568 0.00147 0.00105 0.00095

0.916 0.784 0.796 0.716

(ISSc), indicating that our data was appropriate for phylogenetic analyses. Phylogenetic reconstruction Phylogenetic trees reconstructed by both ML and BI approaches yielded similar topologies. The Shimodaira-Hasegawa test failed to

detect any significant differences (p > 0.05) between the two trees and hence, only the BI tree is presented (Fig. 2). The Iranian clade appears to be strongly-supported (1.0 posterior probability and 100% bootstrap support, Fig. 2), despite including only two samples outside the Iranian distribution range. Moreover, populations from eastern and northeastern Iran, along with the individuals from Iran-Turkmenistan-Afghanistan borderline (hereafter, Clade IR1) were separated from the remaining Iranian populations with rather high posterior probability (0.98) and bootstrap (94.5%) values (yellow clade in Fig. 2). However, the remaining two sister clades, clades 2 and 3 (hereafter, Clade IR2 and Clade IR3, respectively), showed moderate posterior probability (0.82) and bootstrap (79.8%) supports. Clade IR2 includes populations from northeastern and western provinces of Iran (North Khorasan, Razavi Khorasan, Zanjan, and Kurdistan provinces; green clade in Fig. 2), whereas Clade IR3 includes populations from western Iran (Hamedan, Qazvin, and Alborz provinces; purple clade in Fig. 2).

A. Asgharzadeh et al. / Mammalian Biology 98 (2019) 137–145

141

Fig. 3. SplitsTree and TCS haplotype networks based on cyt b. Colorations are concordant with Figs. 1, 2, 4, and 5. (For interpretation of the references to colour in the text, the reader is referred to the web version of this article.)

Table 2 Uncorrected p-distances (below diagonal) with associated standard errors (SE) based on 10,000 bootstrap replicates (bold numbers in parentheses), along with pairwise FST values (above diagonal) among the Yellow ground squirrel phylogroups based on cyt b data (significant value of FST in ␣ = 0.01 is shown by **). IR-TU-AF: Iran-Turkmenistan-Afghanistan borderline populations, NE W Iran: northeastern and western populations, and W Iran: western populations. Clade IR1 (IR-TU-AF) Clade IR1 (IR-TU-AF) Clade IR2 (NE-W Iran) Clade IR3 (W Iran)

0.011 (0.003) 0.009 (0.002)

Clade IR2 (NE-W Iran)

Clade IR3 (W Iran)

0.887**

0.838** 0.871**

0.006 (0.002)

Genetic distances The genetic distances among Iranian Yellow ground squirrel clades (see Fig. 2) varied from 0.6% to 1.1% (Table 2). Clade IR1 (populations from eastern and northeastern Iran along with IranTurkmenistan-Afghanistan border) possessed the greatest genetic distance to the other clades. However, the distance dropped to 0.6% between the two clades of IR2 and IR3.

Haplotype networks Both SplitsTree and TCS haplotype networks (Fig. 3A and B, receptively) were in accordance with the phylogenetic tree and recovered four haplogroups from Russia and Kazakhstan, eastern and northeastern Iran (Clade IR1) to western Iran (clades IR2 and IR3). These four clusters were separated from each other by 4–8 mutational steps. The first haplogroup, corresponding to Clade IR1 (yellow haplogroup in Fig. 3B), comprised seven haplotypes and showed high haplotype diversity (Hd = 0.784; Table 1). We found no common haplotypes between this haplogroup and the others. On the other hand, the other two haplogroups, corresponding to clades IR2 and IR3 (green and purple haplogroups in Fig. 3B, respectively), were separated from each other by 4 mutational steps. They included nine and 10 haplotypes with 0.796 and 0.716 haplotype diversity, respectively (Table 1). Furthermore, the two samples from Russia and Kazakhstan shared one common haplotype (blue haplogroup in Fig. 3B) and were separated from the Iranian haplogroups, IR1 and IR3, by 8 and 5 mutational steps, respectively. It seems that lack of enough samples from high latitudinal countries

142

A. Asgharzadeh et al. / Mammalian Biology 98 (2019) 137–145

has led to small mutational steps between Iran and those countries in our analysis. Spatial genetic structure In line with the phylogenetic tree and haplotype networks, the Bayesian Analysis of Population Structure (BAPS) identified three clusters of S. fulvus in Iran (Fig. 4). These clusters include (i) populations in concordance with Clade IR1 and the relevant haplogroups (from 34.89 N and 60.79 E to 36.70 N and 57.89 E, comprising populations from eastern and northeastern Iran, as well as IranTurkmenistan-Afghanistan borderline (IR-TU-AF; yellow polygons in Fig. 4), (ii) populations from northeastern and western Iran (NEW Iran) (from 36.90 N and 58.85 E to 35.85 N and 47.61 E; green polygons in Fig. 4), exhibiting two large isolated populations, and (iii) populations inhabiting western Iran (W Iran) (from 35.94 N and 50.46 E to 35.24 N and 48.42 E; purple polygons in Fig. 4). Additionally, pairwise fixation index (FST ) confirmed a significant genetic structure among the clusters (Table 2). Divergence time Bayes factor revealed that the exponential relaxed clock model was significantly more adapted to our dataset than the lognormal relaxed (BF = 7.15), or strict (BF = 20.41) clocks. The results of BEAST based on the exponential relaxed clock model provided an estimate of the mean rate of sequence evolution of approximately 0.01 sub-

Fig. 4. Bayesian spatial clustering for groups of individuals of Iranian S. fulvus (based on 1140 cyt b) performed in BAPS. The mixture analysis was set for K = 4. Identical colors represent populations that are clustered together, concordant with Figs. 1, 2, 3, and 5. Two samples from Russia and Kazakhstan shared one polygon, with a large distance between them (blue polygons in high latitudes). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

stitutions per My. Furthermore, effective sample sizes were higher than 500.

Fig. 5. A time-calibrated species tree reconstruction of ground squirrels based on cyt b of 100 sequences (81 S. fulvus, one S. alashanicus, one S. erythrogenys, three S. xanthoprymnus, one S. citellus, one C. leucurus, one C. gunnisoni, one C. parvidens, one C. mexicanus, one C. ludovicianus, one X. perotensis, one X. spilosoma, one X. tereticaudus, one X. mohavensis, one U. parryii, one U. richardsoni, one S. major, and one S. pygmaeus sequences as outgroups) using BEAST 1.8.0 (Drummond et al., 2013). Main node values are indicated in million years before the present time (Mya) and in bold format under the bar. Bayesian posterior probabilities are given in italics above the bar of each branch. Nodes A–C show the divergence of S. fulvus in Iran, where the three stars correspond to calibration points. Colors correspond to each clade as shown in Figs. 1–4.

A. Asgharzadeh et al. / Mammalian Biology 98 (2019) 137–145

143

Table 3 Results of the demographic analysis of S. fulvus phylogroups. n: number of sequences, RI: Harpending’s raggedness index, and SSD: sum of squared deviations. Fu’s Fs and Tajima’s D in bold indicate the rejection of the hypothesis for constant population size (p < 0.05), while non-significant (non-bold numbers) SSD and RI indicate a good fit to the models of population expansion. See Phylogenetic reconstruction section for more details of the three Iranian S. fulvus clades (clades IR1, IR2, and IR3).

All Clade IR1 Clade IR2 Clade IR3

n

Fu’s Fs

Tajima’s D

SSD

RI

Mismatch distribution

79 18 28 33

−4.902 −1.992 −4.625 -5.946

−0.271 1.277 −1.287 −1.696

0.018 0.036 0.019 0.010

0.025 0.155 0.148 0.118

Multimodal Bimodal Unimodal Unimodal

(P < 0.05) can represent the presence of selective sweep or strong population growth signal (Smith and Haigh, 1974; Tajima, 1989; Fu, 1997). However, the sum of the squared deviations (SSD) and Harpending’s raggedness index (RI) failed to reveal any significant differences between the observed and simulated expected distributions in any of the tests (Table 3), indicating that none of the clades deviated from the expected distribution under the model of expansion, even though Clade IR1 exhibited a bimodal mismatch distribution. Discussion Mitochondrial genealogy and genetic structure

Fig. 6. Bimodal (clades IR1) and unimodal (clades IR2 and IR3) mismatch distributions for the Iranian Yellow ground squirrel using cyt b, determined by DnaSP version 5.0 (Librado and Rozas, 2009). Unimodal dotted lines represent populations which have experienced sudden range expansions, while bimodal is associated with a constant population size. Black solid curve shows the distributions expected under a constant population. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Our molecular clock dating based on the mtDNA dataset revealed that the divergence within the S. fulvus started in the mid-Pleistocene, for which three main diversifications could be recognized: (i) the separation of the Iranian S. fulvus from the ancestral pool in higher latitudes that occurred in the mid-Pleistocene (node A, Fig. 5) about 0.427 Mya (95% HPD 0.148–0.854), (ii) the first period of diversification corresponding to the separation of Clade IR1 from the other clades 316 Kya (95% HPD 0.118–0.595) (node B, Fig. 5), (iii) the last divergence of S. fulvus in Iran that occurred between clades IR2 and IR3 at 243 Kya (95% HPD 0.089–0.443) during the mid-Pleistocene (node C, Fig. 5). Demographic history Demographic history was deduced based on the observed mismatch distributions (Fig. 6), the neutrality tests, and the Harpending’s raggedness index (Table 3). Clades IR2 and IR3 depicted a unimodal distribution, whereas Clade IR1 represented a bimodal one. Unimodal and bimodal distributions are usually associated with sudden demographic (or spatial) expansions and populations at demographic equilibrium (Slatkin and Hudson, 1991; Rogers and Harpending, 1992; Excoffier, 2004), respectively. Also, a bimodal distribution may indicate two separate lineages. Fu’s Fs and Tajima’s D statistic tests were negative for all lineages, including the Iranian Yellow ground squirrel as a whole (Table 3), suggesting a recent population expansion (e.g., Rogers and Harpending, 1992; Slatkin and Hudson, 1991). On the other hand, the significant and negative results of Fu’s Fs and Tajima’s D

Our mitochondrial trees reconstructed using ML and BI methods revealed the same three major clades, labeled IR1-IR3 in Fig. 2, for the Iranian Yellow ground squirrel from eastern and northeastern to western Iran. Furthermore, the presence of well-clustered populations (identified by BAPS, Fig. 4) and lack of geographical overlap between haplotypes (Fig. 3) imply allopatric fragmentation (Avise, 2000). Moreover, the significant inter-population genetic differentiation between geographical populations (IR-TU-AF, NE-W Iran, and W Iran), as measured by FST values (Table 2), also supports the subdivision of the Yellow ground squirrel in Iran. Such intraspecific splits were presumably caused by major biogeographical barriers during glacial and interglacial cycles of the Pleistocene, which will be discussed in greater detail later. Both Fu’s Fs and Tajima’s D tests for detection of demographic expansion revealed significantly negative values for clades IR1, IR2, and IR3. In addition, low nucleotide diversity (pi) and high haplotype diversity (Hd) are characteristic of populations that have undergone rapid demographic expansion (Avise, 2000), a pattern also confirmed for clades IR1, IR2, and IR3 through non-significant values of the sum of square deviations and the raggedness index (Harpending, 1994) in the present study. Moreover, the star-shaped structure of our haplotype network is a further proof for this hypothesis (Slatkin and Hudson, 1991). Furthermore, the unimodal shapes of mismatch distributions for clades IR2 and IR3 corroborate rapid range expansion for populations of these lineages. However, a bimodal distribution, as observed for Clade IR1, likely indicates the presence of two distinct lineages (Alvarado-Bremer et al., 2005). It could be assumed that some populations had formerly colonized the areas between clades IR1 and IR2, but apparently became extinct not long ago, as in the case of Clade IR2 with two largely geographically distant groups. Divergence and dispersal during the Pleistocene oscillations The Earth’s climate was greatly influenced through the Last Glacial Maximum (LGM), during which glaciers were allowed to grow at lower latitudes and embrace over 30% of the earth’s surface (Kehl, 2009), leading to significant climatic fluctuations

144

A. Asgharzadeh et al. / Mammalian Biology 98 (2019) 137–145

during the Quaternary period. Indeed, Iran’s climate also experienced alterations (Blanchet et al., 1997), including dramatic falling and moderate rising temperatures during early Würm around 100 Kya (Ehlers, 1980). Moreover, it is known that the climate of Iran was most likely cold in the west and north during glacials and warm during interglacials (Kehl, 2009). Influenced by the LGM, Iran underwent five to eight degrees Celsius reduction in temperature, particularly in its central parts as well as mountains of Zagros and Alborz. Such cold climatic conditions in Iran can be inferred from the effects of glacial/interglacial cycles on Iran’s mountains and the existence of Loess sediments in its central regions (Ehlers, 1980; Yamani, 2002; Moghimi, 2008; Kehl, 2009; Moghimi, 2010). Recently, Stewart et al. (2010) adopted a longitudinal approach to how species respond to glacial/interglacial periods and found that some species might have had continental interglacial refugia. They noted that this longitudinal approach, commonly referred to as the oceanic–continental gradient, is in contrast to the latitudinal approach for northern and southern refugia, and clarifies how steppe species expanded towards the steppe-tundra biome during the Late Pleistocene. Our molecular dating results suggested that during the cold periods (before the Holsteinian temperate period at 430-350 Kya; Lowe and Walker, 1997), the ancestors of the Iranian Yellow ground squirrel expanded their distribution towards eastern Iran through cold and dry steppes of high altitudinal countries. This event was followed by population contractions into interglacial refugia around the Holsteinian temperate period (430-350 Kya; Lowe and Walker, 1997), leading to their isolation from their ancestral pool. Following that, the species dispersed from eastern to northeastern Iran along a longitudinal axis during Fuhne cold period (350-340 Kya; Lowe and Walker, 1997), concurrent with the early Riss period in northern Alpes. Subsequently, the first divergence due to vicariance among Iranian populations occurred between Iran-Turkmenistan-Afghanistan borderline (Clade IR1) and the populations of northeastern Iran (Clade IR2), estimated to have occurred around 316 Kya (95% HPD 595-118 Kya; Fig. 5), apparently at the Domnitz temperate period (340-300 Kya; Lowe and Walker, 1997) when populations were trapped in interglacial refugia. Finally, during the cold climatic conditions of the Drenthe period (300-250 Kya; Lowe and Walker, 1997), populations dispersed from the northeast to western Iran, again along the longitudinal axis, whereas the second divergence between populations (clades IR2 and IR3) occurred some 243 Kya (95% HPD 443-89 Kya; Fig. 5), when populations were confined to interglacial refugia during the Saale/Drenthe temperate period about 250-190 Kya (Lowe and Walker, 1997). These results provide proof for the role climatic oscillations play in isolation of populations. As populations of the species dispersed towards western Iran, it seems their potential routes were obstructed by high elevations of the Zagros Mountains. Thus, it could be assumed that subsequent to their colonization into Iran through eastern and northeastern routes, populations of the Yellow ground squirrel fell under the influence of glacial and interglacial cycles along a longitudinal axis. As discussed earlier, Iranian Yellow ground squirrel populations mainly showed features consistent with the recent demographic expansion in congruence with the rapid expansion of populations during glacial periods, through which populations moved towards eastern and subsequently northeastern and western Iran. This finding is in agreement with Gündüz et al.’s study (2007) on S. xanthoprymnus. On the other hand, no shared haplotypes were observed between the northeastern and western (clades IR2 and IR3) and the eastern clade (Clade IR1), suggesting that the gene flow occurred through a rather unidirectional manner from the northeast to the west of Iran, with no bidirectional gene transfer to the ancestral pool. This complies with the new approach developed by Stewart et al. (2010) which proposed that some Eurasian

mammal species such as ground squirrels Spermophilus spp. and pika Ochotona spp. that currently exhibit continentally-limited Palearctic distributions probably had expansive distributions during certain stages of the last glacial period (Musil, 1985). Our results are relatively new for Iran and in need of caution to postulate the existence of multiple continental interglacial refugia in the steppe biome of Iran during climatic oscillations of the Pleistocene, particularly because our findings merely rely on mtDNA markers. Furthermore, it is not entirely safe to propose these isolated populations as cryptic refugia. As Clade IR2 splits into two largely geographically distant groups, it might have once constituted a continuous population from northeastern to western Iran in the southern part of the Caspian Sea. Over time, central populations disappeared, either as a result of human activities or, as noted earlier, due to the species range expansions and contractions caused by the spread of grasslands and forests, respectively, therefore allowing only the eastern and western populations to persist in the area. Implications for conservation and future work Over the past half-century, poor farming practices have degraded steppe habitats, leading to severe habitat fragmentation and destruction of ground squirrel populations (Kryˇstufek et al., 2009). Also, a decreasing trend in populations of S. fulvus has been reported in some parts of its range, as the species receives no legal protection and is listed as Least Concern (LC) on the IUCN Red List (IUCN, 2019). Additionally, it is not listed as protected under the national regulations of the Department of Environment of Iran, nor are four other squirrel species that occur in Iran. Therefore, it is essential to reconsider their inclusion in the list of nationally protected species and establish practical measures to ensure their conservation. Phylogenetic studies are known to provide the foundation for conservation management of polytypic taxa, as they allow conservationists to direct their efforts towards each clade separately instead of focusing on the species as a single unit (Hung et al., 2012; Nazarizadeh et al., 2016). Such taxa, generally known as Evolutionary Significant Units (ESUs), require distinct conservation attention owing to their independent evolutionary histories (Moritz, 1994; Ryder, 1986). According to the definition by Fraser and Bernatchez, ESUs encompass lineages with extremely limited intraspecific gene flow and are therefore of conservation significance (Fraser and Bernatchez, 2001). Consequently, we proposed the three isolated lineages of S. fulvus identified in Iran as three ESUs. Clade IR1 exhibits a small range of distribution (about 15,000 km2 ), and also forms a separate evolutionary clade, highlighting the importance of its conservation. Furthermore, despite the large distribution of clades IR2 and IR3, their populations are fragmented by urban development and industrialization and are also jeopardized by human-based threats including cultivation of land, killing by farmers, road-kills, and habitat destruction. All in all, we suggest that the Department of Environment of Iran take the following conservation measures into consideration: (i) introduce safe zones for the existing protected areas to protect natural habitats of the Yellow ground squirrel; (ii) prohibit or control legal/illegal hunting of the species populations in rural areas, (iii) add the species to the protected national list of the Department of Environment of Iran. Acknowledgements We would like to thank Masoud Nazarizadeh for providing useful information to carry out the analyses and interpreting the results of this study. We are grateful to the Department of Environment

A. Asgharzadeh et al. / Mammalian Biology 98 (2019) 137–145

of Iran (Razavi Khorasan, North Khorasan, Zanjan, Qazvin, Alborz, Kurdistan, Hamedan, and west Azerbaijan provinces) and Ali Teymuri (Head of Guard Office) for providing sampling permissions and support. We thank Azadeh Hatef (University of Saskatchewan), Ali Bali, Saber Khederzadeh, Leila Mohammad Zaheri, Fatemeh Maleki, Hasan Mohammadian, Mona Samiei, Leila Ezedinloo, and Tahereh Eftekhari (Natural History Museum and Genetic Bank Bureau) for their technical help with the laboratory work. In addition, we wish to acknowledge Ali Khani, Mohsen Jahanpour, Reza Ramazani, Ali Gharcheloo, Alireza Mohammadi, Pouya Rafatpey, Kamran Almasieh, Omid Yusefi, Ali Hatami, Ali jafary, Masoud Lahut, Nikan Norouzzadeh, Arman Menbari, Babak Ghavidel, and kind environmental guards for their assistance with field sampling. Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.mambio.2019. 09.007. References ˜ Alvarado-Bremer, J.R.A., Vinas, J., Mejuto, J., Ely, B., Pla, C., 2005. Comparative phylogeography of Atlantic bluefin tuna and swordfish: the combined effects of vicariance, secondary contact, introgression, and population expansion on the regional phylogenies of two highly migratory pelagic fishes. Mol. Phylogenet. Evol. 36, 169–187. Avise, J.C., 2000. Phylogeography: the History and Formation of Species. Harvard university press. Blanchet, G., Sanlaville, P., Traboulsi, M., 1997. Le Moyen-Orient de 20 000 ans BP à 6 000 ans BP Essai de reconstitution paléoclimatique. Palorient 23, 187–196. Brandley, M.C., Schmitz, A., Reeder, T.W., 2005. Partitioned Bayesian analyses, partition choice, and the phylogenetic relationships of scincid lizards. Syst. Biol. 54, 373–390. Cheng, L., Connor, T.R., Sirén, J., Aanensen, D.M., Corander, J., 2013. Hierarchical and spatially explicit clustering of DNA sequences with BAPS software. Mol. Biol. Evol. 30, 1224–1228. Clement, M., Posada, D., Crandall, K.A., 2000. TCS: a computer program to estimate gene genealogies. Mol. Ecol. 9, 1657–1659. Drummond, A.J., Rambaut, A., Suchard, M., 2013. Beast 1.8.0. Ehlers, E., 1980. Iran: Grundzüge einer geographischen Landeskunde. Broschiert. Excoffier, L., 2004. Patterns of DNA sequence diversity and genetic structure after a range expansion: lessons from the infinite-island model. Mol. Ecol. 13, 853–864. Excoffier, L., Laval, G., Schneider, S., 2005. Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evol. Bioinf. 1, 117693430500100000. Faerman, M., Bar-Gal, G.K., Boaretto, E., Boeskorov, G.G., Dokuchaev, N.E., Ermakov, O.A., Golenishchev, F.N., Gubin, S.V., Mintz, E., Simonov, E., Surin, V.L., Titov, S.V., Zanina, O.G., Formozov, N.A., 2017. DNA analysis of a 30,000-year-old Urocitellus glacialis from northeastern Siberia reveals phylogenetic relationships between ancient and present-day arctic ground squirrels. Sci. Rep. 7, 42639. Fraser, D.J., Bernatchez, L., 2001. Adaptive evolutionary conservation: towards a unified concept for defining conservation units. Mol. Ecol. 10, 2741–2752. Fu, Y.X., 1997. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147, 915–925. Goodwin, H.T., 1995. Pliocene-Pleistocene biogeographic history of prairie dogs, genus Cynomys (Sciuridae). J. Mammal. 76, 100–122. Gündüz, I˙ ., Jaarola, M., Tez, C., Yeniyurt, C., Polly, P.D., Searle, J.B., 2007. Multigenic and morphometric differentiation of ground squirrels (Spermophilus, Scuiridae, Rodentia) in Turkey, with a description of a new species. Mol. Phylogenet. Evol. 43, 916–935. Harpending, H.C., 1994. Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Hum. Biol. 66, 591–600. Harrison, R.G., Bogdanowicz, S.M., Hoffmann, R.S., Yensen, E., Sherman, P.W., 2003. Phylogeny and evolutionary history of the ground squirrels (Rodentia: Marmotinae). J. Mamm. Evol. 10, 249–276. Helgen, K.M., Cole, F.R., Helgen, L.E., Wilson, D.E., 2009. Generic revision in the Holarctic ground squirrel genus Spermophilus. J. Mammal. 90, 270–305. Herron, M.D., Castoe, T.A., Parkinson, C.L., 2004. Sciurid phylogeny and the paraphyly of Holarctic ground squirrels (Spermophilus). Mol. Phylogenet. Evol. 31, 1015–1030. Hung, C., Drovetski, S.V., Zink, R.M., 2012. Multilocus coalescence analyses support a mtDNA-based phylogeographic history for a widespread palearctic passerine bird, Sitta europaea. Evolution 66, 2850–2864. Huson, D.H., Bryant, D., 2006. Application of phylogenetic networks in evolutionary studies. Mol. Biol. Evol. 23, 254–267.

145

IUCN, 2019. The IUCN Red List of Threatened Species. Version 2019-2 (Accessed 08 August 2019) http://www.iucnredlist.org. Kapustina, S.Y., Brandler, O.V., Adiya, Y., 2015. Phylogeny of genus Spermophilus and position of Alashan ground squirrel (Spermophilus alashanicus, Büchner, 1888) on phylogenetic tree of Paleartic short-tailed ground squirrels. Mol. Biol. 49, 391–396. Karami, M., Ghadirian, T., Faizolahi, K., 2016. The Atlas of Mammals of Iran. Jahad daneshgahi, kharazmi Branch. Kehl, M., 2009. Quaternary climate change in Iran – the state of knowledge. Erdkunde 63, 1–17. Kryˇstufek, B., Bryja, J., Buˇzan, E.V., 2009. Mitochondrial phylogeography of the European ground squirrel, Spermophilus citellus, yields evidence on refugia for steppic taxa in the southern Balkans. Heredity 103, 129–135. Kryˇstufek, B., Vohralík, V., 2012. Taxonomic revision of the Palaearctic rodents (Rodentia): Sciuridae: Xerinae 1 (Eutamias and Spermophilus). Lynx, n. s. (Praha) 43 (1–2), 17–111. Lanfear, R., Calcott, B., Ho, S.Y.W., Guindon, S., 2012. PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol. Biol. Evol. 29, 1695–1701. Leigh, J.W., Bryant, D., 2015. POPART: full-feature software for haplotype network construction. Methods Ecol. Evol. 6, 1110–1116. Librado, P., Rozas, J., 2009. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452. Lieberman, B.S., 2002. Phylogenetic biogeography with and without the fossil record: gauging the effects of extinction and paleontological incompleteness. Palaeogeogr. Palaeoclimatol. Palaeoecol. 178, 39–52. Lowe, J.J., Walker, M.J.C., 1997. Reconstructing Quaternary Environments. Longman, Harlow. Moghimi, E., 2008. Climatic Geomorphology Cold and Glacial Territory. University of Tehran Press. Moghimi, E., 2010. Geomorphology of Iran. University of Tehran Press. Moritz, C., 1994. Defining ‘evolutionarily significant units’ for conservation. Trends Ecol. Evol. 9, 373–375. Musil, R., 1985. Paleobiogeography of terrestrial communities in Europe during the Last Glacial. Acta Entomol. Musei Natl. Pragae XLI B 1, 1–83. Nazarizadeh, M., Kaboli, M., Rezaie, H.R., Harisini, J.I., Pasquet, E., 2016. Phylogenetic relationships of Eurasian Nuthatches (Sitta europaea Linnaeus, 1758) from the Alborz and Zagros Mountains, Iran. Zool. Middle East 62, 217–226. Nguyen, L.T., Schmidt, H.A., von Haeseler, A., Minh, B.Q., 2014. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum likelihood phylogenies. Mol. Biol. Evol. 32, 268–274, http://dx.doi.org/10.1093/molbev/msu300. Rambaut, A., Drummond, A.J., Available: 2007. Tracer. v. 1.5. http://beast.bio.ed.ac. uk/Tracer. Rogers, A.R., 1995. Genetic evidence for a Pleistocene population explosion. Evolution 49, 608–615. Rogers, A.R., Harpending, H., 1992. Population growth makes waves in the distribution of pairwise genetic differences. Mol. Biol. Evol. 9, 552–569. Ronquist, F., Huelsenbeck, J.P., 2003. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572–1574. Ryder, O.A., 1986. Species conservation and systematics: the dilemma of subspecies. Trends Ecol. Evol. 1, 9–10. Sambrook, J., Fritsch, E.F., Maniatis, T., 1989. Molecular Cloning: A Laboratory Manual. Cold Spring Harbor Laboratory Press. Sgro, C.M., Lowe, A.J., Hoffmann, A.A., 2011. Building evolutionary resilience for conserving biodiversity under climate change. Evol. Appl. 4, 326–337. Shimodaira, H., Hasegawa, M., 1999. Multiple comparisons of log-likelihoods with applications to phylogenetic inference. Mol. Biol. Evol. 16, 1114–1116. Slatkin, M., Hudson, R.R., 1991. Pairwise comparisons of mitochondrial DNA sequences in stable and exponentially growing populations. Genetics 129, 555–562. Smith, J.M., Haigh, J., 1974. The hitch-hiking effect of a favourable gene. Genet. Res. 23 (1), 23–35. Stewart, J.R., Lister, A.M., Barnes, I., Dalén, L., 2010. Refugia revisited: individualistic responses of species in space and time. Proc. Biol. Sci. 277, 661–671. Swofford, D.L., 2002. PAUP*. Phylogenetic Analysis Using Parsimony (* and Other Methods). Version 4. Sinauer Assoc., Sunderland, MA. Tajima, F., 1989. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123, 585–595. Tamura, K., Stecher, G., Peterson, D., Filipski, A., Kumar, S., 2013. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol. Biol. Evol. 30, 2725–2729. Thompson, J.D., Higgins, D.G., Gibson, T.J., 1994. CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 22, 4673–4680. Thorington, R.W.J., Koprowski, J.L., Steele, M.A., Whatton, J.F., 2012. Squirrels of the World. JHU Press. Xia, X., Lemey, P., 2009. Assessing substitution saturation with DAMBE. Phylogenetic Handb. a Pract. Approach to DNA and Protein Phylogeny, vol. 2., pp. 615–630. Xia, X., Xie, Z., 2001. DAMBE: software package for data analysis in molecular biology and evolution. J. Hered. 92, 371–373. Yamani, M., 2002. The geomorphology of Alamkooh glaciers. Geogr. Res. Q. 34, 1–18.