A total of 126 new samples (See supplementary file 1) were collected in the field, either from road-killed animals or live-trapping or from museum collections (Museum of Natural Science of Barcelona) and stored in 96% ethanol. These specimens come from Tunisia (72 specimens), Algeria (32), Morocco (2 specimens), Spain (15), Balearic (2) and Canary Islands (3). Geographical coordinates for the collected samples were recorded with a Global Positioning System. Species identification in the field was based on morphological characteristics, namely ventral, leg and facial color patterns, as well as ear length (Aulagnier et al. 2008).

Whole genomic DNA was extracted using the QIAamp DNA Mini Kit (QIAGEN), following the manufacturers’ protocols. We amplified a fragment of 456–465 bp of the mitochondrial control region (CR) using the primers developed by Derouiche et al. (2016): CRF: 5’-CATCAACACCCAAAGTTG-3’and Erinaceinae CRR: 5’-TGAAGAAAGAACCAGATG-3’. The amplification conditions consisted of an initial denaturation phase at 94 °C for 3 min, followed by 40 cycles each of denaturation at 94 °C for 40 s, hybridization at 53 °C for 40 s and elongation at 72 °C for 60 s; then a final elongation phase ran at 72 °C for 5 min. The quality and size of the PCR products were then verified by electrophoresis on 1% agarose and ethidium bromide-supplemented gel. The double-stranded PCR products were sent for purification and Sanger sequencing at EUROFINS. All newly determined sequences were submitted to GenBank (accession numbers OL692499–OL692624).

Chromatograms were checked and sequences manually corrected in CodonCode Aligner 8.0.2. In our phylogenetic analyses, we included all the newly sequenced specimens and all specimens of Atelerix and Paraechinus available in the GenBank database for which the control region was sequenced. We also added several references sequences for the hedgehog species Erinaceus europaeus, E. roumanicus and E. concolor. Sequences were aligned with ClustalW implemented in Geneious 4.6.4., and the alignment was corrected by eyes. The final alignment comprised 283 sequences (see supplementary file 1 for specimens list, and see supplementary file 2 for sequence alignment). Sequences of 423 bp were retained for the final analyses as they were available for most specimens.

We estimated evolutionary relationships among sequences by constructing phylogenetic trees using maximum-likelihood (ML) and Bayesian inference (BA). The computer program MrModeltest ver.2 (Nylander 2004) was used to evaluate the fit of 24 nested models of nucleotide substitution to the data. The model (HKY+I+G) chosen by MrModeltest according to the Akaike Information Criterion was then used in ML and BA analyses. The ML analysis was performed with PhyML ver. 3.0 (Guindon et al. 2010) on the PhyML online web server (Guindon et al. 2005). In the PhyML procedure, we used the BIONJ distance-based tree as the starting tree and 100 bootstrap replicates. Bayesian inference was performed with MrBayes ver. 3.2.6 (Ronquist et al. 2012) using 7,000,000 generations (burnin = 25%), sampling every 1,000 steps. Convergence between chains was checked by manual inspection using the ‘sump’ command and by using the mean SD of split frequencies, which was below 0.01 in our runs, as an indicator of convergence.

A haplotype network for each species (See supplementary file 3 for A. algirus sequences alignment, and see supplementary file 4 for P. aethiopicus sequences alignment) was constructed in TCS 1.2.0 (Clement et al. 2000) using statistical parsimony to calculate the minimum number of mutational steps among haplotypes with connections based on 95% confidence interval. Gaps were treated as 5^th state. The web-based program tcsBU was used to draw the network (Murias dos Santos et al. 2016).

Nucleotide diversity, haplotype diversity (Nei 1987) and the mean number of pairwise differences were calculated using DnaSP 6.12 (Librado and Rozas 2009).

Inferences of population expansion were made using two different methods: the Fs statistic (Fu 1997) and the R2 statistic (Ramos-Onsins and Rozas 2002). Ramos-Onsins and Rozas (2002) demonstrated that these statistics have the greatest power to detect population expansion for nonrecombining regions of the genome under a variety of different circumstances, especially when population sample sizes are large or when sample sizes are small. The significance of Fu’s Fs and R2 was obtained by examining the null distribution of 1,000 coalescent simulations of these statistics using DnaSP 6.12 (Rozas et al. 2017). Significantly large negative Fu’s Fs values and significantly positive R2 values were taken as evidence of a population expansion.

To evaluate current and past habitat suitability for the two species of hedgehogs (A. algirus and P. aethiopicus) we produced distribution models under three different climatic conditions (current, mid-Holocene and Last Glacial Maximum [LGM]) using a distribution modeling approach. Species occurrence data were assembled from the literature (Derouiche et al. 2016; Aulagnier et al. 2017; Deef 2020; El-Farhati et al. 2020; de Smet et al. 2021), museum collections (Global Biodiversity Information Facility, GBIF.org (2 December 2020) GBIF Occurrence Download https://doi.org/10.15468/dl.gkz4ya and https://doi.org/10.15468/dl.m8y5zz) and our own records (See supplementary file 5). Datasets of this type often exhibit sampling bias, which can lead to unreliable ecological niche modeling results (Fourcade et al. 2014). To correct sampling bias, a subsample of records regularly distributed in the geographical space was selected. The resolution of the reference grid was set to 2.5 arc min; i.e. the same as that used for the environmental layers. This coarse resolution was chosen to better match the coordinate uncertainties associated with georeferenced, textual localities of museum specimens. After filtering, the final dataset used for modeling consisted of 791 distribution records for A. algirus and 481 distribution records for P. aethiopicus.

As environmental layers, we used available climatic data from the WorldClim database (Hijmans et al. 2005). We removed from our analyses variables BIO2, BIO3, BIO5, BIO7 and BIO15 because those variables show a high level of discrepancy between general circulation models (GCMs) for North Africa (Varela et al. 2015). Colinearity of the initial variables was measured by Pearson’s correlation coefficient. We retained six variables that are considered ecologically important for hedgehogs and for which pairwise correlations were mostly smaller than 0.75: annual mean temperature (BIO1), seasonality of temperature (BIO4), minimal temperature for coldest period (BIO6), mean temperature of the warmest quarter of the year (BIO10), mean temperature of the coldest quarter of the year (BIO11), annual precipitation (BIO12). These climatic variables are ecologically important for hedgehogs because temperature affects hedgehogs’ prey activity, and hedgehogs are known to hibernate in cold winters (Mouhoub-Sayah 2009; Mouhoub Sayah et al. 2009). Climatic variables were downloaded for present conditions, for mid-Holocene (a wetter period, about 6,000 years ago) and for the LGM (a cold and dry period, about 22,000 years ago). Mid-Holocene and LGM climate data were drawn from General Circulation Model (GCM) simulations from two climate models: the Community Climate System Model (CCSM ver. 3) (Collins et al. 2006) and the Model for Interdisciplinary Research on Climate (MIROC ver. 3.2) (Hasumi and Emori 2004). These climate models differ in temperature and precipitation patterns. For a discussion of the uncertainties associated with the climatic data see Schorr et al. (2012) and Varela et al. (2015). The use of these two different climate models enabled us to assess and account for modeling uncertainty due to LGM climate data. To predict the potential distribution of the species in current conditions and in the past, we used Maxent ver. 3.3.3 (Phillips et al. 2006), which outputs a model with relative occurrence probability of a species within the grid cells of the study area. To ensure consistency of model predictions among repeated runs, we performed a 50-fold cross-validation with random seed. To determine whether the predictions for current conditions generated by Maxent were better than random predictions we used the area under the receiver operating characteristic curve (AUC), a commonly used measurement for comparison of model performance (Elith et al. 2006). The AUC ranges from 0 to 1, with greater scores indicating better discrimination ability; an AUC greater than 0.5 indicates that the model discriminates better than random.

The final alignment contained 283 sequences of 423 pb long with 32 sites with alignment gaps. ML and BA analyses resulted in a well resolved tree (Fig. 2) supporting the monophyly of the three included genera (Atelerix, Erinaceus and Paraechinus). All P. aethiopicus specimens cluster in a well-supported lineage. Specimens of A. algirus cluster in two distinct well supported lineages (I and II), which are not sister lineages. The first lineage includes most specimens and is distributed throughout the geographical range of the species. The second lineage includes only two Moroccan specimens. The average number of nucleotide differences between these two lineages is 27.574 with 22 fixed differences. As a matter of comparison, it is 27.000 with 27 fixed differences between A. algirus lineage II and A. albiventris lineage 1, and it is 27.667 with 26 fixed differences between A. algirus lineage II and A. albiventris lineage 2. The two sequences of A. algirus lineage II were previously published by Velo-Antón et al. (2019). Considering the high genetic divergence between these two lineages and the topology of the tree, doubts on the species identification and/or validity of these sequences are raised (see discussion part). Thus, they are not considered as belonging to A. algirus in subsequent analyses.

Figure 2. 

Bayesian phylogenetic tree generated in MrBayes for the hedgohogs. Circles and numbers at nodes represent posterior probabilities (Bayesian analysis) and bootstrap support (Maximum Likelihood analysis). Branch lengths are in substitutions/site.

The 204 specimens of A. algirus lineage I resulted in 58 haplotypes (Table 1). The nucleotide diversity is 0.00745 ± 0.00037, haplotype diversity is 0.924 ± 0.008 and the average number of nucleotide differences is 3.010. No strong phylogeographic structure is recorded within this species (Fig. 3). One abundant haplotype is shared between specimens from Tunisia, Algeria and Morocco, and several haplotypes are shared between specimens from Algeria and Morocco, or from Algeria and Tunisia. One haplotype is shared between specimens from Libya and Tunisia. High genetic diversity is observed within Morocco (51 specimens, 23 haplotypes) with numerous endemic haplotypes. Haplotypes from Morocco can be divided in two groups: the first haplogroup (haplotypes A01 to A13) is composed exclusively of haplotypes endemic to Morocco that differ from haplotypes from other countries by more than three mutations (except one haplotype shared between Morocco and Canary Islands), and the second haplogroup is composed of haplotypes shared (or closely related: one mutation) with specimens from Algeria and Tunisia. The first Moroccan haplogroup is mainly distributed in Atlantic plains of Morocco, and the second one is mainly distributed in southern and eastern Morocco. Both haplogroups occur in Atlas and Rif mountains (Fig. 1). Three haplotypes are found in Spain, and one of them is shared with Balearic Islands. We obtained significant negative values of Fu’s Fs for all haplogroups (Table 1). Positive values of R2 were obtained for all haplogroups, but they were not significant.

Figure 3. 

Statistical parsimony haplotype networks for Atelerix algirus (top) and Paraechinus aethiopicus (bottom) based on the control region. Codes refer to haplotype numbers (see Appendix S1 for the correspondance between haplotypes numbers, specimens numbers and geographic origin).

The 45 specimens of P. aethiopicus resulted in 19 haplotypes. Nucleotide diversity is 0.01086 ± 0.00050, haplotype diversity is 0.905 ± 0.026 and the average number of nucleotide differences is 4.386. A higher number of haplotypes is found in Algeria, but this is probably due to the greater sample size in this country (Fig. 3). No phylogeographic structure is apparent in the network: specimens from Morocco, Mauritania and Senegal can share the same haplotype, and this is also true for specimens from Algeria, Morocco and Mauritania, from Tunisia, Algeria and Mauritania, and from Algeria, Morocco and Tunisia. We obtained a negative value of Fu’s Fs and a positive value of R2, but they were not significant (Table 1).

We built species distribution models for the two African hedgehog species based on the known presence localities. The AUC values were 0.935 and 0.808 for A. algirus and P. aethiopicus, respectively. Thus, they are considered to correspond to a useful predictive model.

Under present-day climatic conditions the model revealed climatic suitability for A. algirus across all the Mediterranean Basin (Fig. 4). This area is larger than the currently known distribution of this species, as it is not present in Eurasian coasts of the Mediterranean basin (except in the Iberian Peninsula). In North Africa, the model shows three patches of high climatic suitability: one in Western Morocco isolated by High Atlas, one in northern part of Eastern Morocco and Northern Algeria and one in Tunisia and coastal part of Libya. At the mid-Holocene, according to both CCSM and MIROC models, areas of high climatic suitability were localized more or less at the same places than at present-day conditions but their size were smaller. During the LGM three areas of high climatic suitability were present: one in south-west Morocco, one in Tunisia and coastal Libya, and one on the coast at the border between Morocco and Algeria.

Figure 4. 

Species distribution modeling of Atelerix algirus as estimated by Maxent for present-day conditions (current), for Holocene and the Last Glacial Maximum (LGM) based on the Model for Interdisciplinary Research on Climate (MIROC) and on the Community Climate System Model (CCSM) paleoclimatic models. Warmer colors show areas with higher probability of presence.

For P. aethiopicus, in North Africa, the model revealed continuous high climatic suitability under present-day bioclimatic conditions in Western part of the Sahara desert from Mauritania to Morocco limited by the Anti-Atlas and the High Atlas (Fig. 5). Continuous high climatic suitability was predicted from Saharan Atlas (Algeria), Southern Tunisia and all coastal parts of Libya and Egypt to the Arabian Peninsula. The distribution predicted for P. aethiopicus during the mid-Holocene is remarkably similar to the present distribution model with more climatic suitability in Algeria and North-Western part of Mauritania according to both CCSM and MIROC data. For LGM period, CCSM and MIROC model show a high climatic suitability in Western part of Sahara desert, Libyan an Egyptian coastal parts.

Figure 5. 

Species distribution modeling of Paraechinus aethiopicus as estimated by Maxent for present-day conditions (current), for Holocene and the Last Glacial Maximum (LGM) based on the Model for Interdisciplinary Research on Climate (MIROC) and on the Community Climate System Model (CCSM) paleoclimatic models. Warmer colors show areas with higher probability of presence.

One of the most remarkable finding from this study is the discrepancy between our results and those published by Velo-Antón et al. (2019) for A. algirus. These authors found this species to be monophyletic and to include three distinct genetic lineages (lineages 1, 2 and 3), while we did not recover the monophyly of this species and found two distantly related lineages within it (lineages I and II). Our lineage II corresponds to the lineage 3 of Velo-Antón et al. (2019). It is composed of only two specimens from Morocco. In our phylogenetic tree this lineage does not cluster with lineage I and the genetic distance between these lineages is comparable to those found between lineage II and the two lineages of A. albiventris. Thus, the question of the species identification of these specimens is raised. Velo-Antón et al. (2019) also published the nuclear RAG1 sequences of these two specimens. They are both identical, with no polymorphic sites, and identical to the most frequent RAG1 haplotype found in A. algirus. They differ by several mutations from A. albiventris haplotypes. This mitonuclear discordance could be the result of an ancient hybridization event between the two species leading to a phenomenon of mitochondrial introgression. Alternatively, the presence of a new cryptic species in Morocco is possible. Furthermore, some species exhibit segmental duplications involving the control region, that, in some cases, leads to the formation of pseudogenes that may co-amplified by PCR (Bensch and Harlid 2000; Cadahia et al. 2009). Multilocus data are necessary to test these alternative hypotheses. If these hypotheses could explain the discrepancy between mitochondrial and nuclear data, we still do not understand why our mitochondrial results differ from those of Velo-Antón et al. (2019), given that they are based on the same DNA fragment.

Our A. algirus lineage I corresponds to the two sister lineages A. algirus 1 and A. algirus 2 of Velo-Antón et al. (2019). In their study, haplotypes from lineages 1 and 2 differ by at least 22 mutations. In our study these two lineages are not recovered in the phylogenetic tree, and they can only be identified in the haplotype network because the minimum number of mutations observed between haplotypes of these two lineages is only of three mutations (between haplotypes A13 and A14, or haplotypes A09 and A49). This discrepancy between the results of these two studies is striking. It is not only due to our greater sample size (allowing the recovery of more haplotypes), but may also reflect alignment differences and/or differences in the way of treating gaps and missing data (our network analyses are based on complete sequences only, the alignment included 7 gaps and gaps were treated as a fifth state; see Appendix S3).

To sum up, our results are congruent with the study of Khaldi et al. (Khaldi et al. 2016) in showing low genetic diversity in A. algirus across its geographic range. However, our more extensive sampling scheme and the use of a more rapidly evolving genetic marker allows us to highlight mtDNA phylogeographic structure within this species: one haplogroup is endemic to Morocco and Canary Islands, and another haplogroup is widely distributed from Morocco to Libya, and through Spain and Balearic Islands. These results suggest that the mitochondrial phylogeographic structure of the species has a recent origin. This low mitochondrial structure at continental scale contrasts with the marked east–west divergence that has been commonly found in previous phylogeographic studies of North African organisms and which identified two common disjunction areas: one around the Algeria-Morocco border and the other one in Kabylia (central Algeria) (Nicolas et al. 2015; Beddek et al. 2018). It is interesting to note that the two haplogroups co-occur in Atlas and Rif mountains. Two hypotheses can be proposed to explain this pattern: (i) vicariance of an eastern and a western haplogroup followed by secondary expansion or (ii) retention of deep mitochondrial DNA diversity in Morocco followed by loss of diversity during eastward expansion of the species out of Morocco. Note that both scenarios imply a population expansion of haplogroup 2, for which, we find evidence only with Fu’s Fs test. Our climatic niche modeling results support the first hypothesis since they show that during the LGM there were several disconnected areas of high climatic suitability for A. algirus: one in South-West Morocco, one at the coastal Moroccan-Algerian border and one in Tunisia-coastal Libya.

Our genetic results confirm that the Algerian hedgehog recently colonized Spain, Balearic and Canary Islands, and that this colonization was probably mediated by humans, as already suggested by several authors (Hutterer 1983; Borg 2005; Morales and Rofes 2008; Khaldi et al. 2016). Our ecological niche modeling analyses suggest that suitable climatic conditions occurred throughout most Southern and Eastern lowland Iberian Peninsula during the last 6000 years. This would have favored the spatial expansion of the Algerian hedgehog after its arrival in Europe. It is interesting to note that the current distribution of this species in Spain is much more reduced than expected based on climatic niche modelling. Species interactions (e.g. mutual competitive exclusion, predation, parasitism, and mutualism) are also important in determining species distribution (Pearman et al. 2008; Sinclair et al. 2010; Giannini et al. 2013). Competition with the European hedgehog, Erinaceus europaeus, a species already well established in Spain before the introduction of the Algerian hedgehog, may explain partially its restricted distribution in Spain.

Our study shows no phylogeograhic structure in North African P. aethiopicus. This is consistent with our climatic niche modeling results showing a more or less continuous high or moderate climatic suitability for this species throughout most of the Saharan desert since the LGM (except in the most central parts). This study significantly increases our knowledge of the distribution of P. aethiopicus expanding currently recognized ranges (IUCN, 2015; Best, 2018). An extension of its distribution range in Atlantic Sahara was already proposed by Aulagnier et al. (2017). Velo-Antón et al. (2019), based on molecular data, proposed its extension in Mauritania, Senegal, Atlantic Sahara and Libya. We show that its distribution range should also be extended in Algeria since we captured it in the M’sila region, up to 35.748°N, confirming, with molecular data, the previous observation of Sellami et al. (1989).

Several subspecies have been described for Atelerix algirus, either from the mainland or islands, Atelerix algirus algirus from North Africa, A. algirus vagans from Spain, Minorca and Balearic Islands and A. algirus girbaensis from Djerba Island (Tunisia). Given that the differential diagnosis is mainly based on pelage variations (hairs, spines), the validity of these subspecies has been regularly questioned (Corbet 1988; Hutterer 2005; Happold 2013). Our results do not show any strong genetic discontinuity between the Algerian, Moroccan, Tunisian, Libyan, Spanish, Balearic and Canarian populations, suggesting that subspecific recognition may be unwarranted. The status of the subspecies from Djerba Island (A. a. girbaensis) still needs to be assessed by genetic data. Several subspecies were also described within P. aethiopicus: P. a. aethiopicus from Soudan and Ethiopia, P. a. pectoralis from Arabian Peninsula, P. a. deserti from North Africa, P. a. albatus from Djerba Island and P. a. ludlowi from Iraq. Our genetic data do not show any genetic discontinuity throughout North Africa, from Senegal to Libya, confirming the results previously obtained by Derouiche et al. (2016) and Velo-Antón et al. (2019). Derouiche et al. (2017), using a combination of nuclear and mitochondrial markers, studied the genetic differentiation between northwest African populations (deserti) and Arabian populations (dorsalis) of Ethiopian hedgehog. They showed a level of genetic differentiation that is compatible with their taxonomic separation postulated from morphology. Additional sampling in Djerba Island, Soudan, Ethiopia, and in the Middle East is needed to test the genetic differentiation of others proposed subspecies.