Genetic characterization of species using phylogeographic approaches represents a basic refer- ence to understand their evolutionary history as well as to identify conservation priorities to protect areas of particular interest regarding evolutionary potential. Even in well-studied regions such in- formation is lacking for the majority of species, including many endemic species with reduced dis- tribution ranges. We investigate the phylogeographic pattern of the Iberian frog Rana iberica,an endemic amphibian restricted to Central and North-Western Iberian Peninsula. Using mitochon- drial sequences, we reconstruct the phylogeographic history of the species to test the effect of Quaternary climate changes on the evolutionary diversiﬁcation of lineages, that is, the differenti- ation of mitochondrial lineages and the formation of genetic diversity melting pots, and integrate phylogeographic evidence for future conservation planning. Our results indicate the existence of 3 main mitochondrial lineages differentiated during the Upper Pleistocene. Both historical demo- graphic analyses and climatic niche modeling show a strong effect of glacial climate changes, sug- gesting recurrent range contractions and expansions. Under such circumstances, differentiation took place most likely by isolation in allopatric interglacial refugia. Secondary lineage admixture in northern Portugal generated a broad mixed zone with highest nucleotide diversity. Given its par- ticular evolutionary potential, its reduced distribution and eventual threats under current climate change scenario, conservation priorities should focus on the isolated lineage from Sierra de Guadalupe. Key words: climate niche modeling, conservation, demography, glacial refugia, phylogeography, Amphibia. Conservation genetics takes advantage of the non-uniform geographic (see definitions by Moritz 2002) and the recognition of hotspots of distribution of genetic diversity to select and design the best-fitting intraspecific genetic diversity, that is geographic areas where genetic areas to conserve a species while preserving its genetic signal. Two of diversity is higher than in the remaining areas (Canestrelli et al. the most used criteria to define these conservation areas are the exist- 2014). While the first criterion has been extensively dealt with (Lesica ence of genetically distinct, but geographically isolated populations and Allendorf 1995; Moritz 2002; Vandergast et al. 2008; V C The Author(s) (2018). Published by Oxford University Press. 1 This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact email@example.com Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zoy010/4827690 by guest on 13 July 2018 2 Current Zoology, 2018, Vol. 0, No. 0 Thomassen et al. 2011), the second is often treated lightly, forgetting this species, only a few studies concerning the population genetics of about the evolutionary potential of these genetically complex areas. R. iberica are available (Arano et al. 1993; Mart´n ı ez-Solano et al. The nature and origin of intraspecific genetically diverse areas are 2005), none covering the whole distribution range of the species and subject to intense study in the fields of ecology and evolutionary biol- characterizing phylogeographical patterns. ogy (centers of origin; climatic sanctuaries; secondary contact and ad- In this study, we investigate the mitochondrial DNA (mtDNA) mixture, etc.), but often overlooked for conservation purposes, since phylogeographic patterns in a large sample of R. iberica comprising delimitation of protected areas is often focused on current patterns of its complete geographical distribution. Specifically, we test the effect organismal and ecosystem-level biodiversity, typically ignoring the of Quaternary climate changes on the evolutionary diversification of evolutionary processes involved in the origin of that biodiversity lineages, that is, the differentiation of mitochondrial lineages and (Vandergast et al. 2008). Neglecting the nature of genetically diverse the formation of genetic diversity melting pots by secondary admix- areas makes difficult to decide on the future conservation value of the ture, and integrate the phylogeographic evidences here presented for chosen area but might also alter the evolutionary potential of a given future conservation planning and immediate actions for this species. species. Some generalities have been observed when comparing phylo- Materials and Methods geographic structure among species within large glacial refugia, such as the Mediterranean peninsulas (Go ´ mez and Lund 2007), although Sampling patterns of genetic diversity are often varied and reflect different, We analyzed 162 specimens of R. iberica from 57 localities across specific responses to the same events (Nieto Feliner 2011). One is northwestern Iberia (Figure 1 and Table 1), covering all the species that the majority of species exhibit a previously unanticipated level range across Portugal, Northern Spain, and all known population of fragmentation, indicating the occurrence of multiple refugia isolates of the species. Tissue samples consisted of tail ends of larvae within the peninsulas. This is in accordance with the hypothesis that or toe clips of adult specimens. All individuals were released in the refugia are characterized by stable ecological conditions that have wild after tissue collection. Sample sizes per locality are indicated in allowed the survival of large populations and, in consequence, the Table 1. persistence of a relatively large part of their genetic legacy (Abella ´n and Svenning 2014). In several cases it has been possible even to DNA amplification and sequencing identify pre-Pleistocene lineages, usually restricted to particular geo- DNA was extracted from alcohol-preserved tissues following the graphical areas that can be defined as sanctuaries of intraspecific di- standard proteinase K protocol (Sambrook and Russell 2001). versity (Recuero and Garcı ´a-Parı ´s 2011). In others, differentiation Amplification of a 1,294-base pair fragment of cytochrome oxidase took place by allopatric fragmentation during the repeated I (Cox1) was performed via PCR using the light chain primer P2F4 Pleistocene climate changes (Knowles 2001; Canestrelli et al. 2014). designed specifically for this work (5 -GGCCACTTTACCCGTGA Every genetic lineage may represent an evolutionary unit with poten- 0 0 TATT-3 ) and the heavy chain primer P3R (5 -ATRATATTTATTA tially unique genetic characteristics and adaptive traits within a spe- TYTGAGAAGC-3 )(San Mauro et al. 2004). Amplification condi- cies, and therefore should be taken into consideration in any specific tions were as in Recuero et al. (2010), with annealing at 48 C. conservation plan (Hampe and Petit 2005). Another assumed gener- Sequences were determined in the ABI PRISM 3700 Genetic alization is that population expansions from multiple refugia gener- Analyzer. Sequences were aligned manually with BioEdit 7.0.9 (Hall ate areas of low genetic diversity but, in some cases, may produce 1999). The nucleotide sequences data reported here are registered in admixture zones between divergent lineages. In the later case, high the GenBank Nucleotide Sequence Database under the accession levels of genetic diversity can be reached in what are called melting numbers MG813567-MG813728. pots of diversity (Petit et al. 2003; Dufresnes et al. 2016), in which genetic diversity is enriched by the immigration of individuals from different genetic lineages. These highly genetically diverse popula- Data analysis tions should also be considered for their conservation relevance, We used PartitionFinder v. 1.0.1 (Lanfear et al. 2012) to select the considering that they should present wider alternatives for coping best-fit model of nucleotide substitution considering 3 possible parti- with both natural and anthropogenic environmental changes (Sgro ` tions for each codon position, using a greedy algorithm and the et al. 2011). Endemic species, with reduced distribution areas and Bayesian information criterion. Accordingly, 3 independent parti- often specific ecological constrains, are particularly sensitive to cur- tions were defined for the first (HKY), second (TrNþG), and third rent levels of anthropogenic habitat alterations and in many cases (K80þI) positions. subject to different levels of protection. Phylogeographic studies for Phylogenetic relationships were inferred using MrBayes v3.2 such taxa should be one of the first steps to efficiently plan conserva- (Huelsenbeck and Ronquist 2001; Ronquist et al. 2012). Analyses tion strategies. However, for most species genetic diversity charac- were run for 10 million generations, sampling every 1,000 gener- terization is incomplete if not totally missing. ations. Convergence was assessed through examination of the This is the case of Rana iberica (Boulenger 1879), an endemic spe- standard deviation of split frequencies, which was well below rec- cies restricted to the northwestern Iberian Peninsula and categorized ommended thresholds. A consensus phylogram was computed after as Near Threatened by the IUCN for its currently decreasing popula- discarding trees reconstructed during the default burn-in period. tion trend (Tejedo et al. 2009). This species preferentially inhabits This analysis was run in the Cipres Science Gateway (Miller et al. swift streams in humid, mountainous areas, where it often co-occurs 2010). Maximum-likelihood (ML) analyses were performed with with other endemics such as Chioglossa lusitanica (Bocage 1864), Garli v2.01 (Zwickl 2006) using the heuristic search algorithm with Lissotriton boscai (Lataste, 1879), and Lacerta schreiberi (Bedriaga model parameters estimated with PartitonFinder. We used non- 1878; Teixeira 2008). In addition, R. iberica also has population iso- parametric bootstrapping (100 pseudoreplicates) to assess the stabil- lates in its southern distribution range, possibly raising relevant con- ity of internal branches. We also explored our mtDNA data by esti- servation issues. In spite of the potential phylogeographical interest of mating relationships between haplotypes using the program TCS Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zoy010/4827690 by guest on 13 July 2018 Teixeira et al. Rana iberica Phylogeography 3 Figure 1. Distribution map of Rana iberica in the Iberian Peninsula (shaded), with indication of the geographical origin of samples analyzed in this study and proportion of lineages by locality: lineage A in gray, lineage B in white (Sistema Central sublineage gridded), lineage C in black. Population codes asin Table 1. version 1.21 (Clement et al. 2000). We used BEAST 1.8.4 expansions a unimodal mismatch distribution is expected (Rogers (Drummond and Rambaut 2007) to estimate phylogenetic relation- and Harpending 1992). Tajima’s D statistic tests for a departure ships, divergence times across clades, and demographic history of from neutrality as measured by the difference between the number the different lineages, under a Bayesian Skyline tree prior of segregating sites (h) and the average number of pairwise nucleo- (Drummond et al. 2005). A lognormal relaxed clock model was im- tide differences (p). Population expansions can cause significant plemented preliminarily to test the fit of the data to a molecular negative departures of Tajima’s D from zero (Tajima 1989). Fu’s Fs clock model. As both ucld.mean and coefficient of variation priors statistic provides a test for population growth as it identifies an ex- tend to zero we cannot reject a strict molecular clock. Further ana- cess of rare alleles in an expanding population when compared with lyses were performed under a strict molecular clock for the whole the number expected in a stationary population (Fu 1997). We car- fragment. Lacking internal calibration points, time estimates were ried out the analysis with 5,000 simulations. based on an uniform prior for the substitution rate, bounded at The nucleotide diversity index (p) and the haplotype diversity 0.012 and 0.014 substitutions/site/My, according to previous mo- (Hd) were used to estimate the genetic diversity for the overall data, lecular studies in other European species of Rana (Canestrelli et al. within distinctive phyloclades, geographical subsets, and popula- 2014). Nucleotide substitution models were unlinked according to tions with enough sample size (n> 4) using DnaSP v. 5.10 (Rozas the PartitionFinder results. Analyses were run for 150 million gener- et al. 2003). Population structure was inferred by the analysis of mo- ations, sampling every 15,000 generations. This analysis was run in lecular variance (AMOVA) (Excoffier et al. 1992) provided in the the Cipres Science Gateway (Miller et al. 2010). We used Tracer v. program package Arlequin V. 3.11 (Excoffier et al. 2005). 1.6 to check that effective sample sizes (ESSs) of parameters were We also used BEAST v.1.8.4 to perform continuous diffusion >200 and convergence of parameter estimates across runs, and to analysis. Geographical coordinates for each sequence were included. build the final Bayesian skyline plots (BSP). A maximum clade cred- A jitter module (window size: 0.01) generating random coordinates ibility tree (MCCT) was reconstructed with TreeAnnotator. was included for samples from the same population. Analyses were In addition, pairwise mismatch distribution analyses (Rogers and run for 20,000 million generations under a strict molecular and a Harpending 1992), Tajima’s D (Tajima 1989), and Fu’s Fs(Fu Yule tree prior. We used a Cauchy model to infer the geographical 1997) statistics were carried out in Arlequin V. 3.11 (Excoffier et al. diffusion process through time across branches of the inferred gene 2005), to test for population expansion in the total data and within tree (Lemey et al. 2010). Convergence was assessed using Tracer. A mitochondrial lineages. Mismatch distributions describe the fre- MCCT was reconstructed with TreeAnnotator and used to generate quency of pairwise substitutional differences among individuals. the reconstruction of the diffusion process using the modules These distributions can be informative about the demographic his- “Continuous Tree” and “Time Slicer” in Spread 1.0.7 (Bielejec et al. tory, as in populations that undergone a recent bottleneck and rapid 2011). Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zoy010/4827690 by guest on 13 July 2018 4 Current Zoology, 2018, Vol. 0, No. 0 Table 1. Sampling localities, population codes, sample size (n), haplotypes, and GeneBank accession numbers for the Cox1 sequences obtained in this study Locality Code n haplotypes GeneBank Accession Nos Portugal: Viana do Castelo VIA 3 H71 MG813567–MG813569 Portugal: Ponte de Lima LIM 4 H19, H43 MG813570–MG813573 Portugal: Paredes de Coura COU 3 H11, H44, H71 MG813574–MG813576 Portugal: Castro Laboreiro LAB 2 H18, H70 MG813577–MG813578 Portugal: Gere ˆ s GER 6 H11, H15, H46, H54, H71 MG813579–MG813584 Portugal: Boticas BOT 2 H23, H41 MG813586–MG813587 Portugal: Montesinho MON 6 H16, H23, H25, H42, H46 MG813588–MG813593 Portugal: Vinhais VNH 3 H23, H27 MG813594–MG813596 Portugal: Santa Comba de Rossas NOG 2 H23, H46 MG813597–MG813598 Portugal: Santo Tirso STS 2 H13, H28 MG813599–MG813600 Portugal: Braga TIB 1 H21 MG813601 Portugal: Vila da Feira FEI 1 H68 MG813602 Portugal: Lousada LSD 2 H14, H40 MG813603–MG813604 Portugal: S. Pedro do Sul SPS 2 H58 MG813605–MG813606 Portugal: Ribadouro RIB 2 H23, H61 MG813607–MG813608 Portugal: Fafe FAF 1 H23 MG813609 Portugal: Montemuro MTM 6 H15, H47, H48, H55, H59, H60 MG813610–MG813615 Portugal: Amarante AMA 3 H23, H38 MG813616–MG813618 Portugal: Ribeira de Pena RPN 1 H69 MG813619 Portugal: Valongo VAL 4 H23, H24, H28, H62 MG813620–MG813623 Portugal: Penacova PNC 1 H45 MG813585 Portugal: Sever do Vouga SVG 2 H66, H67 MG813624–MG813625 Portugal: Vila Pouca de Aguiar VPA 3 H23, H38, H71 MG813626–MG813628 Portugal: Meda MED 1 H57 MG813629 Portugal: Cernache do Bonjardim BNJ 3 H55, H64 MG813630–MG813632 Portugal: Lous~ a LOU 5 H31, H33, H63 MG813633–MG813637 Portugal: Arganil ARG 2 H33 MG813638–MG813639 Portugal: S. Pedro de Moel SPM 7 H55, H65 MG813640–MG813646 Portugal: Pomar, Alvito da Beira MRD 2 H38 MG813648–MG813649 Portugal: Penalva do Castelo PVC 1 H56 MG813647 Portugal: Castelo Branco CBR 2 H34, H36 MG813650–MG813651 Portugal: Fund~ ao FUN 1 H32 MG813652 Portugal: Manteigas MTG 1 H33 MG813653 Portugal: Guarda GRD 1 H35 MG813654 Portugal: Malcata MAL 5 H33, H49 MG813655–MG813659 Portugal: Vila de Rei REI 2 H39, H55 MG813660–MG813661 Portugal: S. Mamede SMA 8 H37 MG813662–MG813669 Spain: Pontevedra PON 6 H4, H5, H6, H11 MG813670–MG813675 Spain: Santiago de Compostela COM 4 H7, H11 MG813676–MG813679 Spain: A Coruna ~ COR 3 H1, H11 MG813680–MG813682 Spain: Lugo LUG 5 H9, H11, H22 MG813683–MG813687 Spain: Donı´s DON 1 H26 MG813688 Spain: Ponferrada PNF 4 H11, H23 MG813689–MG813692 Spain: Ourense OUR 2 H23 MG813693–MG813694 Spain: Monforte de Lemos MNF 2 H8, H12 MG813695–MG813696 Spain: Toreno TOR 2 H23 MG813697–MG813698 Spain: Navia NAV 3 H2, H22 MG813699–MG813701 Spain: Oviedo OVI 3 H3, H11 MG813702–MG813704 Spain: Pena ~ de Francia FRA 3 H50 MG813705–MG813707 Spain: Sierra de Guadalupe, Navezuelas GUA 7 H29, H30 MG813708–MG813714 Spain: Alava ALA 3 H10 MG813715–MG813717 Spain: Madrid MAD 3 H53 MG813718–MG813720 Spain: El Tiemblo, Avila AVI 2 H53 MG813721–MG813722 Spain: Mengamunoz, ~ Avila AVL 1 H53 MG813723 Spain: Plasencia PLA 2 H51, H52 MG813724–MG813725 Spain: Puebla de Sanabria SAN 2 H17, H23 MG813726–MG813727 Spain: Leo ´ n LEO 1 H20 MG813728 Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zoy010/4827690 by guest on 13 July 2018 Teixeira et al. Rana iberica Phylogeography 5 Contemporary and past niche modeling was performed with Historical demographic analyses Maxent v.3.3.3k (Phillips et al. 2006; Phillips and Dud´k ı 2008). We BSP results indicate demographic expansions for the whole species gathered a total of 510 localities from the whole distribution area, as well as for each of the 3 main mitochondrial lineages associated including the complete sampling plus personal observations and some with the last glacial maximum (LGM), and a progressive stabiliza- of the GBIF records (many records were duplicated or based in obvi- tion in the Holocene (Figure 4). These results are generally in ac- ous misidentifications; doi:10.15468/dl.qp7txp). Bioclimatic layers cordance with other historical demographic analyses: the mismatch (30 arcsec), covering the range of the species, were downloaded from analysis produced a bimodal distribution for the total data and uni- the WorldClim database (http://www.worldclim.org) and analyzed modal distributions for all phylogroups and geographic subsets ana- with ENMTools 1.4.3 (Warren et al. 2010) to check for correlations lyzed. Tajima’s D statistic detected significant deviations from among them using the Pearson correlation coefficient. Past climate neutrality in lineages A and B. With the exception of lineage C, Fu’s simulations followed the Community Climate System Model (http:// Fs statistics presented significant negative departures from zero in www2.cesm.ucar.edu/). Taking into consideration the biological rele- all 3 lineages and the total data set. vance of the different bioclimatic variables and the observed correl- Mean pairwise sequence divergence among haplotypes varied ations (pairs of variables with Pearson coefficient value 0.8 were from 0.080% to 1.050% and the overall nucleotide diversity (p) considered highly correlated), we considered Bio4 (Temperature averaged 0.004. Lineage B exhibited the higher genetic diversity Seasonality), Bio10 (Mean Temperature of Warmest Quarter), Bio12 (h ¼ 41, Hd ¼ 0.958), followed by lineage A (h ¼ 28, Hd ¼ 0.887), (Annual Precipitation), Bio15 (Precipitation Seasonality), and Bio17 and away down by lineage C (h ¼ 2, Hd ¼ 0.286). Nucleotide diver- (Precipitation of Driest Quarter) to generate the model, using 70% of sity is higher in populations from intermediate latitudes (Portugal the localities as training data and 30% to test the model. The model north of the Douro River) (Figure 5), while no geographic structure was evaluated using “area under the curve” (AUC) values obtained was observed for Hd. for the test data (Swets 1988; Elith 2000). For the AMOVA analysis, we grouped samples by lineages A–C. The AMOVA revealed significant genetic structuring across all hier- archical levels (all P< 0.0001), with 69.1% of the variation result- Results ing of differences between clades (U ¼ 0.691). Only 18.8% of the CT total variation was attributed to differences among populations mtDNA structure and diversity within these clades (U ¼ 0.609), while the remaining 12.1% was SC A total of 67 Cox1 variable sites and 71 haplotypes were detected due to variation within populations (U ¼ 0.879). ST among the 162 analyzed sequences. Three main matrilineal lineages were recovered by phylogenetic approximations (labeled A–C, Figure 2), with marked phylogeographic structure. With the excep- Glacial and interglacial distribution tion of a single individual from the Serra de Montemuro (MTM), all The models of climatic niche (Figure 6) showed high AUC test val- samples in lineage A were found north of the Douro River, being the ues (0.933, SD ¼ 0.007) and significance for every binomial omis- only one present in all populations from northern Spain, including sion test, supporting a good performance of the model. The 5 those from the Basque Country. Lineage B was found across most of included variables contributed unequally to the model, with 2 of the range of R. iberica, mostly in the samples originated from them (Bio12, Bio17) accounting for 74.2% of relative contributions Central Portugal and the Spanish Sistema Central Mountains. to the Maxent model. According to the generated models, the distri- However, this lineage is also widely present in populations from bution of R. iberica has suffered a series of contractions and expan- Portugal north of the Douro River, in sympatry with lineage A. sions associated to the climatic oscillations of the Pleistocene, with Lineage C is exclusive from the populations of Sierra de Guadalupe wider potential distribution during LGM and more restricted ones (GUA) in the southern limits of the species distribution in Spain in interglacial periods, with relatively stable conditions in the (Figures 1 and 2). North-west of the Iberian Peninsula. All 3 groups are supported as monophyletic, but the sister rela- tionship between lineages A and C is only supported by MrBayes and ML, but not by BEAST analyses (Figure 2). Additionally, haplo- Discussion types of lineage B from the Sistema Central Mountains form a shal- low but well-supported clade differentiated from the Portuguese Evolutionary history and phylogeographic structure populations. This clade is geographically restricted and only 1 indi- The Iberian Peninsula is one of the main biodiversity hotspots in the vidual from the Portuguese population of Malcata, just some tens of Western Palearctic, harboring not only a considerable number of en- kilometers west of Spanish Sistema Central, falls within it. The top- demic species (Me ´dail and Que ´zel 1997), but also a large part of the ology of the haplotype network depicts 2 main star-type groups of intraspecific genetic diversity of widespread species is present in that re- haplotypes, corresponding to lineages A and B (Figure 2), with gion (Hewitt 2011). In many cases, this diversity reflects pre-Pleistocene strong genetic structure within each clade. The lineage C exhibits evolutionary processes and has survived the numerous cycles of very low genetic structure with only 2 haplotypes found. Age esti- Pleistocene climatic oscillations in sanctuaries of diversity (Recuero and mates provide wide intervals for all 3 main lineages and their most Garcı ´a-Parı ´s 2011; Dufresnes et al. 2016) representing type “S” species recent common ancestor (Figure 2), with differentiation likely occur- (sensu Recuero and Garcı ´a-Parı ´s 2011) such as the widespread species ring during the Middle to Upper Pleistocene. Rana temporaria Linnaeus (Vences et al. 2013, 2017). Continuous diffusion analyses inferred an ancestral area located in In other cases, however, species have been severely affected by the current border of north-western Portugal with Spain, with subse- glacial and interglacial changes, presenting a comparatively reduced quent expansions to the North, followed by colonization of marginal genetic diversity as only a small fraction of populations survived in areas in the Sistema Central Mountains, eastern Cantabric Region true refugia (Recuero and Garcı ´a-Parı ´s 2011; Martı ´nez-Freirı ´a et al. and Sierra de Guadalupe, and ending with a progressive expansion in 2015). Rana iberica is an old lineage, sharing a Miocene common Portugal into an admixture zone in northern Portugal (Figure 3). ancestor with a clade including R. temporaria, R. arvalis Nilsson, Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zoy010/4827690 by guest on 13 July 2018 6 Current Zoology, 2018, Vol. 0, No. 0 Figure 2. Mitochondrial diversity in R. iberica. Upper: Maximum clade credibility tree recovered by BEAST, showing time estimates for supported clades. Node bars represent 95% HPD intervals for node ages. Also shown are MrBayes/BEAST/ML posterior probabilities and bootstrap values for main nodes (only values >0.9 in at least one of the analyses or bootstrap values >50 are shown). The scale is in thousands of years. Lower: Maximum parsimony network for Cox1 haplo- types observed in R. iberica showing the 3 haplogroups recovered. Each circle represents a different haplotype (codes as in Table 1) with size proportional to its relative frequency. Small dots represent hypothetical unsampled haplotypes. and R. pyrenaica Serra-Cobo (Yuan et al. 2016). However, the cur- analyses of R. iberica, which show that population expansions of rent mitochondrial phylogeographic pattern observed in R. iberica is every lineage during the last glacial period finished with, or little the consequence of much recent evolutionary processes that took after, the beginning of the current postglacial era. The existence of a place only in the Middle to Upper Pleistocene, in direct association Pleistocene fossil record for this species in Burgos, dated at approxi- with the last 2 or 3 glacial periods. Similar patterns of reduced diver- mately 37,600 years BP (Esteban and Sanchı ´z 1991) and outside the sity, corresponding to type “R” species (sensu Recuero and Garcı ´a- current distribution of R. iberica, indicates that the species reached a Parı ´s 2011) have also been observed in other species of Rana, includ- wider range than its current one with subsequent extinction of per- ing widespread ones such as Rana dalmatina (Vences et al. 2013) ipheral populations. Of course, results presented here correspond to but also in other Mediterranean endemics. For example, in the patterns of mitochondrial diversity and must be complemented with Apennine Peninsula, R. italica Dubois shows 2 mitochondrial lin- informative nuclear marker data to fully reconstruct the species evo- eages with similar ages as in R. iberica (Canestrelli et al. 2008), lutionary history, including recent demographic and dispersal events while another Iberian endemic, R. pyrenaica, presents even higher that may have not leave trace at the mitochondrial level (Zink and diversity reduction with only 2 cytochrome b haplotypes differing in Barrowclough 2008). 1 transversion (Carranza and Arribas 2008). Also, climate niche modeling for R. iberica suggests that its dis- Rana iberica, and probably most of Mediterranean populations tribution was much extended, at least potentially, during glacial of Rana, could be considered as relictual units currently surviving in periods, spreading over a large part of the Iberian Peninsula. interglacial refugia, as suggested also by our historical demographic Potential distribution would be much reduced with interglacial Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zoy010/4827690 by guest on 13 July 2018 Teixeira et al. Rana iberica Phylogeography 7 A B Figure 3. Continuous diffusion phylogeographical reconstruction in R. iberica, representing ancestral distribution and expansion of the identiﬁed mitochondrial lineages. (A) 280 Kyears before present (Kybp); (B) 200 Kybp; (C) 60 Kybp; and (D) 20 Kybp. Figure 4. Bayesian Skyline Plot (BSP) showing historical demographic trends in R. iberica; x-axis: time in thousands of years before present; y-axis: estimated population size [units ¼ Net, the product of effective population size and generation length in years (log transformed)]. Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zoy010/4827690 by guest on 13 July 2018 8 Current Zoology, 2018, Vol. 0, No. 0 0.0050 0.0045 0.0040 0.0035 0.0030 0.0025 0.0020 0.0015 0.0010 0.0005 Latitude Figure 5. Latitudinal distribution of nucleotide diversity in R. iberica along a South–North transect along the species main distribution area (Central Portugal to NW Spain). conditions, confined into a range similar to that observed today. Successive events of range contractions, associated with strong demographic oscillations, provide an excellent scenario to promote the geographical isolation of more or less reduced population groups. Under such circumstances, processes of allopatric differenti- ation gave rise to the main lineages identified by our phylogenetic analyses. Even if more recent, this could be also the case of the Sistema Central sub-lineage. In this case, the species would have been isolated in a small refugium between Spain and Portugal, from where it expanded eastward, as suggested by the decrease of genetic diversity in that direction (Martı ´nez-Solano et al. 2005), while west- ward expansion of this mitochondrial lineage would have been pre- vented by competitive exclusion (Waters 2011). Alternatively, isolation by the presence of geological barriers to gene flow has been proposed in other species with overlapping dis- tributions. The Douro River has been suggested to have acted as an Figure 6. Climate niche models for R. iberica showing potential distribution of important barrier for other amphibian species such as L. boscai and the species at (A) the present, (B) during the Last Glacial Maximum (LGM, Alytes obstetricans (Laurenti 1768), with different phylloclades 21 Ka), and (C) during the last interglacial period (LIG, up to ca. 120–140 Ka). being observed northward and southward of this river (Fonseca et al. 2003; Martı ´nez-Solano et al. 2006; Gonc ¸alves et al. 2015). Additionally, Alexandrino et al. (2000) also indicate this river as a Conservation implications major barrier for C. lusitanica during its post-glacial expansion, re- Conservation of relictual species with relatively large distribution sulting in depleted levels of genetic diversity for this species in north- ranges is problematic. The case of R. iberica exemplifies the situ- ern Portugal and in the Spanish regions of Galicia and Asturias. ation. Despite of being an old lineage, intraspecific diversification in However, the extensive presence of R. iberica lineage B north and R. iberica is very limited and recent, suggesting that almost all units south of the Douro River contradicts this hypothesis. resulting from pre-Pleistocenic diversification processes along the Lineages A and B are present in several populations in northern evolutionary history of the lineage became extinct. The shallow gen- Portugal. This lineage combination results in the highest population etic differentiation obtained for R. iberica is corroborated by previ- nucleotide diversity across the entire R. iberica range. The most ous allozyme and microsatellite studies, which detected low genetic likely explanation for this is the formation of a diversity melting pot variability among populations (D ¼ 0.000–0.020, Arano et al. NEI generated by admixture of the 2 lineages after range expansions, 1993; F ¼ 0.145–0.250, Martı ´nez-Solano et al. 2005). Evidence of ST and this is supported by continuous diffusion analysis. This region is a recent geographic distribution contraction, based on recent fossil characterized by climatic niche modeling as one of the most stable records (Esteban and Sanchı ´z 1991), combined with a reduction in areas for R. iberica, so the admixture zone may have formed early climatic suitability for amphibians and reptiles in the Iberian after lineage differentiation and maintained ever since without much Peninsula (Arau ´ jo et al. 2006) makes the situation complicated. change, allowing the accumulation and maintenance of intraspecific Teixeira and Arntzen (2002) predicted a reduction of favorable area diversity, while its geographic limits were probably favored by com- for the distribution of the co-distributed and ecologically similar petitive exclusion of mitochondrial genomes (Waters 2011). C. lusitanica up to 35% until the end of this century. The isolated Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zoy010/4827690 by guest on 13 July 2018 Nucleotide diversity 39.30 39.75 40.10 40.28 40.98 41.18 41.76 41.84 41.89 42.49 42.49 42.93 43.03 Teixeira et al. Rana iberica Phylogeography 9 populations of R. iberica at Serra de S~ ao Mamede (SMA), S~ ao Pedro Funding de Moel (SPM), and Basque Country (ALA) exhibit very low genetic This work was supported by the Portuguese Foundation for Science and differentiation from the nearby populations, which suggest that they Technology (FCT) through post doc grants (SFRH/BPD/27173/2006 to J.T.) are very recent and that the species distribution is already retreating and (SFRH/BPD/26555/2006 to H.G.). Currently, H.G. is supported by a from its glacial optimum. To complicate matters, most of the mito- postdoctoral Grant from FCT (SFRH/BPD/102966/2014). Additional support chondrial genetic diversity of R. iberica is located in the southern was provided from project CGL2015-66571-P by Spanish Ministerio de Economı´a y Competitividad (MINECO/FEDER). part of its range, a scenario also described for several co-distributed species (Paulo et al. 2001; Alexandrino et al. 2002, 2007; Martı ´nez- Solano et al. 2006; Teixeira et al. 2015). Lineage B exhibits by far References the higher mitochondrial diversity of the species. Lineage A follows with a moderate diversity and lineage C exhibits a highly depauper- Abella ´ n P, Svenning JC, 2014. Refugia within refugia-patterns in endemism ate genetic diversity. and genetic divergence are linked to Late Quaternary climate stability in the In terms of conservation efforts, it is clear that the southernmost Iberian Peninsula. Biol J Linn Soc 113:13–28. population from the Sierra de Guadalupe deserves a special atten- Alexandrino J, Froufe E, Arntzen JW, Ferrand N, 2000. Genetic subdivision, glacial refugia and postglacial recolonization in the golden-striped salaman- tion due to its very restricted distribution in the species rear edge, der, Chioglossa lusitanica (Amphibia: urodela). Mol Ecol 9:771–81. that is, marginal, isolated populations that have persisted through Alexandrino J, Arntzen JW, Ferrand N, 2002. Nested clade analysis and the genetic some of the Quaternary climate changes (Hampe and Petit 2005). In evidence for population expansion in the phylogeography of the golden-striped this respect, Hampe and Petit (2005) enhanced the potential critical salamander Chioglossa lusitanica (Amphibia: urodela). Heredity 88:66–74. importance of these populations as long-term stores of species’ gen- Alexandrino J, Teixeira J, Arntzen JW, Ferrand N, 2007. The historical biogeog- etic diversity and focuses of evolutionary and speciation potential. raphy and conservation of the Golden-striped salamander Chioglossa lusitan- In addition, the tendency for climatic warming observed during the ica: bringing together ecological and genetic data. In: Weiss S, Ferrand N, last decades (Hughes 2000) is expected to have a strong effect on editors. Phylogeography of European Refugia. New York: Springer, 189–205. amphibians (Pounds 2001) and to produce a poleward range shift of Arano B, Esteban M, Herrero P, 1993. Evolutionary divergence of the species organisms if dispersal is unlimited (Parmesan and Yohe 2003; of the Rana temporaria group. Ann Sci Nat 14:49–57. Arau ´ jo MB, Thuiller W, Pearson RG, 2006. Climate warming and the decline Arau ´ jo et al. 2006). Therefore, because of its peripheral of amphibians and reptiles in Europe. J Biogeogr 33:1712–28. Mediterranean location and its geographic isolation, the Sierra de Bedriaga J, 1878. Herpetologische studien. Arch Naturgesch 44:259–320. Guadalupe population is prone to a higher susceptibility to climate Bielejec F, Rambaut A, Suchard MA, Lemey P, 2011. SPREAD: spatial phylo- warming and an increased risk of extinction. In this respect, the genetic reconstruction of evolutionary dynamics. Bioinformatics 27:2910–12. population of the Sierra de San Pedro raises an intriguing problem. Bocage JVB, 1864. Note sur un nouveau batracien du Portugal, Chioglossa Given the potential phylogeographical interest of this population lusitanica, et sur une grenouille nouvelle de l’Afrique occidentale. Rev Mag isolate, at mid-distance between the population isolates of Sierra de Zool Ser 2 16:248–53. Guadalupe (GUA) and SMA in Portugal, we realized several field ´ Boulenger GA, 1879. Etude sur les grenouilles rousses, Ranae temporariae et trips to that area, but no larval or adult specimens were found. As description d’espe ` ces nouvelles ou me ´ connues. Bull Soc Zool Fr 4:158–93. so, we believe that this species does not currently occur in that Canestrelli D, Bisconti R, Sacco F, Nascetti G, 2014. What triggers the rising of an intraspeciﬁc biodiversity hotspot? Hints from the agile frog. Sci Rep 4:5042. mountain area, which is marked with a “?” in Figure 1. This fact Canestrelli D, Cimmaruta R, Nascetti G, 2008. Population genetic structure can be explained either by an unlikely erroneous location of the spe- and diversity of the Apennine endemic stream frog Rana italic: insights on cies attributed for that area by Meijide (1985), or by an alarming the Pleistocene evolutionary history of the Italian peninsular biota. Mol evidence of local extinction of southern rear edge populations. Ecol 17:3857–72. Considering this second hypothesis, as this location constitutes the Carranza S, Arribas O, 2008. Genetic uniformity of Rana pyrenaica smallest population isolate with the lowest altitude amplitude from Serra-Cobo, 1993 across its distribution range: a preliminary study with all rear edge populations, it would be the most vulnerable to local mtDNA sequences. Amphib Reptil 29:579–82. extinction events, and it highlights the risks for the rest of isolated, Clement M, Posada D, Crandall KA, 2000. TCS: a computer program to esti- southern populations and particularly for the endemic mitochon- mate gene genealogies. Mol Ecol 9:1657–60. drial lineage of Sierra de Guadalupe. Phylogeographic studies are a Drummond AJ, Rambaut A, 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol 7:214. necessary tool to identify the existence of isolated lineages, with in- Drummond AJ, Rambaut A, Shaphiro B, Pybus OG, 2005. Bayesian coales- dependent evolutionary trajectories from the main population core cent inference of past population dynamics from molecular sequences. Mol of otherwise genetically homogeneous taxa. Direct conservation Biol Evol 22:1185–92. measures are necessary to avoid losing these southernmost rear edge Dufresnes C, Litvinchuk SN, Leuenberger J, Ghali K, Zinenko O et al., 2016. populations and unique lineages, because they might represent not Evolutionary melting pots: a biodiversity hotspot shaped by ring diversiﬁca- only a critical fraction of the total amount of genetic diversity within tions around the Black Sea in the Eastern tree frog Hyla orientalis. Mol Ecol the species, as it is the case for R. iberica, but also an important 25:4285–300. element of the species evolutionary potential. Elith J, 2000. Quantitative methods for modeling species habitat: comparative performance and an application to Australian plants. In: Ferson S, Burgman M, editors. Quantitative Methods for Modeling Species Habitat. New York: Acknowledgments Springer, 39–58. Esteban M, Sanchı ´z B, 1991. On the presence of Rana iberica in the Collecting permits were issued by ICN for Portugal and by the “Juntas de Pleistocene of Burgos, Spain. Rev Esp Herp 5:93–9. Medio Ambiente” of Galicia, Asturias, Extremadura, Castilla y Leo ´ n, and Excofﬁer L, Smouse PE, Quattro JM, 1992. Analysis of molecular variance Paı´s Vasco for Spain. We are grateful to Alberto Gosa ´ , Ana Fonseca, inferred from metric distances among DNA haplotypes: application to Armando Loureiro, Bruno Ribeiro, Clau ´ dia Soares, Fernando Queiro´s, human mitochondrial DNA restriction data. Genetics 131:479–91. Fernando Sequeira, Inigo ~ Mart´ne ı ´ z-Solano, Neftal´ı Silero, Sara Rocha, Excofﬁer L, Laval G, Schneider S, 2005. Arlequin ver. 3.0: an integrated software Raquel Ribeiro, and Ricardo Pereira for their help during ﬁeldwork and/or package for population genetics data analysis. Evol Bioinform Online 1:47–50. for providing samples. Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zoy010/4827690 by guest on 13 July 2018 10 Current Zoology, 2018, Vol. 0, No. 0 Fonseca A, Arntzen JW, Crespo EG, Ferrand N, 2003. Regional differentiation Phillips SJ, Dudı´k M, 2008. Modeling of species distribution with MAXENT: in the common midwife toad Alytes obstetricans in Portugal: a picture from new extensions and a comprehensive evaluation. Ecography 31:161–75. mitochondrial DNA. Z Feldherpetol 10:83–9. Phillips S, Anderson R, Schapire R, 2006. Maximum entropy modeling of spe- Fu YX, 1997. Statistical tests of neutrality of mutations against population cies geographic distributions. Ecol Modell 190:231–59. growth, hitchhiking and background selection. Genetics 147:915–25. Pounds JA, 2001. Climate and amphibian declines. Nature 410:639–40. Go ´ mez A, Lunt DH, 2007. Refugia within refugia: patterns of phylogeographic Recuero E, Garcı´a-Parı´s M, 2011. Evolutionary history of Lissotriton helveti- concordance in the Iberian Peninsula. In: Weiss S, Ferrand N, editors. cus: multilocus assessment of ancestral vs. recent colonization of the Iberian Phylogeography of European Refugia. Amsterdam: Springer, 155–88. Peninsula. Mol Phylogenet Evol 60:170–82. Gonc ¸alves H, Maia-Carvalho B, Sousa-Neves T, Garcı ´a-Parı´s M, Sequeira F Recuero E, Cruzado-Cortes J, Parra-Olea G, Zamudio KR, 2010. Urban aquatic et al., 2015. Multilocus phylogeography of the common midwife toad Alytes habitats and conservation of highly endangered species: the case of Ambystoma obstetricans (Anura, Alytidae), contrasting patterns of lineage diversiﬁcation mexicanum (Caudata, Ambystomatidae). Ann Zool Fenn 47:223–38. and genetic structure in the Iberian refugium. Mol Phylogenet Evol 93:363–79. Rogers AR, Harpending HC, 1992. Population growth makes waves in the dis- Hall TA, 1999. BioEdit: a user-friendly biological sequence alignment editor tribution of pairwise genetic differences. Mol Biol Evol 9:552–69. and analysis program for Windows 95/98/NT. Nucleic Acids Symp Ser Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A et al., 2012. (Oxf) 41:95–8. MrBayes 3.2: efﬁcient Bayesian phylogenetic inference and model choice Hampe A, Petit RJ, 2005. Conserving biodiversity under climate change: the across a large model space. Syst Biol 61:539–42. rear edge matters. Ecol Lett 8:461–7. Rozas J, Sanchez-Del Barrio JC, Messeguer X, Rozas R, 2003. DnaSP, DNA Hewitt GM, 2011. Mediterranean peninsulas: the evolution of hotspots. In: Zachos polymorphism analyses by the coalescent and other methods. FE, Habel JC, editors. Biodiversity Hotspots. Heidelberg: Springer, 123–47. Bioinformatics 19:2496–7. rd Huelsenbeck JP, Ronquist F, 2001. MRBAYES: Bayesian inference of phyl- Sambrook J, Russell D, 2001. Molecular Cloning: A Laboratory Manual.3 ogeny. Bioinformatics 17:754–5. edn. New York: Cold Spring Harbor Laboratory Press. Hughes L, 2000. Biological consequences of global warming: is the signal al- San Mauro D, Gower DJ, Oommen OV, Wilkinson M, Zardoya R, 2004. ready apparent? Trends Ecol Evol 15:56–61. Phylogeny of caecilian amphibians (Gymnophiona) based on complete mito- Knowles LL, 2001. Did the Pleistocene glaciations promote divergence? Test chondrial genomes and nuclear RAG1. Mol Phylogenet Evol 33:413–27. of explicit refugial models in montane grasshoppers. Mol Ecol 10:691–701. Sgro ` CM, Lowe AJ, Hoffman AA, 2011. Building evolutionary resilience for Lanfear R, Calcott B, Ho SYW, Guindon S, 2012. PartitionFinder: combined conserving biodiversity under climate change. Evol Appl 4:326–37. selection of partitioning schemes and substitution models for phylogenetic Swets J, 1988. Measuring the accuracy of diagnostic systems. Science 240:1285–93. analyses. Mol Biol Evol 29:1695–701. Tajima F, 1989. Statistical method for testing the neutral mutation hypothesis. Lataste MF, 1879. Pelonectes boscai. In: Tourneville A, editor. Description Genetics 123:585–95. d’une Nouvelle Espece de Batracien Urode ` le d’Espagne (Pelonectes boscai Teixeira J, 2008. Rana iberica. In: Loureiro A, Ferrand N, Carretero MA, Lataste). Bull Soc Zool Fr 4:69–87. Paulo OS, editors. Atlas dos Anfı´bios e Re´pteis de Portugal. Lisboa: Laurenti JN, 1768. Specimen Medicum, Exhibens Synopsin Reptilium Instituto da Conservac¸~ ao da Natureza, 124–25. Emendatum cum Experimentis Circa Venena et Antidota Reptilium Teixeira J, Arntzen JW, 2002. Potential impact of climate warming on the dis- Austriacorum. Wien: Joan Thom Nob de Trattnem. tribution of the golden-striped salamander Chioglossa lusitanica in Iberian Lemey P, Rambaut A, Welch JJ, Suchard MA, 2010. Phylogeography takes a Peninsula. Biodivers Conserv 11:2167–76. relaxed random walk in continuous space and time. Mol Biol Evol 27:1877–85. Teixeira J, Martı´nez-Solano I, Buckley D, Tarroso P, Garcı´a-Parı´s M et al., Lesica P, Allendorf FW, 1995. When are peripheral populations valuable for 2015. Genealogy of the nuclear b-ﬁbrinogen intron 7 in Lissotriton boscai conservation? Conserv Biol 9:753–60. (Caudata, Salamandridae): concordance with mtDNA and implications for Martı´nez-Freirı´a F, Velo-Anto ´ n G, Brito JC, 2015. Trapped by climate: inter- phylogeography and speciation. Contr Zool 84:193–215. glacial refuge and recent population expansion in the endemic Iberian adder Tejedo M, Bosch J, Martinez Solano I, Salvador A, Garcı´a Parı´s M et al., Vipera seoanei. Divers Distrib 21:331–44. 2009. Rana iberica. The IUCN Red List of Threatened Species 2009: Martı´nez-Solano I, Rey I, Garcı´a-Parı´s M, 2005. The impact of historical and e.T58622A86641528. [cited 2017 August 2] Available from: http://dx.doi. recent factors on genetic variability in a mountain frog: the case of Rana org/10.2305/IUCN.UK.2009.RLTS.T58622A11813789.en. iberica (Anura: ranidae). Anim Conserv 8:1–11. Thomassen HA, Fuller T, Buermann W, Mila ´ B, Kieswetter CM et al., 2011. Martı´nez-Solano I, Teixeira J, Buckley D, Garcı´a-Parı´s M, 2006. Mt-DNA Mapping evolutionary process: a multi-taxa approach to conservation pri- phylogeography of Lissotriton boscai (Caudata, Salamandridae): evidence oritization. Evol Appl 4:397–413. for old, multiple refugia in an Iberian endemic. Mol Ecol 15:3375–88. Vandergast AG, Bohonak AJ, Hathaway SA, Boys J, Fisher RN, 2008. Are Me ´ dail F, Que ´ zel P, 1997. Hot-spot analysis for conservation of plant bio- hotspots of evolutionary potential adequately protected in southern diversity in the Mediterranean Basin. Ann Mo Bot Gard 84:112–27. California? Biol Conserv 141:1648–64. Meijide MW, 1985. Localidades nuevas o poco conocidas de anﬁbios y reptiles Vences M, Hauswaldt JS, Steinfartz S, Rupp O, Goesmann A et al., 2013. de la Espana ~ continental. Donana ~ Acta Vert 12:318–23. Radically different phylogeographies and patterns of genetic variation in Miller MA, Pfeiffer W, Schwartz T, 2010. Creating the CIPRES Science two European brown frogs, genus Rana. Mol Phylogenet Evol 68:657–70. Gateway for inference of large phylogenetic trees. In Proceedings of the Vences M, Sarasola-Puente V, Sanchez E, Amat F, Hauswaldt JS, 2017. 2010 Gateway Computing Environments Workshop (GCE). New Orleans: Diversity and distribution of deep mitochondrial lineages of the common Institute of Electrical and Electronics Engineers, 1–8. frog, Rana temporaria, in northern Spain. Salamandra 53:25–33. Moritz C, 2002. Strategies to protect biological diversity and the evolutionary Warren DL, Glor RE, Turelli M, 2010. ENMTOOLS: a toolbox for compara- processes that sustain it. Syst Biol 51:238–54. tive studies of environmental niche models. Ecography 33:607–11. Nieto Feliner G, 2011. Southern European glacial refugia: a tale of tales. Waters JM, 2011. Competitive exclusion: phylogeography’s ‘elephant in the Taxon 60:365–72. room’? Mol Ecol 20:4388–94. Parmesan C, Yohe G, 2003. A globally coherent ﬁngerprint of climate change Yuan ZY, Zhou WW, Chen X, Poyarkov NA, Chen HM et al., 2016. impacts across natural systems. Nature 421:37–42. Spatiotemporal diversiﬁcation of the true frogs (Genus Rana): a historical frame- Paulo OS, Dias C, Bruford MW, Jordan WC, Nichols RA, 2001. The persistence work for a widely studied group of model organisms. Syst Biol 65:824–42. of Pliocene populations through the Pleistocene climatic cycles: evidence from Zink RM, Barrowclough GH, 2008. Mitochondrial DNA under siege in avian the phylogeography of an Iberian lizard. Proc R Soc B 268:1625–30. phylogeography. Mol Ecol 17:2107–21. Petit R, Aguinagalde I, de Beaulieu JL, Bittkau C, Brewer S et al., 2003. Zwickl D, 2006. Genetic Algorithm Approaches for the Phylogenetic Analysis Glacial refugia: hotspots but not melting pots of genetic diversity. Science of Large Biological Sequence Datasets under the Maximum Likelihood 300:1563–5. Criterion. Austin: School of Biological Sciences, University of Texas. Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zoy010/4827690 by guest on 13 July 2018
Current Zoology – Oxford University Press
Published: Jan 27, 2018
It’s your single place to instantly
discover and read the research
that matters to you.
Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.
All for just $49/month
Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly
Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.
Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.
Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.
All the latest content is available, no embargo periods.
“Hi guys, I cannot tell you how much I love this resource. Incredible. I really believe you've hit the nail on the head with this site in regards to solving the research-purchase issue.”Daniel C.
“Whoa! It’s like Spotify but for academic articles.”@Phil_Robichaud
“I must say, @deepdyve is a fabulous solution to the independent researcher's problem of #access to #information.”@deepthiw
“My last article couldn't be possible without the platform @deepdyve that makes journal papers cheaper.”@JoseServera