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 14515775, 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 Article Information Received 03 August 2019 Revised 25 December 2019 Accepted 14 February 2020 Available online 06 July 2020


INTRODUCTION
A s 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 O n l i n e

F i r s t A r t i c l e
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.

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 phenolchloroform 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 O n l i n e

F i r s t A r t i c l e
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)

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 H d = 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).

F i r s t A r t i c l e
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.

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 O n l i n e

F i r s t A r t i c l e
Genetic Diversity, Population Structure and Demographic History 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.

O n l i n e F i r s t A r t i c l e O n l i n e F i r s t A r t i c l e
Genetic Diversity, Population Structure and Demographic History of Hedgehogs