Aude Lalis
email: lalis@mnhn.fr

Raphael Leblois

Emmanuelle Stoetzel

Touria Benazzou

Karim Souttou

Christiane Denys

Violaine Nicolas

Phylogeography and demographic history of Shaw's Jird (Meriones shawii complex) in North Africa

Keywords:

Palaeoenvironmental data and climatic reconstructions show that the Mediterranean ecoregion of North Africa underwent drastic ecological changes during the Pleistocene. Given its rich palaeontological record, North Africa is a pertinent region for documenting the role of climate change and human mediated-habitat changes on the demography and genetic structure of faunal species. In the present study, we collected data from this species in Morocco, Algeria, and Tunisia, and we combined molecular (mitochondrial and nuclear DNA sequences, microsatellites), fossil, palaeoenvironmental, and human context data to propose an explanation for the fluctuations of populations belonging to the Meriones shawii complex in the past. Genetic and fossil data both indicate a strong bottleneck in Moroccan populations at the Middle Holocene (last interglacial optimum) compared to the Late Pleistocene. Our mitochondrial DNA data suggest a diversification event within Morocco corresponding to the 130-125 kya interglacial optimum. Given that (1) major demographic changes in the M. shawii complex coincide with the interglacial optimums, and (2) the impact of human activities on the landscape and faunal communities was moderate during the Middle Holocene (beginnings of the Neolithic culture), our results demonstrate that climate, rather than anthropogenic influences, likely explains the M. shawii complex population decline in the Holocene.

INTRODUCTION

The Mediterranean ecoregion of North-western Africa is the biogeographical area extending westerly to the Atlantic Ocean, north to the Mediterranean Sea, and east and south to the Sahara desert. Palaeoenvironmental data and climatic reconstructions show that this region underwent drastic ecological changes during the Pleistocene, mainly related to the alternation of expansion/reduction in the size of the Saharan desert, reduction/expansion of the Mediterranean vegetation, and reduction/development of water ponds, lakes, and rivers [START_REF] Brun | Micoflores et pal eov eg etations en Afrique du Nord depuis 30 000 ans[END_REF][START_REF] Brun | R eflexions sur les Pluviaux et Arides au Pl eistoc ene sup erieur et a l'Holoc ene en Tunisie[END_REF][START_REF] Hooghiemstra | Vegetational and climatic changes at the northern fringe of the Sahara 250,000-5000 years BP: evidence from 4 marine pollen records located between Portugal and the Canary Islands[END_REF][START_REF] Jolly | Simulated climate and biomes of Africa during the late Quaternary: comparison with pollen and lake status data[END_REF][START_REF] Trauth | Trends, rhythms and events in Plio-Pleistocene African climate[END_REF][START_REF] Stoetzel | Taphonomy and palaeoecology of the Late Pleistocene to Middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF]. These environmental changes could have led to allopatric differentiation in separate refuges for many species. However, the impacts of climate change on the phylogeographical patterns of species remain poorly understood in the region [START_REF] Harris | Phylogeography and genetic relationships of North African Bufo mauritanicus Schlegel, 1841 estimated from mitochondrial DNA sequences[END_REF][START_REF] Faleh | Phylogeography of two cryptic species of african desert jerboas (Dipodidae: Jaculus)[END_REF][START_REF] Boratynski | The origin of two cryptic species of African desert jerboas (Dipodidae: Jaculus)[END_REF][START_REF] Husemann | Palaearctic biogeography revisited: evidence for the existence of a North African refugium for Western Palaearctic biota[END_REF][START_REF] Nicolas | Phylogeography of the North African Dipodil (Rodentia: Muridae) based on cytochrome b sequences[END_REF]. North-western Africa has yielded a rich palaeontological record allowing a better understanding of the evolution of small mammal species and communities over time in relation with Quaternary climatic changes [START_REF] Stoetzel | Late Cenozoic micromammal biochronology of northwestern Africa[END_REF]. Moreover, this region is now recognized as a major area for the emergence and dispersal of anatomically modern humans from at least 160 kya [START_REF] Mcbrearty | The revolution that wasn't: a new interpretation of the origin of modern human behavior[END_REF][START_REF] Manica | The effect of ancient population bottlenecks on human phenotypic variation[END_REF][START_REF] Smith | Earliest evidence of modern human life history in north african early Homo sapiens[END_REF][START_REF] Garcea | Successes and failures of human dispersals from North Africa[END_REF]. North-western Africa 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. Understanding how fluctuations in climate and human-mediated habitat change influenced the demography and genetic structure of historic species is not only interesting for fundamental research, but also can be employed to predict future demographic and genetic parameters of species. This is particularly true in the Mediterranean region of North Africa where species and landscape diversity is currently threatened by global climate change and increasing human activities, resulting in the progressive disappearance of certain biotopes (e.g. Mediterranean dry woodlands and steppe ecoregions; [START_REF] Olson | Terrestrial ecoregions of the world: a new map of life on earth a new global map of terrestrial ecoregions provides an innovative tool for conserving biodiversity[END_REF].

Rodents are good models for historical reconstructions of environment influences on biota [START_REF] Wang | Intraspecific molecular phylogeny and phylogeography of the Meriones meridianus (Rodentia: Cricetidae) complex in northern China reflect the processes of desertification and the Tianshan mountains uplift[END_REF]. This is a result of their short generation time, rapid mitochondrial (mt)DNA substitution rate, relatively limited dispersal ability, and strong associations with particular habitats. This leads to informative contemporary patterns of genetic variation [START_REF] Tolley | Biogeographic patterns and phylogeography of dwarf chameleons (Bradypodion) in an african biodiversity hotspot[END_REF][START_REF] Fedorov | Comparative phylogeography and demographic history of the wood lemming (Myopus schisticolor): implications for late quaternary history of the taiga species in Eurasia[END_REF][START_REF] Nicolas | Comparative phylogeography of two sibling species of forest-dwelling rodent (Praomys rostratus and P. tullbergi) in west Africa: different reactions to past forest fragmentation[END_REF][START_REF] Bryja | Plio-pleistocene history of west african sudanian savanna and the phylogeography of the Praomys daltoni complex (rodentia): the environment/geography/genetic interplay[END_REF][START_REF] Russo | Bioregion heterogeneity correlates with extensive mitochondrial DNA diversity in the namaqua rock mouse, Micaelamys namaquensis (Rodentia: Muridae) from southern africaevidence for a species complex[END_REF][START_REF] Wang | Intraspecific molecular phylogeny and phylogeography of the Meriones meridianus (Rodentia: Cricetidae) complex in northern China reflect the processes of desertification and the Tianshan mountains uplift[END_REF][START_REF] Boratynski | Molecular and ecological signs of mitochondrial adaptation: consequences for introgression?[END_REF]. The Jird Meriones is an appropriate model for investigating the roles of humans and climate change in shaping faunal genetic diversity and distribution. The Meriones shawii complex belongs to the Muridae family and Gerbillinae subfamily. This complex is widely distributed from Morocco to Algeria, Tunisia, Libya, and Egypt, and along the west side of the Nile [START_REF] Aulagnier | Meriones shawi[END_REF]in IUCN 2012). It is not found in mesic environment, such as forests, grasslands, wetlands, lakes, and rivers [START_REF] Aulagnier | Meriones shawi[END_REF]in IUCN 2012), and it avoids rocky basins. The M. shawii complex remains poorly known, according to its systematics, ecology, and geographical distribution. Based on external and cranial measurements, some studies recognize two valid species in this group: M. shawii and Meriones grandis [START_REF] Cabrera | Meriones grandis. IN: Algunos roedores nuevos de Marruecos[END_REF][START_REF] Pavlinov | Contribution to craniometric variation and taxonomy of jirds from the group 'shawi-grandis' of the genus Meriones (Gerbillidae)[END_REF] but with a large biometric overlapping and without clear geographical distribution. Other studies consider them only as subspecies [START_REF] Petter | R epartition et ecologie des rongeurs d esertiques (du Sahara occidental a l'Iran oriental)[END_REF][START_REF] Aulagnier | Catalogue des mammif eres sauvages du maroc[END_REF], with the question remaining unresolved.

The Meriones genus originated in North Africa at the end of the Middle Pleistocene, and has remained well represented in fossil assemblages subsequent to the beginning of the Late Pleistocene (approximately 130 kya) [START_REF] Ouahbi | Recent Quaternary fossil mammals of Chrafate and Ez Zarka. The origin of modern fauna in the Northern Rif (NW Morocco, Northern Africa)[END_REF][START_REF] Reed | A preliminary account of the rodents from Pleistocene levels at Grotte des Contrebandiers (Smuggler's Cave), Morocco[END_REF][START_REF] Stoetzel | Pleistocene and Holocene small vertebrates of El Harhoura 2 Cave (Rabat-Temara, Morocco): an annotated preliminary taxonomic list[END_REF][START_REF] Lopez-Garcia | The small mammals from the Holocene site of Guenfouda (Jerada, Eastern Morocco): chronological and paleoecological implications[END_REF][START_REF] Stoetzel | Late Cenozoic micromammal biochronology of northwestern Africa[END_REF]. All fossil Meriones remains of 'modern' morphology are assigned to M. shawii, without any mention of M. grandis. Exceptional fossil records covering the last 120 000 years are available from the El Harhoura 2 cave, in the region of T emara on the north-Atlantic coast of Morocco, a few kilometres to the south of Rabat [START_REF] Stoetzel | Pleistocene and Holocene small vertebrates of El Harhoura 2 Cave (Rabat-Temara, Morocco): an annotated preliminary taxonomic list[END_REF][START_REF] Stoetzel | Taphonomy and palaeoecology of the Late Pleistocene to Middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF][START_REF] Stoetzel | Context of modern human occupations in North Africa: contribution of the T emara caves data[END_REF]. This site has recorded the entire last climatic cycle along 11 archaeological levels, from approximately 120 kya to 5.8 kya [START_REF] Jacobs | Singlegrain OSL chronologies for Middle Palaeolithic deposits at El Mnasra and El Harhoura 2, Morocco: implications for Late Pleistocene human-environment interactions along the Atlantic coast of northwest Africa[END_REF][START_REF] Stoetzel | Context of modern human occupations in North Africa: contribution of the T emara caves data[END_REF]. All Meriones remains belong to the M. shawii complex, and represent the most abundant small vertebrate taxon in the El Harhoura 2 assemblage, whereas it is now poorly represented in the septentrional Atlantic plains, including the Rabat-T emara region [START_REF] Aulagnier | Catalogue des mammif eres sauvages du maroc[END_REF][START_REF] Aulagnier | Zoog eographie des Mammif eres du Maroc: de l'analyse sp ecifique a la typologie de peuplement a l' echelle r egionale[END_REF]. The human context in which this settlement occurred is well known: several cultures succeeded during Late Pleistocene and Holocene [START_REF] Nespoulet | Palaeolithic and neolithic occupations in the temara region (Rabat, Morocco): recent data on hominin contexts and behavior[END_REF][START_REF] Jacobs | Singlegrain OSL chronologies for Middle Palaeolithic deposits at El Mnasra and El Harhoura 2, Morocco: implications for Late Pleistocene human-environment interactions along the Atlantic coast of northwest Africa[END_REF][START_REF] Stoetzel | Context of modern human occupations in North Africa: contribution of the T emara caves data[END_REF]. Moreover, the palaeoenvironmental context is also well known: palaeoecological analysis has shown an alternation of arid and more humid periods during the Late Pleistocene, ending with a humid period during the Middle Holocene, corresponding to the last climatic optimum [START_REF] Stoetzel | Taphonomy and palaeoecology of the Late Pleistocene to Middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF][START_REF] Stoetzel | Context of modern human occupations in North Africa: contribution of the T emara caves data[END_REF]. These environmental changes are accompanied by differences in the relative proportion of small vertebrate species between the El Harhoura 2 levels. Despite these palaeoenvironmental changes, Meriones are dominant in relatively stable proportions all along the Late Pleistocene record. However, a significant decrease in the proportion of Meriones is observed in the assemblage of level 1 (Holocene: approximately 5.8 kya BP) compared to Pleistocene levels [START_REF] Stoetzel | Taphonomy and palaeoecology of the Late Pleistocene to Middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF]. It is important to note that the taphonomic study has shown no or only very low biases or perturbations, indicating that the relative proportion of species along the archaeological sequence is neither related to the type of predators at the origin of the fossil accumulations, nor other taphonomic agents [START_REF] Stoetzel | Taphonomy and palaeoecology of the Late Pleistocene to Middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF].

In the present study, we combined mitochondrial and nuclear data to assess the taxonomic status of the M. shawii complex, and also whether significant demographic changes occurred during the last 120 000 years in this taxa. The combination of molecular, fossil, palaeoenvironmmental, and human context data should allow us to propose a scenario to explain the historic fluctuations of populations within this taxon in North Africa.

MATERIAL AND METHODS

SAMPLING AND LOCALITIES

Samples of M. shawii complex were collected between January 2010 and April 2012 in Morocco and Algeria (Fig. 1; see also Supporting information, Table S1) with permission from the 'Haut Commissariat aux Eaux et forêts et a la Lutte contre la d esertification' (autorization no. 15 HCEFLCD/ DLCDPN/DPRN/CFF) in Morocco and the Ministry of Forestry in Algeria. Animals were live-captured using Sherman traps and handled in accordance with the guidelines of the American Society of Mammalogists (http://www. mammalogy.org/committees/ index.asp; Animal Care and Use Committee, 2011) and also in accordance with standard procedures for BSL3 work in the field (Federal Guidelines for Field Work, CDC 1997). 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 isofluorane, followed by cervical dislocation. The protocol was approved by Comit e Cuvier (permission no. 68.009). Moreover, four specimens from Tunisia housed in the collections of the National Museum of Natural history of Paris (France) were included in the present study (MNHN). In total, we used 178 samples representing 14 localities (details of sampling localities are provided in Table 1). Total genomic DNA was extracted from ethanol-fixed intercostal muscle tissues using NucleoSpinR 96 Tissues (Macherey-Nagel) in accordance with the manufacturer's instructions. The carcasses were fixed in formalin for later preparation as skin and skull specimens. Specimens from Morocco are temporarily housed at the MNHN and will be deposited at the Institut Scientifique de Rabat (ISR, Morocco). Specimens from Algeria are housed at the Institut of Agronomy of El Harrach (Alger, Algeria).

DNA AMPLIFICATION SEQUENCING AND GENOTYPING

We amplified and sequenced 178 individuals (1-33 per population) (Table 1) for the cytochrome b gene (cytb) fragment of the mitochondrial DNA (mtDNA) genome using primers L14723 [START_REF] Ducroz | An assessment of the systematics of Arvicanthine rodents using mitochondrial DNA sequences: evolutionary and biogeographical implications[END_REF] and H6 [START_REF] Mongelard | Molecular systematics of the Scirugnathi (Rodentia): the mithochondrial cytochrome b and the 12S rRNA genes support the Anomaluroidea (Petetidae and Anomaluridae)[END_REF]. Intron 7 of the nuclear fragment, the b-fibrinogen We constructed two median-joining haplotype networks for the cytb and BFIBR datasets using NET-WORK, version 4.500 [START_REF] Bandelt | Median-joining networks for inferring intraspecific phylogenies[END_REF]. Prior to this analysis, the existence of heterozygous positions for the nuclear gene fragment was investigated in accordance with the procedure described in [START_REF] Nicolas | A remarkable case of micro-endemism in Laonastes aenigmamus (Diatomyidae, Rodentia) revealed by nuclear and mitochondrial DNA sequence data[END_REF]. The results obtained from four individuals (MA97, MA127, MA848, MA853) were unclear (i.e. several possible haplotype pairs for each individual with similar probabilities). They were thus removed from all subsequent analyses. The Bayesian Markov chain Monte Carlo (MCMC) approach implemented in BEAST, version 1.8 [START_REF] Drummond | BEAST: Bayesian evolutionary analysis by sampling trees[END_REF] was employed to simultaneously estimate phylogenetic relationships and time to the most recent common ancestor (TMRCA). Divergence times and their credibility intervals were estimated using a relaxed clock model with branch rates drawn from an uncorrelated log-normal distribution. We used the model of sequence evolution retained by MRMO-DELTEST, version 3.04 [START_REF] Nylander | MrModeltest v2. Program distributed by the author[END_REF]) and a coalescent model with varying population sizes (skyline model). Two independent runs of 100 million iterations with burn-ins of 25% were performed. The results were inspected visually using TRACER, version 1.8 to ensure proper mixing of the MCMC. A consensus chronogram with median age estimates and 95% higher posterior density intervals was generated and visualized with TREEANNOTATOR (BEAST package, version 1.8) and FIGTREE, version 1.3.1 [START_REF] Rambaut | FigTree version 1.3.1[END_REF]. Four fossil calibrations were used to calibrate the chronogram. All calibrations were applied as log-normal prior distributions, and the means and SDs of these distributions were chosen to construct 95% confidence intervals (CIs) spanning the 90-95% Marshal indices [START_REF] Marshall | Confidence intervals on stratigraphic ranges: partial relaxation of the assumption of randomly distributed fossil horizons[END_REF] reported by the Paleobiology Database (PDB) [START_REF] Jaeger | The age of the Mus-Rattus divergence: paleontological data compared with the molecular clock[END_REF][START_REF] Pdb | The Paleobiology Database[END_REF]. These represent the estimated 95% CIs for the actual origination of a taxon based on first occurrences and stratigraphic sampling. The two first calibrations applied in the present study have been used in previous studies [START_REF] Schenk | Ecological opportunity and incumbency in the diversification of repeated continental colonizations by muroid rodents[END_REF] and the other two are derived from the PDB: (1) The Gerbillinae-Deomyinae split based on the first occurrence of Gerbillinae in the Lower Miocene fauna of Saudi Arabia (offset = 15.868, range = 16.000-23.700); (2) the Lophuromys-Acomys-Deomys split, hence the origin of Acomys, based on the earliest known Acomys fossil from Kenya (offset = 5.258, range = 5.3-29.050); (3) the Meriones-Psammomys-Rhombomys split, hence the origin of Meriones, based on the earliest known Meriones fossil from Kazakstan in the Pliocene (offset = 2.6, range = 2.6-5.3); and (4) the Gerbillus-Sekeetamys split based on the first occurrence of the genus Gerbillus in Armenia in the Pliocene (offset = 2.6; range = 2.6-5.3). The species used as outgroup were: Mus musculus (Genbank AB819920), Acomys airensis (AJ012021), Acomys cahirinus (AJ233953), Acomys chudeaui (FJ415538), Acomys cilicicus (AJ233957), Acomys dimidiatus (AJ233959), Acomys ignitus (Z96064), Acomys johannis (HM635823), Acomys minous (GU046553), Acomys nesiotes (AJ233952), Acomys percivali (EF187818), Acomys russatus (FJ415485), Acomys spinosissimus (AM409396), Acomys subspinosus (JN247673), Acomys wilsoni (EF187799), Deomys ferrugineus (FJ415478), Lophuromys flavopunctatus (EU349754), L. sikapusi (AJ012023), Desmodillus auricularis (AJ851272), Gerbilliscus robustus (AM409374), Gerbilliscus guinea (AJ430562), Gerbillurus paeba (AJ430557), Gerbilliscus tytonis (AJ430559), Sekeetamys calurus (AJ851276), Gerbillus campestris (AJ851271), Gerbilliscus gerbillus (AJ851269), Gerbilliscus henleyi (JQ753050), Gerbilliscus nanus (JQ753051), Gerbilliscus poecilops (JQ753064), Gerbilliscus simoni (GU356577), Gerbilliscus tarabuli (GU356573), Desmodilliscus braueri (AJ851273), Taterillus arenarius (AJ851261), Psammomys obesus (AJ851275), Rhombomys opimus (AJ430556), Meriones vinogradovi (VV1989001), Meriones chengi (AB381900), Meriones crassus (AJ851267), Meriones libycus (JQ927411), Meriones meridianus (AJ851268), Meriones rex (AJ851265), Meriones unguiculatus (AF119264), Meriones tristrami (KU189331), and Meriones persicus (KT949958). These phylogenetic and divergence time analyses were only performed on the cytb dataset as a result of the lack of outgroup sequences for the BFIBR gene.

Morocco Algeria

Tunisia Spain Portugal

GENETIC DIVERSITY AND POPULATION STRUCTURE

For sequence data, the nucleotide diversity, haplotype diversity, and mean number of nucleotide differences were calculated using DNASP, version 5.0 [START_REF] Librado | Dnasp V5: a software for comprehensive analysis of DNA polymorphism data[END_REF]. For microsatellite loci, the number of alleles (N A ) and the expected and observed heterozygosities (H E and H O ) were calculated using the R package ADEGENET, version 1.2-7 [START_REF] Jombart | Adegenet: a R package for the multivariate analysis of genetic markers[END_REF]. Allelic richness (R) was calculated using FSTAT. All loci were tested for Hardy-Weinberg equilibrium (HWE) and the presence of linkage disequilibrium using GENEPOP, version 4.1.3 [START_REF] Rousset | Genepop'007: a complete re-implementation of the genepop software for windows and linux[END_REF]. MICRO-CHECKER, version 2.2.1 [START_REF] Van Oosterhout | Micro-Checker: software for identifying and correcting genotyping errors in microsatellite data[END_REF] was used to test for the presence of null alleles, large allele dropouts, and scoring errors. For the microsatellite data, we employed clustering analyses to describe the genetic structure of our sample and to estimate the most likely number of genetically homogeneous clusters (K) using STRUCTURE, version 2.3.3 [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF][START_REF] Falush | Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies[END_REF]. STRUC-TURE was run using the admixture model, with localities as priors and assuming correlated allelic frequencies. Final outputs were obtained for 20 independent runs, testing K = 1 to K = 6, each with a total of 250 000 iterations and a burn-in of 150 000. The number of contributing populations was tested using the ad-hoc Evanno statistic Delta K [START_REF] Evanno | Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study[END_REF]. Membership probabilities (i.e. Q-values) of the 20 runs for K = 2 and 3 were averaged using CLUMPP, version 1.2 [START_REF] Jakobsson | CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure[END_REF] and displayed using DISTRUCT, version 1.1 [START_REF] Rosenberg | Distruct: a program for the graphical display of population structure[END_REF]. Population differentiation was inferred from the mtDNA and microsatellite data sets using F ST estimates. Values for mtDNA data were computed with ARLEQUIN, version 3.5.1 [START_REF] Excoffier | Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows[END_REF] using populations with over 12 individuals. Population differentiation was also inferred from the microsatellite data set by comparing F ST estimates [START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF] between all population pairs and among all populations, at each locus and over all loci, as computed by FSTAT [START_REF] Goudet | FSTAT: a computer program to calculate fstatistics[END_REF][START_REF] Goudet | Testing differenciation in diploid populations[END_REF] and GENEPOP, version 4.1.3. We tested whether F ST estimates were significantly > 0 by permuting multilocus genotypes among samples with FSTAT. Genotypic differentiation between populations was investigated by exact tests using Markov chain algorithms implemented in GENEPOP, version 4.1.3. [START_REF] Rousset | Genepop'007: a complete re-implementation of the genepop software for windows and linux[END_REF]. Pairwise F ST values were used to build Neighbour-joining trees using POPULATIONS, version 1.2.30b (http://bioinformatics.org/~tryphon/ populations/).

We also used spatial analysis of shared alleles (SaShA) using SaShA, version 2.0 (http://sasha.stanford.edu) to test for genetic subdivision across the geographical space for cytb and BFIBR datasets. This analysis uses spatial and haplotypic information to detect nonrandom allele distribution against an expectation of panmixia [START_REF] Kelly | A method for detecting population genetic structure in diverse, high gene-flow species[END_REF]. The test statistic describes the observed mean distance between alleles (OM). When OM is less than the expected mean (EM), alleles are considered to be aggregated. When OM is larger than EM, alleles are considered randomly distributed. A jackknife procedure identifies which alleles are strongly influencing the observed distribution. Isolation-by-distance patterns were tested with a Mantel's test [START_REF] Mantel | The detection of disease clustering and a generalized regression approach[END_REF] with 30 000 permutations. For sequence data, the test was performed using ARLEQUIN, version 3.5.1 [START_REF] Excoffier | Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows[END_REF] on the relationship between the mean number of pairwise nucleotide differences and geographical distances between sampling localities. For microsatellite data, the test was performed using GENEPOP, version 4.1.3, by regressing F ST /(1 À F ST ) between populations over the logarithm of geographical distances [START_REF] Rousset | Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance[END_REF].

DEMOGRAPHIC HISTORY

Demographic history was explored using the MIGRAINE (http://kimura.univ-montp2.fr/~rousset/ Migraine.htm) and the model with historic variations in population size [START_REF] Leblois | Maximum likelihood inference of population size contractions from microsatellite data[END_REF]. This model consists of a single, isolated panmictic population that undergoes continuous exponential change in population size starting at time T generations in the past and continuing until the moment of sampling (i.e. present). MIGRAINE uses the class of importance sampling algorithms developed by de Iorio & Griffiths (2004a, b), de [START_REF] De Iorio | Stepwise mutation likelihood computation by sequential importance sampling in subdivided population models[END_REF] and extended in [START_REF] Leblois | Maximum likelihood inference of population size contractions from microsatellite data[END_REF]. MIGRAINE was used to test for past change in population size and to estimate current and ancestral scaled population size (h = 2Nl and h anc = 2N anc l, where N and N anc are the current and ancestral haploid population size, respectively, and l is the mutation rate per generation of the whole locus) and D, the time when the demographic change starts, scaled by population size (i.e. D = T/2N). MIGRAINE was first applied on the mitochondrial and nuclear sequence data sets separately, and only for the Moroccan samples. Prior to this analysis, alignments of nuclear genes were pruned to exclude stretches with missing data at the beginning and the end of some sequences. Also, because MIGRAINE is based in the infinitely many sites model (ISM) for analysis of sequence data, different datasets were produced for both the mtDNA cytb region and for the nuclear BFIBR region to fit this model. There are two reasons why a sequence data set may not fit the ISM: sites can show more than two nucleotidic states or pairwise comparisons of sites may not comply to the four gamete test [START_REF] Hudson | Statistical properties of the number of recombination events in the history of a sample of DNA sequences[END_REF]. For one data set, we chose to systematically remove incompatible sites for all individuals. For the second, we removed haplotypes with incompatible sites. For the MIGRAINE analysis, mtDNA and nuclear DNA (nDNA) sequence data were pooled for all individuals from Moroccan populations (136 individuals). This initial data set contained 138 and 26 segregating sites, with 42 and 12 unique haplotypes in 136 and 76 individuals, for cytb and BFIBR respectively. Fitting the data sets to the ISM resulted in four modified data sets (Table 3) distinguished by the remaining number of sites, haplotypes, and individuals. All runs with MIGRAINE consisted of 1 000 000 trees, 2400 points, and two iterations. The MIGRAINE software was also used on the microsatellite data set to infer past changes in population size on the pooled Moroccan dataset. A benefit of using MIGRAINE over MSVAR [START_REF] Beaumont | Detecting population expansion and decline using microsatellites[END_REF][START_REF] Storz | Testing for genetic evidence of population expansion and contraction: an empirical analysis of microsatellite DNA variation using a hierarchical Bayesian model[END_REF], which does similar analyses, is that it allows departure from the strict stepwise mutation model by use of a generalized stepwise mutation model. All MIGRAINE analyses for microsatellites were only run on the Moroccan samples and used 20 000-200 000 trees, 2400 points, and three iterations (see the MIGRAINE manual for details concerning these settings). To convert scaled population sizes (h and h anc ) into effective population sizes (N and N anc ), and scaled times (D) into times in years (T), we used: (1) a generation time of 1 year; (2) mutation rates of 5 9 10 À4 mutations per locus per generation for all microsatellite loci [START_REF] Dib | A comprehensive map of the human genome based on 5,264 microsatellites[END_REF][START_REF] Ellegren | Microsatellite mutations in the germline: implications for evolutionary inference[END_REF][START_REF] Sun | A direct characterization of human mutation based on microsatellites[END_REF] and 10 À7 and 10 À8 mutations per site per generation for the mitochondrial and nuclear locus, respectively. Those values are similar to those found in our BEAST analyses, and also similar to other values found in the literature for rodents [START_REF] G€ Und€ Uz | Mitochondrial DNA variation in the western house mouse (Mus musculus domesticus) close to its site of origin: studies in Turkey[END_REF][START_REF] Nabholz | Strong variations of mitochondrial mutation rate across mammals, the longevity hypothesis[END_REF]; (3) the length in pb of the two sequenced loci; and (4) a correction factor for the loss of polymorphic sites that do not fit the ISM for the sequenced loci. The F-statistic analysis yields a highly significant level of population structure among all localities (F ST = 0.71, P < 0.001). The F ST value between Morocco and Algeria is high and significantly different from 0 (F ST = 0.90, P < 0.001). Pairwise F ST values between Moroccan and Algerian populations are also high (0.96 > F ST > 0.87) and significantly different from 0 (Table 2). Pairwise F ST values between Moroccan populations are lower (0.49 < F ST < 0.17) but still significantly different from 0. Lastly, the F ST value between the two Algerian populations is low and nonsignificant (F ST = 0.04, P > 0.05). The SAShA analysis shows that M. shawii complex haplotypes are significantly aggregated (OM = 13 km, expected = 508 km; P = 0.001). This is also true when only Moroccan specimens of either clade A + B (OM = 7 km, expected = 236 km; P = 0.001) or clade A (OM = 7 km, expected = 225 km; P = 0.001) are considered. For the three analyses, the haplotypeby-haplotype analysis shows that all common haplotypes are significantly aggregated and the jacknife analyses indicate the robustness of the overall results, which remain qualitatively the same when any single haplotype is removed. When all specimens are considered, a significant positive relationship between geographical and genetic distances is found (Mantel test, P < 0.001, slope of 0.736). However, no significant correlation between geographical and genetic distances is observed within clade A (P = 0.071, slope of 0.451), which is the clade with the highest sample size and the best geographical coverage. Concerning the BFIBR, the alignment of the 90 phased haplotypes required the addition of 11 gaps. According to our network analyses, two main groups of haplotypes differing by 19 mutations can be identified (Fig. 3): clade 1 groups all but two Moroccan individuals; clade 2 groups all Algerian and individuals, as well as two individuals from the Moroccan locality of Guenfouda (MA205 and MA210). The mean percentage of nucleotide differences between clades 1 and 2 is 1.9%. For a given specimen the two phased haplotypes always cluster within the same clade. There is congruence between the cytb and BFIBR data: the cytb monophyletic group A + B corresponds to BFIBR clade 1, and cytb clade C corresponds to BFIBR clade 2. The only exception concerns specimens from Guenfouda: two specimens of clade A cluster within clade 2 and one specimen of clade C clusters within clade 1. Branch lengths are proportional to the number of mutations between haplotypes (except between clades A and C for the mtDNA dataset). Colours represent the country of sampling (red = Morocco, blue = Algeria, green = Tunisia). Two-or three-letter codes refer to the locality of sampling (for localities codes, see Table 1). The analysis of 137 M. shawii complex genotypes reveals a variable degree of polymorphism across the six locations studied. A total of 249 alleles are detected across all loci, ranging from a maximum of nine alleles detected in population for locus MS-10 to a minimum of one allele detected in populations GF and IF for locus MS-3, MS and OB for locus MS-2, and OB for locus MS-4 (see Supporting information, Table S2). Considering the 10 loci together, the population from OB shows the highest allele mean number (4.6), whereas the IF population has the lowest allele mean number (3.8). Observed heterozygosity ranges from 0.355 (OB) to 0.528 (GF), whereas the expected heterozygosity ranges from 0.411 (OB) to 0.558 (MS) (Table 1). No significant linkage disequilibrium is detected after adjusting for multiple comparisons. Significant deviations from HWE are observed in 18 out of 60 locus/population combinations (see Supporting information, Table S2). The highest posterior probability of the model given by STRUCTURE is obtained for K = 4, and the highest posterior probability of Delta K is obtained for K = 2 (Fig. 4). One cluster is comprised of the Algerian population, whereas all of the Moroccan populations are part of the second cluster. Some Moroccan individuals show a high proportion of their genome belonging to the Algerian cluster; however, they appear to belong to the Moroccan clusters for K = 3 and for K = 4 (Fig. 4). The F ST analysis between Moroccan and Algerian populations is highly significant (F ST = 0.19, P < 0.001). Indeed, pairwise F ST values between Moroccan and Algerian populations are non-negligable (0.17 > F ST > 0.24) and the differentiation is also always significant (Table 2). By contrast, pairwise F ST values between Moroccan populations only are low (0.03 < F ST < 0.15) and differentiation is not always significant. We did not find any significant correlation between F ST /(1 À F ST ) values and the logarithm of the geographical distance in separating Moroccan populations (P = 0.582, slope of 0.065).

RESULTS

DEMOGRAPHIC HISTORY

All MIGRAINE results are presented in Table 3. Analyses of modified data sets (mtDNA and nDNA sequence) fit to the ISM give nonhomogeneous results, suggesting that modifications to fit the ISM may have been too sharp to retain most of the information within the original sequences. This could be a result of multiple recurrent mutations or recombination events that occurred on those two loci. However, for both genes, a significant degree of past contraction is detected by one of the two analyses, although parameter estimation differs slightly for the two markers. First, the nuclear BFIBR locus provides more information on the time when the past contraction began and suggests that it was a relatively recent contraction [i.e. D 95% CI = (0.01-1.4)], whereas the mtDNA cytb sequence analysis shows very wide CIs for this estimation [D 95% CI = (0.008-∞)]. Second, point estimates and the CI of the scaled mutation rates (h and h anc ) show larger values for the mtDNA data set than for the BFIBR one. However, given that the effective population size in the number of genes is four times higher for the nuclear locus but the mutation rate for the mtDNA cytb locus may be much higher than that for the nuclear BFIBR locus, those differences may not be unrealistic. MIGRAINE results from the microsatellite data show a highly significant signal of past contraction with parameter estimations that are not always concordant with those obtained with the mtDNA and nDNA data sets. As expected, because of the greater number of loci, estimation precision is much better than for the DNA sequence analyses. The magnitude of past contraction is better inferred with a point estimate for the ratio of current over ancestral population sizes of 0.034 (0.0021-0.042, 95% CI).

The conversion of our estimates of scaled parameters into unscaled demographic parameters is shown in Table 3. Population contraction was substantial (from several hundreds of thousands individuals to a few hundred or thousands) and likely occurred a few thousand years ago. Interestingly, all analyses leads to very similar results in terms of ancestral population size, the parameter for which there is the more information in the data.

DISCUSSION

SPATIAL GENETIC STRUCTURE

With both mtDNA and nDNA, we recover two main genetic clades: (1) an eastern clade, which includes all individuals from Tunisia and Algeria, as well as some individuals from the most eastern Moroccan locality (Guenfouda), and (2) a western clade, comprising most specimens from Morocco. The cytb K2P genetic difference between these two clades is 7.8%, which often corresponds to the genetic distance observed between sister species in rodents [START_REF] Baker | Speciation in mammals and the genetic species concept[END_REF][START_REF] Boratynski | The origin of two cryptic species of African desert jerboas (Dipodidae: Jaculus)[END_REF][START_REF] Ndiaye | Evolutionary systematics and biogeography of endemic gerbils (Rodentia, Muridae) from Morocco: an integrative approach[END_REF]). It is interesting to note that the specimens from Guenfouda, which group with either the eastern or the western clades, differ in the two datasets. For cytb, MA234 specimen clusters with the eastern clade, whereas all other specimens from this locality cluster with the western clade (Figs 2 and 3), but, for BFIBR, MA205 and MA210 specimens cluster with the eastern clade (Fig. 3), whereas all other specimens cluster with the western clade. For a given specimen, the two phased BFIBR haplotypes always cluster within the same clade. In animal species, nuclear and mitochondrial markers differ in effective size, presence or the absence of recombination and biparental vs. maternal inheritance. Thus, it is not surprising to find different results for mitochondrial and nuclear data. The pattern observed in the locality of Guenfouda could be explained by mtDNA introgression as a result of past hybridization followed by back-crosses with paternal lineages. mtDNA introgression is not rare in nature, particularly for rodents [START_REF] Bryja | Plio-pleistocene history of west african sudanian savanna and the phylogeography of the Praomys daltoni complex (rodentia): the environment/geography/genetic interplay[END_REF][START_REF] Boratynski | Molecular and ecological signs of mitochondrial adaptation: consequences for introgression?[END_REF]. Of the 10 microsatellite loci studied, three loci (MS-6 MS-7 and MS-8), have specific fixed alleles in the Algerian population (not shown) that can be considered as diagnostic alleles of the Algerian populations. Four Moroccan specimens are associated with the Algerian population (specimens MA924 and MA928 from Ouled Boughadi and specimens MA960 and MA961 from Ifrane; Fig. 4) because of the presence of some of these specific alleles (M924 and M928: homozygous on locus MS-7; M960 and M961: heterozygous on locus MS-6). This result is similar to the one found in North African sympatric mice species (Mus musculus domesticus and Mus spretus; [START_REF] Orth | Hybridation naturelle entre deux esp eces sympatriques de souris Mus musculus domesticus L. et Mus spretus Lataste[END_REF] and is the result of a single locus effect of the 10 loci analysis. Our results tend to favour the hypothesis of two sister species in North Africa with possible hybridization between them. Additional analyses with more specimens and more genetic markers are necessary to confirm this hypothesis.

The suggested two species would have diverged in the early Pleistocene. This type of genetic structuring, with one or several eastern clades (Algeria/Tunisia) and one or several western clade (Morocco), has been found in many other vertebrate species [START_REF] Arano | PHYLOGEOGRAPHICAL STRUCTURE OF THE MERIONES SHAWII COMPLEX 275 molecular and morphological evidence[END_REF][START_REF] Recuero | Mitochondrial differentiation and biogeography of Hyla meridionalis (Anura: Hylidae): an unusual phylogeographical pattern[END_REF][START_REF] Barata | Comparative phylogeography of northwest African Natrix maura (Serpentes: Colubridae) inferred from mtDNA sequences[END_REF][START_REF] Fritz | Mitochondrial phylogeography of Testudo graeca in the Western Mediterranean: old complex divergence in North Africa and recent arrival in Europe[END_REF][START_REF] Nicolas | Phylogeography of the North African Dipodil (Rodentia: Muridae) based on cytochrome b sequences[END_REF][START_REF] Vences | New insights on phylogeography and distribution of painted frogs (Discoglossus) in northern Africa and the Iberian peninsula[END_REF]. In most of these cases, the Moulouya River and/or the wide arid valley extending along much of the river (except close to the estuary) was referred to as the major geographical and climatic barrier responsible for this structuring. However, [START_REF] Fritz | Mitochondrial phylogeography of Testudo graeca in the Western Mediterranean: old complex divergence in North Africa and recent arrival in Europe[END_REF] suggested that this bipartite east-west differentiation may be too simplistic, reflecting incomplete sampling rather than true geographical differentiation. Indeed, most studies lack proper sampling in Algeria.

General phylogeographical patterns in Morocco are diverse. Although several species show a lack of phylogeographical structure [START_REF] Batista | Genetic variation within Bufo viridis: evidence from mitochondrial 12S and 16S RRNA DNA sequences[END_REF][START_REF] Harris | Phylogeography and genetic relationships of North African Bufo mauritanicus Schlegel, 1841 estimated from mitochondrial DNA sequences[END_REF][START_REF] De Pous | Integrating mtDNA analyses and ecological niche modelling to infer the evolutionary history of Alytes maurus (Amphibia; Alytidae) from Morocco[END_REF], others show the opposite [START_REF] Pinho | Contrasting patterns of population subdivision and historical demography in three western Mediterranean lizard species inferred from mitochondrial DNA variation[END_REF][START_REF] Recuero | Mitochondrial differentiation and biogeography of Hyla meridionalis (Anura: Hylidae): an unusual phylogeographical pattern[END_REF][START_REF] St€ Ock | Mitochondrial and nuclear phylogeny of circum-Mediterranean tree frogs from the Hyla arborea group[END_REF][START_REF] Fritz | Mitochondrial phylogeography of Testudo graeca in the Western Mediterranean: old complex divergence in North Africa and recent arrival in Europe[END_REF][START_REF] Beukema | Biogeography and contemporary climatic differentation among Moroccan Salamandra algira[END_REF][START_REF] Nicolas | Phylogeography of the North African Dipodil (Rodentia: Muridae) based on cytochrome b sequences[END_REF]. The strong phylogeographical pattern within Morocco can be largely explained by climate changes during the Plio-Pleistocene, which may have created multiple refugia. In North Africa, Morocco has a high physiographical complexity, with several large mountain ranges oriented primarily east-west (Anti-Atlas, High Atlas, Middle Atlas, Rif). These mountains may have aided the survival of populations through altitudinal shifts, allowing them to track suitable microclimates as the general climate fluctuated. Moreover, as a result of its geographical location, Morocco is influenced by both the Atlantic and the Mediterranean sea, resulting in a wide range of climates. These characteristics make it unlikely that Morocco offered a single homogeneous and continuous refuge area throughout the Pleistocene. Instead, the variable distribution and fragmented nature of suitable habitats would have favoured the occurrence of multiple refugia isolated from one another. Our mtDNA data show the existence of two genetic clades within Morocco that likely diverged in the Middle Pleistocene (414 kya; range: 240-617 kya). However, presently, these clades are broadly sympatric, rendering the localization of the two possible refugia impossible to determine. Within Morocco, our mtDNA results show some spatial clustering of haplotypes, although this is less evident in the nuclear (BFIBR and microsatellite) datasets. Isolation-by-distance patterns and geographical aggregation of haplotypes suggest spatially limited dispersal and thus geographically restricted gene flow. Significant mtDNA spatial patterns and nonsignificant nuclear ones may indicate sex-biased dispersal with males dispersing more frequently than females. This pattern is common in mammals (Lawson [START_REF] Handley | Advances in our understanding of mammalian sex-biased dispersal[END_REF] and has been observed in M. unguiculatus [START_REF] Liu | Effects of supplemental food on the social organization of Mongolian gerbils during the breeding season[END_REF].

DEMOGRAPHIC HISTORY OF THE SPECIES BASED ON PALEONTOLOGICAL AND GENETIC DATA

According to El Harhoura 2 cave fossil data, Meriones are dominant in relatively stable proportions all along the Late Pleistocene record, although a significant decrease in the proportion of Meriones is observed in the assemblage of level 1 (Holocene approximately 5.8 kya) compared to Pleistocene levels [START_REF] Stoetzel | Taphonomy and palaeoecology of the Late Pleistocene to Middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF]. Our genetic data confirm a bottleneck in Moroccan M. shawii complex at least during the Holocene, and possibly slightly earlier. However, our dating estimate that marks the beginning of this bottleneck has a large uncertainty. Our data suggest that the pattern observed in the fossil record represents demographic changes. Two hypotheses can be proposed to explain M. shawii complex population size changes over time:

(1) environmental changes caused by humans, and (2) climatic changes unrelated to human activities. Despite the fact that the Holocene level of El Harhoura 2 is attributed to the Neolithic culture, few faunal remains represent domestic taxa. Furthermore, neither commensal species (Mus domesticus, Rattus rattus, Rattus norvegicus), nor the opening of habitat under anthropic intervention (agriculture) occurred [START_REF] Stoetzel | Context of modern human occupations in North Africa: contribution of the T emara caves data[END_REF]. Thus, at this period in the region of T emara, the impact of human activities on the landscape and the faunal communities should have been moderate. The difference in species proportions observed in recent times could therefore be a result of environmental factors. Indeed, the Holocene level of El Harhoura 2 (approximately 5.8 kya) is characterized by a more humid climate than previous periods [START_REF] Stoetzel | Taphonomy and palaeoecology of the Late Pleistocene to Middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF] corresponding to the last climatic optimum (warm and humid climate, development of Mediterranean vegetation in North Africa, and reduction of the desert areas) and extant Meriones species avoid semi-humid and humid regions [START_REF] Aulagnier | Catalogue des mammif eres sauvages du maroc[END_REF][START_REF] Aulagnier | Zoog eographie des Mammif eres du Maroc: de l'analyse sp ecifique a la typologie de peuplement a l' echelle r egionale[END_REF]. Thus, we could easily argue a climatic context to explain the M. shawii complex population size reduction in the Holocene. Despite pronounced climate changes during the period approximately 20-120 kya [START_REF] Stoetzel | Taphonomy and palaeoecology of the Late Pleistocene to Middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF][START_REF] Stoetzel | Context of modern human occupations in North Africa: contribution of the T emara caves data[END_REF], no significant change in the proportion of M. shawii complex in the small mammal assemblage of El Harhoura 2 was observed [START_REF] Stoetzel | Taphonomy and palaeoecology of the Late Pleistocene to Middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF]. Unfortunately, the previous glacial optimum (Eemian, approximately 125-130 kya) is not recorded at El Harhoura 2 (or the corresponding level was not yet reached), preventing any comparison between these two interglacials (Eemian and Holocene). Between them, during the Late Pleistocene glacial period, the relatively humid interstadials were far less humid and their impact would have been milder than the last climatic optimum. According to our phylogenetic tree analyses, a splitting event occurred within this species in Morocco at the end of the Middle Pleistocene (TMRCA of the two clades: 0.292 Mya, with large CI: 0.120-0.559 Mya). This may correspond to the isolation of M. shawii complex populations in two distinct arid refugia during the 125-130 kya glacial optimum; however, the large CIs hinder the precision of any estimates.

CONCLUSIONS

The combination of three different data sets (fossil remains, molecular data, and palaeoenvironmmental data) points towards a climatic rather than anthropogenic influence with respect to explaining M. shawii complex population size reduction in the Holocene. The ability of the M. shawii complex to thrive in cultivated crops and irrigation structures (to the point where it is considered to be a crop pest) [START_REF] Aulagnier | Meriones shawi[END_REF]in IUCN 2012) must have appeared very recently. The present study demonstrates the advantage of using integrative and multidisciplinary approaches for investigating species abundance and demography. We also provide relevant information on the genetic structure of M. shawii complex, structured geographically in two major clades, which can be of importance to those concerned with pest control and human disease because this species is considered to represent an economically significant pest [START_REF] Aulagnier | Meriones shawi[END_REF]in IUCN 2012).

Figure 1 .

 1 Figure 1. Map of the sampled geographical localities.

Figure 2 .

 2 Figure 2. Time calibrated phylogeny inferred using BEAST, depicting phylogenetic relationships within the species Meriones shawii. Numbers at nodes represent clade posterior probabilities. To improve clarity, outgroup taxa are not shown.

Figure 3 .

 3 Figure 3. Minimum spanning network of Meriones shawii cytb [mitochondrial (mt)DNA, left] and BFIBR (nuclear DNA, right) haplotypes. Circle sizes are proportional to the number of similar haplotypes observed in the data set.Branch lengths are proportional to the number of mutations between haplotypes (except between clades A and C for the mtDNA dataset). Colours represent the country of sampling (red = Morocco, blue = Algeria, green = Tunisia). Two-or three-letter codes refer to the locality of sampling (for localities codes, see Table1).

Figure 4 .

 4 Figure 4. Meriones shawii complex populations clustering based on STRUCTURE Bayesian inference (K = 2-4); the graph illustrates the log posterior probabilities of the microsatellite data [Ln P(K)] for each number of genetic groups (K) and the number of contributing populations was tested using the ad-hoc Evanno statistic (Delta K) for K = 1-6. Each individual is represented by a thin vertical line, which is partitioned into K-coloured segments that indicate an individual's estimated membership fraction in K clusters.

Table 1 .

 1 Meriones shawii sampling locality details and genetic diversity showing the number of individuals analyzed (N), averaged values of observed heterozy-

	E ), mean number of alleles (N A ), allelic richness (A) within-population coefficient of inbreeding (FIS), nucleotide dive-	nDNA diversity:	Microsatellite diversity mtDNA diversity: cytb BFIBR	Longitude N H O A F IS H E N A N Pi H D H D Pi N		À7.97 29 0.488 0.462 4.20 3.30 À0.038 28 0.01397 0.918 0.883 0.00404 8	À2.01 18 0.528 0.556 4.20 3.59 0.080 20 0.00866 0.632 0.835 0.01061 7	À5.04 11 0.518 0.493 3.80 3.54 À0.004 12 0.00484 0.318 0.827 0.00393 11	À7.01 33 0.509 0.464 4.00 3.03 À0.081 33 0.00289 0.845 0.844 0.00216 5	À6.72 33 0.355 0.411 4.60 3.20 0.048 33 0.00406 0.748 0.530 0.00222 6	À9.77 À 2 À	À9.71 À 2 À	À7.97 À 1 À	À8.92 À 3 À	À9.00 À 2 1		4.5 13 0.43 0.558 4.10 3.74 0.009 14 0.00275 0.912 0.964 0.01001 4	2.97 -25 0.00446 0.923 -		9.6 -2 1	9.5 -1 2	137 178 45
	gosity (H O ) and expected heterozygosity (H	D ) rsity (Pi), and haplotypic diversity (H		Latitude Locality Code	Morocco	32.19 BenGuerir BG	MM Guenfouda GF	33.42 Ifrane IF	31.84 Sour El Az $EA	33.07 Ouled Boughadi OB	29.83 Aglou AG	31.32 2okm $ Essaoulra E$	31.83 20 km N Manakech MA	3040 Qued Souss OS	31.53 Sidl El Moctar $EM	Algeria	35.7 MSila MS	34.3 Hadjeb El Djmel HEA	Tunisia	35.07 Garat An Njila GAN	35.03 Sidi Bouzid SB	Total/mean	mtDNA, mitochondrial DNA.
																				262-279

© 2015 The Linnean Society of London, Biological Journal of the Linnean Society, 2016, 118,

Table 2 .

 2 Mitochondrial DNA sequence (above diagonal) and microsatellite (below diagonal) pairwise F ST statistics between Meriones shawii populations. (P < 0.001)

	Population	BG	GF	IF	SEA	OB	MS	HEA
	BG		0.3309*	0.2153*	0.2233*	0.2634*	0.8713*	0.8803*
	GF	0.0549		0.4651*	0.4932*	0.4711*	0.9139*	0.9131*
	IF	0.0966	0.1135*		0.2983*	0.2344*	0.9509*	0.9399*
	SEA	0.0389	0.067	0.1178*		0.1678*	0.9625*	0.9535*
	OB	0.0675	0.1113*	0.1542*	0.038		0.9514*	0.9445*
	MS	0.1877*	0.1767*	0.2420*	0.1976*	0.2050*		0.0401
	*Significant values.

© 2015 The Linnean Society of London, Biological Journal of the Linnean Society, 2016, 118, 262-279 PHYLOGEOGRAPHICAL STRUCTURE OF THE MERIONES SHAWII COMPLEX 263

© 2015 The Linnean Society of London, Biological Journal of the Linnean Society, 2016, 118, 262-279

ACKNOWLEDGEMENTS

The present study was supported by ANR 6 eme extinction ANR-09-PEXT-004, MOHMIE, 'Modern Human installation in Morocco, Influence on the small terrestrial vertebrate biodiversity and Evolution' (CD) and the project CMEP TASSILI MDU 09MDU755. RL has been partially funded by two ANR projects (EMILE 09-blan-0145-01 and IM-Model@CORAL.FISH 2010-BLAN-1726-01) and by the INRA (Project INRA Starting Group 'IGGiPop'). ES was funded by both the ANR-09-PEXT-004, MOH-MIE and the LabEx BCDiv. Fieldwork was made possible through collaboration with the 'Institut Scientifique de Rabat' (A. El Hassani, M. Fekhaoui) and the 'Haut Commissariat aux Eaux et Forêts et Lutte contre la D esertification'. We are grateful to all collectors, particularly Abderrahmane Mataame, Hicham El Brini, Arnaud Delapr e, L ea Bourg, Loubna Tifarouine, and Laurent Granjon. Molecular analyses were financially supported by the 'ATM MNHN: Taxonomie mol eculaire, DNA Barcode & Gestion Durable des Collections', the 'Service de Syst ematique Mol eculaire' of the MNHN (UMS 2700, Paris, France), and the network 'Biblioth eque du Vivant' funded by the CNRS, MNHN, INRA, and CEA (Genoscope). Part of this work was carried out using resources from the MNHN UMS 2700, INRA GENO-TOUL, and MIGALE bioinformatics platforms and the computing grid of the CBGP laboratory. We thank K. Gavrilchuk for helpful comments on an earlier version of the manuscript and three anonymous reviewers for their helpful comments.

SUPPORTING INFORMATION

Additional Supporting Information may be found online in the supporting information tab for this article: Table S1. List of specimens included in the molecular analyses. Table S2. Genetic variability of Meriones shawii populations at microsatellite loci. Genetic variability of M. shawii at six locations and 10 loci. Sample size; number of alleles; observed heterozigosity (H O ) and expected heterozigosity (H E ); and probability of significant deviation from Hardy-Weinberg equilibrium (P, Markov chain procedure, a=0.05) are reported. Significant P values are shown in bold.