Mitochondrial DNA sequence analysis reveals multiple Pleistocene glacial refugia for the Yellow‐spotted mountain newt, Neurergus derjugini (Caudata: Salamandridae) in the mid‐Zagros range in Iran and Iraq

Abstract Phylogeography is often used to investigate the effects of glacial cycles on current genetic structure of various plant and animal species. This approach can also identify the number and location of glacial refugia as well as the recolonization routes from those refugia to the current locations. To identify the location of glacial refugia of the Yellow‐spotted mountain newt, Neurergus derjugini, we employed phylogeography patterns and genetic variability of this species by analyzing partial ND4 sequences (867 bp) of 67 specimens from 15 sampling localities from the whole species range in Iran and Iraq. Phylogenetic trees concordant with haplotype networks showed a clear genetic structure among populations as three groups corresponding to the populations in the north, center, and south. Evolutionary ages of clades north and south ranging from 0.15 to 0.17 Myr, while the oldest clade is the central clade, corresponding to 0.32 Myr. Bayesian skyline plots of population size change through time show a relatively slight increase until about 25 kyr (around the last glacial maximum) and a decline of population size about 2.5 kyr. The presence of geographically structured clades in north, center, and south sections of the species range signifies the disjunct populations that have emerged in three different refugium. This study illustrates the importance of the effect of previous glacial cycles in shaping the genetic structure of mountain species in the Zagros range. These areas are important in terms of long‐term species persistence and therefore valuable areas for conservation of biodiversity.


| INTRODUC TI ON
Quaternary glacial-interglacial periods had a significant impact on the current pattern of species distribution and interpopulation gene flow (Fedorov & Stenseth, 2002;Ikeda et al., 2017). Generally, glacial plates at higher latitude (and elevation) made habitats inappropriate to inhabit, while climate at the same location at lower latitude (and elevation) was habitable (Ohlemüller, Huntley, Normand, & Svenning, 2012). Populations at higher latitude may cope through adaptation or phenotypic plasticity, or alternatively they may track suitable habitat (Hoffmann & Sgro, 2011). However, a more likely outcome is that such populations go extinct (Stewart, Lister, Barnes, & Dalén, 2009).
At lower latitudes, populations can endure glacial periods relatively unaffected in so called glacial refugia (Bennett & Provan, 2008).
Glacial refugium refers to a location where the mountain species can survived in the glacial periods, regardless of the geographical position or spatial extent of the location (Holderegger & Thiel-Egenter, 2009). Biogeographic patterns are one approaches for identifying and describing refugia and suggesting that refugia existed in an area at some stage in the past (Keppel et al., 2012). Numerous studies are now revealing that within each of the main refugial areas, intraspecific lineages exhibited differential responses to climatic changes, often surviving in distinct, allopatric local refugia and exhibiting population fragmentation, allopatri differentiation, demographic expansion, and population admixture (Branco, Monnerot, Ferrand, & Templeton, 2002;García-París, Alcobendas, Buckley, & Wake, 2003;Mattoccia, Marta, Romano, & Sbordoni, 2011). On a broader scale, in the case of widespread species a detailed knowledge of phylogeographic patterns within refugial areas is crucial for deriving robust inferences about the evolutionary history of species at a continental scale (Gómez & Lunt, 2007). On a finer scale, this approach depicts the existence of multiple refugial areas and the succession of events of population fragmentation, allopatric differentiation, demographic expansion, and population admixture (Sequeira, Alexandrino, Rocha, Arntzen, & Ferrand, 2005).
In higher latitudes in North American and Northern Europe, the influence of Quaternary climatic fluctuations on plant and animal species is well documented (Hewitt, 2011); however, such influence of climate change on range dynamics of species in other parts of the world is much less reported (Jandzik et al., 2018). Various geomorphic studies in Iran, including lake sediments, deserts, glacial moraines and periglacial features, salt domes, alluvial sediments, pediments and alluvial fans, fluvial and marine terraces, and loesssoil sequences have shown that Iran was affected by Pleistocene and Holocene climate fluctuations (Kehl, 2009). The available paleoecological and palynological records in Iran indicate that the last glacial maximum (LGM) was characterized by a cooler and more arid climate compared to the Holocene (Djamali et al., 2012;Kehl, 2009). At the present time, deep valleys of several high altitude mountain areas in the Iranian Plateau are covered by glaciers. In the LGM, the expansion of such climatic conditions encouraged the spread of steppe habitats in lower altitudes, while high altitudes in the mountain areas were covered by glaciers (Ghahroudi Tali, Naeimi, & Gharnaie, 2016). In Zagros range, there are reports of isolated populations of various temperate evergreen plant species such as Olea europaea and Myrtus communis (Sharifi, Najafi, Yossefshahi, & Hemmati, 2001) and pockets of plant communities such as mires and patterned mires (Wasylikowa & Walanus, 2004), floating meadow and highland peatlands (Sharifi, Rezaii, Hosseini, & Raji, 2004). These relict individuals and plant communities are assumed to persist in enclaves of environmental conditions within an inhospitable regional climate and left behind during range shifts caused by Pleistocene climatic fluctuations (Hampe & Jump, 2011;Wright, McAndrews, & Van Zeist, 1967). According to Djamali et al. (2008), the vegetation composition during the LGM in the Lake Urmia, in northwestern Iran, demonstrates lower winter temperatures than today, and higher July temperatures about 11-12°C (Djamali et al., 2008).
Neurergus derjugini is a critically endangered species occurring in 42 highland streams of the Zagros Mountain range in western Iran and eastern Iraq (Afroosheh, Akmali, Esmaeili-Rineh, & Sharifi, 2016 (Afroosheh et al., 2016). The extent of occurrence of N. derjugini as indicated by a minimum convex polygon is 6,366 km 2 for the 42 known localities (Afroosheh et al., 2016). This polygon is positioned along the western edge of the Zagros range with elevations ranging from 630 to 2057 masl. This area is covered by an oak open woodland. The breeding habitat of N. derjugini in the Zagros range has been degraded recently by water pollution, water extraction, and severe droughts, which have led to the extirpation of some populations (Afroosheh et al., 2016).
The present study analyzed spatial variation of mtDNA sequences, aiming to investigate the genetic structure within species, as well as attempting to trace the recent demographic and phylogeographical history of the N. derjugini in the Zagros range. We tested the multiple glacial refugia hypothesis for N. derjugini using molecular data in the whole specie range in mid-Zagros range at the border of Iran and Iraq. We therefore aim to: (a) examine the extent of in-

| Population sampling and sequencing
We were sampled from 67 individuals of Neurergus derjugini from 15 localities in Iran and Iraq during 2012-2014 (Table 1, Figure 1).
Tissue samples were collected from each individual by removing a small section of tail or toe using sterile equipment. Tissue samples were stored in 95% ethanol at −20°C until DNA was extracted.
Genomic DNA was extracted using tissue kit GenNetBio TM . We sequenced an 867 bp fragment of a mtDNA generation consisting of 726 bp fragment of NADH subunit 4, the whole tRNA-His and tRNA-Ser and 21 bp from the "5′" end of tRNA-Leu using the primers ND4-Leu: ND4, CAC CTA TGA CTA CCA AAA GCT CAT GTA GAA GC and Leu, CAT TAC TTT TAC TTG GAT TTG CAC CA (Arevalo, Davis, & Sites, 1994). Amplification of mtDNA was conducted using denaturation at 94°C for 2 min, 58°C for 45 s, 72°C for 2 min followed by 35 cycles at 94°C for 30 s, 58°C for 45 s, 72°C for 60 s for annealing and extension at 72°C for 3 min. Sequencing was performed by Microsynth Switzerland Laboratories.
We used ArcMap 10.3 to process variables. The matrix of environmental distances was computed by SPSS version 16.0. We computed geographical distances using DIVA-GIS v 7.5.0 (Hijmans, Guarino, & Mathur, 2004). We performed a Mantel test to evaluate correlations between genetic, geographical, and environmental distances using Arlequin 3.0 (Excoffier et al., 2005). This analysis runs with 10,000 random permutations. In addition, a three-way Mantel test was applied between the matrix of pairwise genetic differentiation and the matrix of environmental distances and geographical distances among populations. We also performed principal components analysis (PCA) to investigate ecological differentiation within N. derjugini distribution range.

| Demographic analysis
Mismatch distribution analysis was performed using DnaSP v 5.10.01 software (Rozas et al., 2010)

| Divergence time estimate
We estimate times of divergence between lineages of N. derjugini using BEAST v 2.4.5 (Bouckaert et al., 2014). We ran for 30,000,000 generations sampled every 1,000 steps, with the first 3,000,000 generations regarded as burn-in 95% highest posterior density (95% HPD) were considered significant support. The maximum clade credibility tree summarized by Tree Annotator v1.8.4 program (Drummond & Rambaut, 2007

| Genetic variation
We   The statistical parsimony haplotype network (Figure 4) suggested three haplotype subnetworks, which were in agreement with the topology described in the phylogenetic trees. The haplotype network pattern suggestive of little or no gene flow between regions because northern, central, and southern haplogroups were not intermingled.

| Population analysis
The AMOVA analysis showed most genetic variation was significantly explained by differences among regions (74.87% genetic variation, F CT = 0.74, p < .01), and genetic difference was detected among populations within regions (19.07 genetic variation, F SC = 0.75, p < .001; Table 4). The genetic differences were significant among populations without the grouping. Most of the diversity was observed among populations (92.16%) and among regions (74.78%), while low percentage of variance was detected within populations (6.06%-7.84%; Table 4 (Table 5).

| Demographic analysis
The mismatch distribution for the N. derjugini appeared smooth and unimodal, consistent with a model of population expansion ( Figure 5).
The Harpending's Raggedness Index (r = .08, p(r) = .1) and sum of squared deviations (SSD) (r = .03, p(r) = .1) were not significant, indicating population expansion. Results of Tajima's D test was positive and not significant (0.60, P(r) = .10). However, negative value (−0.63, P Fs = 0.04) for Fu's Fs, which is more sensitive to recent population expansions, suggested the recent demographic expansion. Bayesian skyline plots suggest that population size was increased very slightly 25 kyr (around LGM), (Figure 6). BSP s showed a decline of population size about 2.5 kyr. Thus, the BSP s suggest that populations went through a bottleneck in the cold period around 2.5 kyr.

| Environmental data
The analyzed populations revealed significant levels of isolation by distance (r = .66, p < .0001; Figure 7a) and correlation between genetic divergence and environmental distance (r = .68, p < .0001; Figure 7b). Mantel test detected a signal of isolation by environment had significant influence on genetic distance. Each of two variables in isolation had a significant effect on genetic differentiation. In the three-way Mantel test, the correlation between genetic and geographical (r = .36, p = .0068) and environmental distances remained significant (r = .30, p = .0069). TA B L E 3 K2P genetic distances values among 15 populations of Neurergus derjugini. Numbers are representative of localities as indicated in Table 1 No The first two principal components (PC1 and PC2) explained 72.51% and 22.17% of the total variation and did not separate the northern, central, and southern localities along temperature and precipitation gradients (Table 6, Figure 8). The occurrence sites were not divided into separated environmental spaces in the Cartesian coordinates formed by the first two principal components (Figure 8).

| Divergence time estimate
Dating of major mtDNA cladogenetic events within N. derjugini haplogroups is presented in Figure 9.

| D ISCUSS I ON
The Yellow-spotted mountain newt exhibited a complex phylogeographic pattern with multiple divergent mtDNA clades across its relatively small range (6,366 km 2 ). A total of 3 allopatric and 11 parapatric clades were identified for ND4, most of them with a very  (Schneider & Schneider, 2011). However, use of this feature in taxonomy of the Yellow-spotted newts in Middle East has caused uncertainty because more than one species of newt in Iran and neighboring Iraq and Turkey have yellow spots (Schneider & Schneider, 2011;Sharifi, Naderi, & Hashemi, 2013). Although a taxonomic revision is beyond the scope of this study, our results indicate that the current intraspecific taxonomy  Table 1 No 1  of "southern genetic richness and northern purity" (Fritz et al., 2012(Fritz et al., , 2008(Fritz et al., , 2009). Contrarily to this configuration, the pattern of distribution of genetic diversity in N. derjugini is an exception. The presence of several monophyletic groups indicate the existence of multiple refugia and has been used for explanation of genetic structure in many amphibian species (Giovannotti, Nisi-Cerioni, & Caputo, 2010). The pattern of haplogroups in present study is, also, confirmed by AMOVA analysis, showing a pronounced geographical partition of N. derjugini genetic variation.
One reason for a highly structured genetic diversity in populations occurring in close distances may be high fidelity to small home ranges in breeding streams reported for many species of

ACK N OWLED G M ENTS
This work was supported by the Razi University authorities, Kermanshah as a part of a PhD research project.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no competing interests.
F I G U R E 9 Chronogram of diversification implemented in BEAST based on ND4 gene for N. derjugini. The divergence times with 95% highest posterior density (95% HPD) are presented below the branches

AUTH O R S ' CO NTR I B UTI O N S
SV and MSH conceived the main idea with input from MM and SV collected the data; and MM and SV analyzed the data with support from MSH. Sampling was done by SV and MSH. The manuscript was written by SV, MSH, and MM. All authors made substantial contributions to the interpretation of results and the editing of the manuscript.

E TH I C S A PPROVA L A N D CO N S E NT TO PA RTI CI PATE
This study complied with the appropriate institutional, national, and international guidelines.

DATA AVA I L A B I L I T Y S TAT E M E N T
MtDNA haplotypes used in this study were deposited in the NCBI