Genetic Diversity, Population Structure and Demographic History of Brandt and Long Eared Hedgehogs in Iran

Ehsan Kashani1, Hamid Reza Rezaei2*, Morteza Naderi3 and Nematolah Khorasani4

1Department of Environmental Sciences, Faculty of Natural Resources and Environment, Science and Research Branch, Islamic Azad University P.O. Box 14515-775, Tehran, Iran

2Department of Environment, Faculty of Fisheries and Environmental Science, Gorgan University of Agricultural Sciences and Natural Resources, Iran

3Department of Environmental Science, Faculty of Agriculture and Natural Resources, Arak University, Arak, Iran

ABSTRACT

Knowledge about the genetic diversity of the mammalian species and their plasticity in ever-changing environments is emphasized by scientists and conservation biologists. Central deserts and steppes of Iran host a different kind of insectivores like hedgehogs; however, their ecology and genetic properties are rarely documented. For species of hedgehogs that are included in the subfamily Erinaceinae, at least four known species are listed for Iran. Brandt’s (Paraechinus hypomelas) and Long-eared (Hemiechinus auritus) hedgehogs were two candidates in which we tried to investigate their molecular peculiarities in central steppes and deserts of Iran. Mitochondrial evidence indicates rapid population expansion and bottleneck event for both species during the past. Phylogenetic analysis successfully grouped the specimens in two distinct clades without any common haplotypes. Long-eared specimens grouped in two distinct phylogenetic groups which belong to two different geographic areas while for P. hypomelas no significant divergence happened in the species life history.

Article Information

Received 03 August 2019

Revised 25 December 2019

Accepted 14 February 2020

Available online 22 February 2021

Authors’ Contribution

NK formulated the hypothesis and designed the field projects. Field sampling was done by EK. Manuscript was drafted and edited by HRR and MN. All authors equally contributed in bioinformatics and data analysis.

Key words

Conservation genetics, Geographical distribution, Hemiechinus auritus, Paraechinus hypomelas, Phylogenetic, Phyloegeography

DOI: https://dx.doi.org/10.17582/journal.pjz/20190803220849

* Corresponding author: [email protected]

0030-9923/2021/0002-0643 $ 9.00/0

Copyright 2021 Zoological Society of Pakistan

INTRODUCTION

As the oldest known living placental mammals, hedgehogs split from their sister groups around 70 million years ago (He et al., 2012) which diversified and dispersed to the areas mostly in Eurasia and Africa (Corbet, 1988). Their activity pattern usually depends on air temperature and food availability (Feldhamer et al., 2003). Subfamily Erinaceinae members constitute an enigmatic group that their genetic related peculiarities are under debates by the evolutionary biologists. For the species which are included in this subfamily, at least four known species (Corbet, 1988; Frost et al., 1991; He et al., 2012) are listed for Iran with no documented molecular investigation focusing on these species.

Brandt’s (Paraechinus hypomelas) and long-eared (Hemiechinus auritus) hedgehogs occupy vast areas of the western Palearctic biogeographic region where substantial parts of Iran are also located (Fig. 1). Brandt’s hedgehog has been reported only from Tunb Island in southern Iran and south eastern parts of the country (Lay, 1967). Southern White-breasted hedgehog (Erinaceus concolor) has been reported from the north and north western area of Iran (Lay, 1967). These parapatric species have a distribution range overlap in some parts of the country, primarily in the central, eastern, and southeastern parts. The comparison of the genetic diversity and phylogenetic stats of these species where they can be found in the same habitats can yield interesting results.

 

Few genetic studies that have been carried out so far are limited to hedgehog species in Central and Eastern Europe (Seddon et al., 2001; Berggren et al., 2005), and there is little information on the genetic diversity and structure of these species in Asia, especially in Iran. In this research we aimed to evaluate and compare the species’ genetic diversity and structure based on cytb marker in their contact zone in Iran.

MATERIAL AND METHODS

Sampling and laboratory procedures

The eastern and southeastern desert steppes of Iran host both of our target species (Fig. 2) Major climatic regime of the region is arid to semi-arid, and the highest rainfall occurs in the winter and autumn seasons. A total of 76 samples belong to both species (34 samples of P. hypomelas, and 42 samples from Hemiechinus auritus) were collected using noninvasive methods like sampling from hairs, ear pad or saliva. The geographic coordinates of the sampling points were recorded using GPS.

 

The DNA was extracted using standard phenol-chloroform method. The following primers were used to sequence Cytb: TCGTGACCCATGATATGAAAAATC and H-eri: GAATATCAACTTTG GGTGTTGATGGTA (Bannikova et al., 2014). The amplification of genetic loci and polymerase chain reaction (PCR) was performed using a 2X Taq PerMix kit in a volume of 20 μl. The PCR conditions included 10 μl of single-polymerase enzyme, 1μl of DNA, 0.4 μM of mix of five pM primer, 200 μM of dNTP, 1 μl of MgCl2 and doubled-distilled water to reach 20-μl volume. The PCR program of 5 minutes at 95 °C, followed by 35 cycles, including 45 seconds at 94 °C, 45 seconds at 58 °C, one minute at 72 °C, and the final expansion with 72 °C for 10 minutes was followed. The observation of amplified parts was performed using 2% agarose gel electrophoresis. Finally, the PCR product was sequenced using the 3130xl Genetic Analyzer device.

Data analysis

Sequences were edited for correction with SeqScape v2.6 (Applied Biosystems). All new sequences have been submitted to GenBank then the BioEdit 7.2.5 and the Clustalw algorithm was used to align sequences (Thompson, et al., 1994). The genetic diversity indices, number of haplotypes (H), polymorphic sites (S), haplotype diversity (Hd) (Nei, 1987), and nucleotide diversity (π) (Lynch and Crease, 1990) were calculated using DNAsp 5 (Librado and Rozas, 2009). Totally, 37 sequences belonging to the Erinaceinae subfamily were downloaded from NCBI data base (Bannikova, et al., 2014; He et al., 2012). The nucleotide composition and genetic distances between different hedgehog species were measured using unmatched genetic distances with 1000 bootstrap in MEGA 5.

Demographic history and phylogenetic analysis

Demographic events like bottleneck and population expansion were evaluated using various methods. In the first step, Tajima (1989D), Fu’s Fs (Fu, 1997), and Fu and Li’s F * indices were calculated in DnaSP 5.10 using 1000 permutations (Rozas et al., 2003). The Negative values of these indices indicate the expansion of the species population in the past. In addition, the raggedness index (rg) and the population expansion parameter tau (t) were calculated in Arlequin 3.0 (Excoffier et al., 2005). Demographic history and probable population expansion of each species in the past were estimated using the distribution of the pairwise differences in the sequences by using the mismatch distribution test (Schneider and Excoffier, 1999) in Arlequin 3.0. The quality of fit between the observed and nucleotide differences in the obtained sequences was evaluated using 1000 replications, bootstrap test and the sum of squared deviations (SSD).

Furthermore, to estimate the population size more accurately, Coalescent-based Bayesian skyline plots were plotted in BEAST 1.4.8 to consider the molecular clock (Drummond and Rambaut, 2007). The evolutionary model of the sequences was determined using the software J Model Test 2.1.5 and the Akaike index (Darriba et al., 2012). This method uses Marco Chain (MCMC) to estimate the posterior distribution of the effective population size. The Marco process was replicated three times for each species with 5000000 replications while saving each 1000 replications. The results were then combined in the Log Combiner environment. The phylogenetic relationship of the haplotypes was investigated using a neighbor net, an uncorrected patristic distance and bootstrap analysis with 1000 replications in SPLITSTREE v4.6 (Huson and Bryant, 2005).

The best nucleotide substitution model was determined based on AIC criteria using the J Model Test 2.1.5 software (Darriba et al., 2012). Then, according to the selected model, Mr Bayes 3.2.2 software was used to draw the Bayesian tree (Huelsenbeck and Ronquist, 2001). Bayesian analyses were performed based on the Markov Chain Monte Carlo with five million repetitions, burn-in = 500,000, and a sampling frequency of 500. Finally, the plotted trees were drawn by a consensus majority-rule. The concept of the “barcoding gap” was used to determine the threshold for the separation of Erinaceinae subfamily species as well as to evaluate the separation level of the two genetic groups obtained from H. auritus. The genetic distance between H. auritus and other species of this subfamily was calculated using uncorrected p-distance in MEGA v.5 (Tamura et al., 2013). Using Structure software, we compared our data with those of deposited in Gene bank for species including Atelerix algirus, A. albiventris, A. frontalis, Paraechinus aethiopicus, P. hypomelas, Mesechinus hughi, M. dauuricus, Erinaceus roumanicus, E. europaeus, E. concolor, E. amurensis.

RESULTS

Genetic diversity

Haplotype diversity of P. hypomelas specimens was higher than H. auritus. Evaluation of genetic variation showed that among the 34 sequences with the length of 1121bp obtained, there were 24 haplotypes with a mean genetic distance of 0.115 ± 0.010 and 39 polymorphic sites, of which 17 were the most parsimonious sites. Haplotype and nucleotide diversity in H. auritus (17 haplotypes) were calculated as Hd = 0.922 ± 021 and Pi = 0.00497 ± 0.0043, respectively, indicating a high level of nucleotide diversity in this species compared to P. hypomelas. There were also 30 polymorphic sites in H. auritus sequences that included the 15 most parsimonious sites. Most of the haplotypes identified in the P. hypomelas had a lower frequency in the samples, and the frequency of haplotypes in this species varied from 2.94 to 8.2. We found that 16 of 24 haplotypes identified in the P. hypomelas were unique. By contrast, the frequency of haplotypes identified in H. auritus was greater, ranging from 2.38 to 14.29, of which 10 were identified as unique haplotypes.

The evaluation of the genetic structure of H. auritus and P. hypomelas based on the Bayesian algorithm showed the sequences of the two species are separated in two genetically distinct groups with no overlap (Fig. 2a). Our results showed that samples belonged to H. auritus were clustered in two distinct groups and coincide with the haplotype network and phylogenetic tree (Fig. 2b). These two groups consisted of samples from southeast (group G1), center and east of the country (group G2). Also, the results of the evaluation of the genetic structure of P. hypomelas samples showed that there were four structures in this species in the studied region (Fig. 2c). Despite the expectation of the existence of four genetic clusters existence, these genetic groups did not show a specific geographical distribution pattern. On the other hand, many samples of this species belong to more than one group (Qi <0.8). Therefore, it can be concluded that the amount of admixture between the genetic structures of these species is high and this causes a mosaic geographical pattern in the genetic groups identified in the studied area.

The results of the assignment of samples based on the Geneland model confirmed the clustering results of the Structure software. The Geneland model represents three mosaic-based genetic groups in H. hypomelas (Fig. 3a). As with STRUCTURE results, the distribution of these groups did not show a specific geographic pattern. H. auritus samples were divided into two distinct geographic groups which constitute samples from southeastern parts of the country (group G1, gray color in Fig. 3b) and samples from central and eastern regions (G2 group, green color in Fig. 3b).

 

Demographic history and population expansion

Negative values were found in neutrality tests (TD, Fu’s Fs, Fu and Li’s F*) for H. auritus (P>0.05) and H. hypomelas (P<0.05). The population expansion index of H. auritus (t = 6.705) was greater than H. hypomelas (t = 3.742) (Table I). Our results indicated that the occurrence of sudden expansion of the population in the past among these species is possible which also confirmed by mismatch distribution test for P. hypomelas (SSD= 0.0033, r= 0.021, P> = 0.05) and H. auritus (SSD= 0.0038, r= 0.012, P> = 0.05). On the other hand, the unimodal form of the mean observed and expected mismatch between any two sequences, derived from the rapid expansion model in the past, confirms the occurrence of the population expansion in the both species (Fig. 4). Bayesian skyline plots regarding the demographic history (Fig. 5 a, b) of the species confirm population expansion in both species, especially H. auritus.

 

Phylogenetic analysis

For phylogenetic analysis, we used to show that the TrN+I+G model based on BIC output as the best nucleotide evolution model for the obtained sequences. Mean genetic distances between two clusters of H. auritus and other subfamily members is shown in Table II. Haplotype network indicates that our target species grouped in different clusters by considerable amount of mutational steps (Fig. 6a). Bayesian analysis divided H. auritus specimens to two distinct sister clusters which exactly coincide with their isolated geographic regions including southeastern (Iran_G1), eastern and center Iran (Iran_G2; Figs. 6b and 7). Phylogenetic tree indicates that Turkmenia Long-eared hedgehogs are closets group to Iranian samples.

 

 

Table I. Genetic diversity and population expansion indices for two species of H. auritus and P. hypomelas based on the sequences of cytochrome b.

	N	Nh	t	rg	Fu and Li's F*	Fu's Fs	TD	h ± s.d.	± s.d.π
H. auritus	42	17	6.705	0.012	-1.86812	-2.965	-0.80961	0.913 ±0.024	0.00481 ±0.00046
P. hypomelas	34	24	3.742	0.021	-2.58197*	-15.737*	-1.88228*	0.979 ±0.012	0.00423 ±0.00045

 

D, Fu and Li’s F*; Fu’s Fs; h, haplotype diversity; N, number of individuals; Nh, number of haplotypes; p, nucleotide diversity; rg, raggedness index; s.d., standard deviation; t, time from expansion in mutation units TD, Tajima’s D.

 

Table II. Uncorrected P_distance for the cytochrome b between the target species. The bold numbers show the standard deviation. G1 and G2 show the two distinct genetic structures of H. auritus in the Structure software.

	A. alg-irus	A. albiv-entris	A. fron-talis	M. hu-ghi	M. dauu-ricus	E. rouma-nicus	E. europ-aeus	E. conc-olor	E. amur-ensis	H. aur-itus	H. auritus (Iran-_G1)	P._hypo-melas (Iran)
A. algirus		0.01	0.01	0.01	0.01	0.01	0.01	0.01	0.01	0.01	0.01	0.01
A. albiventris	0.11		0.01	0.01	0.01	0.01	0.01	0.01	0.01	0.01	0.01	0.01
A. frontalis	0.13	0.13		0.01	0.01	0.01	0.01	0.01	0.01	0.01	0.01	0.01
P. aethiopicus	0.21	0.19	0.20	0.01	0.01	0.01	0.01	0.01	0.01	0.01	0.01	0.01
M. hughi	0.20	0.20	0.19		0.01	0.01	0.01	0.01	0.01	0.01	0.01	0.01
M. dauuricus	0.21	0.20	0.19	0.07		0.01	0.01	0.01	0.01	0.01	0.01	0.01
E. roumanicus	0.20	0.19	0.21	0.19	0.19		0.00	0.01	0.01	0.01	0.01	0.01
E. europaeus	0.21	0.20	0.21	0.20	0.19	0.07		0.01	0.01	0.01	0.01	0.01
E. concolor	0.21	0.21	0.21	0.18	0.18	0.10	0.13		0.01	0.01	0.01	0.01
E. amurensis	0.21	0.19	0.19	0.18	0.18	0.12	0.11	0.13		0.01	0.01	0.01
H. auritus	0.21	0.21	0.21	0.11	0.11	0.19	0.19	0.18	0.19		0.00	0.01
P. hypomelas	0.19	0.21	0.20	0.19	0.18	0.19	0.19	0.19	0.19	0.18	0.01	0.00
H. auritus (Iran_G2)	0.21	0.21	0.20	0.11	0.11	0.19	0.19	0.19	0.18	0.03	0.00	0.01
H. auritus (Iran_G1)	0.21	0.21	0.20	0.11	0.11	0.19	0.19	0.19	0.18	0.03		0.01
P. hypomelas (Iran)	0.19	0.21	0.20	0.19	0.18	0.19	0.19	0.19	0.19	0.18	0.18

 

DISCUSSION

Knowledge of genetic diversity, population structure, and demographic history of species has an important role in assessing adaptation, population survival, establishment of management (MUs), significant evolutionary units (ESUs) and the development of effective management plans. Topographic variation due to the existence of the Zagros Mountains created diverse habitats and landscapes as well as barriers that block gene flow among the animals with low dispersal abilities. Such pattern of genetic isolation and variation among the mammalian species has been already documented (Naderi et al., 2014, Kashani et al., 2019). The analysis of cytochrome b in two species of H. auritus and P. hypomelas shows (1) high haplotype diversity in both species, especially in P. hypomelas, (2) genetic structure in both species, especially in the H. auritus in the west to east of the study area, and (3) the occurrence of sudden demographic expansion in both species, especially in P. hypomelas. Results confirmed the importance of future studies on the determination of taxonomic units and species boundaries in this landscape. Regarding the obtained genetic distance between long eared hedgehogs clusters which is around 1 percent, more data can be useful in the future.

Genetic variation and structure in H. auritus and P. hypomelas

Genetic diversity is affected by historical factors, and low rates of evolution in the mitochondrial genome (Nabholz et al., 2008). Nucleotide diversity is one of the most important indicators of genetic diversity in a population, and the low level of this index can be indicative of a low effective population size. Although the genetic variation in the two studied species was partly similar, the haplotype diversity of H. auritus was lower than P. hypomelas

 

despite the larger number of the former one and 17 haplotypes were seen in 42 individuals (compared to P. hypomelas which had 24 haplotypes in 34 samples). In contrast, greater nucleotide diversity was observed in H. auritus (0.00481±0.00046) in comparison to P. hypomelas (0.00423±0.00045). The observed haplotype diversity in two species of H. auritus (Hd= 0.922 ± 02) and P. hypomelas (Hd= 0.979 ± 012) was more than the variation observed in other studies on hedgehogs in the world (Bolfíková and Hulva, 2012; Stefanović, et al., 2016; Djan et al., 2017; Derouiche et al., 2017). Despite the population expansion, the high genetic diversity in the studied species can be due to the presence of small populations and the mixing of these populations after the spread of populations. Lower nucleotide diversity in P. hypomelas, as compared to H. auritus, may also be due to a more severe bottleneck phenomenon, resulting in a greater effect of the founder effect on this species. The rapid rate of population expansion has a greater impact on founder effect, resulting in a greater loss of genetic diversity in this species.

The genetic structure in the two studied species is not surprising, given the low mobility of hedgehogs, the demographical events occurred in these species, as well as the shape of the land changes in this part of Iran. Because the genetic structure of the populations is not necessarily based on the degree of geographic proximity, any predefined number of populations can lead to more or less accurate estimates of the actual amount of genetic structures. Therefore, in this study, there was no prediction regarding the pre-grouping of samples. Genetic structure in P. hypomelas represents three mosaic-like genetic structures in this species and no significant geographic pattern was observed among the genetic structures. In contrast, the genetic structure of H. auritus showed a geographical pattern that strengthens the hypothesis of limiting the gene flow between the populations as well as the time for development of the genetic structure in this species. The absence of shared haplotypes between the two genetic groups in H. auritus is indicative of a lack of gene flow between the two populations for long periods.

The existence of distinct genetic structures can be attributed to their geographic isolation with relatively no gene flow between them. In addition to geographical distance between the populations of H. auritus in the west of the study area with eastern and southeastern populations, human disturbances in their habitat, land use changes (like urban development and road network) and lack of appropriate habitat patches among identified populations can be other reasons for the genetic structuring in this species. Besides, the occurrence of bottleneck events in the populations of this species in the past may have caused great differences in the frequency of haplotypes in different populations, which also affects on their genetic structure. Colonization of H. auritus from the east to the west can also result in genetic structuring and the existence of different allelic abundances between the populations of this species after population expansion. On the other hand, a more distinct genetic structure in H. auritus may be indicative of sex-based dispersal patterns in this species compared to P. hypomelas as reported before for hedgehogs. For example, studies on the European hedgehogs (E. europaeus) showed that their males usually look for a larger territory and more displacement, especially in mating seasons (Reeve, 1982). Considering the limited dispersal ability of hedgehogs (Reeve, 1982), as well as the fragmentation of their habitat in the studied area, especially in areas affected by human activities, the genetic structuring and occurrence of the Wahlund event (Baker et al., 2017) is probable in this species. Low mobility, as well as relatively high birth rates in hedgehogs can accelerate the effects of contemporary changes in their habitat (such as habitat destruction and fragmentation) leading to genetic diversity and structuring of these species.

The haplotype network and the phylogenic tree grouped H. auritus sequences in two distinct clades, which corresponded to the Bayesian genetic structure. Also, none of the sequences obtained in this study showed an unexpected phylogeny position, and the separation between different species of the Ericaniane was well established and confirmed the results of previous studies of the separation between the species of this subfamily (Corbet, 1988; Grenyer and Purvis, 2003; Bannikova et al., 2014).

The results of demographic analysis indicated instability and increased population size in both species, especially in P. hypomelas. The population expansion has been reported in many species of hedgehogs such as E. roumanicus (Bolfíková and Hulva, 2012; Djan et al., 2017). First, the unimodal distribution of the pairwise differences between sequences indicates bottleneck and population expansion in the studied species, especially in P. hypomelas. On the other hand, a lower nucleotide diversity in P. hypomelas compared to H. auritus, as well as the negative and insignificant Fu’s Fs and Fu and Li’s F in P. hypomelas can be a reason for the occurrence of bottleneck in the past in this species leading to population expansion with small initial effective population size (Avise, 2000; Low et al., 2014; Tajima, 1989; Zhou et al., 2013). The sudden population expansion may increase the number of rare haplotypes due to the high probability of new mutations in haplotypes compared to the removal of haplotypes as a result of genetic drift (Stopar et al., 2010). The small number of mutations between the haplotypes identified in P. hypomelas on one hand and the discrepancy between the genetic structures observed in this species and the geographical regions on the other hand, indicate a recent divergence between these structures and the same origin of all genetic structures identified in this species. The interspecies relations in the sympatric zones have led to the acceleration of speciation processes in the secondary contact zone of two species. This is a consequence of the development of adaptation in order to prevent the occurrence of mating between related species (Bolfíková and Hulva, 2012).

Consequently, it can be concluded that genetic structure in two species of H. auritus and P. hypomelas can be attributed to bottleneck, patterns of colonization, and sex-based dispersal. The high sensitivity of hedgehogs to climate changes and their low motility (Fowler and Racey, 1987) makes it necessary to protect their genetic diversity and gene flow among populations of these species. On the other hand, severe habitat degradation, road network and climate changes can affect the gene flow among populations leading to contemporary micro-scale genetic structuring in these species sensitive to climate change and human interference (Bannikova et al., 2014). Habitat fragmentation in the studied area causes more isolation between the two species and increases the genetic drift in the remaining populations. Therefore, a thorough understanding of the impact of these factors on spatial patterns of genetic variation requires a wider sampling of larger geographic scale and use of other molecular markers.

CONCLUSION

In the present study, population scenarios were only performed on two species of hedgehogs based on mitochondrial sequencing. The mitochondrial marker is very sensitive to the effective population size and provides only information on the maternal line. Hence, further genetic studies are needed by using molecular markers containing useful information such as microsatellites to conduct genetic analysis on a small spatial scale. Such studies play an important role in assessing effects of human interferences as well as habitat fragmentation on genetic variation, gene flow, and sex-based dispersal.

Statement of conflict of interest

The authors have declared no conflicts of interest.

REFERENCES

Avise, J., 2000. Phylogeography: The history and formation of species. Harvard University Press, Cambridge, Massachusetts and London.

Baker, K.H., Gray, H.W.I., Ramovs, V., Mertzanidou, D., Akın Pekşen, Ç., Bilgin, C.C., Sykes, N. and Hoelzel, A.R., 2017. Strong population structure in a species manipulated by humans since the Neolithic: The European fallow deer (Dama dama dama). Heredity, 119: 16–26. https://doi.org/10.1038/hdy.2017.11

Bannikova, A.A., Lebedev, V.S., Abramov, A.V. and Rozhnov, V.V., 2014. Contrasting evolutionary history of hedgehogs and gymnures (Mammalia: Erinaceomorpha) as inferred from a multigene study. Biol. J. Linn Soc., 112: 499-519. https://doi.org/10.1111/bij.12299

Berggren, K., Ellegren, H., Hewitt, G. and Seddon, J., 2005. Understanding the phylogeographic patterns of European hedgehogs, Erinaceus concolor and E. europaeus using the MHC. Heredity, 95: 84-90. https://doi.org/10.1038/sj.hdy.6800694

Bolfíková, B. and Hulva, P., 2012. Microevolution of sympatry: landscape genetics of hedgehogs Erinaceus europaeus and E. roumanicus in Central Europe. Heredity, 108: 248-255. https://doi.org/10.1038/hdy.2011.67

Corbet, G.B., 1988. The family Erinaceidae: a synthesis of its taxonomy, phylogeny, ecology and zoogeography. Mammal Rev., 18: 117-172. https://doi.org/10.1111/j.1365-2907.1988.tb00082.x

Darriba, D., Taboada, G.L., Doallo, R. and Posada, D., 2012. jModelTest 2: more models, new heuristics and parallel computing. Nat. Methods, 9: 772-780. https://doi.org/10.1038/nmeth.2109

Derouiche, L., Vercammen, P., Bouhadad, R. and Fernandes, C., 2017. Genetic evidence supporting the taxonomic separation of the Arabian and Northwest African subspecies of the desert hedgehog (Paraechinus aethiopicus). Gene, 620: 54-65. https://doi.org/10.1016/j.gene.2017.04.009

Djan, M., Stefanovic, M., Velickovic, N., Cosic, N. and Cirovic, D., 2017. Mitochondrial genetic diversity and structuring of northern white-breasted hedgehogs from the Central Balkans. Turk. J. Zool., 41: 774-782. https://doi.org/10.3906/zoo-1603-23

Drummond, A.J. and Rambaut, A., 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol., 7: 214-220. https://doi.org/10.1186/1471-2148-7-214

Excoffier, L., Laval, G. and Schneider, S., 2005. Arlequin (version 3.0): An integrated software package for population genetics data analysis. Evol. Bioinform., 1: 47-50. https://doi.org/10.1177/117693430500100003

Feldhamer, G.A., Thompson, B.C. and Chapman, J.A., 2003. Wild mammals of North America: Biology, management and conservation. The Johns Hopkins University Press, 2nd edition

Fowler, P.A. and Racey, P.A., 1987. Relationship between body and testis temperatures in the European hedgehog, Erinaceus europaeus, during hibernation and sexual reactivation. J. Reprod. Fert., 81: 567-573. https://doi.org/10.1530/jrf.0.0810567

Frost, D.R., Wozencraft, W.C. and Hoffmann, R.S., 1991. Phylogenetic relationships of hedgehogs and gymnures (Mammalia, Insectivora, Erinaceidae). Smithson. Contrib. Zool., 1: 1-69. https://doi.org/10.5479/si.00810282.518

Fu, YX., 1997. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics, 147: 915-925.

Grenyer, R. and Purvis, A., 2003. A composite species-level phylogeny of the ‘Insectivora’ (Mammalia: Order Lipotyphla Haeckel, 1866). J. Zool., 260: 245-257. https://doi.org/10.1017/S0952836903003716

He, K., Chen, J.H., Gould, G.C., Yamaguchi, N., Ai, HS., Wang, Y.X., Zhang, Y.P. and Jiang, X.L., 2012. An estimation of Erinaceidae phylogeny: A combined analysis approach. PLoS One, 7: e39304. https://doi.org/10.1371/journal.pone.0039304

Huelsenbeck, J.P. and Ronquist, F., 2001. MRBAYES: Bayesian inference of phylogenetic trees. Bioinformatics, 17: 754-755. https://doi.org/10.1093/bioinformatics/17.8.754

Huson, D.H. and Bryant, D., 2005. Application of phylogenetic networks in evolutionary studies. Mol. Biol. Evol., 23: 254-267. https://doi.org/10.1093/molbev/msj030

Kashani, E., Rezaei, H., Naderi, M. and Khorasani, N., 2019. Population genetics of hedgehogs in Kerman and Sistan and Baluchistan Provinces. PhD thesis, Islamic Azad University, Science and Research branch, Tehran, Iran.

Lay, D.M., 1967. A study of the mammals of Iran: Resulting from the street expedition of 1962-63. Fieldiana Zoology, 54, Chicago, Field Museum of Natural History. https://doi.org/10.5962/bhl.title.2961

Librado, P., Rozas, J. 2009. DnaSP v5: A software for comprehensive analysis of DNA polymorphism data. Bioinformatics, 25: 1451-1452. https://doi.org/10.1093/bioinformatics/btp187

Low, V.L., Adler, P.H., Takaoka, H., Yacob, Z., Lim, P.E., Tan, T.K., Lim, Y.A., Chen, C.D., Norma-Rashid, Y. and Sofian-Azirun, M., 2014. Mitochondrial DNA markers reveal high genetic diversity but low genetic differentiation in the black fly Simulium tani Takaoka and Davies along an elevational gradient in Malaysia. PLoS One, 9: e100512. https://doi.org/10.1371/journal.pone.0100512

Lynch, M. and Crease, T., 1990. The analysis of population survey data on DNA sequence variation. Mol. Biol. Evol., 7: 377-394.

Nabholz, B., Mauffrey, J.F., Bazin, E., Galtier, N. and Glemin, S., 2008. Determination of mitochondrial genetic diversity in mammals. Genetics, 178: 351–361. https://doi.org/10.1534/genetics.107.073346

Naderi, M., Kaboli, M., Krystufek, B., Koren, T. and Zupan, S., 2014. Mitochondrial evidence uncovers a refugium for the fat dormouse (Glis glis Linnaeus, 1766) in Hyrcanian forests of northern Iran. Mammal. Biol., 79: 202-207. https://doi.org/10.1016/j.mambio.2013.12.001

Nei, M., 1987. Molecular evolutionary genetics. Columbia University Press, New York. https://doi.org/10.7312/nei-92038

Reeve, N., 1982. The home range of the hedgehog as revealed by a radio tracking study. Symp. Zool. Soc. Lond., 49: 207-230.

Rosenberg, N.A., 2004. DISTRUCT: A program for the graphical display of population structure. Mol. Ecol. Notes, 4: 137-138. https://doi.org/10.1046/j.1471-8286.2003.00566.x

Rozas, J., Sánchez-DelBarrio, J.C., Messeguer, X. and Rozas, R., 2003. DnaSP, DNA polymorphism analyses by the coalescent and other methods. Bioinformatics, 19: 2496-2497. https://doi.org/10.1093/bioinformatics/btg359

Schneider, S. and Excoffier, L., 1999. Estimation of past demographic parameters from the distribution of pairwise differences when the mutation rates vary among sites: Application to human mitochondrial DNA. Genetics, 152: 1079-1089.

Seddon, J., Santucci, F., Reeve, N. and Hewitt, G., 2001. DNA, footprints, of European hedgehogs, Erinaceus europaeus and E. concolor: Pleistocene refugia, postglacial expansion and colonization routes. Mol. Ecol., 10: 2187-2198. https://doi.org/10.1046/j.0962-1083.2001.01357.x

Stefanović, M., Veličković, N., Ćosić, N., Ćirović, D. and Đan, M., 2016. Genetic Variability of Northern White-Breasted Hedgehogs (Erinaceus roumanicus) From Serbia. Balkan J. Wildl. Res., 3: 1-6. https://doi.org/10.15679/bjwr.v3i1.35

Stopar, K., Ramšak, A., Trontelj, P. and Malej, A., 2010. Lack of genetic structure in the jellyfish Pelagia noctiluca (Cnidaria: Scyphozoa: Semaeostomeae) across European seas. Mol. Phylogen. Evol., 57: 417-428. https://doi.org/10.1016/j.ympev.2010.07.004

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. and Kumar, S., 2013. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol. Biol. Evol., 30: 2725-2729. https://doi.org/10.1093/molbev/mst197

Thompson, J.D., Higgins, D.G. and 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. Nucl. Acids Res., 22: 4673-4680. https://doi.org/10.1093/nar/22.22.4673

Zhou, X., Xie, Y., Zhang, Z.-h., Wang, C.-d., Sun, Y., Gu, X.-b., Wang S.-x., Peng X.R. and Yang, G.Y., 2013. Analysis of the genetic diversity of the nematode parasite Baylisascaris schroederi from wild giant pandas in different mountain ranges in China. Parasit. Vect., 6: 233-240. https://doi.org/10.1186/1756-3305-6-233