Comparative phylogeography of amphibians and reptiles in Algeria suggests common causes for the east-west phylogeographic breaks in the Maghreb

A series of phylogeographic studies in the Maghreb identified a repeated pattern of deep genetic divergence between an eastern (Tunisia) and western (Morocco) lineage for several taxa but lack of sampling in Algeria made it difficult to know if the range limits between the eastern and western lineages were shared among taxa or not. To address this question, we designed a comparative phylogeographic study using 8 reptile and 3 amphibian species with wide distribution in the Maghreb as models. We selected species where previous studies had identified an East-West phylogeographic divide and collected sampled in Algeria to 1) examine whether the simple East-West divergence pattern still holds after filling the sampling gap in Algeria or if more complex diversity patterns emerge; 2) if the E-W pattern still holds, test whether the limits between the E and W clades are shared between species, suggesting that common historical process caused the E-W divergences; 3) if E-W limits are shared between species, use information on the age of the divergence to identify possible geological or climatic events that could have triggered these E-W differentiations. We found that the E-W pattern was generally maintained after additional sampling in Algeria and identified two common disjunction areas, one around the Algeria-Morocco border, the other one in Kabylia (central Algeria), suggesting that common historical mechanisms caused the E-W divergences in the Maghreb. Our estimates for the times to most common recent ancestors to the E and W clades span a wide range between the Messinian salinity crisis and the Plio-Pleistocene limit (except for one older split), suggesting different origins for the initial divergences and subsequent preservation of the E and W lineages in common climatic refugia in the west and the east of the Maghreb.


Introduction
While we still lack a proper comparative phylogeography approach in the Maghreb, the existing phylogeographic studies identify a seemingly repeated geographical pattern of deep lineages divergence between Tunisia in the east and Morocco in the west of the region. This E-W divergence is verified in the majority of the studied taxa: reptiles (Natrix maura, [34,35], Chalcides ocellatus, [36], Trogonophis wiegmanni [37], Ptyodactylus oudrii, [38,39], Timon tangitanus / pater, [40,41], amphibians (Hyla meridionalis, [42,43], Pelophylax saharicus, [44]), scorpions [45], rodents (Jaculus orientalis, [46], Meriones shawii complex, [47]), molluscs: Cornu aspersum [48]. In birds the pattern does not seem to be as general but it has been recovered in Dupont's Lark [49] and Crested Lark [50] at least. This common pattern is currently difficult to interpret because in most cases there are few or no sample from Algeria. As a consequence, the locations of the suture zones between the E and W clades are unknown. The reasons for the lack of sampling in Algeria are mainly the current political unrest and safety issues, resulting in a large sampling gap that currently precludes our understanding of the origin of the E-W genetic divide in Maghreban lineages.
In order to improve our knowledge on the origin of current biodiversity patterns in the Maghreb, we designed a comparative phylogeography approach in the Maghreb, specifically aiming at reducing the sampling gap in Algeria. We propose that reptiles and amphibians constitute a relevant model for this study, firstly because previous studies allowed us to select several taxa showing deep E-W genetic divergence, secondly because their relatively lower dispersal ability relative to birds or large mammals make them more sensitive to landscape barriers or climatic oscillations. The aim of this study were thus: 1) to examine whether the simple East-West divergence pattern that divides the populations of many species into two main clades distributed east and west of the Maghreb still holds after filling the sampling gap in Algeria or if more complex diversity patterns with different geographical distributions emerge; 2) if the E-W pattern still holds, to test whether the limits between the E and W clades are shared between species, suggesting that common historical process caused the E-W divergences; 3) if E-W limits are shared between species, to use information on the age of the divergence to identify possible geological or climatic events that could have triggered these E-W differentiations.
To carry out this study, we selected species that fulfil two criteria: 1) they occur throughout most of the Maghreb and 2) they have been the subject of previous phylogeographic studies and exhibit deep divergence between an eastern and western clade. We conducted targeted sampling in Algeria to produce original sequence data that were merged with sequence available in Genbank. We use mitochondrial DNA (mtDNA) as a proxy for evolutionary units because it allowed us to investigate a larger number of species (11, see material and methods) with our available resources. We are fully aware of the risks of single-markers studies [51] and of the benefits of multi-locus approaches but multi-locus analyses would have constrained us to use a reduced number of species. Moreover, most of former studies were mainly or exclusively based on mtDNA, so mtDNA is the only way to merge all existing data together. Even if several evolutionary mechanisms can cause discordance in divergence patterns between mtDNA and the nuclear genome [52], the most common pattern is by far general congruence between the two types of markers [53,54], making us confident that our general conclusion will not be affected by this problem.

Ethic statement
Specimens from Morocco were collected under the authorisation n˚15 HCEFLCD / DLCDPN / DPRN / CFF delivered by the « Haut-Commissariat aux Eaux et Forêts et à la Lutte Contre la Désertification ». The list of protected species of reptiles and amphibians in Algeria is has been determined by the Executive Decree n˚12-235 05/24/2012. Among the species we have used, the following are protected under Algerian law: Acanthodactylus savignyi, A. blanci, Chalcides ocellatus, Timon pater, T. tangitanus. However, the legislative authority supposed to deliver capture licenses for protected species according to the Law n˚14-07 relative to Natural Resources is not identified in any official document and Algerian biologists therefore have to collect tissue samples or voucher specimens without legal permits. Most tissue samples used for this study have been obtained by non-invasive processes (tail tips or muscle samples from road-killed specimens). Voucher collection has followed practice of Algerian herpetologists. Field work in protected areas has been approved, and often supported in the field, by local authorities (see list in Acknowledgments).

Sampling
A total of 323 tissue samples (mainly muscle samples, sometimes buccal swabs) mostly collected from the field in Northern Algeria from 2013 to 2015 and deposited in the CEFE-CNRS tissues collection in Montpellier (BEV collection) were used in this study (S1 Table). The samples belong to 11 taxa, 8 reptiles and 3 amphibians. All these taxa had been subjected to previous phylogeographic studies revealing an E-W divergence in mtDNA: Acanthodactylus erythrurus [3], Chalcides ocellatus [36], Hemorrhois hippocrepis [55], Natrix maura [34], Podarcis vaucheri [56][57][58], Ptyodactylus oudrii, [38,39], Timon pater / tangitanus [40,41]. The samples origin and their GenBank accession codes are detailed in S1 Table. DNA extraction, amplification and sequencing DNA was extracted with the Qiagen DNeasy Blood and Tissues kit following the manufacturer's instructions. Standard protocols were used for Polymerase Chain Reaction (PCR) amplifications with the Sigma kit. The choice of the molecular marker sequenced for each species was dependant on previous published studies for this species. We used sequences of ND4 and adjacent tRNA for Acanthodactylus erythrurus, Podarcis vaucheri and Natrix maura, CytB for Chalcides ocellatus, Hemorrhois hippocrepis and Discoglossus pictus, CO1 for Hyla meridionalis and Pelophylax saharicus, 16S RNA for Timon pater / tangitanus and 12S RNA for Ptyodactylus oudrii, allowing us to merge our sequences with the sequences already published. PCR conditions and primers are detailed in S2 Table. Purification and sequencing of PCR products were carried out by Eurofins MWG Germany. The chromatograms were checked using Codon code aligner version 4.2.7 for Windows (CodonCode Corporation, Dedham, MA, USA). The alignment of the sequences was carried out with KlustalW implemented in Mega version 7.0 [59] and verified by eye.

Phylogenetic analyses
Rooted trees were produced with two approaches. Firstly, RaxML [60] with partitioning and the online platform PhyML [61] with automatic AIC criterion based substitution model selection and NNI option for tree improvement were used for maximum likelihood reconstructions. The support of the nodes was evaluated by bootstrapping with 1000 replicates. As GTR is the only nucleotide substitution model available in RAxML, GTR + I + G was applied to all partitions. Secondly, MrBayes 3.2 [62] was used for Bayesian inference (BI). Two separate runs of 6 million generations were conducted with a sampling every 500 trees. Four Monte Carlo Markov chains were initiated for which the starting topology was randomly selected. The convergence of the runs was visualized using Tracer 1.4 [63]. The first 10% trees were eliminated as the burn-in period. The majority consensus tree was computed from the remaining trees. For BI methods, best evolution models and partitioning were calculated using Partitionfinder [64]. Fifty-six different models were tested and the Bayesian information criterion (BIC) was used for the final choice of model. Both reconstruction methods gave similar topologies for all taxa except Timon and we retained the BI topology for inference and mapping of the deepest lineages (see below). For Timon we did not manage to obtain the same topology with and without our short sequences so we inferred the topology using long sequences from Genbank only and placed the new Algerian lineage (only represented by our short sequences) in this tree in a second step. For P. oudrii and T. wiegmanni, we added longer Genbank sequences to the main alignment comprising our sequences (12 S for P. oudrii, 16 S for T. wiegmanni, see S1 Table) to improve the nodal robustness.

Time calibration
The time to most recent common ancestor (TMRCA) of the eastern and western mitochondrial clades were estimated as follows: time calibrated trees were obtained with Beast v1.8.2 [65] under an uncorrelated lognormal relaxed clock; we used a Yule prior on rates of cladogenesis for almost all taxa since the trees included outgroup species but the constant size prior was used for P. saharicus; four runs were conducted for each taxon (parameters of chain's lengths and substitution rates are given in Table 1). For each time calibration tree, we used published substitution rates: [66] for A. erythrurus, medium evolution rate (see [67]) for H. meridionalis and D. pictus, [39] for P. oudrii, [68] for P. vaucheri, [69] for T. pater and T. tangitanus, [70] for P. saharicus, except N. maura for which we used published time calibrated nodes [71]. The four runs were first analysed in Tracer 1.4 [63] to check the convergence of the chains and then combined in LogCombiner v1.8.2 [72]. TreeAnnotator v1.7.0 [73] was used to produce the chronograms. We performed the dating analysis only on 8 of the 11 taxa (see Table 1). For C. ocellatus the analysis was already done by [36]. For H. hippocrepis and T. wiegmanni because of non-convergence of the runs we did not manage to have accurate results even after several runs (>4).

Lineage delimitation and genetic diversity
To estimate genetic diversity for each species we used Dnasp [74] to compute the following indices: S: number of segregating sites, H: number of haplotypes, Hd: haplotypic diversity. To assess if an E-W divide is the main process shaping mitochondrial diversity in our taxa, we used the deepest split in the phylogeny of every taxon to identify the two main clades and mapped the distribution of these two clades. To explore further significant structure in mitochondrial diversity, we used the automatic barcode gap discovery (ABGD) method [75]. Although this program is used for species delimitation, we use it to objectively describe distinct evolutionary units regardless to their taxonomic rank. The ABGD approach is based on the principle that divergence between individuals belonging to the same "species" (here lineage) is smaller than divergence between individuals from different "species". This approach is designed to work with short DNA sequences, hence it is appropriate for our short sequence data. The analyses were run with default settings. However, as recommended by the program's author when the run with the default gap width gives only one partition, the relative gap width (X = 1.5 by default) was decreased to 1.2 for Trogonophis wiegmanni and to 1.0 for Acanthodactylus erythrurus. The most stable distribution value of the initial partition was used to decide the most likely number of distinct lineages. We spatially displayed the results using ARC map10 to get an explicit map of lineage distribution, then, the geographical limits between linages were plotted using the mid distance between the nearest individuals of adjacent lineages.

Mitochondrial diversity and geographical distribution of major clades
Delimitation of mtDNA clades by ABGD was concordant with previous studies but our additional sampling in Algeria revealed previously overlooked lineages in several taxa (   Sequence divergence between clades and diversity within clades are given in Table 2. In spite of the diversity of phylogeographic patterns uncovered by sampling in Algeria, the previously identified general pattern of an E-W divergence corresponding with the major break in mtDNA diversity (deepest node within North Africa) still holds for most taxa (see . In addition, a more complex situation emerges alongside the E-W general pattern in Acanthodactylus (an endemic lineage in the Atlas mountain of Morocco is of uncertain position and might be basal), Podarcis (a deeply divergent lineage from the Atlas Mountain of Morocco groups with a high support with the E clade, see REF) and Ptyodactylus (a basal lineage found in the pre-Sahara of Algeria is of uncertain position).

Location of E-W suture zones
Algerian populations belong to the E or W main lineages depending on the taxon (Fig 14). For most species the limits between the E and W lineages span a narrow longitudinal range but for Timon, specimens of the E lineages occupy the north of Algeria while populations of the W lineage are found along the Saharan Atlas in Algeria, resulting in a north-south phylogeographic divide within Algeria that contrast with the situation with the other taxa.
Two repeated patterns are identified for the location of the E-W phylogeographic breaks: the Moulouya suture zone around the Algeria-Morocco border and the Kabylian suture zone in the eastern part of Algeria (Fig 14). Six taxa out of 11 obey the "Moulouya pattern", with suture zones located between 3.4 and 0.1˚W around the Algerian-Moroccan boundaries; four species have suture zones in the Kabylia region between 3.6 and 5.3˚E. Hyla meridionalis is the only species that deviates from these patterns as its disjunction area is located further east around 7.5˚E. As can be seen from Fig 13, the limits in the Kabylian suture zone are more tightly spaced than around the Moulouya valley.  Table 3. Three main diversification periods can be identified ( Table 3). The first period, represented by P. oudrii, N. maura, and Timon.spp., with an old diversification dated to the Miocene, long before the Messinian crisis. The second diversification event corresponds with the Messinian salinity crisis (A. erythrurus and C. ocellatus). The last period corresponds with more recent divergence within the Maghreb after the Messinian crisis between 2 and 3.6 MYA (D. pictus, H. meridionalis, P. saharicus, P. vaucheri)

Geographical distribution of internal clades and secondary phylogeographic suture zones
As can be seen from Figs 1 and 2, the secondary suture zones between the internal clades delimited by ABGD do not follow a general trend as internal clade limits are spread over the whole of Algeria. The only common pattern is that most sutures are E-W breaks, reinforcing Comparative phylogeography of amphibians and reptiles in the Maghreb the idea that for Mediterranean or mesic species, longitudinal fragmentation of ancestral areas was the main force shaping divergence in mtDNA.

How significant is the persistence of the E-W pattern after inclusion of Algerian samples?
The fact that we identified an E-W split as the most basal divergence within North Africa in most species (and that this E-W split structures most of the mitochondrial diversity in all of them) should not be interpreted as a support for the generality of the E-W divergence pattern in North Africa. We have chosen species where the E-W pattern had been previously recovered and we wanted to know if this pattern was disrupted by inclusion of Algerian samples. However, the persistence of this E-W divergence was not the only possible outcome of our study. In theory, new Algerian lineages could place at the basis of the mtDNA tree, leaving Moroccan and Tunisian samples in a monophyletic clade, if parts of Algeria harboured lineages that are older than the separation of Tunisian and Moroccan populations. The fact that we never recovered this pattern supports the idea that allopatry in E and W refugia is an important factor A rapid survey of published phylogeography for Maghreban amphibians and reptiles uncovers one other general patterns in addition to the main E-W divergence: the deepest divergence is located in Morocco, with the rest of the Maghreb colonized by one of the Moroccan lineages (Mauremys leprosa: [76]; Daboia mauritanica: [77]; Bufotes boulengeri: [78]. The other taxa either follow the general E-W divergence pattern: Testudo graeca: [79], Salamandra algira [80][81], or their phylogeography within North Africa was too poorly resolved to assign them to one of these situations: Sclerophrys mauritanica: [82]; Psammodromus algirus: [83].
In addition, more complex patterns have also emerged aside the main E-W divergence in North Africa. For example, in Timon spp, the disjunction corresponds with the western break area in the northern part of the range but isolated populations along the Saharan atlas to eastcentral Algeria belong to the western clade. These populations are geographically separated from the northern ones in Algeria carrying the eastern mtDNA haplotypes by the arid areas of the High Plateaux. Isolated populations genetically assigned to the Eastern clade are sometimes found within the range of western populations, for example in Podarcis. Last, basal clades

Concordance of suture zones between E and W lineages in several taxa suggest common historical processes
The limits of the E and W clades are concentrated in two main areas: a suture zone around the boundaries of Algeria and Morocco, broadly corresponding with the previously suggested suture zone of the Moulouya valley and one newly identified disjunction zone in the Kabylia region (see [34,44]). Both areas are also suture zones between ABGD lineages for several species. This repeated pattern, supported by results of other studies, indicate that common historical processes shaped the current divergence around these two areas. Note that there is no apparent difference in TMRCA of E and W lineages between the two main disjunction zones (see Table 3).
The Kabylian break is located around the Soumam Valley in the centre-east of Algeria. It had not been identified before this study, even if Podarcis vaucheri [56] and Cornu aspersum [84] seemed to conform to this pattern. We show here that this area also constitutes a major divide for Chalcides ocellatus, Hemorrhois hippocrepis (Fig 1) and Pelophylax saharicus (Fig 2), hence for four out of 11 taxa. The Kabylian break is also the second main disjunction area for Natrix maura and Timon pater (Fig 1).
Six of our 11 taxa conform more or less strictly to the western break pattern (see Fig 1), while this area is also a secondary break for Hyla meridonalis and Pelophylax saharicus (Fig 2). This suture zone has been identified by several previous studies [34,44,70,85]. These studies suggested that this suture zone was a direct consequence of the Moulouya River valley and its arid corridor, acting as a barrier for Mediterranean species (but see below).

Geographical barriers and glacial refugia
Our results suggest that the transition zone is wider than the Moulouya River valley itself (see Figs 1 and 2), suggesting that it cannot be the direct consequence of this valley. In addition, the W suture zone affects species with very diverse ecology, some of them occur in the Moulouya River valley itself (for ex. Pelophylax, Natrix maura, Discoglossus, Acanthodactylus erythrurus). It is thus obvious that the western suture zone cannot be explained by the Moulouya River itself acting as a barrier to dispersal. The Moulouya River is located in the Rifan corridor, a former marine corridor joining the Mediterranean east of the Rif Mountains with the Atlantic near present-day Rabat [86,87]. The width of this marine corridor was at least 60 km. The corridor opened between~8.0 and 7.6 Ma, and the maximum sea level was reached in the late-Tortonian and decreased at the Tortonian-Messinian limit (7.2 MYA) but the connexion between the Mediterranean and the Atlantic continued during the Messinian. The final closure The Rifan corridor was undoubtedly a major geographic barrier that split the Maghreb for a long period (Miocene to mid Pliocene) and several E-W divergences that we have identified can be dated to this period (see Table 3). The fact that the western break zone coincides with the area of the Rifan Corridor and that the molecular dating of these splits corresponds broadly to this period suggests that the distribution of many species were fragmented on either side of this corridor. This major past geographical barrier could explain East-west divergence in the North and North-South divergence for other species in the southern part of the corridor (Middle Atlas).
The Kabylian break is located around the Soumam valley. Although, like the Moulouya, the Soumam is a relatively dry valley compared to adjacent areas, several species that are currently affected by the Kabylian break inhabit the valley itself. We thus suggest that the Kabylian break is not a direct consequence of the presence of the Soumam River either. Unfortunately, the palaeogeography of the area is poorly documented, but no major barrier comparable to the Rifan Corridor has been identified in Kabylia.
As mentioned above (introduction), the role of the glacial / interglacial cycles in generating biodiversity below or around the species level has been thoroughly documented in Europe and North America. Nothing comparable exists for Northern Africa but we show here with a Comparative phylogeography of amphibians and reptiles in the Maghreb comparative phylogeography approach a good congruence between the Kabylia suture zones and glacial refugia identified by [20] on the basis of intra-specific phylogeographies of plant species.
The major paleo barrier of the Rifan Corridor could thus explain a fair proportion of ancient divergence events in Maghreban species but the origin of many cladogenesis events remain elusive. However, the conservation of these old divergent lineages requires a long term Comparative phylogeography of amphibians and reptiles in the Maghreb separation or a repeated expansion / contraction of the ranges of these lineages. Given the congruence in Kabylia between the location of suture zones and the location of plant refugia, the hypothesis of repeated contraction to glacial refugia seem the most likely scenario to explain the current diversity at species and lineage scales.

Evolutionary significance of mitochondrial diversity
The lack of contemporary geographical barrier and the narrowness of the contact zones suggests the possibility of competitive exclusion between lineages. Such reciprocal exclusion can be neutral (drift in local demes) or selective (if the mtDNA lineages do not freely interbred and exclude each other). We need a multilocus approach to examine these hypotheses. Firstly, we need to verify that the mtDNA lineages that we have identified correspond with real evolutionary units, as profound mitochondrial divergence can sometime persist without nuclear divergence (see [91]). We already know that some of the mitochondrial divergence we have studied are mirrored in nuclear genes: for Pelophylax, the central and W lineages differ in allozymes [92], (no sample from the E lineage were included) and for Podarcis E and W lineages are highly divergent in allozymes [67]. Secondly, we need to assess whether these evolutionary units persist in the face of gene flow (implying restriction to dispersal) or are associated with speciation events. Indeed, several of the "species" that we have examined are probably made of several species, implying that current taxonomy is inadequate.
A better understanding of the role of speciation events in generating the pattern of diversity that we have exposed here would also have important consequences in conservation. If some Comparative phylogeography of amphibians and reptiles in the Maghreb of the lineages with restricted distribution that we identified in Algeria correspond to valid biological species, the level of endemism and conservation status of the Algerian herpetofauna would be significantly changed. Indeed, the western and eastern lineages of Timon are already widely recognized as two species (T. tangitanus in the west and T. pater in the east) even if this might require confirmation. A recent detailed investigation on Tarentola mauritanica in Morocco [93] identified several unrecognized species and Ptyodactylus oudrii has been suggested to constitute a species complex [39]. Many similar cases probably await discovery in Algeria, where the amount of detailed studies on small vertebrates is much lower and ABGD identified multiple lineages in several "species" that we investigated. While some of these lineages undoubtedly correspond with truly intraspecific diversity, a multilocus study is needed in their contact zones to figure out the level of reproductive isolation between these lineages.

Impact of poorly resolved phylogenies
The mitochondrial phylogenies we used were sometimes poorly supported. While in several taxa, it was a consequence of us using short sequences and the original studies reported better supported phylogenies, in Trogonophis, the obtained phylogeny lacks strong nodal supports (posterior probability = 0.55), however a recent study resulted in the same overall topology with highly supported nodes [94]. In the case of Timon spp, the original study [42] reported a Comparative phylogeography of amphibians and reptiles in the Maghreb well-supported phylogeny but we were unable to recover the same topology with our own analyses of the same dataset whatever the method used. These cases highlight the need for improved phylogenetic resolution through longer mitochondrial sequences and addition of multilocus data.
Supporting information S1 Table. The list of used tissue samples, their geographic origin and GenBank accession numbers. Latitudes (LAT) and Longitudes (LON) are given in WGS84 geodetic system. The locus is specified only when there is concatenation of 2 loci for the phylogenetic