Phylogenetic and phylogeographic patterns in Mediterranean worm lizards - cryptic diversity in Blanus and Trogonophis

Phylogenetic and phylogeographic patterns in Mediterranean worm lizards Cryptic diversity in Blanus and Trogonophis Filipa Leão Sampaio 2012 Phylogenetic and phylogeographic patterns in Mediterranean worm lizards Cryptic diversity in Blanus and Trogonophis Filipa Leão Sampaio Mestrado em Biodiversidade, Genética e Evolução Departamento de Biologia (FCUP) CIBIO-UP – Centro de Investigação em Biodiversidade e Recursos Genéticos 2012 Orientador D. James Harris Coorientadores Daniele Salvi Ana Perera Phylogenetic and phylogeographic patterns in Mediterranean worm lizards Cryptic diversity in Blanus and Trogonophis Filipa Leão Sampaio Dissertação de Mestrado em Biodiversidade, Genética e Evolução apresentada à Faculdade de Ciências da Universidade do Porto Todas as correções determinadas pelo júri, e só essas, foram efetuadas. O Presidente do Júri, Porto, ______/______/_________ AGRADECIMENTOS À Fundação para a Ciência e Tecnologia (FCT) pelo financiamento do trabalho: PTDC/BIA-BDE/74349/2006 (JH) e BII/CIBIO/INTBIO/2009 (FS). Ao James pela proposta do tema e pela oportunidade de fazer trabalho de campo em Marrocos – que me possibilitou ver mais do que um pedaçinho da cauda dos bichos com que estive a trabalhar! Para além disso, todos os comentários foram importantes para a concretização desta tese. Ao Daniele pela boa disposição e por todo o empenho e ajuda no laboratório, os comentários e ensinamentos ao longo deste ano. À Anna toda a paciência, disponibilidade e simpatia com que posso sempre contar. Parte do trabalho de laboratório foi realizado no Institut de Biologia Evolutiva (CSICUPF) em Barcelona, no âmbito do Programa de Mobilidade Erasmus Estágio, onde estive durante três meses (de 02-11-2011 a 03-02-2012). Fiquei sob a supervisão do Salvador Carranza, a quem tenho que agradecer toda a disponibilidade, ajuda e simpatia com que me recebeu. E também as amostras que cedeu para este estudo. Para além disso, durante a minha estadia foi bom contar com o grupo que me recebeu em Barcelona, por toda a ajuda e boa disposição no laboratório, o que tornou tudo mais fácil. À Sara Costa, com quem fiz o estágio de licenciatura no CIBIO, que começou por fazer uma parte do trabalho de laboratório com amostras usadas neste trabalho. Aos meus colegas do CIBIO, em especial ao João Maia, Joana Mendes e Duarte pela ajuda no laboratório e nas minhas dúvidas. E claro a todos, pelos momentos de descontracção e patetices nas longas viagens de metro e nos intervalos de trabalho. Aos meus amigos, em especial à Catarina, Fernando e Sandra, por partilharem comigo a atribulada tarefa de acabar uma tese. Ao Joel pelo incentivo, apoio e amizade. À minha família, por todo o apoio e paciência ao verem-me passar tantas horas em frente ao computador. VI V. FINAL REMARKS AND FUTURE PERSPECTIVES .............................................................................................. 57 1. CONCLUSION ............................................................................................................................................................. 59 2. FUTURE WORK .......................................................................................................................................................... 60 VI. REFERENCES ...................................................................................................................................................................... 63 VII LIST OF FIGURES Figure 1 Higher-level phylogeny of amphisbaenids based on Vidal et al. (2008), with additional data from Kearney and Stuart (2004). African amphisbaenids are paraphyletic. (Vidal and Hedges, 2009) .................................................................................................................................. 6 Figure 2 Blanus specimens. ............................................................................................................................. 9 Figure 3 Map with the region where Blanus species occur. Distribution range of each species according to IUCN (2012). ................................................................................................................ 9 Figure 4 (A) Maximum parsimony (MP) tree inferred from Blanus ND4 mitochondrial sequences. Bootstrap support for MP and maximum likelihood analyses are indicated above and below the nodes. (B) Map showing sampled localities. Figures adapted from Vasconcelos et al. (2006). ................................................................................................................................ 11 Figure 5 (A) Maximum likelihood (ML) tree based on Iberian Blanus ND4 and 16S mitochondrial haplotypes. Numbers on each node represent from top to bottom, ML bootstrap values, Bayesian posterior probabilities, maximum parsimony (MP) and minimum evolution bootstrap values. Nodes with either ML or MP bootstrap values above 70% are shown, otherwise are collapsed. (B) Map of the Iberian Peninsula showing sampling localities of Blanus cinereus. Grey dots represent populations of the southwestern clade, and black dots represent populations of the central clade. Figures adapted from Albert et al. (2007) and Albert and Fernández (2009). ......................................... 12 Figure 6 Trogonophis wiegmanni specimens. ........................................................................................ 13 Figure 7 Map with the region where Trogonophis wiegmanni occurs, with its distribution range according to IUCN (2012). .................................................................................................................. 13 Figure 8 (A) Maximum likelihood tree (ML) based on 12S and 16S rRNA sequences from Mendonça and Harris (2007). All analyses produced identical relationships to the one shown. Near the nodes, Bayesian posterior probabilities and ML bootstrap values. For both analyses, only bootstrap values above 50% are represented. (B) Map showing sampling locations of T. wiegmanni specimens analysed. Figures adapted from Mendonça and Harris (2007). .............................................................................................................................................. 15 Figure 9 Map of the study region - the Mediterranean Basin. ......................................................... 16 Figure 10 Maps of the study areas with the geographic location of the Blanus samples analysed. Samples are represented by circles coloured according to phylogenetic analyses results. Samples with (*) after code are from a published study (Albert et al., 2007). (Samples B33, B34 and Bs* are not displayed on the map since their exact locations are unknown.) Dashed areas in the central map represent the distribution range of each species according to IUCN (2012). .............................................................................................................. 37 Figure 11 Phylogenetic tree representing the relationships within Blanus, derived from the Bayesian analysis of the combined markers for a total fragment of 2303 bp [16S (474 bp), ND4 (802 bp), POMC (461 bp) and MC1R (566 bp)]. The topology is similar in ML VIII analysis (data not shown). For the major clades and lineages, Bayesian posterior probabilities (BPP) and ML bootstrap support (BP) are given above and below nodes, respectively. The tree was rooted using Diplometopon zarudnyi and Trogonophis wiegmanni. ............................................................................................................................................................. 39 Figure 12 Statistical parsimony networks representing relationships among Blanus. Haplotypes inferred from 1294 bp mitochondrial sequences (16S and ND4 – 493 bp and 801 bp, respectively). Circles represent different haplotypes with size proportional to sample frequency. Small grey circles represent missing or extinct haplotypes. ...................... 40 Figure 13 Haplotype networks recovered by statistical parsimony analyses representing relationships among Blanus, inferred from nuclear markers POMC (461 bp) and MC1R (566 bp). Circles represent different haplotypes with size proportional to sample frequency. Small grey circles represent missing or extinct haplotypes. Allelles from heterozygotes samples are represented by ‘i’ or ‘ii’ after the sample code. ............................... 42 Figure 14 Map of the study area with the geographic location of the Trogonophis samples analysed. Samples are represented by circles coloured according to phylogenetic analyses results. Dashed area in the bottom right map represents the distribution range of this genus according to IUCN (2012). ................................................................................................................. 42 Figure 15 Phylogenetic tree representing the relationships within Trogonophis, derived from the Bayesian analysis of the combined markers [12S (353 bp), 16S (465 bp), CYTB (284 bp), CMOS (350 bp), POMC (467 bp) and RAG2 (798 bp)] for a total fragment of 2717 bp. The topology is similar in ML analysis (data not shown). For the major clades and lineages, Bayesian posterior probabilities (BPP) and ML bootstrap support (BP) are given above and below nodes, respectively. The tree was rooted using Diplometopon zarudnyi and Blanus mettetali. ......................................................................................................................................... 44 Figure 16 Haplotype networks recovered by statistical parsimony analysis representing relationships among Trogonophis wiegmanni, inferred from 828 bp mitochondrial sequences (12S and 16S – 361 bp and 472 bp, respectively). Circles represent different haplotypes with size proportional to sample frequency. Small grey circles represent missing or extinct haplotypes. ....................................................................................................................... 45 Figure 17 Haplotype networks recovered by statistical parsimony analysis representing relationships among T. wiegmanni main clades, inferred from nuclear markers CMOS (350 bp), POMC (467 bp) and RAG2 (798 bp). Circles represent different haplotypes with size proportional to sample frequency. Small grey circles represent missing or extinct haplotypes. Heterozygotes were included as independent samples and are represented by ‘i’ or ‘ii’ after sample code. ............................................................................................................................... 47 IX LIST OF TABLES Table 1 Blanus samples used in this study – sample code, species, country, sampling locality and GenBank accession number for sequences from Albert et al. (2007). ................. 25 Table 2 Trogonophis samples used in this study, sample code, subspecies, country, sampling locality and GenBank accession number for sequences from Mendonça and Harris (2007). ....................................................................................................................................................... 27 Table 3 Gene, primer name, sequence, source and PCR conditions (temperature, time and number of cycles). Amplification of RAG2 fragments was done with a nested PCR with two sets of primers - first PCR using external primers 31FN.Venk and Lung.460R and a second PCR with internal primers Lung.35F and Lung.320R. The latter ones were used for sequencing. ............................................................................................................................................................ 30 Table 4 Net distance values between Blanus species and uncorrected p-distance values within Blanus species. Lower left values correspond to 16S and top right to ND4 net distance values between species sequences. Shaded diagonal values correspond to uncorrected p-distance values within species haplotypes. ............................................................... 41 Table 5 Net distance between Trogonophis clades and uncorrected p-distances within Trogonophis clades. Lower left values correspond to 12S and top right to 16S values between clades sequences. Shaded diagonal values correspond to uncorrected p-distance values within clades haplotypes. .................................................................................................................. 45 0 1 I. INTRODUCTION 2 3 1. SPECIES AND CRYPTIC DIVERSITY Up until this day the definition of ‘species’ still remains a controversial topic in evolutionary biology. However, nowadays it is broadly accepted that this term is not fixed and contemporary biological diversity has discontinuities along morphological, genetic and ecological axes (Niemiller et al., 2012). For a long time morphology has been used as the main tool for taxonomists to describe and identify species (Mayr, 1949). Nevertheless, the development of molecular techniques has uncovered unexpectedly high levels of genetic diversity. Even though biodiversity is mostly discussed at species or even higher taxonomical level, intraspecific genetic variation is also an integral part of biodiversity (Taberlet et al., 1998). Therefore, the level of genetic diversity is an important tool to be used as a complementary strategy to more traditional conservation approaches, when prioritizing populations for protection purposes (Bonin et al., 2007). Indeed, molecular data have changed the way biodiversity is perceived. For instance, it has enabled the discovery of cryptic species, i.e., species which are morphologically similar but genetically well differentiated and reproductively isolated (Bickford et al., 2007). In some cases, morphological similarity may be due to convergence and genetic markers can provide powerful information for disentangling evolutionary relationships. Even though the concept of cryptic species is not new (Winker, 2005), in the last decades molecular tools have been effective in identifying cryptic diversity, and the number of described cryptic species has been significantly increasing (Bickford et al., 2007 and references therein). Particularly interesting cases are phylogenetic studies which have found discordance between morphological and genetic differentiation in subterranean taxa. Fossorial organisms have evolved in extreme environmental settings, limiting possible adaptive responses in which organisms are able to adapt (Nevo, 2001). Thus, morphological changes associated with speciation may often be reduced or non-existant (Bickford et al., 2007), or there may be morphological convergence of adaptive characters. Therefore in subterranean taxa, species delimitation based on morphology is particularly difficult due to the possible occurrence of morphological convergence (Niemiller et al., 2012). In these cases, biodiversity assessment based only on morphological traits could be strongly biased (Lefébure et al., 2006). Amphisbaenians are a group of squamate burrowing reptiles which have been recently the focus of molecular analyses with mitochondrial markers. Recent studies have recovered highly complex phylogenetic relationships, revealing the existence of cryptic diversity, with high genetic variation in morphologically indistinguishable taxa. However, 4 these results are based mainly on mitochondrial data. Even though it used to be a common practice in phylogenetic studies, setbacks to the use of only this type of molecular marker to infer diversity and phylogenetic patterns are currently known. 2. MOLECULAR MARKERS IN THE ASSESSMENT OF CRYPTIC DIVERSITY In the last decades, many phylogenetic and phylogeographic studies of herpetofauna have relied only on mitochondrial DNA (mtDNA) sequence variation. In fact, mtDNA analyses have been considered to be adequate to assess cryptic species with similar morphology (Slade and Moritz, 1998). Mitochondrial DNA has several properties that make it suitable to infer evolutionary relationships (Avise, 2000; Avise et al., 1987). It can be easily obtained due to high copy number in cells; it has a small genome size and simple sequence organization; it is transmitted maternally, having a non-recombining mode of inheritance; evolves rapidly in animal populations and has extensive intraspecific polymorphism (Avise, 2000). Even though mtDNA may be extremely useful to address phylogenetic studies because of its characteristics, the application of a single locus approach may lead to misleading interpretations. The representation of a single locus may not necessarily reflect species evolutionary history, which could lead to a biased interpretation of the results. Also, it may not detect introgression and incomplete lineage sorting phenomena (reviewed in Ballard and Whitlock, 2004; Bazin et al., 2006; Zhang and Hewitt, 2003). Due to the downsides of the use of only mtDNA sequences, recently it has been a more common approach the inference of phylogenies based on multilocus datasets, constituted by a combination of multiple mitochondrial and nuclear markers. These are more valuable to perform more robust phylogenies and give further insights on the evolutionary relationships of the taxa under study. 3. PHYLOGENETIC INFERENCE ANALYSES The progress of laboratory procedures to amplify DNA fragments has led to an increase of DNA data in the last decades. Also, an interest on the estimation of taxa evolutionary history has arisen and improved methodologies and computational procedures have been designed to make phylogenetic inferences based on DNA sequences. There are currently several methods available to reconstruct phylogenetic trees such as maximum parsimony, maximum likelihood and Bayesian inference. These methods are used for depicting relationships on a deeper level, among species. 5 Maximum likelihood (ML) is a method of statistical inference to estimate an evolutionary tree from DNA sequences (Felsenstein, 1981) based on a chosen model of sequence evolution. It assigns quantitative probabilities to mutational events to compare possible phylogenetic trees and find the evolutionary tree which best predicts the observed data (Felsenstein, 1981; Makarenkov et al., 2006). There are various software available to perform ML analysis – PAUP (Swofford, 2003), GARLI (Zwickl, 2006), PHYML (Guindon and Gascuel, 2003) and RAXML (Stamatakis, 2006). In a Bayesian inference (BI) analysis, the inference on phylogeny is based on the posterior probabilities of a tree (Huelsenbeck and Bollback, 2001b; Huelsenbeck and Ronquist, 2001). MRBAYES (Huelsenbeck and Ronquist, 2001) implements this analysis with a Markov chain Monte Carlo (MCMC) approach to approximate the posterior probabilities of a tree distribution of topologies. Most ML and BI analyses employ models of DNA substitution, so that the appropriate model for each alignment is used in the phylogenetic analysis. To select the appropriate model, there is software available, such as MODELTEST (Posada and Crandall, 1998) or JMODELTEST (Posada, 2008). For more shallow phylogenies, at the intraspecific level, phylogenetic relationships are better represented by networks, because they offer more resolution of the relationships among haplotypes, than phylogenetic trees. This is also an important tool for dealing with genes at population level (Makarenkov et al., 2006). There are many methods and software available that produce networks, being often used in phylogenetic and population genetics studies. For instance, median-joining in NETWORK (Fluxus Technology; Bandelt et al., 1999), or statistical parsimony approach implemented in TCS (Clement et al., 2000). The method developed by Templeton et al. (1992) (TCS) collapses identical sequences into haplotypes and calculates the haplotypes’ frequency. Missing intermediates are also estimated. The statistical parsimony algorithm estimates the maximum number of differences among haplotypes which are caused by single substitution events with a 95% parsimony connection limit (by default the limit of parsimony is 95%, but it can it can be used a cut-off between 90 to 99%). This translates into the maximum number of single nucleotide mutations that can be connected in a single haplotype network; haplotypes separated by more mutational steps remain disconnected. The statistical parsimony method implemented in TCS connects haplotypes with small differences, displaying the similarities rather than the dissimilarities between the haplotypes (Clement et al., 2000; Makarenkov et al., 2006). 12 In fact, another study by Albert et al. (2007) (Figure 5) using 16S and ND4 mitochondrial sequences and one anonymous nuclear marker revealed similar results. B. strauchi was recovered as a sister group to the remaining Blanus species; North African Blanus formed a monophyletic group, which is the sister group of another clade including all Iberian haplotypes. Within the Iberian clade, two distinct monophyletic groups were recovered with genetic distance (10.5–12.4% uncorrected p-distance for ND4) – as high has those found between B. mettetali and B. tingitanus (12.3%) (Albert et al., 2007). These results led to description of a new species for the south western region of the Iberian Peninsula – B. mariae – based on molecular and morphological evidence by Albert and Fernández (2009). While this was tentatively accepted in a recent checklist of European herpetofauna (Speybroeck et al., 2010), it was noted that it is “impossible to really evaluate the degree of concordance between nuclear and mtDNA data”, and further state that it is not clear if morphological differences are maintained near contact zones. Figure 5 (A) Maximum likelihood (ML) tree based on Iberian Blanus ND4 and 16S mitochondrial haplotypes. Numbers on each node represent from top to bottom, ML bootstrap values, Bayesian posterior probabilities, maximum parsimony (MP) and minimum evolution bootstrap values. Nodes with either ML or MP bootstrap values above 70% are shown, otherwise are collapsed. (B) Map of the Iberian Peninsula showing sampling localities of Blanus cinereus. Grey dots represent populations of the southwestern clade, and black dots represent populations of the central clade. Figures adapted from Albert et al. (2007) and Albert and Fernández (2009). 13 4.1.2. THE CHECKERBOARD WORM LIZARDS – TROGONOPHIS Trogonophis wiegmanni Kaup, 1830 (Figure 6) is the only representative of the Trogonophidae family in North Africa. This species is endemic to the Maghreb, ranging from southwest Morocco to northeast Tunisia, within a Mediterranean biome (Bons and Geniez, 1996) (Figure 7). T. wiegmanni is the only species of the genus, having two currently recognized subspecies: T. wiegmanni wiegmanni Kaup, 1830 and T. wiegmanni elegans (Gervais, 1835) (Bons and Geniez, 1996; Schleich et al., 1996). T. w. wiegmanni is distributed in central and eastern Morocco, northern Algeria and western Tunisia (Gans, 2005). This subspecies inhabits relatively dry regions, being found in altitudes up to 1600 meters (m) (Schleich et al., 1996). T. w. elegans is endemic to western Morocco, ranging from the Rif to Souss Valley in southwest Morocco, except in the highest mountains, being seldom found higher than 900 m. It occupies relatively moist regions influenced by the temperate Atlantic climate (Bons and Geniez, 1996; Schleich et al., 1996). The apparently very different ecological demands of both forms suggest a considerable step towards speciation (Schleich et al., 1996). Figure 7 Map with the region where Trogonophis wiegmanni occurs, with its distribution range according to IUCN (2012). Figure 6 Trogonophis wiegmanni specimens. 14 Morphologically, this eyeless worm-like lizard has scales arranged in annuli, a short and conical tail, and a chessboard pattern (Schleich et al., 1996), hence its common name – checkerboard worm lizard. Colouration is different in the two subspecies – T. w. wiegmanni has a ground yellow colour and T. w. elegans has a malve or pink colouration, particularly visible in juveniles (Bons and Geniez, 1996; Schleich et al., 1996). However, the yellow pigmentation of T. w. wiegmanni tends to disappear in preserved specimens, making it is impossible to distinguish morphologically between the two subspecies (Schleich et al., 1996). It has been observed specimens from Algeria (La Chiffaand Biskra) which were fuliginous grey and may represent a third subspecies (Schleich et al., 1996). Mendonça and Harris (2007) studied levels of genetic variation based on mitochondrial sequences of 12S and 16S rRNA (Figure 8). Phylogenetic analyses showed two monophyletic clades corresponding to T. w. elegans and T. w. wiegmanni from Morocco, separated by the Atlas Mountains, with high genetic distance (3.8% uncorrected p-distance for 16S). This result combined with different morphology and the seemingly distinct ecological needs of these two forms indicate that they could possibly be elevated to species status. Additionally, it was included in the analyses a T. w. wiegmanni sample from Tunisia, arising as a separate lineage, with a high level of genetic distance to all Moroccan samples (4.8% uncorrected p-distance for 16S). This Tunisian sample appeared to be more closely related to T. w. elegans (4.4%) than to Moroccan T. w. wiegmanni (4.8%) – however the analyses of this study did not support the inclusion of the Tunisian specimen in either Trogonophis clades. The authors further suggested that this Tunisian form may be a different subspecies or even species, given that some authors already recognize elegans as a full species (e.g. Gans, 2005). This study shows another exceptional example of an amphisbenid species complex in North Africa. Further studies including a larger sampling coverage – especially from Algeria and Tunisia – performing molecular analyses with nuclear markers, as well as assessing morphological variation in this genus, particularly between the Tunisian and Moroccan T. w. wiegmanni forms, are necessary in order to clarify the taxonomic status within this genus. 15 5. BIOGEOGRAPHY OF THE STUDY REGION The Mediterranean Basin (Figure 9) is considered to be one of the world’s biodiversity hotspots with many areas presenting significant high levels of endemism (Myers et al., 2000). This region has suffered a series of events that have helped shaping the existing biodiversity and taxa genetic diversity. The European Southern Peninsulas, Anatolia and the Maghreb are fascinating areas to study phylogeographic patterns due to their complex geological and climatic histories, heterogeneous landscapes, diversity of habitats, well defined barriers and the known age for some geological events. Figure 8 (A) Maximum likelihood tree (ML) based on 12S and 16S rRNA sequences from Mendonça and Harris (2007). All analyses produced identical relationships to the one shown. Near the nodes, Bayesian posterior probabilities and ML bootstrap values. For both analyses, only bootstrap values above 50% are represented. (B) Map showing sampling locations of T. wiegmanni specimens analysed. Figures adapted from Mendonça and Harris (2007). 16 5.1. HISTORICAL EVENTS – GEOLOGICAL AND CLIMATIC CHANGES The Mediterranean Basin between southern Europe and North Africa has had a complex geological history resulting from the movement to the north of the African plate towards western Eurasia. The eastern Mediterranean closed 15-19 MYA in the midMiocene, allowing biotic dispersal between Eurasian and African taxa (Carranza et al., 2004 and references therein). During the Tortonian stage it occurred the uplift of the Atlas Mountains. At that time there was an archipelago not colonisable by land, between the Iberian mainland and Northwest Africa, which nowadays corresponds to the Betic cordillera in southeastern Iberia and the Rif Mountains in northern Morocco (De Jong, 1998). Then, at the end of the Tortonian (7.2 MYA), the Betic strait closed (Duggen et al., 2003), which led to the Tortonian salinity crisis, coinciding with a global aridification of the climate (Krijgsman et al., 2000). The connection between the Atlantic Ocean and the Mediterranean Sea ceased 5.96 MYA, leading to the Messinian salinity crisis, during which parts of the Mediterranean dried. At that time, North Africa and the Iberian Peninsula were connected, with extensive land-bridges allowing the dispersal of fauna all over the basin. This period ended 5.33 MYA with the opening of the Strait of Gibraltar (Duggen et al., 2003; Hsü, 1973; Krijgsman et al., 1999). Even though it was first proposed that the reopening of the Strait of Gibraltar acted as a barrier promoting vicariant events in several organisms, recent phylogeographic studies, mainly on amphibians and reptiles, have advanced that the Strait may have acted as a dispersal channel rather than an impermeable barrier for these taxa (Santos et al., 2012 and references therein). Figure 9 Map of the study region - the Mediterranean Basin. 17 During the Quaternary, which occurred approximately 2.4 MYA, it took place a series of Ice Ages with successive cooling and warming of the Earth’s climate (Hewitt, 2000). These worldwide climatic changes have played a role on organisms’ phylogeographic patterns, being responsible for a major change of organisms’ demographic structure and distribution all over the world and providing opportunities for adaptation to occur (Hewitt, 2000; Hewitt, 2004). The interpretation of historical factors that have led to the present geographic distribution of evolutionary units, at the intra and interspecific level, has allowed the identification of glacial refugia, post-glacial colonization routes and secondary contact zones (Hewitt, 1996, 1999, 2000, 2001; Taberlet et al., 1998). Southern European peninsulas – Iberian, Apennine and Balkan peninsulas – Anatolia and the Maghreb have been considered important refugia areas of genetic diversity during the Pleistocenic climatic oscillations (Hewitt, 1996, 1999; Schmitt, 2007). These regions harbour large numbers of endemic taxa (Hewitt, 1996, 1999, 2000; Taberlet et al., 1998), including high amounts of herpetofauna diversity. In Africa, glaciations corresponded to humid periods (Prentice and Jolly, 2000), allowing species dependent of humid habitats to range widely, while interglacials would have restricted them to the mountains, more similar to the present distribution. In Europe, during glacial cycles, advancing ice and tundra in higher latitudes forced organisms to retreat to southern peninsular refugia, less affected and with a more stable climate. This was followed by expansion of persistent species from these refugia throughout interglacial warming periods to fastly recolonise newly deglaciated areas in northern and central Europe (Hewitt, 1996, 1999, 2000). Intrinsic to this perspective is the fact that southern refugia harbour higher levels of species richness, serving as a resource for later demographic expansions as well as evolutionary radiations (Hewitt, 1996, 2001). Populations in the northern region of the refugia would spread out over long distances, when the climate ameliorated, and colonise large, suitable and available habitats, forming areas of secondary contact after the expansion of divergent lineages. Interglacial periods represented a chance for interaction and introgression between lineages. High levels of genetic diversity have also been found in suture zones outside of refugia, where species or evolutionary lineages originating from different refugia met after the last glaciations (Petit et al., 2003; Taberlet et al., 1998). 5.2. NORTH AFRICA The Maghreb region is constituted by Northwest African countries, including Morocco, Algeria and Tunisia (Bons and Geniez, 1996). This region is rich in herpetofauna diversity due to a complexity of geological and climatic factors (Schleich et al., 1996) – it has 18 Mediterranean climate, landscape and vegetation, and fauna is constituted by a mixture of African, Saharian and Mediterranean elements (Bons and Geniez, 1996). During the Miocene (23-5 MYA), geological changes have produced a massive impact on North Africa’s floral and faunal diversity, such as the uplift of the Atlas Mountains in the midHolocene. Most of the recent studies regarding Maghrebian flora and fauna revealed unexpected patterns of genetic diversity. Morocco is an interesting area to study genetic diversity in reptiles because it is one of the most species-rich countries in North Africa (Bons and Geniez, 1996), having the highest percentage of endemic reptile species in the Mediterranean region (Pleguezuelos et al., 2010). This richness is due to several characteristics of the country. First of all, Morocco has a big area – 458 730 square kilometres (km²). It is delimited to the west by the Atlantic Ocean, to the north by the Mediterranean Sea, to the south by the Sahara desert, an arid zone and to the east by the Moulouya river basin (Bons and Geniez, 1996). In this country there are several geographical units, such as the Rif and Atlas Mountain systems, and the Moulouya river basin. The Rif Mountains, located in the north of the country, are oriented northwest to the northeast. The Atlas Mountains in Morocco, with three sub-systems – Anti-Atlas, High Atlas and Middle Atlas – are located in the centre and south of the country, and are oriented from the northeast to southwest. The Atlas mountainous barrier divides Morocco into two different bioclimatic regions: north and east Morocco have Mediterranean climate and South of the Atlas have Saharan climate (Bons and Geniez, 1996). Finally, the Moulouya River basin, located in eastern Morocco, may also act as a barrier to a great number of Moroccan and Algerian taxa (Bons and Geniez, 1996). Phylogeographic studies show patterns that reflect the influence that geological factors may have had on the biogeography of the region. For example, vicariance phenomena associated to the uplift of the Atlas Mountains and the existence of the Moulouya River have been proposed to be the causes of the phylogeographical patterns observed in Agama impalearis (Brown et al., 2002) and Testudo graeca (Álvarez et al., 2000), respectively. 5.3. IBERIAN PENINSULA The Iberian Peninsula is a species rich area, which may have worked as one of the foremost Pleistocene glacial refugia in Southern Europe (Hewitt, 1999, 2001). Its high habitat diversity and complex geological history may have made this region an ideal survival refuge during the Pleistocene. The Iberian Peninsula is formed by several mountain ranges and river systems, most of them being east-west orientated. These act as apparent barriers to gene flow to north- 19 south dispersal for many species and at the same time provide a high variety of microclimates, allowing populations to move in altitude in response to the climate fluctuations (Hewitt, 1996). Because of its location, the Iberian Peninsula is under the influence of both the North Atlantic and Mediterranean, having a wide variety of climates – desert, Mediterranean, Alpine and Atlantic. At the same time, it has a large area of 580 000 km2, which makes it improbable that it may have worked as a single homogenous and continual refugial area during the ice ages (Gómez and Lunt, 2007). The wide variety of climate and geological structures might have created a differential distribution and fragmented suitable habitats across Iberia that may have favoured the existence of multiple glacial refugia distant from each where species persisted during Pleistocene climatic cycles (Gómez and Lunt, 2007; Martínez-Solano et al., 2006). Two areas of deciduous forest became acknowledged as the Lusitanian and Andalucian refugia and have been considered as refugia for several species in the Iberian Peninsula (Paulo et al., 2001). Phylogeographic studies in the Iberian Peninsula showed extremely divergent lineages with strong correspondence with geography, due to deep population fragmentation related with isolation in glacial refugia. This observation led to a ‘refugia-within-refugia’ paradigm in this peninsula postulated by Gómez and Lunt (2007), suggesting that most species would have persisted across the Pleistocene Ice Ages in different and isolated locations in Iberia. Southwestern and southeastern Iberia revealed to be important refugial areas. Several Iberian endemic species show undeniable phylogeographic concordance with an alleged refugium located in or near the southern Betic Ranges (Gómez and Lunt, 2007). This is supported, for instance, by the high diversity found within Mauremys leprosa (Fritz et al., 2006), Lacerta schreiberi (Godinho et al., 2008; Paulo et al., 2001), Alytes cisternasii (Gonçalves et al., 2009) and Lissotriton boscai (Martínez-Solano et al., 2006) in southern Portugal populations. 5.4. ANATOLIA Anatolia, also known as Asia Minor, is a region in western Asia, with an area of 755 688 km2 and about 7 000 km2 additional coastal islands, including territories in Turkey, Greece, Asia Aegean and Mediterranean coast. It is limited by the Aegean, the Mediterranean and the Black Sea to the west, south and north, respectively, while to the northeast and the east by the Caucasus and the Armenian highlands (González, 2012). This region has a complex climatic and geological history which has shaped endemic taxa biogeographic patterns, being characterized by a rich species and taxonomic diversity of reptiles (Sindaco et al., 2000). During Pliocene and Pleistocene climatic oscillations, this region served as a major refuge, which led to successive vicariance and dispersal events in 20 interglacial periods and subsequent extinctions. Relict refugial populations established centres of endemism (Veith et al., 2003; Wilke et al., 2007). Due to its position and long palaeogeographic and palaeoclimatic history, in the past Anatolia acted as a bridge and as a barrier for species dispersal between Asia and Europe (Kornilios et al., 2011; Sindaco et al., 2000). It is situated in a mixture of European, Asian and African biomes, from where it experienced repeated invasions since the Late Oligocene and several present-day Anatolian endemics are relics of such invasion processes (Veith et al., 2003). It is mostly a mountainous region, with true low lands confined to coastal fringes. Four main relief regions can be identified: the northern and southern mountains, the central massif and the Arabian platform (Sindaco et al., 2000). The relief also affects the climate in this region, which is often much harder than might be expected for that latitude, producing different climatic regions (Sindaco et al., 2000). Therefore, Anatolian mountains may have played an important role in the speciation and definition of biogeographical subregions and have been defined as biodiversity hotspots (Kornilios et al., 2011). 6. OBJECTIVES The high levels of mtDNA divergence within Iberian Blanus and North African Trogonophis lineages and the apparent occurrence of cryptic species (Albert and Fernández, 2009; Albert et al., 2007; Mendonça and Harris, 2007; Vasconcelos et al., 2006) may suggest similarly complex patterns in other worm lizards taxa. Particularly, North African and eastern Mediterranean Blanus await further phylogenetic studies to verify the occurrence of potential cryptic diversity in those taxa. Additionally, mitochondrial phylogenetic inferences ought to be further reassessed with nuclear markers, to corroborate previous results. The aim of this study was to estimate levels of genetic diversity and phylogenetic relationships within Blanus and Trogonophis wiegmanni, by analysing multiple mitochondrial and nuclear sequence data. In Blanus, by including nuclear markers and additional specimens for each one the five Blanus species it was further tested: i) whether nuclear DNA supports the phylogenetic relationships of the genus and the occurrence of cryptic species within the Iberian Blanus suggested by previous studies mostly based on mtDNA (Albert and Fernández, 2009; Albert et al., 2007; Vasconcelos et al., 2006); ii) whether the northern and southern clades of B. tingitanus, identified by Vasconcelos et al. (2006) based on mtDNA, constitute two allopatric clades both in mitochondrial and nuclear genealogies, thus suggesting that they may represent independent taxa; and iii) whether also B. mettetali and B. strauchi show high intraspecific differentiation and cryptic diversity. 21 Regarding Trogonophis, by adding more samples, particularly from the poorly studied eastern distribution of Trogonophis, and including nuclear sequences, it was intended to: i) corroborate the mitochondrial phylogenetic relationships within T. wiegmanni forms, inferred in a previous study (Mendonça and Harris, 2007); ii) analyse the differentiation among the Moroccan T. wiegmanni subspecies iii) test the variation between the Tunisian T. w. wiegmanni form and the remaining Moroccan Trogonophis and its inclusion in either clade by adding samples from Algeria and Tunisia. Additionally in this study, it was also briefly discussed the utility of the nuclear markers employed to uncover the levels of differentiation with are observed at the mitochondrial level. 28 2. MOLECULAR MARKERS SELECTION For this study it was assembled a multilocus dataset, with mitochondrial and nuclear gene fragments. Mitochondrial markers included 12S ribosomal RNA (12S), 16S ribosomal RNA (16S), cytochrome b (CYTB) and NADH dehydrogenase subunit 4 (ND4). Nuclear markers included melanocortin-1 receptor (MC1R), proopiomelanocortin (POMC), oocyte maturation factor (CMOS) and recombination activation gene 2 (RAG2). Mitochondrial markers such as 12S and 16S rRNA and protein coding genes CYTB and ND4 are fast evolving genes that have been broadly used to infer phylogenetic relationships in several herpetofauna taxa (Perera and Harris, 2010; Salvi et al., 2010) as well in amphisbaenians (Mott and Vieites, 2009; Mulvaney et al., 2005; Vidal et al., 2008). Besides their characteristics – universal primers, useful and informative for an initial phylogenetic assessment, hence being widely used – most mitochondrial genes were chosen for this study because they had already been used for Blanus and Trogonophis in previous studies (Albert et al., 2007; Mendonça and Harris, 2007; Vasconcelos et al., 2006), with sequences available on GenBank. Nuclear markers were selected after trials with an extended set of primers. Selected nuclear markers are protein-coding single-copy genes, which vary in degree of conservation. CMOS is a slow-evolving proto-oncogene that encodes a protein that regulates meiotic maturation (Saint et al., 1998). It is a useful maker to test relationships within and among squamate families (Harris et al., 1999; Saint et al., 1998). It has been found to be appropriate at both deep and shallow divergences (Saint et al. 1998), being informative among taxa that diverged up to 400 MYA (Graybeal, 1994). MC1R is a critical regulator of melanin synthesis. When developed by Pinho et al. (2010) it showed high levels of polymorphism, being useful in population genetics and phylogenetic analyses for a variety of taxa. Even though this marker has not been used in amphisbaenians phylogenetic assessments, it has been used to infer other squamata phylogenies (Barata et al., 2012; Gonçalves et al., 2012). POMC is a polypeptide hormone precursor only found in vertebrates, which undergoes post-translational modification to produce multiple hormones. It has both conserved and variable regions, suggesting it may be a useful phylogenetic marker (reviewed in Becker et al., 2011). Even though it has not been used to estimate amphisbaenians phylogenies, it has been selected based on its proven phylogenetic utility for other groups of squamate reptiles (Vieites et al., 2011). RAG2 is a slow-evolving and highly conserved nuclear marker in vertebrates. It encodes components of the recombinase involved in recombination of immunoglobin and T-cell receptor genes (reviewed in Lovejoy et al., 2001). It has been used both in the inference of amphisbaenian 29 relationships (Vidal et al., 2008), as well of other squamate reptiles (Crottini et al., 2009; Rato et al., 2010). 3. LABORATORY PROCEDURES – DNA EXTRACTION, AMPLIFICATION AND SEQUENCING Genomic DNA was extracted from tails preserved in ethanol using a standard saline method (Sambrook et al., 1989). Small amounts of tissue were digested in extraction buffer with proteinase K solution, occurring lysis of cell membranes. Then, ammonium acetate was added to precipitate proteins, which were removed. To the supernatant, icecold isopropanol was added to precipitate DNA, forming a pellet after centrifugation. DNA pellets were washed with ethanol and then air dried and hydrated in 70µL of ultra-pure water. The obtained DNA was then ready to be used as template in Polymerase Chain Reaction (PCR). However, if the amount of tissue were too small, DNA was extracted using a commercial kit, QIAmp® DNeasy Blood & Tissue kit (Quiagen, Valencia, California), following the manufacturer’s instructions. The advantage of this method is that it allows the rapid isolation of highly pure genomic DNA from a small piece of tissue. Nevertheless, it is more expensive than the standard saline method often used in DNA extraction. Two mitochondrial markers – 16S and ND4 were amplified for all 49 Blanus samples. Additionally, two nuclear markers – POMC and MC1R – were amplified for 26 and 21 samples, respectively. After genetically analysing mitochondrial fragments, samples were selected from each lineage to be amplified for the nuclear fragments POMC and MC1R. The primers used for MC1R failed to amplify that fragment for B. strauchi and outgroups. For Trogonophis, three mitochondrial – 12S, 16S and CYTB – and three nuclear markers – CMOS, POMC and RAG2 were amplified. Twenty eight T. wiegmanni samples were amplified only for 12S and 16S and 22 samples were amplified for all genes. For both Blanus and Trogonophis samples, amplifications were performed in final volumes of 25µL with a thermocycler, containing 5µL 5x reaction buffer, 2-3mM MgCl2, 0.2-0.4µM each dNTP, 0.2µM each primer, 2U of Taq polymerase and 0.5-1µL DNA template. When necessary, annealing temperatures and/or the amount of magnesium ions were adjusted to increase amplification yield and specificity on a case by case basis. For further detail on primer names, sequences and references, and specific PCR conditions temperatures for each marker consult Table 3. Amplified products were visualized with 2% agarose gel electrophoresis, in order to confirm if the PCR reactions were successful. Amplified products were sequenced by a commercial sequencing facility (Macrogen Inc.) using the same primers used for amplification. 30 Table 3 Gene, primer name, sequence, source and PCR conditions (temperature, time and number of cycles). Amplification of RAG2 fragments was done with a nested PCR with two sets of primers - first PCR using external primers 31FN.Venk and Lung.460R and a second PCR with internal primers Lung.35F and Lung.320R. The latter ones were used for sequencing. Gene Primer Sequence (5’ → 3’) Reference PCR condition 12S 12Sa AAACTGGGATTAGATACCCCACTAT Kocher et al. (1989) 92ºC(2m), [30x 92ºC(30s), 48ºC(40s), 72ºC(45s)], 72ºC(5m) 12Sb GAGGGTGACGGGCGGTGTGT Kocher et al. (1989) 16S 16SL CGCCTGTTTATCAAAAACAT Palumbi et al. (1996) 92ºC(2m), [30x 92ºC(30s), 48ºC(40s), 72ºC(45s)], 72ºC(5m) 16SH CCGGTCTGAACTCAGATCACGT Palumbi et al. (1996) ND4 ND4 CACCTATGACTACCAAAAGCTCATGTAGAAGC Arévalo et al. (1994) 94ºC(3m), [40x 94ºC(30s), 50ºC(30s), 72ºC(45s)], 72ºC(4m) LEU CATTACTTTTACTTGGATTTGCACCA Arévalo et al. (1994) CYTB CYTB1 CCATCCAACATCTCAGCATGATGAAA Kocher et al. (1989) 94ºC(3m), [30x 94ºC(30s), 48ºC(30s), 72ºC(1s)], 72ºC(5m) CYTB2 CCCTCAGAATGATATTTGTCCTCA Kocher et al. (1989) CMOS CMOS G73 GCGGTAAAGCAGGTGAAGAAA Saint et al. (1998) 94ºC(3m), [35x 94ºC(45s), 48ºC(45s), 72ºC(1m30s)], 72ºC(5m) G74 TGAGCATCCAAAGTCTCCAATC Saint et al. (1998) RAG2 31FN.Venk TTYGGICARAARGGITGGCC Venkatesh et al. (2001) 94ºC(5m), [35x 94ºC(30s), 50ºC(50s), 68ºC(1m30s)], 68ºC(5m) Lung.460R GCATYGRGCATGGACCCARTGCC Brinkmann et al. (2004) Lung.35F GGCCAAAGAGRTCYTGTCCIACTGG Hoegg et al. (2004) 94ºC(5m), [35x 94ºC(30s), 50ºC(50s), 68ºC(1m30s)], 68ºC(5m) Lung.320R AYCACCCATATYRCTACCAAACC Hoegg et al. (2004) MC1R MC1RF GGCNGCCATYGTCAAGAACCGGAACC Pinho et al. (2010) 94ºC(3m), [30x 94ºC(30s), 50ºC(30s), 72ºC(1m)], 72ºC(5m) MC1RR CTCCGRAAGGCRTAAATGATGGGGTCCAC Pinho et al. (2010) POMC POMCF ATATGTCATGASCCAYTTYCGCTGGAA Vieites et al. (2007) 94ºC(3m), [30x 94ºC(30s), 50ºC(45s), 72ºC(1m)], 72ºC(5m) POMCR GGCRTTYTTGAAWAGAGTCATTAGWGG Vieites et al. (2007) 30 31 4. PHYLOGENETIC INFERENCE ANALYSES 4.1. DATA ANALYSES DNA sequences’ chromatographs were checked and sequences were edited in GENEIOUS v5.3.6 (Drummond et al., 2010). For nuclear sequences, nucleotide ambiguities with similar peak size in chromatograms were considered heterozygous positions. Previously published sequences of 16S and ND4 for Blanus (Albert et al., 2007) and 12S and 16S for Trogonophis (Mendonça and Harris, 2007) were added to the analyses. DNA sequences for each gene independently were aligned using MAFFT v6.814b (Katoh et al., 2002) with default parameters (gap open penalty=1.53, gap extension=0.0). In case sequences were shorter than the rest of the alignment, the initial or end gaps were substituted by “N”, meaning that there is an equal probability for any nucleotide to be present in those positions. Before carrying out phylogenetic analyses, 12S and 16S alignments were analysed with GBLOCKS (online version 0.91b; Castresana, 2000) using a less stringent selection to remove regions that could not be unambiguously aligned. Summary statistics for all markers were calculated in DNASP v5 (Librado and Rozas, 2009). For Blanus, networks were produced using 16S and ND4 mitochondrial sequences concatenated, including some sequences, mainly of Iberian species, from a previous study (Albert et al., 2007) to better represent the haplotype diversity and distribution coverage of these species. Representatives from the main lineages were selected and sequenced for nuclear genes MC1R and POMC, and further ML and BI analyses were conducted with mitochondrial and nuclear alignments combined. Then, nuclear sequences were used to produce individual gene networks. For Trogonophis, not all previously published samples were available to amplify for new markers. As a result, it was produced a mitochondrial network using 12S and 16S sequences combined, including sequences from Mendonça and Harris (2007), representing a larger distribution coverage dataset. Then, samples available for this study were sequenced for CYTB, CMOS, POMC and RAG2 and combined for further ML and BI inference analyses. Also, individual nuclear alignments were used to produce networks. 32 4.2. PHYLOGENETIC ANALYSES Phylogenetic analyses were performed using maximum likelihood (ML) and Bayesian Inference (BI) methods to produce gene trees. Nuclear alignments had unphased sequences, with heterozygous positions coded as ambiguities (IUPAC codes). First, analyses on each gene were done independently to identify potential incongruence between partitions (results not presented) and then analyses were performed on concatenated datasets. Different sets of outgroups were tested – for Blanus phylogenetic analyses were conducted with D. zarudnyi and/or T. wiegmanni; for Trogonophis were tested D. zarudnyi and/or B. mettetali. For each, it was chosen and it is further presented the set of outgroups which provided a better resolution of the phylogenetic relationships with the best support values. For both Blanus and Trogonophis, BI and ML phylogenies were conducted with mitochondrial and nuclear alignments concatenated. Analyses were carried out as partitioned analyses of molecular data, since data from different DNA regions were combined. Gene by gene partitions were used for all concatenated analyses. JMODELTEST v0.1 (Posada, 2008) was used to select the best fitting models of nucleotide substitution for each gene for BI analyses, based on likelihood scores for 88 different models under the Akaike Information Criterion corrected for small sample sizes (AICc). ML analyses were performed using a graphical user interface (GUI) for RAXML (Stamatakis, 2006) – RAXML GUI v1.2 (Silvestro and Michalak, 2010) – in individual sequences for the partitioned concatenated dataset under the GTR+G+I model and perpartition branch lengths. It was carried out a ML search and thorough bootstrapping, with 1000 replications to evaluate the stability of nodes of the phylogenetic tree (BP) (Felsenstein, 1985). All partitioned Bayesian analyses were performed with MRBAYES v3.1.2 (Huelsenbeck and Ronquist, 2001), using the selected model of sequence evolution and model parameters for each partition, in individual sequences. Bayesian posterior probability (BPP) values were estimated using a Metropolis-coupled, Markov chain Monte Carlo (MCMCMC) sampling approach. Bayesian analyses started with randomly generated trees and ran for 2x106 generations, using four incrementally heated Markov chains with default heating values. Markov chains were sampled at intervals of 100 generations, producing 20 000 trees. All analyses ended with the standard deviation of split frequencies less than 0.01. Stabilization and convergence between runs were assessed in TRACER v1.5 (Rambaut and Drummond, 2009), in terms of likelihood scores and parameters. The log- 33 likelihood values of the 20 000 trees in each analysis were plotted against the generation time, using the “sump” command generated in MRBAYES. Burn-in data sampled from generations preceding the stationarity of the Markov Chain were discarded. Runs became stationary after 100 000 – 200 000 generations, and the corresponding first trees were discarded as burn-in to assess posterior probabilities for nodal support (BPP). Remaining trees were combined and 50% majority-rule consensus trees were generated. Two independent replicates were carried out to check that analyses were not trapped at local optima (Huelsenbeck and Bollback, 2001a). Nodes were considered strongly supported if they received BP≥70% and BBP≥0.95. Mitochondrial sequence variation was analysed by producing haplotype networks. For Blanus and Trogonophis combined mitochondrial alignments (16S and ND4, and 12S and 16S, respectively), statistical parsimony haplotype networks were carried out in TCS v1.2.1 (Clement et al., 2000), under the 95% probability criterion. Net distances based on a p-distance method (the proportion of nucleotide sites at which two sequences being compared are different) were calculated between mtDNA clades sequences, and uncorrected p-distances were calculated within mtDNA clades haplotypes in MEGA v5 (Tamura et al., 2011). Haplotype diversity and structure for each nuclear marker was represented by haplotype networks to compare the phylogenetic signal with the phylogenies produced. Nuclear markers sequences were computationally phased using a coalescent-based Bayesian method in PHASE v2.1.1 (Stephens and Donnelly, 2003; Stephens et al., 2001) using default parameters (thresholds: p=q=90%), as implemented in DNASP v5 (Librado and Rozas, 2009). Three runs were carried for each dataset to check for consistency of results. Output files were used to construct haplotype networks under the statistical parsimony approach implemented in TCS v1.2.1 (Clement et al., 2000), under the 95% probability criterion. 34 35 III. RESULTS 36 37 1. BLANUS Figure 10 Maps of the study areas with the geographic location of the Blanus samples analysed. Samples are represented by circles coloured according to phylogenetic analyses results. Samples with (*) after code are from a published study (Albert et al., 2007). (Samples B33, B34 and Bs* are not displayed on the map since their exact locations are unknown.) Dashed areas in the central map represent the distribution range of each species according to IUCN (2012). 44 12S) (Table 5). The separation of these lineages has geographic concordance, with T. w. elegans occurring in western Morocco and T. w. wiegmanni occurring in the eastern region of the country (Figure 14). The basal lineage corresponding to the eastern distribution of T. w. wiegmanni in Algeria and Tunisia is well supported (BPP=1.00/BP=96). Figure 15 Phylogenetic tree representing the relationships within Trogonophis, derived from the Bayesian analysis of the combined markers [12S (353 bp), 16S (465 bp), CYTB (284 bp), CMOS (350 bp), POMC (467 bp) and RAG2 (798 bp)] for a total fragment of 2717 bp. The topology is similar in ML analysis (data not shown). For the major clades and lineages, Bayesian posterior probabilities (BPP) and ML bootstrap support (BP) are given above and below nodes, respectively. The tree was rooted using Diplometopon zarudnyi and Blanus mettetali. 45 Table 5 Net distance between Trogonophis clades and uncorrected p-distances within Trogonophis clades. Lower left values correspond to 12S and top right to 16S values between clades sequences. Shaded diagonal values correspond to uncorrected p-distance values within clades haplotypes. 16S 12S T. w. elegans T. w. wiegmanni Morocco T. w. wiegmanni Algeria-Tunisia T. w. elegans 0.015 0.010 0.025 0.030 T. w. wiegmanni Morocco 0.031 0.013 0.015 0.028 T. w. wiegmanni Algeria-Tunisia 0.031 0.042 0.022 0.006 The mitochondrial statistical parsimony network resulted in 22 haplotypes, with a connection limit of 12 mutational steps (Figure 16). The networks supported the combined mitochondrial and nuclear phylogenetic relationships. They also showed that T. w. elegans has high levels of intraspecific diversity. Within T. w. wiegmanni in Morocco, the mitochondrial analysis – like the combined phylogeny – showed sub-structuring into two unconnected networks, revealing high distances between them. Among T. w. wiegmanni Figure 16 Haplotype networks recovered by statistical parsimony analysis representing relationships among Trogonophis wiegmanni, inferred from 828 bp mitochondrial sequences (12S and 16S – 361 bp and 472 bp, respectively). Circles represent different haplotypes with size proportional to sample frequency. Small grey circles represent missing or extinct haplotypes. 46 eastern samples from Algeria and Tunisia, there is a relatively high level of genetic divergence, with 2.2% uncorrected p-distance for 16S, among the three individuals analysed (Table 5). Both the combined phylogeny (Figure 15) and the mitochondrial haplotype networks (Figure 16) showed divergence between the Algerian and the Tunisian samples. Statistical parsimony analyses (Figure 17) revealed that for CMOS and RAG2 nuclear haplotype networks, the three different clades found in the phylogenetic tree analyses (Figure 15) shared haplotypes. In the POMC network, a clear structure separating the three lineages is evident, with the exception of the T. w. elegans sample T7, which shares the same haplotype with T. w. wiegmanni samples from Morocco. 47 Figure 17 Haplotype networks recovered by statistical parsimony analysis representing relationships among T. wiegmanni main clades, inferred from nuclear markers CMOS (350 bp), POMC (467 bp) and RAG2 (798 bp). Circles represent different haplotypes with size proportional to sample frequency. Small grey circles represent missing or extinct haplotypes. Heterozygotes were included as independent samples and are represented by ‘i’ or ‘ii’ after sample code. 48 49 IV. DISCUSSION 50 51 Cryptic genetic diversity is evident in the amphisbaenids occurring in the Mediterranean region, Blanus and Trogonophis. In these taxa, genetic analyses using mitochondrial and nuclear gene fragments revealed high levels of differentiation among apparently morphologically similar forms. Even though previous studies had already analysed genetic variation and phylogenetic relationships in these forms, they were based mainly on mitochondrial markers. New nuclear data revealed similar patterns to the ones previously recovered, with some markers found to be more informative than others. Moreover, the addition of samples from new localities revealed possible contact zones areas, previously unknown and that might be worthy to further investigate in the future. 1. BLANUS The phylogenetic analyses of both mitochondrial and nuclear gene fragments support the monophyly of all currently described Blanus species. Furthermore, the analyses show the existence of three main clades grouping the Iberian, North African and Anatolian species, with the latter being basal to the other species of this genus. The monophyly of the North African clade had been questioned by Albert et al. (2007), stating the possibility of ancient paraphyletic lineages but, with the addition of more sampled Moroccan localities, this study found it to be a monophyletic clade. The phylogenetic relationships between species recovered in this study corroborate molecular findings by Vasconcelos et al. (2006) and Albert et al. (2007). The present study provides a more complete phylogenetic inference in terms of sampling and also number of analysed markers, thus providing new insights into the genetic diversity and differentiation among Blanus, which are discussed in the next sections. 1.1. IBERIAN BLANUS CRYPTIC SPECIES The combined mitochondrial and nuclear phylogeny recovered two distinct Blanus clades in the Iberian Peninsula, which is concordant with the previous results by Vasconcelos et al. (2006) and Albert et al. (2007) based mainly on mitochondrial data. Nevertheless, the description of B. mariae by Albert and Fernández (2009) made it difficult to make direct comparisons of patterns between mitochondrial and nuclear data, due to the use of an anonymous nuclear marker by Albert et al. (2007) (Speybroeck et al., 2010). Networks results from two nuclear markers (Figure 13) support the genetic distinction between B. mariae and B. cinereus, with no haplotype sharing between the two Iberian species. Moreover, this study results indicate that the distribution of B. mariae is wider than previously described by Albert and Fernández (2009). The authors suggested that 52 the limit to the distribution in western Iberian Peninsula was in the lower third of Portugal, from the Algarve up to Elvas. However, the new data analysed in this study extended the distribution range of the recently described B. mariae to central Portugal, in Carvalhão (samples B12-13) (Figure 10). This clearly indicates the need for further sampling and field observations, particularly in centre and north Portugal, in order to establish an accurate distribution of this species, of great relevance to establish adequate conservation policies. Also, the distribution of the two Iberian species may partially overlap in central Portugal – B. cinereus B19 from São Mamede and B. mariae B12-13 from Carvalhão (Figure 10). Detailed analyses in contact areas, using both mitochondrial and nuclear markers and morphology would be of great to evaluate the occurrence of introgression phenomena between the two forms, and understand the evolutionary history of the Iberian Blanus. 1.2. INTRASPECIFIC DIVERSITY WITHIN MOROCCAN B. TINGITANUS This study indicates that within B. tingitanus exists considerable genetic differentiation between northern and southern lineages, with 2.6% distance (net distance for ND4), although this is less than between currently accepted species. The finding of these two lineages with a combined mitochondrial and nuclear phylogeny corroborates results by Vasconcelos et al. (2006). Comparatively, this new study analyses a larger sampling coverage in Morocco, particularly for the western distribution of the species, providing a more accurate distribution of these two lineages. In Kenitra, a locality in northwestern Morocco, both lineages are in simpatry (samples B1 and B7) (Figure 10 and Figure 11) and nuclear networks reveal that these samples share haplotypes with southern samples (Figure 13). Also, in this area B. mettetali has also been reported (Bons and Geniez, 1996). Therefore, this ought to be an interesting region to be further investigated. The distribution of the B. tingitanus lineages has geographic concordance, with the northern lineages including specimens restricted to the Rif region, and the southern lineages with samples restricted to south of the Rif and north of the Middle Atlas Mountains. The phylogeographic break observed in B. tingitanus may indicate that the Moroccan mountain systems, such as the Rif or the Atlas Mountains may have played a role in shaping the genetic diversity in this species. Indeed, the Atlas Mountains have been proposed as the cause for phylogeographical division for other reptile species, such as Agama impalearis (Brown et al., 2002) and Mauremys leprosa (Fritz et al., 2005). 53 1.3. PRELIMINARY ASSESSMENT OF GENETIC DIVERSITY WITHIN B. METTETALI AND B. STRAUCHI Intraspecific genetic variation is high within B. mettetali, as revealed by high levels of mitochondrial divergence (4.9% uncorrected p-distance for ND4), despite only five samples were available to study across a wide range (Figure 10 and Figure 12). Within B. strauchi, variation was exceptionally high between the sample from Turkey and the remaining samples from Greece, as revealed by the mitochondrial network (Figure 12) and a 13% net distance for ND4, a level of divergence similar or even higher to that observed between accepted species (e.g. 10% between B. cinereus and B. mariae) – but again few samples were available. Nevertheless, it seems possible that undescribed Blanus species may occur in the Anatolian region. This region is a crossroad between Palearctic, Oriental and Afrotropic ecozones, and was a climatic refugia during the quaternary climatic fluctuations (Hewitt, 2001; Kornilios et al., 2011; Sindaco et al., 2000). For these reasons, this still under-studied region is considered a hotspot of biodiversity. In effect, several recent studies confirm high levels of diversity in plants (Ansell et al., 2011). turtles (Fritz et al., 2009) or mammals (Gündüz et al., 2007). A recent genetic study on another Anatolian reptile, the burrowing snakes Typhlops vermicularis, also revealed high variation (up to 8.4% for ND2) (Kornilios et al., 2011). 2. TROGONOPHIS As first suggested by Mendonça and Harris (2007) and confirmed in the present study, the phylogeny of T. wiegmanni is composed of three monophyletic clades in North Africa. This study inferences based on mitochondrial and nuclear sequences found two monophyletic clades in Morocco, correspondent to T. w. elegans in western Morocco and T. w. wiegmanni in eastern Morocco, and a third clade clustering T. w. wiegmanni samples from Algeria and Tunisia, forming a basal clade (Figure 14). This suggests that T. w. wiegmanni is a paraphyletic subspecies. 2.1. MOROCCAN T. WIEGMANNI FORMS Within the subspecies T. w. wiegmanni in Morroco, the mitochondrial analysis showed a sub-structuring into two lineages (Figure 16Figure 14). One of them seemed to be restricted to areas with lower altitude while the other one was found in mountainous areas (Figure 14). This suggests that within this subspecies there are two forms which may have different ecological requirements, but this will require further assessment. 60 includes genetic diversity, and with the assistance of molecular tools it is possible to determine which areas deserve higher conservation priorities. Hopefully this study will be useful as a further stepping stone in the understanding of cryptic diversity in amphisbaenid taxa. 2. FUTURE WORK Regarding Blanus and T. wiegmanni worm lizards, further sampling is necessary throughout the taxa distribution in order to obtain an accurate distribution of the different lineages and species. At this respect, sampling should be more intense on B. mettetali in southern Morocco, B. strauchi in Anatolia, and T. w. wiegmanni in Algeria and Tunisia. It would be particularly interesting to conduct more sampling around possible contact zones between the different lineages found in this study, in order to better understand the evolutionary processes occurring in areas of sympatry (if such areas exist). The addition of fieldwork on the presumable contact zone between B. mariae and B. cinereus in central Portugal, between B. tingitanus lineages and B. mettetali near Kenitra region in Morocco, and in central Morocco where both T. w. elegans and T. w. wiegmanni occur will be essential to understand the processes involved in the isolation and maintenance of the lineages in contact. This will be valuable not just for assessing the status of these taxa, but to draw comparisons with other similar species in this region. It is also planned to further investigate the correspondence between genetic and morphological variation in Trogonophis, particularly between the Tunisian and Moroccan T. w. wiegmanni forms, in order to clarify their taxonomic status. To do so, we have already contacted several museums to lend specimens to be morphologically analysed. To better understand phylogeographic patterns in these amphisbaenids, it can also be useful to estimate divergence times between lineages. This can be done by using mitochondrial or nuclear genes and several calibration points. Calibration can be obtained from fossil record – which is scarce for Blanus and Trogonophis and therefore may lead to unreliable results – or by biogeographical events, such as the formation of mountain ranges or volcanic islands (Weir and Schluter, 2008). Finally, the two genus of worm lizards studied here represent a challenge for delimiting species due to their cryptic nature. The actual lack of diagnosable morphological characters makes species diagnosis difficult. The use of new approaches in taxonomy, such as an “integrative taxonomy” (Padial et al., 2010) might be useful in this case. Along with genetic studies using a multilocus approach, and a morphological reassessment, it would be also valuable to implement an ecological niche modelling approach to assess which climatic and physical factors, such as temperature, precipitation 61 or altitude, affect the distribution of genetic lineages driving to lineages divergence and also aid in species delimitation (Rissler and Apodaca, 2007). All these evidences could be integrated for taxonomic purposes. 62 63 VI. REFERENCES 64 65 Albert, E.M., Fernández, A., 2009. Evidence of cryptic speciation in a fossorial reptile: description of a new species of Blanus (Squamata: Amphisbaenia: Blanidae) from the Iberian Peninsula. Zootaxa 2234, 56-68. Albert, E.M., Zardoya, R., García-París, M., 2007. Phylogeographical and speciation patterns in subterranean worm lizards of the genus Blanus (Amphisbaenia: Blanidae). Molecular Ecolology 16, 1519-1531. Alexander, A.A., 1966. Taxonomy and variation of Blanus strauchi (Amphisbaenia, Reptilia). Copeia 1966, 205-224. Álvarez, Y., Mateo, J.A., Andreu, A.C., Díaz-Paniagua, C., Diez, A., Bautista, J.M., 2000. Mitochondrial DNA haplotyping of Testudo graeca on both continental sides of the Straits of Gibraltar. Journal of Heredity 91, 39. Ansell, S.W., Stenøien, H.K., Grundmann, M., Russell, S.J., Koch, M.A., Schneider, H., Vogel, J.C., 2011. The importance of Anatolian mountains as the cradle of global diversity in Arabis alpina, a key arctic–alpine species. Annals of botany 108, 241-252. Arévalo, E., Davis, S.K., Sites, J.W., 1994. Mitochondrial DNA sequence divergence and phylogenetic relationships among eight chromosome races of the Sceloporus grammicus complex (Phrynosomatidae) in central Mexico. Syst. Biol. 43, 387. Avise, J., 2004. Molecular markers, natural history, and evolution Sinauer. Sunderland, MA. Avise, J.C., 2000. Phylogeography: the history and formation of species. Harvard University Press. Avise, J.C., Arnold, J., Ball, R.M., Bermingham, E., Lamb, T., Neigel, J.E., Reeb, C.A., Saunders, N.C., 1987. Intraspecific phylogeography: the mitochondrial DNA bridge between population genetics and systematics. Annual Review of Ecology and Systematics 18, 489522. Ballard, J.W.O., Chernoff, B., James, A.C., 2002. Divergence of mitochondrial DNA is not corroborated by nuclear DNA, morphology, or behavior in Drosophila simulans. Evolution 56, 527-545. Ballard, J.W.O., Whitlock, M.C., 2004. The incomplete natural history of mitochondria. Molecular Ecology 13, 729-744. Bandelt, H.J., Forster, P., Röhl, A., 1999. Median-joining networks for inferring intraspecific phylogenies. Molecular Biology and Evolution 16, 37-48. Barata, M., Perera, A., Martínez-Freiría, F., Harris, D.J., 2012. Cryptic diversity within the Moroccan endemic day geckos Quedenfeldtia (Squamata: Gekkonidae): a multidisciplinary approach using genetic, morphological and ecological data. Biological Journal of the Linnean Society. Bazin, E., Glémin, S., Galtier, N., 2006. Population size does not influence mitochondrial genetic diversity in animals. Science 312, 570-572. Becker, R.E., Valverde, R.A., Crother, B.I., 2011. Proopiomelanocortin (POMC) and testing the phylogenetic position of turtles (Testudines). Journal of Zoological Systematics and Evolutionary Research 49, 148-159. Bezy, R.L., Gorman, G.C., Kim, Y.J., Wright, J.W., 1977. Chromosomal and genetic divergence in the fossorial lizards of the family Anniellidae. Systematic Biology 26, 57-71. Bickford, D., Lohman, D.J., Sodhi, N.S., Ng, P.K.L., Meier, R., Winker, K., Ingram, K.K., Das, I., 2007. Cryptic species as a window on diversity and conservation. Trends in Ecology & Evolution 22, 148-155. 66 Bonin, A., Nicole, F., Pompanon, F.C.O., Miaud, C., Taberlet, P., 2007. Population adaptive index: a new method to help measure intraspecific genetic diversity and prioritize populations for conservation. Conservation Biology 21, 697-708. Bons, J., Geniez, P., 1996. Amphibians and Reptiles of Morocco. Asociación Herpetológica Española, Barcelona. Brinkmann, H., Venkatesh, B., Brenner, S., Meyer, A., 2004. Nuclear protein-coding genes support lungfish and not the coelacanth as the closest living relatives of land vertebrates. Proceedings of the National Academy of Sciences of the United States of America 101, 4900-4905. Brown, R.P., Suárez, N.M., Pestano, J., 2002. The Atlas mountains as a biogeographical divide in North-West Africa: evidence from mtDNA evolution in the Agamid lizard Agama impalearis. Molecular Phylogenetics and Evolution 24, 324-332. Busack, S.D., 1988. Biochemical and morphological differentiation in Spanish and Moroccan populations of Blanus and the description of a new species from Northern Morocco (Reptilia, Amphisbaenia, Amphisbaenidae). Copeia, 101-109. Carranza, S., Arnold, E., Wade, E., Fahd, S., 2004. Phylogeography of the false smooth snakes, Macroprotodon (Serpentes, Colubridae): mitochondrial DNA sequences show European populations arrived recently from Northwest Africa. Molecular Phylogenetics and Evolution 33, 523-532. Castresana, J., 2000. Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Molecular Biology and Evolution 17, 540. Clement, M., Posada, D., Crandall, K.A., 2000. TCS: a computer program to estimate gene genealogies. Molecular Ecology 9, 1657-1659. Crottini, A., Dordel, J., Köhler, J., Glaw, F., Schmitz, A., Vences, M., 2009. A multilocus phylogeny of Malagasy scincid lizards elucidates the relationships of the fossorial genera Androngo and Cryptoscincus. Molecular Phylogenetics and Evolution 53, 345-350. De Jong, H., 1998. In search of historical biogeographic patterns in the western Mediterranean terrestrial fauna. Biological Journal of the Linnean Society 65, 99-164. Drummond, A.J., Ashton, B., Buxton, S., Cheung, M., Cooper, A., Heled, J., Kearse, M., Moir, R., Stones-Havas, S., Sturrock, S., 2010. Geneious version 5.3. Available at http://www.geneious.com. Duggen, S., Hoernle, K., van den Bogaard, P., Rüpke, L., Morgan, J.P., 2003. Deep roots of the Messinian salinity crisis. Nature 422, 602-606. Felsenstein, J., 1981. Evolutionary trees from DNA sequences: a maximum likelihood approach. Journal of Molecular Evolution 17, 368-376. Felsenstein, J., 1985. Confidence limits on phylogenies: an approach using the bootstrap. Evolution, 783-791. Fritz, U., Ayaz, D., Hundsdörfer, A.K., Kotenko, T., Guicking, D., Wink, M., Tok, C.V., Çiçek, K., Buschbom, J., 2009. Mitochondrial diversity of European pond turtles (Emys orbicularis) in Anatolia and the Ponto-Caspian Region: Multiple old refuges, hotspot of extant diversification and critically endangered endemics. Organisms Diversity & Evolution 9, 100-114. Fritz, U., Barata, M., Busack, S.D., Fritzsch, G., Castilho, R., 2006. Impact of mountain chains, sea straits and peripheral populations on genetic and taxonomic structure of a freshwater turtle, Mauremys leprosa (Reptilia, Testudines, Geoemydidae). Zoologica Scripta 35, 97-108. 67 Fritz, U., Fritzsch, G., Lehr, E., Ducotterd, J.M., Müller, A., 2005. The Atlas Mountains, not the Strait of Gibraltar, as a biogeographic barrier for Mauremys leprosa (Reptilia: Testudines). Salamandra 41, 97-106. Gans, C., 1974. Biomechanics: an approach to vertebrate biology. Philadelphia, Lippincott. Gans, C., 1978. The characteristics and affinities of the Amphisbaenia. The Transactions of the Zoological Society of London 34, 347-416. Gans, C., 1990. Patterns in amphisbaenian biogeography: a preliminary analysis. Vertebrates in the Tropics. Alexander Koenig Zoological Research Institute and Zoological Museum, Bonn, 133-143. Gans, C., 2005. Checklist and bibliography of the amphisbaenia of the world. Bulletin of the American Museum of Natural History, 1-130. Godinho, R., Crespo, E.G., Ferrand, N., 2008. The limits of mtDNA phylogeography: complex patterns of population history in a highly structured Iberian lizard are only revealed by the use of nuclear markers. Molecular Ecology 17, 4670-4683. Gómez, A., Lunt, D.H., 2007. Refugia within refugia: patterns of phylogeographic concordance in the Iberian Peninsula. In: Weiss, S., Ferrand, N. (Eds.), Phylogeography of southern European refugia. Springer, pp. 155-188. Gonçalves, D.V., Brito, J.C., Crochet, P.A., Geniez, P., Padial, J.M., James Harris, D., 2012. Phylogeny of North African Agama lizards (Reptilia: Agamidae) and the role of the Sahara desert in vertebrate speciation. Molecular Phylogenetics and Evolution 64, 582– 591. Gonçalves, H., Martínez-Solano, I., Pereira, R.J., Carvalho, B., García-París, M., Ferrand, N., 2009. High levels of population subdivision in a morphologically conserved Mediterranean toad (Alytes cisternasii) result from recent, multiple refugia: evidence from mtDNA, microsatellites and nuclear genealogies. Molecular Ecology 18, 5143-5160. González, M., 2012. Determining the difference in the geographic overlap of the potential distribtuion of the green and ocellated lizards at continental and regional extents in the Mediterranean Basin. (MSc Thesis) University of Twente. Graybeal, A., 1994. Evaluating the phylogenetic utility of genes: a search for genes informative about deep divergences among vertebrates. Systematic Biology 43, 174-193. Gündüz, İ., Jaarola, M., Tez, C., Yeniyurt, C., Polly, P.D., Searle, J.B., 2007. Multigenic and morphometric differentiation of ground squirrels (Spermophilus, Scuiridae, Rodentia) in Turkey, with a description of a new species. Molecular Phylogenetics and Evolution 43, 916-935. Guindon, S., Gascuel, O., 2003. A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Systematic Biology 52, 696-704. Harris, D.J., Sinclair, E.A., Mercader, N.L., Marshall, J.C., Crandall, K.A., 1999. Squamate relationships based on C-mos nuclear DNA sequences. Herpetological Journal 9, 147-152. Hembree, D.I., 2006. Amphisbaenian paleobiogeography: evidence of vicariance and geodispersal patterns. Palaeogeography, Palaeoclimatology, Palaeoecology 235, 340-354. Hewitt, G.M., 1996. Some genetic consequences of ice ages, and their role, in divergence and speciation. Biological Journal of the Linnean Society 58, 247-276. Hewitt, G.M., 1999. Post glacial re-colonization of European biota. Biological Journal of the Linnean Society 68, 87-112. Hewitt, G.M., 2000. The genetic legacy of the Quaternary ice ages. Nature 405, 907-913. 68 Hewitt, G.M., 2001. Speciation, hybrid zones and phylogeography - or seeing genes in space and time. Molecular Ecology 10, 537-549. Hewitt, G.M., 2004. Genetic consequences of climatic oscillations in the Quaternary. Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences 359, 183. Hoegg, S., Vences, M., Brinkmann, H., Meyer, A., 2004. Phylogeny and comparative substitution rates of frogs inferred from sequences of three nuclear genes. Molecular Biology and Evolution 21, 1188-1200. Hsü, K.J., 1973. Late Miocene desiccation of the Mediterranean. Nature 242, 240-244. Huelsenbeck, J.P., Bollback, J.P., 2001a. Empirical and hierarchical Bayesian estimation of ancestral states. Systematic Biology 50, 351-366. Huelsenbeck, J.P., Bollback, J.P., 2001b. Empirical and hierarchical Bayesian estimation of ancestral states. Systematic Biology 50, 351. Huelsenbeck, J.P., Ronquist, F., 2001. MRBAYES: Bayesian inference of phylogenetic trees. Bioinformatics 17, 754-755. IUCN, 2012. IUCN Red List of Threatened Species Version 2012.1. www.iucnredlist.org. Downloaded on 15 September 2012. Katoh, K., Misawa, K., Kuma, K., Miyata, T., 2002. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic acids research 30, 30593066. Kearney, M., 2003. Systematics of the Amphisbaenia (Lepidosauria: Squamata) based on morphological evidence from recent and fossil forms. Herpetological Monographs, 1-74. Kearney, M., Stuart, B.L., 2004. Repeated evolution of limblessness and digging heads in worm lizards revealed by DNA from old bones. Proceedings of the Royal Society B: Biological Sciences 271, 1677-1687. Kocher, T.D., Thomas, W.K., Meyer, A., Edwards, S.V., Paabo, S., Villablanca, F.X., Wilson, A.C., 1989. Dynamics of mitochondrial DNA evolution in animals: Amplification and sequencing with conserved primers. Proceedings of the National Academy of Sciences of the United States of America 86, 6196-6200. Kornilios, P., Ilgaz, C., Kumlutas, Y., Giokas, S., Fraguedakis-Tsolis, S., Chondropoulos, B., 2011. The role of Anatolian refugia in herpetofaunal diversity: an mtDNA analysis of Typhlops vermicularis Merrem, 1820 (Squamata, Typhlopidae). Amphibia-Reptilia 32, 351-363. Krijgsman, W., Garcés, M., Agustí, J., Raffi, I., Taberner, C., Zachariasse, W.J., 2000. The Tortonian salinity crisis' of the eastern Betics (Spain). Earth and Planetary Science Letters 181, 497-511. Krijgsman, W., Hilgen, F.J., Raffi, I., Sierro, F.J., Wilson, D.S., 1999. Chronology, causes and progression of the Messinian salinity crisis. Nature 400, 652-655. Lefébure, T., Douady, C., Gouy, M., Trontelj, P., Briolay, J., Gibert, J., 2006. Phylogeography of a subterranean amphipod reveals cryptic diversity and dynamic evolution in extreme environments. Molecular Ecology 15, 1797-1806. Librado, P., Rozas, J., 2009. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451. Lovejoy, N.R., Collette, B.B., McEachran, J., 2001. Phylogenetic relationships of New World needlefishes (Teleostei: Belonidae) and the biogeography of transitions between marine and freshwater habitats. Copeia 2001, 324-338. 69 Macey, J.R., Papenfuss, T.J., Kuehl, J.V., Fourcade, H.M., Boore, J.L., 2004. Phylogenetic relationships among amphisbaenian reptiles based on complete mitochondrial genomic sequences. Molecular Phylogenetics and Evolution 33, 22-31. Makarenkov, V., Kevorkov, D., Legendre, P., 2006. Phylogenetic network construction approaches. Applied mycology and biotechnology 6, 61-97. Martínez-Solano, I., Teixeira, J., Buckley, D., García-París, M., 2006. Mitochondrial DNA phylogeography of Lissotriton boscai (Caudata, Salamandridae): evidence for old, multiple refugia in an Iberian endemic. Molecular Ecology 15, 3375-3388. Mayr, E., 1949. Speciation and selection. Proceedings of the American Philosophical Society 93, 514-519. Mendonça, B., Harris, D.J., 2007. Genetic variation within Trogonophis wiegmanni Kaup 1830. Belgian Journal of Zoology 137, 239-242. Mott, T., Vieites, D.R., 2009. Molecular phylogenetics reveals extreme morphological homoplasy in Brazilian worm lizards challenging current taxonomy. Molecular Phylogenetics and Evolution 51, 190-200. Müller, J., Hipsley, C.A., Head, J.J., Kardjilov, N., Hilger, A., Wuttke, M., Reisz, R.R., 2011. Eocene lizard from Germany reveals amphisbaenian origins. Nature 473, 364-367. Mulvaney, A., Castoe, T.A., Ashton, K.G., Krysko, K.L., Parkinson, C.L., 2005. Evidence of population genetic structure within the Florida worm lizard, Rhineura floridana (Amphisbaenia: Rhineuridae). J. Herpetol. 39, 118-124. Myers, N., Mittermeier, R.A., Mittermeier, C.G., da Fonseca, G.A.B., Kent, J., 2000. Biodiversity hotspots for conservation priorities. Nature 403, 853-858. Navas, C.A., Antoniazzi, M.M., Carvalho, J.E., Chaui-Berlink, J.G., James, R.S., Jared, C., Kohlsdorf, T., Dal Pai-Silva, M., Wilson, R.S., 2004. Morphological and physiological specialization for digging in amphisbaenians, an ancient lineage of fossorial vertebrates. Journal of Experimental Biology 207, 2433-2441. Nevo, E., 2001. Evolution of genome–phenome diversity under environmental stress. Proceedings of the National Academy of Sciences 98, 6233. Niemiller, M.L., Near, T.J., Fitzpatrick, B.M., 2012. Delimiting species using multilocus data: diagnosing cryptic diversity in the southern cavefish, Typhlichthys subterraneus (Teleostei: Amblyopsidae). Evolution. Padial, J.M., Miralles, A., De la Riva, I., Vences, M., 2010. The integrative future of taxonomy. Frontiers in Zoology 7. Palumbi, S.R., 1996. Nucleic acids II: the polymerase chain reaction. In: Hillis, D.M., Moritz, C., Mable, B.K. (Eds.), Molecular Systematics. Sinauer Associates, Inc., Sunderland, MA., pp. 205-247. Paulo, O.S., Dias, C., Bruford, M.W., Jordan, W.C., Nichols, R.A., 2001. The persistence of Pliocene populations through the Pleistocene climatic cycles: evidence from the phylogeography of an Iberian lizard. Proceedings of the Royal Society of London. Series B: Biological Sciences 268, 1625. Pearse, D.E., Pogson, G.H., 2000. Parallel evolution of the melanic form of the California legless lizard, Anniella pulchra, inferred from mitochondrial DNA sequence variation. Evolution 54, 1041-1046. Perera, A., Harris, D.J., 2010. Genetic variability within the Oudri's fan-footed gecko Ptyodactylus oudrii in North Africa assessed using mitochondrial and nuclear DNA sequences. Molecular Phylogenetics and Evolution 54, 634-639.