Introduction

The Mediterranean Basin, between southern Europe and North Africa, is particularly suited to study the role of climate versus human-mediated changes on the evolutionary history of species for several reasons. First, this area has yielded a rich paleontological record allowing a better understanding of the morphological evolution of small mammal species and of community assemblages over time in relation to Quaternary climatic changes (Stoetzel 2013). Second, North Africa is now recognized as a major area for the emergence and dispersal of anatomically modern humans from at least 315 kya (Garcea 2012; Hublin et al. 2017; Richter et al. 2017). It thus constitutes a region of interest for investigating how the arrival of early representatives of Homo sapiens and their recent impact on the environment influenced faunal genetic diversity and distribution. During the Neolithic and after, numerous contacts between the two shores of the Mediterranean Sea occurred, due to the development of navigation (Mulazzani et al. 2010; Souville 1998; Zilhão 2014). This led to the translocation (Dobson 1998) and an almost complete turnover of mammal species in the Mediterranean islands (Vigne 1992).

The Algerian mouse, Mus spretus Lataste, is an emblematic species endemic to the Mediterranean Basin. It occurs in North Africa, from Morocco through northern Algeria and Tunisia to northern Libya, and in south-west Europe (southern France, most of the Iberian Peninsula and the Balearic Islands; Palomo et al. 2009). The oldest fossil of M. spretus is dated from the Middle Pleistocene, about 600 kyr BP, in Tunisia (Jebel Ressas 2, 3; Stoetzel 2013). This species is also known from several sites in Morocco in the Middle-Late Pleistocene, early Holocene and middle-late Holocene (Stoetzel 2013), while it is reported from Europe only in the Neolithic (Palomo et al. 2009). Fossil evidence thus suggests a North African origin of the species and recent spreading to southwestern Europe via Neolithic navigators (Dobson 1998). Based on restriction fragment length polymorphism (RFLP) data, Boursot et al. (1985) showed a strong genetic isolation between European and North African populations, but they could not characterize the tempo and mode of such divergence. Based on few D-loop mitochondrial DNA sequences (mtDNA), Desmarais (1989) showed the absence of actual gene flow between European and African populations, and a probable recent colonization event of Europe from North Africa with a strong founder effect. He also proposed that during an unfavorable climatic period, allopatric diversification occurred in two refugia in North Africa (Tunisia and Morocco).

The Late Pleistocene to middle Holocene succession of El Harhoura 2 cave (Morocco, region of Témara), yields an exceptionally rich small vertebrate assemblage and showed a continuous fossil record for M. spretus from 107 ± 7 kyr BP to 5800 years BP, with nonsignificant change in abundance through time, despite pronounced environmental changes (glacial–interglacial cycles, Stoetzel et al. 2011b). The morphology of the molars of Algerian mouse through this fossil record was also surprisingly stable (Stoetzel et al. 2013). The authors proposed that such evolutionary stability despite environmental changes could be explained by the fact that preferred habitats (grasslands, dry shrubland and open woodland) persisted in the surroundings of El Harhoura 2 throughout the record as a “mosaic” of habitats despite variations in their extension. In addition, the coastal location of this cave probably reduced the impact of global climatic variation compared with further inland. This long-term morphological stability contrasts with the molar size decrease observed between the last El Harhoura 2 Holocene level (5800 years BP) and modern Moroccan mouse populations (Stoetzel et al. 2013). Geographical morphological variation have been documented across the modern distribution range of this mouse (Darviche and Orsini 1982; Darviche et al. 2006; Khammes-El Homsi and Aulagnier 2010; Palomo et al. 1983; Palomo et al. 1985; Stoetzel et al. 2013; Valenzuela-Lamas et al. 2011). The modern geographical variation was hypothesized to be due to relative isolation between populations leading to divergence and possibly local adaptation. The differentiation of Algerian mice from a pool corresponding to Morocco and Southern Europe (Boursot et al. 1985; Khammes-El Homsi and Aulagnier 2010) may have been triggered by North African mice populations being fragmented by places where the desert reaches the sea or in altitude, isolating favorable places (“continental islands”) for the mice (Stoetzel et al. 2013). However, even at a local scale, great morphological variation is observed between present day Moroccan populations that cannot be attributed to such large-scale biogeographical structure (Stoetzel et al. 2013). The three last millennia are known as a period of important environmental perturbations, with an increase of aridity (Hassan 1997; Le Houérou 1997) enhanced by the progressive increase in anthropogenic activities (Brun 1989; Brun 1991; Stambouli-Essassi et al. 2007) and concomitant introduction of commensal rodents such as the house mouse, and the black and brown rats (Audoin-Rouzeau and Vigne 1994; Bonhomme et al. 2011; Cucchi and Vigne 2006). The competition and possible hybridization with the house mouse Mus musculus domesticus, together with anthropogenic changes in landscape, may have particularly impacted the local evolution of M. spretus by restricting its habitat, triggering divergence of isolated populations, which were previously interconnected. The anthropogenic action may also have modified the environment in a more global way. By the fact, among other effects, it was observed a significant size decrease of several wild animals during historical time (Sheridan and Bickford 2011). The size decrease of the Western Mediterranean mouse during historical times may be the discrete signature of a perturbed ecosystem, pushing the species away from its previous stasis as it was evidenced in islands for other species (Sondaar 2000).

In the present study, we combined mitochondrial and nuclear data to infer the evolutionary history of the Algerian mouse. Specifically, we wanted to answer the following questions: (1) When and how did the Algerian mouse arrive in Europe? (2) Can we detect a recent population bottleneck in North African populations linked to the development of anthropogenic activities and/or the arrival of the house mouse? (3) Is the high morphological variability observed in present day populations due to isolation? Moreover, we evaluated these genetic data with climatic niche modeling to reveal the evolutionary history of M. spretus in response to climate change during the Quaternary.

The combination of molecular, fossil, palaeoenvironmental and human context data should allow us to reconstruct the historic fluctuations of populations within this taxon in North Africa and Europe.

Materials and methods

Ethics statement

Animals were live-trapped and handled under the guidelines of the American Society of Mammalogists (Sikes and Mammalogists atACaUCotASo 2016). The protocol was approved by “Comité Cuvier” (permission no. 68.009). All manipulations of animals were made in Morocco in agreement with the global law 11-03 relative to the protection and the development of the environment. Alive animals were euthanized by the injection of a lethal dose of isoflurane, followed by cervical dislocation. Capture permits were obtained through the “Haut Commissariat aux Eaux et forêts et la Lutte contre la desertification” (Authorisation no. 15 HCEFLCD/DLCDPN/DPRN/CFF) in Morocco, and through the Ministry of Forestry in Algeria. No capture permit is necessary in Europe, because this species is under no specific status.

Sampling and DNA extraction

A total of 635 newly collected individuals were sampled in 51 localities (81 specimens from 4 Algerian localities, 448 specimens from 28 Moroccan localities, 7 specimens from 4 Tunisian localities, 61 specimens from 7 French localities, 1 specimen from one locality in Portugal and 37 specimens from 7 Spanish localities; Table 1, Table S1 and Fig. 1). Genomic DNA was extracted from tissues preserved in 95% ethanol using the NucleoSpin Tissue Core kit (MACHEREY-NAGEL, Hoerdt, France).

Table 1 Collection localities of Mus spretus specimens included in this study, with the number of specimens sequenced per DNA region and per locality
Full size table
Fig. 1
figure 1

Map of sampling localities and the names of the main geographical features discussed in the text. For locality names, see Table 1. Gray areas indicate mountains higher than 1000 m

Full size image

DNA amplification and sequencing

Three mitochondrial genes (Cytochrome b—Cytb-, Cytochrome oxidase I—CO1- and control region—CR-) were amplified following (Ducroz et al. 2001; Rajabi-Maham et al. 2008; Robins et al. 2007). All newly determined sequences were submitted to GenBank (accession numbers MH241060 to MH242499 and MH253776 to MH253817).

We genotyped 536 individuals from 19 localities (Table 1) for 18 microsatellite loci using primers and genotyping protocols described in Lalis and Orth (2014).

Sequence analyses

Median-joining networks

In our phylogenetic analyses, we included 471 specimens for Cytb (coming from Algeria, Morocco, Spain, France and Portugal; Table S1), 499 specimens for CO1 (Algeria, Morocco, Tunisia, Spain, France and Portugal) and 518 specimens for CR (Algeria, Morocco, Tunisia, Spain, France and Portugal). Phylogenetic relationships between haplotypes were inferred by constructing a network using the median-joining (MJ) method available in NETWORK v4.500 (Bandelt et al. 1999).

Sequences of 1085 bp for Cytb, 655 bp for CO1 and 840 bp for CR were retained for the network analysis to minimize the number of incomplete sequence as adding ambiguous data in MJ trees is problematic. Analyses were performed for each gene separately and for the three mtDNA genes combined. For the concatenate analysis, only the 425 specimens for which the three mtDNA regions were sequenced were included in the analysis.

Genetic diversity and population differentiation

Only populations with at least 12 individuals were retained in our analyses. Nucleotide diversity, haplotype diversity (Nei 1987) and the mean number of pairwise differences were calculated using DnaSP 5.10 (Librado and Rozas 2009). Population differentiation was inferred from the mtDNA dataset using Fst estimates, computed with ARLEQUIN version 3.5.1 (Excoffier and Lischer 2010).

Demographic history and divergence time estimates

We tested for signatures of population expansion through the Fu’s Fs (Fu 1997), Tajima’s D (Tajima 1989) and Ramos-Onsins and Rozas R2 (Ramos-Onsins and Rozas 2002) tests of neutrality. These analyses were performed with Arlequin 3.5 (Excoffier and Lischer 2010) and DnaSP 5.10 (Librado and Rozas 2009), and the significance of the neutrality statistics was assessed using 1000 coalescent simulations under a model of constant population size using empirical estimates of theta.

We further analyzed the demography of the main populations using the Extended Bayesian Skyline Plot (EBSP) implemented in BEAST v.1.8 (Drummond and Rambaut 2007). EBSP is a non-parametric model that does not specify any prior hypothesis on the tempo and mode of changes of the effective population size through time (Heled and Drummond 2008). EBSP considers populations as fully isolated, meaning that they have evolved separately from any others up to the time to the most recent comment ancestor of the locus(i) under investigation. This hypothesis is clearly violated in our case. However, EBSP can provide important insights if correctly interpreted in the light of the presence of population structure (Maisano Delser et al. 2016; Mazet et al. 2016). We ran 100 million Markov chain Monte Carlo (MCMC) iterations with a thinning interval of 10,000. We checked convergence by visually inspecting the trace and computing the effective sample size (which was always higher than 100) for each parameter in two independent runs using the program TRACER 1.5 (http://tree.bio.ed.ac.uk/software/tracer/). EBSP were plotted using R version 3.3.2 (R Development Core Team 2016) setting, the burn-in to 10%. We set a strict clock with rate of 3 × 10-7 substitutions per site per year and the Hasegawa–Kishino–Yano plus Gamma model of molecular evolution. It is now widely recognized that intraspecific and pedigree-based estimates of substitution rates are generally higher than interspecific phylogenetically calibrated rates (Henn et al. 2009; Ho et al. 2007b; Loogväli et al. 2009). Estimating intraspecific substitution rate is not trivial (Bonhomme et al. 2011). A rate of 4 × 10-7 substitutions per site per year was proposed for the CR of the congeneric rodent species Mus musculus domesticus (Rajabi-Maham et al. 2008). A Cytb clock rate of 3 × 10-7 substitutions per site per year was estimated using radiocarbon dated post-glacial remains in another rodent, Microtus arvalis (Martinkova et al. 2013). Similar clock rates have been obtained with dated ancient DNA from a wide range of mammals and birds (Ho et al. 2007a).

Microsatellite analyses

Genetic diversity and differentiation

Genetic variability of the microsatellite markers was measured for each population by the number of alleles (Na), gene diversity (He, expected heterozygosity) and observed heterozygosity (Ho) using R package ADEGENET V3.3.1 (Jombart 2008). The departures from Hardy–Weinberg equilibrium (HWE) and the linkage disequilibrium between pairs of loci, with significance set at P < 0.05, were assessed with GENEPOP v. 4.1 (Raymond and Rousset 1995) and were later adjusted using a sequential Bonferroni correction (Rice 1989).

The program BOTTLENECK was used to test for past decrease in effective population size (Piry et al. 1999). Most microsatellites fit a two-phase mutation model (TPM) better than a strict stepwise mutation model (SMM) or an infinite alleles model (IAM). However, microsatellites with imperfect repeat patterns are more likely to evolve under the IAM model (Cornuet and Luikart 1996). As we found no strong reason to select or reject one model over another, we ran Bottleneck under all three mutational models. For the TPM model, the proportion of SMM was set to 70% and the variance to 30 (default values). The significance was assessed by performing 10,000 replicates.

Population structure, assignment and admixture analyses

STRUCTURE v.2.3.353 (Pritchard et al. 2000) was used to infer geographical genetic structure. We assumed an admixture model with correlated allele frequencies (Falush et al. 2003). For each value of the number (K) of clusters set between 1 and 19 (number of sampled sites), we performed 10 independent runs of 106 MCMC iterations with a burn-in of 105. The optimal K value was inferred using the ad-hoc Evanno statistic Delta K (Evanno et al. 2005). For the selected values of K, the software CLUMPP v1.2 (Jakobsson and Rosenberg 2007) was used to align cluster membership coefficients from the 10 replicate cluster analyses using the Greedy algorithm with 10,000 random input orders. The results were graphically displayed with DISTRUCT v 1.1 (Rosenberg 2004).

We also performed a discriminant analysis of principal components (DAPC), using R package ADEGENET version 1.2-7 (Jombart 2008) for the analysis of genetically structured populations. Dataset was also depicted in an unrooted neighbor-joining (NJ) tree (Takezaki and Nei 1996) based on the pairwise Reynold’s distances computed with POPULATION 1.2.32 (Langella 2010). The robustness of nodes was evaluated by carrying out 1000 bootstrap replicates. Population differentiation was further analyzed by computing estimates of FST (Weir and Cockerham 1984) between all population pairs using GENEPOP V4.1.3, whose significance was assessed through 10,000 permutations. P-values were adjusted using sequential Bonferroni correction.

Genetic inference of colonization pathways and demographic history of the three main populations (Africa, Spain, France)

We devised four models to decipher the colonization pathways and the demographic history of M. spretus by means of approximate Bayesian computation (ABC) (Bertorelle et al. 2010). These models were designed to address two main issues: (i) the timing and colonization dynamics of Europe; (ii) the response of M. spretus to recent environmental changes in North Africa. The four models and the associated parameters are shown in Figure S1 with the prior distributions presented in Table S2. Briefly, models 1 and 2 represent an independent colonization of France and Spain from North Africa, whereas models 3 and 4 represent a sequential colonization of Europe, with M. spretus invading first Spain and then France. Models 1 and 3 assumed an instantaneous change in effective population size in North Africa followed by an exponential change in effective population size until present, whereas models 2 and 4 assumed a constant effective population size through time. In all the four models, we assumed an instantaneous bottleneck during the colonization process of both France and Spain followed by an exponential change in effective population size until present. The coefficient of the exponential change was computed as follows: first, we extract a time in generations when the change occurred, than a new effective population size (i.e., an instantaneous change), and last, a new modern effective population size. Taking as example MOD1 (or MOD3), this would correspond to extract, respectively, “T change Africa”, “N change Africa” and “N Africa”. The variation of the Ne from “N change Africa” to “N Africa” follows an exponential distribution whose coefficient is determined by these three parameters. Depending on the prior extraction, the modern population would be constant (if “N change Africa” = “N Africa”), expanding (if “N change Africa” < “N Africa”) or contracting (if “N change Africa” > “N Africa”). The model is therefore free to vary, in order to detect if modern population is stable or is it changing through time. We note that MOD2 and MOD4 are simply nested in MOD1 and MOD3 (i.e., “N change Africa” = “N Africa”). Exponential size change in France and Spain are modeled following the same scheme. Each of the three populations is a pool of the sampled populations coming from that geographic area. This was done to lower the number of parameters. We generated 500,000 simulations for each demographic model using fastsimcoal2 version 2.5.1 (Excoffier et al. 2013). Random values from the prior distribution of each parameter were extracted using an R script. Model posterior probabilities were calculated by a weighted multinomial logistic regression (Beaumont 2008) for which we retained the best 60,000 simulations. The demographic parameters within each model were estimated from the 5000 simulations closest to the observed dataset using a local linear regression (Beaumont et al. 2002). Analyses were performed in the R environment with the library abc (Csillery et al. 2012).

We run the four models with two independent set of markers: (i) 18 microsatellites (“MSAT”); (ii) mtDNA, combining d-loop, Cytb and CO1 (“MTDNA”). For each set, we considered a generation time of 0.5 year (Palomo et al. 2009; Vargas et al. 1991). In MSAT, we used per generation per locus mutation rate of 10-4 (Dib et al. 1996) and we computed the following summary statistics: mean and variance over loci for each population of the number of alleles, Garza–Williamson statistic, allelic range and the pairwise Rst. In MTDNA, we used per generation per site mutation rate of 3×10-7 (Martinkova et al. 2013) and we compute the following summary statistics: number of haplotypes, number of segregating sites, Tajima’s D, mean pairwise difference in each population and the pairwise molecular Fst. All summary statistics were computed with arlsumstat (Excoffier and Lischer 2010).

To test the validity of model assumptions we performed a principal component analysis (PCA) on the reference table of summary statistics simulated from prior distributions adding the vector computed from real data. We visually inspect the resulting two first dimensions of the PCA plot to check if the observed vector fell within the clouds of simulated points. The four models passed this test for both datasets (Figure S2).

A cross-validation procedure was implemented to evaluate the model selection performance using the same threshold as in the real data. A similar cross-validation procedure was implemented to evaluate the precision and accuracy of ABC parameter estimates: pseudo-observed datasets (pods) were simulated from prior distributions and re-analyzed using the reference table and the same parameters of the ABC algorithm. We computed the 95% coverage, the scaled mean error (SME) and the scaled root mean square error (SRMSE) calculated as in Walther and Moore (2005) following Maisano Delser et al. (2016). Finally, we performed a Bayesian posterior predictive test (Gelman et al. 2004) using the global Rst in MSAT or Fst in MTDNA as predictive statistics (Figure S3). We chose these two statistics as they were not used in the estimation process as suggested by Bertorelle et al. (2010).

Finally, preliminary results (see below) suggested the presence of significant population structure in the metapopulation from North Africa. Pooling the populations without taking structure into account may bias the demographic estimates, but the extent of this bias is unknown. To remove the influence of population structure, we followed the sampling theory proposed by previous works (Wakeley 1999; Stadler et al. 2009; Mona et al. 2014) by sampling few lineages per deme and pooling them. This scattered sample (sensu Wakeley 1999) eliminates the influence of local coalescence and describes the behavior of the whole metapopulation. To this end, we created 1000 resampled dataset by randomly sampling 25 and 20 African individuals from the MSAT (RES-MSAT) and MTDNA (RES-MTDNA) dataset, respectively. We performed 500,000 simulations for each of the four models as before and applied the same ABC framework for each resampled dataset (both model selection and parameter estimation). Parameters estimations are based on the average of the 1000 posterior distributions.

Species distribution modeling

We modeled the climatic niche of M. spretus to approximate its current, last glacial maximum (LGM) and mid-Holocene distributions. We applied ecological niche modeling methods, where environmental data are extracted from current occurrence records, and habitat suitability is evaluated across the landscape using program-specific algorithms (Elith et al. 2006). The present day models were then projected on the climatic reconstructions of the LGM (a cold and dry period, about 22,000 years ago) and mid-Holocene (a wetter period, about 6000 years ago) under the assumption that the climatic niche of the species remained conserved through time (Elith et al. 2010).

Distribution records of Algerian mouse were assembled from literature (Aulagnier et al. 2017; Brahmi et al. 2010; Brahmi et al. 2012; Khammes-El Homsi and Aulagnier 2010), museum collections (Global Biodiversity Information Facility, GBIF.org (7th September 2016) GBIF Occurrence Download http://doi.org/10.15468/dl.pcgrno) and our own unpublished databases. 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 using ENMtools ver. 1.3 (Warren et al. 2010). The resolution of the reference grid was set to 2.5 arc min (nearly 5 × 5 km), which is 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 702 distribution records. However, this final dataset still exhibited sampling bias: most of these records came from Europe and few from Africa (Figure S4). To overcome this bias it is possible to use a bias file: MAXENT can be fed with a bias grid that is a sampling probability surface. Creating this bias file is difficult because we need to know the actual sampling intensity across the study area. We used a different approach: we considered that sampling intensity was different in Europe and North Africa, and performed two distinct analyses. One analysis with only the distribution records from Europe, which was used to estimate present and past habitat suitability for M. spretus in Europe, and one analysis with only the distribution records from Africa, which was used to estimate habitat suitability in Africa.

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 in ENMtools v1.3 (Warren et al. 2010). A total of five variables, all of which had a Pearson correlation coefficient lower than 0.75, were retained. The final set of environmental predictor variables used for the species distribution models (SDM) consisted of: annual temperature (BIO1), precipitation of wettest quarter (BIO16), precipitation of driest quarter (BIO17), precipitation of warmest quarter (BIO18) and precipitation of coldest quarter (BIO19).

Climatic variables were used for present conditions and for the LGM and mid-Holocene. Paleoclimate data were drawn from the GCM simulations based on two climate models: the Community Climate System Model (CCSM, version 3) (Collins et al. 2006) and the Model for Interdisciplinary Research on Climate (MIROC, version 3.2) (Hasumi and Emori 2004). 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 LGM, we used MAXENT ver. 3.3.3 (Phillips et al. 2006), which returns 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 10-fold cross-validation with random seed. We set the regularization multiplier to 2 (Radosavljevic and Anderson 2014). 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 >0.5 indicates that the model discriminates better than random.

Results

MtDNA results

The three mtDNA genes were sequenced for 425 specimens, resulting in 263 haplotypes and 299 polymorphic sites. The average number of nucleotide differences is 17.725.

Results obtained for the three mtDNA genes were congruent (see Figure S5 for the results obtained for each gene), thus only the results obtained for the combined dataset (Cytb + CO1 + CR) are shown on Fig. 2. Several haplogroups can be identified on the network. Nearly all specimens (67/69) from Europe cluster together in the network; the only exceptions being one specimen from Grenada (Spain), which clusters with African haplotypes of several localities and one specimen from Barcelona (Spain) who shares its haplotype with one specimen from Tetouan (Morocco). Within this European haplogroup, two sub-groups can be identified: one grouping all specimens from France, and one grouping most specimens from the Iberian Peninsula. In both France and the Iberian Peninsula, a clear starlike pattern is observed, suggesting a recent demographic expansion. Several haplogroups grouping together specimens from Morocco and Algeria can be identified, none of which showed a starlike pattern. In continental Europe, branches lengths are short (always <4 mutations between haplotypes), whereas in North Africa branches lengths tend to be longer (many branch lengths >6 mutations, with a maximum of 13 mutations).

Fig. 2
figure 2

Minimum spanning networks of M. spretus obtained for the three mtDNA regions concatenated. Circle sizes are proportional to the number of similar haplotypes observed in the dataset. Branch lengths are proportional to the number of mutations between haplotypes. Colors represent the country of sampling (red = Portugal; blue = France; pink = Spain; yellow = Morocco; green = Algeria)

Full size image

Haplotype diversity is high in all localities (Table 2), with values ranging from  0.874 ± 0.048 to 0.999 ± 0.004. Nucleotide diversity is low in the two European localities (Madrid in Spain = 0.00146 ± 0.00013 and Montbazin in France = 0.00093 ± 0.00015). In North African localities, nucleotide diversity is higher, varying from 0.00468 ± 0.00071 to 0.00742 ± 0.00064. The high haplotype and low nucleotide diversities observed in European localities suggest a rapid population growth from an ancestral population with small Ne, provided the time was sufficient for recovery of haplotype variation via mutation yet too short for an accumulation of large sequence differences. The high haplotype and nucleotide diversities observed in North African populations suggest stable populations with large long-term Ne or admixed sample of individuals from historically sundered populations.

Table 2 Genetic diversity estimates and neutrality tests for the Algerian mouse based on the mitochondrial DNA data
Full size table

We tested for signatures of population size change for each population with neutrality tests (Table 2) and EBSP (Figure S6). In all North African populations, neutrality tests are nonsignificant, except in Berrechid. In Berrechid, a signal of demographic expansion is revealed by the Fs test and the EBSP (started ca 7000 ya). In the two European localities (Madrid and Montbazin), the neutrality tests and EBSP analyses reveal a signal of demographic expansion, which probably started ca 1000 ya.

The results of the differentiation analyses using FST are shown in Table S3. In total, 87% of pairwise comparisons are significant. MerjaZerga, Oued Nfifikh, BenGuerir, Ifrane and Taza are the localities that showed fewer significantly different comparisons with other populations, whereas BenSlimane, DayetErroumi, Montbazin and SidiBoughaba showed significant differentiation in all cases.

µsat results

All the 18 microsatellites were polymorphic in the genotyped 536 individuals. Tests for genotypic linkage disequilibrium confirmed the null hypothesis of the independence of loci pairs. According to locus and population, the expected heterozygosity (He) values vary between 0.587 and 0.744 (Table S4). After sequential Bonferroni corrections were applied, all populations were in HWE with FIS estimate ranging from 0.034 to 0.167 (Table S4).

Three models were used to test for recent genetic bottlenecks based on allele frequency data. Although tests based on the IAM mutation model suggested potential signatures of past genetic bottlenecks in samples collected in Algeria (Meftah) and Morocco (BenGuerir, BeniHadifa, BenSlimane, Berrechid, Cherouda, Clementine, Ifrane, Oued Nfifikh, SidiBoughaba and Taza), those based on the most realistic SMM mutation model were only significant for five Moroccan localities (Berrechid, Cherouda, Clementine, Oued Nfifikh and SidiBoughaba; Table S4). In contrast, no significant values were observed using the TPM model (P > 0.05, Wilcoxon test).

The results of the Bayesian clustering analysis (using STRUCTURE) and Evanno’s method supported K = 4 as the most likely number of clusters, which corresponded to a “Moroccan cluster” (yellow) including all populations from Morocco excepted Cherouda, Tetouan and Ifrane, an “Algerian cluster” (green) consisting of Algerian populations and Ifrane, a “Spanish cluster” (pink) with the population from Madrid and the Moroccan Cherouda and Tetouan populations, and a “French cluster” (blue) with the two French populations (Fig. 3a). Assuming K = 3, the two French populations, the Spanish population and the Moroccan populations of Cherouda and Tetouan were regrouped in the same cluster (blue, Figure S7). A second cluster includes the two Algerian populations and the Moroccan Ifrane population (green). And finally, the last cluster includes all other Moroccan populations (yellow). At K varying from 5 to 6, the Moroccan populations were gradually separated into distinct clusters. The Moroccan population from Ifrane (orange) is isolated from other Moroccan populations, and hence, the Algerian populations form an isolated cluster (green). In contrast, the Spanish population and the Moroccan Cherouda and Tetouan populations remained regrouped in the same cluster, suggesting that genetic divergence among Moroccan populations was larger than between Spanish and Moroccan Cherouda—Tetouan populations (Figure S7).

Fig. 3
figure 3

a Bayesian clustering analyses of M. spretus populations genotyped with 18 microsatellite loci. Clustering was performed in STRUCTURE. Each vertical line represents an individual, and each color represents a cluster (in blue: France; in pink: Madrid—Cherouda—Tetouan; in green: Algerian populations and Ifrane, and in yellow: the others Moroccan populations). Individuals are grouped by sampling location (locality number as in Table 1). b Discriminant analysis of principal components (axis 1 and 2). The axes represent the first two discriminant functions (DA) respectively. Individuals (dots) and populations (colored ellipses) are plotted within the orthogonal space defined by the first two 2 DA eigenvalues (inserts). c Neighbor-joining tree based on Reynold’s distance. Bootstrap values are calculated over 1000 replicates. For b, c, each color represents a cluster: in blue: France; in pink: Madrid; in purple: Moroccan populations Cherouda and Tetouan; in green: Algeria and in yellow: others Moroccan populations

Full size image

On the DAPCs (Fig. 3b), the first two axes explained >97% of the total inertia. The first axis allows to discriminate three groups of populations: the French populations (Montbazin and La Gardiole), the Spanish population (Madrid) and all African populations. On axis 2, the populations from Cherouda is distinct from all other populations. Results of the NJ tree (Fig. 3c) are congruent with those of the DAPC and confirms the close affinity between the two French populations and the distinctiveness of the Spanish population. All African populations are grouped together, and the populations from Cherouda and Tetouan are the closest to the Spanish population.

Despite the geographical distances involved, pairwise FST estimates based on microsatellite data remained relatively moderated (FST ≤ 0.152; Table 3).

Table 3 Pairwise Fst values between M. spretus populations samples (microsatellites data)
Full size table

When considering all the populations, 128 of the 171 pairwise comparisons tested were significantly different from zero (Table 3). FST was not significant between: the Moroccan population Ifrane versus the two Algerian populations; the Moroccan populations Tetouan and Cherouda versus the Spanish population Madrid; between the two Algerian populations Cap Djinet versus Meftah; between the two French populations La Gardiole versus Montbazin (Table 3), which appeared homogenous, with no significant differentiation following the use of Bonferroni correction, suggesting the occurrence of extensive gene flow among these populations. Un-significant genetic differentiation was also observed among most remaining Moroccan populations (37/91 comparisons). Significant results were mainly due to the populations of BenSlimane, MerjaZerga and SidiBoughaba. Interestingly, BenSlimane seems to be an isolated population as it was significantly differentiated from the other populations except Aounate and Clementine.

Genetic inference of colonization pathways and demographic history of the three main populations (Africa, Spain, France)

We tested the potential routes of colonization of M. spretus from North Africa to Europe using ABC methods. Results of model selection for the two datasets are shown in Table 4, whereas cross-validation experiments are presented in Table S5. First, cross-validation suggests that the four models are not easily distinguishable. If we sum the posterior probability of models supporting independent (model 1 + model 2) versus serial (model 3 + model 4) colonization of Europe, MSAT was not able to distinguish between the two hypotheses (both having approximately P = 0.5), whereas MTDNA attained ~0.8 probability of a serial colonization (Table 4). If we consider the results from the cross-validation experiment, it is clear that such a high probability value obtained with MTDNA cannot occur if the independent colonization was the true model. Moreover, both RES-MSAT and RES-MTDNA favored the serial colonization (P = 0.6 and P = 0.87, respectively). Considering only the serial models (models 3 and 4), MSAT and RES-MSAT supported a colonization of Spain and France very close in time: ~5000 and ~3000 years BP for MSAT, respectively; and ~13,000 and ~11,000 years BP for RES-MSAT, respectively (Table 5). For MTDNA and RES-MTDNA, the posterior distribution of this parameter is flat, implying this parameter is not well estimated under such model (Table S6).

Table 4 Posterior probabilities of the four ABC models under microsatellite data and mtDNA data both for the pooled African populations (MSAT and MTDNA, respectively) and the resampling strategy (RES-MSAT and RES-MTDNA, respectively)
Full size table
Table 5 Parameters estimation for the four models under microsatellite data and mtDNA data both for the pooled African populations (MSAT and MTDNA, respectively) and the resampling strategy (RES-MSAT and RES-MTDNA, respectively)
Full size table

We found a higher probability of having a change in the effective population size in Africa under MSAT (model 1 + model 3 have p ~ 0.63), whereas a signature of constant population size under MTDNA (model 2 + model 4 have p ~ 1). When looking at the posterior distribution of the effective size in Africa under MSAT, we note that we have an ancestral population of ~400,000 individuals that instantaneously expand to ~2,800,000 ~8000 years BP and go down again to ~550,000 individuals. However, effective size 95% credible intervals are largely overlapping (model 1 and model 3, Table 5). The SRMSE of the ancestral population size, the time when the change occurred and the instantaneous size change suggest low accuracy in parameter estimation, supporting the idea that if a change occurred it did not leave a strong signature in the pattern of genetic variability (Table S6). RES-MSAT favor a constant population in Africa (P = 0.7). RES-MTDNA was not able to distinguish between the two hypotheses (both having approximately P = 0.5), but the scenario with the highest probability (MOD3, P = 0.46) corresponds to a scenario with population size change in Africa. Values obtained for RES-MTDNA suggest that if a population size change occurred (MOD3), it would have started from an ancestral population of ~29,000 individuals, followed by an instantaneous population expansion to ~531,000 individuals ~4000 years BP and then a population reduction to ~120,000 individuals. Finally, we note that the EBSP highlight a complex demography in Africa with a strong population expansion ~5000 years BP and a following recent bottleneck (~600 years BP; Fig. 4).

Fig. 4
figure 4

Extended Bayesian skyline plot (mtDNA) computed in Africa, France and Spain. X axis: times in years; y axis: effective population size. Median values are shown in black, 95% high posterior density in red

Full size image

Independently of the model and of the markers, modern effective population size in Africa is ~3 to ~60 times higher than in Europe (Table 5), consistently with the EBSP on mtDNA (Fig. 4). EBSP and ABC (under MSAT, MTDNA, RES-MSAT and RES-MTDNA) both highlighted a recent expansion in Spain and France starting from a low effective population size (few thousands individuals, independently on the model considered, Table 5).

Species distribution modeling

We built a SDM for M. spretus based on the known presence localities of this species (Fig. 5). The AUC value in both analyses (Europe and Africa) is higher than 0.985, which is considered to correspond to a useful predictive model. Under present day climatic conditions, the model reveals high climatic suitability for this species across most the Iberian Peninsula (except the western and northern parts, and the Pyrenean Mountains) and coastal Mediterranean part of France, in lowland northern Morocco ( < 1000 m) from the Atlantic and Mediterranean coasts to the Atlas mountains, and in northern Algeria and Tunisia on a small coastal band of approximately 50–100 km large.

Fig. 5
figure 5

Species distribution modeling of M. spretus as estimated by Maxent for present day conditions (actual), for the last glacial maximum (LGM) and for the mid-Holocene. For paleoclimatic data, estimates are 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

Full size image

For the LGM, the CCSM and MIROC models both inferred patches of climatically suitable habitat for M. spretus on the Atlantic coast of Morocco, and in the southern part of the Atlas Mountains in Morocco. According to the CCSM model high habitat suitability could also have occurred, as small patches, on the Mediterranean coast at the Moroccan-Algerian border, and the Algerian-Tunisian border. Climatically suitable habitat for this species also occurred at this period in the Iberian Peninsula (except the western and northern parts, and the Pyrenean Mountains) and the coastal Mediterranean part of France.

For the mid-Holocene, the distribution of suitable habitat for M. spretus is predicted to be more or less similar to present day conditions, except lower habitat suitability in the Rif and High Atlas Mountains, and in the coastal Mediterranean part of France.

Discussion

When and how did the Algerian mouse arrive in Europe?

Currently, for a small terrestrial species, the Strait of Gibraltar is an important barrier to dispersal, being 14 km wide at its narrowest point and exceeding 200 m in depth. However, the fossil evidence suggests a North African origin of M. spretus (Stoetzel 2013) and recent spreading to southwestern Europe (Dobson 1998; Palomo et al. 2009). Previous molecular data (few allozyme markers, RFLP and CR sequences) suggested strong genetic differentiation between European and North African populations (Boursot et al. 1985; Desmarais 1989), as well as a likely recent colonization of Europe from North Africa with a strong founding effect (Desmarais 1989; Jacquart 1986). However, the timing and colonization pathways were not determined.

Here we confirmed the strong genetic isolation between African and European populations for both mitochondrial and nuclear datasets (Figs. 2 and 3, Table 3). By sampling individuals from both France and Spain, we have further shown a significant genetic divergence between these two areas with mitochondrial and nuclear data.

To gain more insights into the evolutionary history of M. spretus, we reconstructed the pathways of colonization of Europe from North Africa by means of ABC. Our mtDNA supported a serial colonization of Europe (i.e., first Spain and then France), whereas nuclear data conferred equal posterior probability to the serial and independent colonization scenario. We argue that the serial colonization scenario is more likely than the independent colonization model because the cross-validation experiment suggest that the high posterior probability value obtained with mtDNA is unlikely to occur if the independent colonization was the true model. Consistently, France has less genetic diversity at both markers compared with Spain. Moreover, coalescent argument can explain why nuclear data were not able to discriminate between the two scenarios under these demographic settings: the colonization of Europe has been recent and it started from a large ancestral population in North Africa with a relative mild bottleneck (see below). Both factors limited the influence of drift particularly at the nuclear level, where the effective population size is four times higher than at the mtDNA level (considering a sex ratio ~1). This explains why mtDNA had more power to distinguish between models in our particular case. Another explanation of the difficulty for ABC analyses to distinguish the two scenarios (serial versus independent colonization scenario) could be due to inadequate sampling in Spain, as almost all samples used are from Madrid. Madrid is located on a high plateau in Central Spain, and does not match the typical climatic niche of M. spretus (Fig. 5). Thus, it was probably colonized recently. It could have been easier to discriminate the two scenarios in our ABC analyses if we had included more specimens from the South of Spain, and particularly from coastal regions. Finally, we highlight that modeling the African metapopulation as panmictic, although simplifying the ABC analyses, may come at some risks. Previous works have shown that the shape of the gene genealogy of lineages sampled from a single deme, a pool of demes, or a pool of single chromosome per deme can vary as a function of the extent of differentiation among demes from the metapopulation (Wakeley 1999; Stadler et al. 2009; Mona et al. 2013, 2014). Several biases can arise if metapopulation structure and the exact sampling scheme are not correctly modeled (Maisano Delser et al. 2016; Mazet et al. 2016), but the influence of this choice on the outcome of the model selection procedure is difficult to evaluate. For this reason, we adopted a resampling strategy to mimic a scatter sample sensu Wakeley (1999) hence removing the influence of population structure on our modeling approach. Results of both RES-MSAT and RES-MTDNA confirmed that the serial colonization model is largely favored. It would be interested to test using simulations, which of the two approaches leads to more accurate results (pooling of demes vs scatter sample), but this is beyond the scope of this work.

Based on our ABC analyses, the microsatellites data support a colonization of Spain and France very close in time, between ~5000 and ~3000 years BP for MSAT, respectively, and between ~13,000 and ~11,000 years BP for RES-MSAT. Credible intervals are quite large and more data are needed to refine the estimates. These parameters were poorly estimated with ABC from mtDNA. All our analyses (neutrality tests, EBSP and ABC) highlight a recent demographic expansion in Spain and France. Our ecological niche modeling analyses suggest that suitable climatic conditions occurred throughout most central and eastern Iberian Peninsula during the last 6000 years (Fig. 5). This would have favored the demographic and spatial expansion of the Algerian mice after its arrival in Europe. It would have first spread throughout Spain and then in the South of France. The lower genetic diversity observed in French M. spretus populations compared with Spanish ones is consistent with this hypothesis, as well as the branching pattern observed in the NJ tree (Fig. 3c) and the haplotype network obtained for the Cyt b gene (Figure S5A). This hypothesis is further strengthen by the study of Jacquart (1986) who observed, based on allozyme data, a gradual diminution of the enzymatic polymorphism in M. spretus from southern Spain to France, which could be explained by founder events.

Our genetic and ecological niche modeling data are consistent with fossil data (Palomo et al. 2009) and favor an anthropogenic translocation from North Africa to the Iberian Peninsula via Neolithic navigators. During the Neolithic and after, numerous contacts between the two shores of the Mediterranean Sea occurred (seafaring evidence attested at ca 7000–9000 ya), due to the development of navigation (Mulazzani et al. 2010; Souville 1998; Zilhão 2014). Moreover, our genetic data show that the ancestral population size (before expansion) comprised a high number of individuals, reinforcing the idea that the Algerian mouse did not colonize the Iberian Peninsula through rafting of a few individuals, but via recurrent temporal anthropogenic translocations. These translocations are likely to have been accidental, mice being transported with agricultural and timber products (Brahmi et al. 2010; Dobson 1998). Several other mammals were also probably translocated by humans from one shore of the Mediterranean Sea to the other during the Neolithic period (Dobson 1998). Moreover from the beginning of the Neolithic to the present day, the turnover of mammal species is almost complete in the Mediterranean islands: most of the mammal fauna has been exterminated and replaced by a diverse assemblage of introduced species (Vigne 1992).

At first sight, it could be surprising that the Algerian mice colonized Europe via anthropogenic translocation several thousand years ago and that its dispersal between the two continents did not continue until today since maritime trade has increased. A similar phenomenon was previously observed in another rodent species, the wood mouse (Lalis et al. 2016; Stoetzel 2017). One hypothesis could be that at the beginning of the Holocene these two species were more commensal than today. Molecular and zooarchaeological data showed that the commensal species Mus musculus only reached Western Europe and North Africa during the first millennium BC and onwards, related to the generalization of maritime trade (Bonhomme et al. 2011). The arrival of the house mouse at this time could have led to a shift in the degree of commensality of M. spretus. Our mtDNA network analyses show that nearly all specimens of M. spretus from Europe cluster together in the network, with two exceptions (two specimens from Spain cluster with African specimens). This result suggests that one main event of colonization of Europe occurred several thousand years ago, but that several limited exchanges between Africa and Spain could also have occurred more recently.

Our microsatellite results show a strong affinity between the Spanish population from Madrid (locality 49 on Fig. 1) and the Moroccan Cherouda (locality 11) and Tetouan (locality 32) populations. These two Moroccan populations are located in the Tingitane Peninsula, close to the Gibraltar Strait, where the sea distance between Africa and Iberia is the lowest. Even if the navigation is not easy in this area due to the presence of strong ocean currents and counter currents, human navigation probably occurred across the Strait of Gibraltar at the Neolithic (Souville 1998; Zilhão 2014). Thus, it could be hypothesized that the Algerian mice colonized Iberia from the Tingitane Peninsula.

Can we detect a recent population bottleneck in North African populations linked to the development of anthropogenic activities and/or the arrival of the house mouse?

The El Harhoura 2 cave succession (Morocco, region of Témara) showed a continuous fossil record for M. spretus from 107 ± 7 kyr BP to 5800 years BP, with no significant change in abundance through time (Stoetzel et al. 2011a). Our climatic niche modeling analyses show that, according to the MIROC model, climatically suitable habitats (high or medium probability of occurrence) probably occurred around Témara throughout the Pleistocene (Fig. 4), which could explain the demographic stability observed in El Harhoura 2 cave.

An important molar size decrease was observed between the El Harhoura 2 Pleistocene and Holocene levels, including the most recent ones (5800 years BP) on one hand, and modern Moroccan mouse populations on the other hand, which contrasts with the stability observed through the long-term fossil record of the El Harhoura 2 cave. It was hypothesized that this could be the result of the competition with the house mouse M. m. domesticus, which only reached northwestern Africa and the Iberian Peninsula in the last millennium BC (Bonhomme et al. 2011), together with anthropogenic changes in landscape (Stoetzel et al. 2013). Important differences in the physiological, spatial and social behavior of the two species exist. It was shown that when food resources are scarce or variable in time, M. m. domesticus is favored and is the only of the two species being able to colonize partly covered habitats (Cassaing and Croset 1985). Moreover, in wild enclosure where both species co-occur only M. m. domesticus is able to reproduce normally (Orsini et al. 1982). Thus, we could expect to find a recent bottleneck (after 3000 ya) in Algerian mice African populations, after the arrival of the house mouse. According to our ABC microsatellite analyses, we found a higher probability of having a change in the effective population size in Africa (p ~ 0.65), with an instantaneous expansion followed by a bottleneck. However, when looking at the posterior distribution of the effective sizes, we note that 95% credible intervals are largely overlapping. Results of our tests for recent genetic bottlenecks based on allele frequency data are also inconclusive: tests based on the IAM mutation model suggested potential signatures of past genetic bottlenecks in samples collected in Algeria and Morocco, whereas those based on the most realistic SMM or TPM mutation models were significant for few or no localities. Based on our ABC analyses, a signature of constant population size was observed for mtDNA based on all specimens (p ~ 1). However, the EBSP highlight a more complex demography in Africa with a strong population expansion ~5000 years BP (in agreement with our climatic niche modeling analyses) and a following recent bottleneck (~600 years BP; Fig. 4 and Figure S6). RES-MSAT favor a constant population size (p ~ 0.7), and RES-MTDNA tend to favor a population size change. Results on the recent demography of M. spretus in North Africa are therefore not conclusive. Our ABC framework could have failed to retrieve the recent bottleneck because we modeled only a change of effective size in North Africa: we detected the oldest expansion as it had more time to shape the genetic diversity. Another explanation is that likelihood approaches have more power to detect more recent events than approximate methods (Girod et al. 2011): bottlenecks are particularly difficult to identify when the Ne before the decrease is large and the intensity of the bottleneck is mild, as in our case (Table 5). If we assume that this bottleneck occurred at the mtDNA level, it remains to be understood why nuclear data do not support the same conclusion. One possibility is that nuclear loci having a larger Ne than mtDNA, it renders harder to detect recent changes of effective size. The power to detect subtle and/or recent changes in the effective size increases with the number of loci studied: genome-wide approaches (i.e., whole-genome sequencing, Rad-seq, transcriptomics) will help us to shed light on this issue in the very next future.

Is the high morphological variability observed in present day populations due to isolation?

Geographical differences in body size, molar or skull morphology have been documented across the modern distribution range of M. spretus (Darviche and Orsini 1982; Darviche et al. 2006; Khammes-El Homsi and Aulagnier 2010; Palomo et al. 1985; Valenzuela-Lamas et al. 2011). The observed morphological variability could be hypothesized to be due to relative isolation between populations, leading to divergence (founder effect, drift) and possibly local adaptation. Our genetic results suggest that these phenomena could easily explain the morphological differentiation observed between the Spanish, French and North African populations.

Our climatic niche modeling analyses reveals that the LGM cold and dry period was unfavorable for M. spretus: climatically suitable habitat occurred mostly on the Atlantic coast of Morocco, and possibly also as small patches on the Mediterranean coast at the Moroccan-Algerian border, and at the Algerian-Tunisian border. Based on these results, one could expect to find genetically distinct populations in Algeria and Atlantic Morocco. However, our genetic data do not indicate a strong phylogeographic structure in North Africa. To further test this hypothesis, we run an additional model under the ABC framework (data not shown): we simulated an ancestral population in North Africa that around LGM splits into two. After ~15,000 years of separation, the two populations merged ~3000 years BP. For each model (1–4), we run the corresponding split and merge alternative. All pairwise comparison for both MSAT and MTDNA returned a posterior probability ~0.50, suggesting that either we do not have enough discriminatory power with the present data, or the two hypotheses do not lead to large differences in the predicted diversity. Further test and more genomic data will help clarify this issue, but our present data do not allow us to conclude that the geographical variability in morphology observed in North Africa (Khammes-El Homsi and Aulagnier 2010; Stoetzel et al. 2013) can be explained by long-term isolation. Numerous environmental factors such as diet, habitat, elevation, temperature, precipitation, urbanization and pollution can affect morphometric variation in rodents (Caumul and Polly 2005; Khammes-El Homsi and Aulagnier 2010; Nunes et al. 2001). More detailed analyses incorporating these factors are necessary to fully understand factors triggering phenotypic variability in the Algerian mice.

To conclude, our genetic and ecological niche modeling data are consistent with fossil and human context data, and favor an anthropogenic translocation of the Algerian mouse from the Tingitane Peninsula to the Iberian Peninsula via Neolithic navigators. More genetic data should be added to strengthen this finding. Once arrived in Spain, suitable climatic conditions would then have favored the dispersion of the Algerian mice to France. Finally, the morphological differentiation observed between the Spanish, French and North African populations could be due to founder effect and possibly local adaptation. This article highlights the importance of both climate and human-mediated changes on the evolutionary history of mammal species in the Mediterranean Basin.

Data archiving

Sequence data available from GenBank: accession numbers MH241060 to MH241525 for Cytb, MH241526 to MH242042 for CR, MH242043 to MH242499 and MH253776 to MH253817 for CO1 sequences.

Microsatellite data available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.3sv4634.