What triggers the rising of an intraspecific biodiversity hotspot? Hints from the agile frog

Hotspots of genetic diversity are regions of utmost importance for species survival and conservation, and their intimate link with the geographic location of glacial refugia has been well established. Nonetheless, the microevolutionary processes underlying the generation of hotspots in such regions have only recently become a fervent field of research. We investigated the phylogeographic and population genetic structure of the agile frog, Rana dalmatina, within its putative refugium in peninsular Italy. We found this region to harbour far more diversity, phylogeographic structure, and lineages of ancient origin than that by the rest of the species' range in Europe. This pattern appeared to be well explained by climate-driven microevolutionary processes that occurred during both glacial and interglacial epochs. Therefore, the inferred evolutionary history of R. dalmatina in Italy supports a view of glacial refugia as ‘factories' rather than as repositories of genetic diversity, with significant implications for conservation strategies for hotspots.

The agile frog, Rana dalmatina Fitzinger in Bonaparte, 1839, is a brown frog that is widely distributed in western, central, and southern Europe, where it inhabits a variety of landscapes, mostly within mixed deciduous forests from sea level to 2000 m a.s.l. R. dalmatina spends most of its life on land, breeding in a variety of natural and anthropogenic water bodies, including flooded meadows, ponds, lakes, temporary pools, and slow-running waters 24 . Although this species is not considered to be under imminent conservation risk, population declines and disappearances have recently been observed 25 . In a recent survey of the genetic variation throughout the species' range, R. dalmatina populations have appeared substantially homogeneous, with the single exception of a strongly differentiated population in southern Italy, which was considered as evidence for a glacial refugium of this species being located within this area 26 .
In this study, we further investigate the population genetic structure and evolutionary history of R. dalmatina, focusing on its putative refugial range in peninsular Italy. This geographic area has recently emerged as a hotspot of genetic diversity and as a site of multiple refugia for many temperate animal species, including amphibians 6,21,[27][28][29][30][31] . Here we adopted a dense sampling strategy within this area, with the following aims: i) to assess the phylogeographic and population genetic structure, and the historical demographic trends of R. dalmatina within its prospected glacial refugium, ii) to investigate the microevolutionary processes that have moulded its genetic diversity within this area; and iii) to evaluate the main implications of our results for the conservation of this species.

Results
Mitochondrial DNA. The final dataset included 741 bp of the mitochondrial cytochrome oxidase subunit I gene (COI; 46 variable positions, 37 parsimony informative) and 653 bp for the mitochondrial cytochrome B gene (CYTB; 55 variable positions, 42 parsimony informative), for all the 244 individuals analysed. Fortynine composite haplotypes were found (see Table 1), showing a sequence divergence (p-distance) ranging between 0.1% and 4.5%.
The statistical parsimony analysis did not connect all haplotypes into a single phylogenetic network ( Figure 1A). Instead, 2 main groups of haplotypes were observed (haplogroups N and S) having an ML-corrected net sequence divergence of 4.0% (0.5% S.E.), each showing a discrete geographic distribution. Haplogroup N was found in north-central and northern Italy, Croatia, and Slovenia, (samples [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19], whereas haplogroup S was distributed in southern Italy (samples . Both of these haplogroups were further divided into 2 sub-haplogroups, the distribution of which is presented in Figure 1C. While haplogroups N and S were never found in syntopy, the respective sub-haplogroups (NI-NII and SI-SII) showed overlapping distributions mostly limited to geographically intermediate populations. The haplotype genealogy yielded by HaplotypeViewer showed the same groupings as the statistical parsimony networks, with minor topological differences that did not affect the general pattern ( Figure 1B). Within this genealogy, haplogroups N and S were separated by 49 mutational steps.
Estimates of genetic diversity for the 4 haplogroups are shown in Table 2. Particularly low values of h and p were observed for the northernmost haplogroup (NI), whereas high values of both indices were found for the southernmost haplogroup (SII), and to a lesser extent, for SI.
The analysis of the divergence time dating with BEAST yielded effective sample size (ESS) values that largely exceeded 200 for all parameters of interest and full convergence for parameter estimates between runs. The divergence time between haplogroups N and S was estimated at 2.109 My (95% HPD: 1.433-2.929 My). On the other hand, the divergence between haplogroups NI-NII and SI-SII was estimated to have occurred at 0.381 My (95% HPD: 0.172-0.594 My) and 0.346 My (95% HPD: 0.197-0.519 My), respectively.  Table 1. The map was drawn using the software Canvas 11 (ACD Systems of America, Inc.). www.nature.com/scientificreports The historical demographic trends for all the main haplogroups were investigated by means of a bayesian skyline plot (BSP) analysis, except for NII because of its low sample size (n 5 10). The BSPs estimated for the remaining haplogroups are shown in Figure 2. The time to the most recent common ancestors (TMRCAs) was estimated as 23 000 years (95% HPD: 8300-45 600), 84 500 years (95% HPD: 27 700-158 800), and 155 000 years (95% HPD: 63 200-266 400) for haplogroups NI, SI, and SII, respectively. For haplogroup NI, a marked demographic expansion was inferred that began shortly after the TMRCA and lasted until recent times. In contrast, the BSPs for haplogroups SI and SII showed no marked demographic changes; the median estimates of effective population size for these groups never fell outside the 95% HPDs estimated along the entire time axis.
The estimated values of the neutrality test statistics are given in Table 2. In agreement with the BSP analysis, large negative values of F S and D and a small positive value of R 2 (all suggesting recent and marked population growth) were only observed for haplogroup NI (all P , 0.05). On the other hand, none of the estimated values of the statistics for the haplogroups NII, SI, and SII was significant.
Microsatellites. A total of 170 individuals from 38 locations were genotyped at 6 microsatellite loci (see Table 1). Among these individuals, none was scored at less than 5 loci, and the overall proportion of missing data in the final dataset was 0.2%.
Evidence of Hardy-Weinberg or linkage disequilibria was not found across samples or loci. For all populations, the observed heterozygosity ranged from 0.69 (H1 locus) to 0.89 (G12 locus), the expected heterozygosity ranged from 0.84 (H1 locus) to 0.96 (E8 locus), and the number of alleles observed at each locus was as follows: F5 5 19, G12 5 28, G11 5 16, E8 5 38, H1 5 16, and B5 Consequently, the total number of alleles found within the study area (n 5 139) was 12.9% higher than that previously reported for the same loci throughout the species range 26 .
The analysis of population genetic structure by using TESS indicated K 5 4 as the best grouping option for our data, since only a minor decrease in DIC values was observed at higher K values ( Figure 3A). As previously observed for the mtDNA data, 2 groups were found in south-central and southern Italy (samples 21-40; see Figure 3B). However, unlike the mtDNA data, the 2 population groups were largely admixed in most of the areas (Figures 3B and 3C), with some exceptions in marginal samples. Two groups were observed in the north-central and northern portion of the study area. However, the geographic distribution of the 2 groups did not conform to the geographic pattern observed in the mtDNA data. Indeed, the area of admixture between the 2 groups was not located in central Italy (samples 17-20), but on the north-eastern side of the Italian Peninsula (samples 8-10, 14). Finally, some instances of unidir-ectional gene exchange from south to north were observed among samples located in central Italy (samples [17][18][19][20], which otherwise belonged to the north-central group ( Figure 3C). Finally, the mean expected heterozygosity estimated for the 4 groups was 0.90 (0.01 SE), 0.92 (0.01 SE), 0.88 (0.03 SE), 0.88 (0.02 SE) for the southern, south-central, north-central and northern groups respectively.

Discussion
Our analysis of genetic variation in the agile frog in peninsular Italy allowed us to identify this area as not only a site of the most divergent lineage but also as the major hotspot of genetic diversity within this species. Indeed, populations from this area harbour far more genetic variation than that found in the central, eastern, and southern European populations combined 26 . Further, most of the genetic variation was accounted for by differences between 4 genetically divergent and geographically well-defined lineages, suggesting a long-lasting and multi-episode evolutionary history for R. dalmatina within this area. In turn, the occurrence of 4 lineages with relatively ancient divergence clearly indicates that peninsular Italy did not act as a single refugium for the species, but rather as a mosaic of refugia, at least during the last Pleistocene glacial-interglacial cycles.
The most ancient event inferred in the history of R. dalmatina populations was the divergence between the N and S mtDNA haplogroups, which has been roughly dated at 2.1 million years (MY). The precise dating of this divergence was prevented by the large 95% CI for the time estimate (1.4-2.9 My), our use of an externally calibrated substitution rate, along with the many uncertainties of the molecular clock 32,33 . Nevertheless, it is worth noting that the phylogeographic discontinuity between the N and S clades is located close to the northern edge of a mid-peninsular region that was repeatedly subjected to marine-flooding during the Plio-Pleistocene, and where phylogeographic discontinuities have been repeatedly reported 6,28,29,34,35 . Thus, it is plausible that a marine transgression that occurred in the early Pleistocene might have primed the observed divergence. This divergence is conspicuous (sequence divergence was 4.0% in the mtDNA data), syntopy was not observed between the 2 mtDNA clades, and the admixture of nuclear genomes appeared both geographically and quantitatively limited (see Figure 3C). Therefore, further study is required to locate the secondary contact zone between both clades and to investigate what microevolutionary processes occurred following the formation of this zone.
The transition zone between the NI and NII lineages in the mtDNA data ( Figure 2) was not reflected by the transition zone between the north-central and the northern population groups in the microsatellite data ( Figure 3). Although not unusual, it may be difficult to reconcile the discordance between nuclear and mitochondrial data to infer the history of populations, given the diversity of microevolutionary processes that are potentially involved 36 . However, while the northern Apennine is a well known suture zone between both intra-and interspecific lineages 27,37-40 , a contact zone in the area close to our sites 17-19 would not have any analogues in the phylogeographic literature of co-distributed temperate species. In addition, no environmental features (present or past) were present that could explain the observed pattern. This observation led us to speculate that, after the secondary contact between the 2 lineages in an area close to the northern Apennines, the NI mitochondrial lineage incurred massive introgression into the NII range, resulting in the former replacing the latter throughout most of central Italy. On the other hand, unfavourable climatic conditions along the central Apennines chain during the glacial stages 41,42 might provide a viable explanation for the divergence between the 2 lineages. During these stages, the NI and NII lineages might have persisted in temperate lowland areas along the eastern Adriatic and the central Tyrrhenian coasts, respectively 42 .
The phylogeographic discontinuity observed between the SI and SII mtDNA clades was located in central Calabria. Although a more admixed picture emerged at the nuclear microsatellites, central Calabria was also the main area of transition between the southcentral and the southern groups of populations (marked as purple and green, respectively, in Figure 3). Interestingly, several phylogeographic breaks and secondary contact zones were clustered in this area 6,21,[27][28][29][30] , which was cyclically inundated by marine transgressions until the late-Middle/Late Pleistocene 43 . Therefore, it seems plausible that these transgressions also primed the divergence between the 2 R. dalmatina groups found in south-central and southern Italy.
Despite the relative geographic contiguity of the alleged past ranges of the divergent R. dalmatina lineages found, the lineages showed substantially different demographic trends during the Late Pleistocene. The star-like structure of the haplotype genealogy for lineage NI, its BSP reconstruction, and the estimated values of the neutrality test statistics jointly indicated that, after a marked demographic bottleneck occurred at the last glacial maximum (LGM; as revealed by the estimated TMRCA around 23 000 years ago), NI underwent a substantial and abrupt demographic expansion. Since NI corresponds to the widespread lineage found by Vences and colleagues 26 (2013; data not shown) throughout most of R. dalmatina's range in Europe, and in light of its Middle-Pleistocene divergence from NII (0.38 [0.17-0.59] My), we suggest that lineage NI re-colonized most of its current range after the LGM from a glacial refugium located somewhere between the western Balkans and the wide coastal plain that occupied what is now the Adriatic Sea 42,44 . In contrast, lineages SI and SII showed the genetic imprints of demographic stability, lasting for most of the Late Pleistocene (SI), or before (SII). Similar trends of prolonged demographic stability have also been hypothesized for intraspecific lineages of other temperate species studied to date in southern Italy 6,29 . Thus, contrary to what is frequently observed farther north, the most striking effect of glacialinterglacial cycles on the history of populations from south-central and southern Italy seems to be driven not by climate forcing on range/demographic size variations, but by glacio-eustatic oscillations in the sea level and the resulting allopatric fragmentation of populations imposed by the intervening seaways. Consequently, this area could be viewed as a 'glacial refugium' ('southern refugium' sensu 45 ) from the worsening climate during the ice ages and as a set of 'interglacial refugia' from the rising sea levels during interglacial epochs.
Among wide-ranging species in Western Europe, geographic patterns of genetic diversity as observed for R. dalmatina (more diversity to the south and less diversity to the north, i.e. the 'southern-richness, northern-purity' pattern 16 ) are not unusual 4,7,46 . In most of these species (including R. dalmatina 26 ), the 'southern-richness' component of this pattern, i.e. the existence and location of hotspots, has been interpreted as the outcome of the long-term persistence of ancestral populations within southern peninsulas 2,3 . In contrast with this view for wide-ranging species and in parallel with recent observations in species endemic to these peninsulas 5,6,8,47 , our results emphasize the key role played by the 'diversity-enhancing' effects of the cycles of allopatric divergence within these areas. Since a sampling scheme specifically focused on the putative refugial ranges is crucial for disentangling the above scenarios, we anticipate that as more studies adopt such sampling schemes, our view of the diversity of processes involved in the formation of hotspots of intraspecific biodiversity will continue to change.
Finally, our results have at least 3 major implications for the conservation of R. dalmatina populations, both within our study area and at the whole species' range. First, in order to preserve the species' resources for coping with future natural and anthropogenic environmental changes, peninsular Italy should receive special attention since most of R. dalmatina's genetic diversity is apportioned there, and in light of the historical role played by this area in the species' survival. This appears especially important when considering the declining trend recently observed for R. dalmatina populations from several portions of its range, including the Italian peninsula 25 . Second, since the deep genetic structure of the populations indicates the occurrence of at least two ESUs (sensu 23 ) within this area, specific conservation efforts should be developed to optimize the amount of genetic variation preserved. Finally, a study of the genetic structure of the presumed secondary contact zone between the N and S lineages in central Italy should be carried out to define the patterns of past and present interaction among them, their most appropriate taxonomic arrangement, and the best strategy for their conservation.

Methods
Sampling and laboratory procedures. We collected tissue samples from 244 individuals of R. dalmatina from 40 localities, 33 of which were located along the Italian peninsula south of the Alps, 3 from Slovenia, and 4 from Croatia. The geographical position and number of individuals sampled at each location are given in Table 1 and Figure 1. Tissue samples were obtained through toe-clipping, and all individuals were released at the place of collection immediately after tissue removal. Samples were transported to the laboratory in liquid nitrogen and were stored at 280uC until subsequent analyses. We extracted genomic DNA by using the standard phenol-chloroform protocol 48 . Two partial fragments of the mitochondrial DNA (mtDNA) genes encoding CYTB and COI were amplified and sequenced. The same PCR protocol was used for both mtDNA fragments. Amplifications were performed in a final volume of 25 ml containing MgCl 2 (2.5 mM), reaction buffer (13; Promega), 4 dNTPs (0.2 mM each), 2 primers (0.2 mM each), Taq polymerase (2-unit; Go-Taq Promega), and 2 ml of template DNA. PCR cycles were as follow: 94uC for 5 min followed by 36 cycles of 94uC for 45 s, annealing temperature (50uC for CYTB and 52uC for COI) for 45 s, 72uC for 90 s, and a final step at 72uC for 10 min. The primers used for amplification and sequencing of the CYTB fragment were 494 49 and MVZ16 50 ; primers for the COI fragment were Amp-P3r and Amp-P3f 51 . Sequencing reactions were carried out by Macrogen Inc. (www.macrogen.com) by using an ABI PRISM 377 DNA sequencer (PE Applied Biosystems) following the ABI PRISM BigDye Terminator Cycle Sequencing protocol. All individuals analysed were double sequenced.
Fragments were sized using the GeneScan 400HD internal size standard, and alleles were scored using PEAKSCANNER 1.0 (Applied Biosystems).
We used the method implemented in the software MICROCHECKER 53 to check for the presence of null alleles and large allele dropout; however, there was no evidence of such problems.
Mitochondrial DNA data analysis. Sequence electropherograms were checked by eye, by using ChromasPro (Technelysium Ltd.), and multiple sequence alignments were carried out using ClustalX 54 . Sequences of the CYTB and COI fragments were concatenated using Concatenator 1.1.0 55 , after a partition homogeneity test 56 conducted in Paup*4.0b10 did not reject the hypothesis of homogeneity in the phylogenetic signal carried by both fragments. The optimal model of sequence evolution for both fragments separately and the concatenated dataset were selected among 88 possible models using Jmodeltest 0.1.1 57 . The best-fit models selected using the Akaike information criterion were TrN 1 C 58 (gamma distribution shape parameter 5 0.032) for COI, TrN 1 I (proportion of invariable sites 5 0.073) for CYTB, and TrN 1 I (proportion of invariable sites 5 0.824) for the concatenated dataset. Therefore, the latter was used in all downstream analyses on the unpartitioned dataset.
Sequences were collapsed into unique haplotypes and their nucleotide variation was evaluated with the software DnaSP 5 59 . The same software was used to estimate the haplotype (h) and nucleotide (p) diversity 60 among major groups of haplotypes as identified by phylogenetic analyses; Mega 5 61 was then used to compute the net sequence divergence among haplotype groups.
The genealogical relationships among the identified haplotypes were investigated using 2 approaches. First, we estimated phylogenetic networks using the statistical parsimony procedure of Templeton 62 , implemented in the software TCS 1.21 63 using the 95% limit for a parsimonious connection. Second, to account for recent concerns about the accuracy of network-building methods 64 , we further investigated the genealogical relationships among haplotypes by using the approach suggested by 64 . A phylogenetic tree topology was estimated using the maximum likelihood algorithm as implemented in PhyML 65 . We used the default setting in PhyML for all parameters except the substitution model (TrN 1 I) and the type of tree improvement (SPR & NNI 66 ). The estimated tree topology was then converted into a haplotype genealogy by using the software HaplotypeViewer 65 .
Estimates of the divergence time between identified haplogroups were obtained using BEAST 1.7.5 67 using the Bayesian Skyline as coalescent tree prior, which imposes minimal assumptions on the data 67 . Preliminary analyses were run using a relaxed molecular clock model, assuming an uncorrelated lognormal distribution. However, since the standard deviation of the uncorrelated lognormal relaxed clock (the parameter ucld.stdev) was close to zero, subsequent analyses were run enforcing a strict clock model. Since we had neither internal calibration points nor previous estimates of the substitution rates for the 2 mtDNA fragments studied in R. dalmatina, this and subsequent demographic analyses were carried out using a uniform prior for the substitution rate, bounded at 0.012 and 0.014 substitutions/site/My, to reflect the range of substitution rates previously estimated for CYTB and COI fragments in other European brown frogs 68,69 . The final analysis was run in triplicate to check for consistency among runs, each with a Markov chain Monte Carlo (MCMC) run for 20 million generations and sampled every 2000 generations. Further, satisfactorily high ESS values for all parameters in the analysis were obtained using these settings. Convergence, stationarity, and the appropriate number of steps to be discarded as burn-in were assessed using TRACER 1.6.
A hierarchical analysis of molecular variance (AMOVA 70 ) was carried out by using ARLEQUIN 3.1 71 , in order to estimate the amount of variation attributable to differences among population groups, among populations within groups, and within populations. This analysis was run using the TrN model, with population groups defined according to the geographical distribution of the main haplogroups.
The historical demography of the main haplogroups was inferred using the Bayesian skyline plot method (BSP) implemented in BEAST 72,73 . Five preliminary replicates of the BSP analysis were carried out to fine-tune the parameters of the MCMC using the auto-optimization option. Then, we ran the final analyses in triplicates for each haplogroup to check for consistency. The number of groups (m) was set to 10 (but setting different values in preliminary analyses did not affect the resulting demographic trends); the MCMCs were run for 30 million generations and sampled every 1000 generations, with the first 10% discarded as burn-in. BSP analyses were run using the TrN 1 I, a strict clock model, and a uniform prior for the substitution rate, bounded at 0.012 and 0.014 substitutions/site/My. The results were analysed and summarised using TRACER 1.6. Finally, the occurrence of past demographic changes was evaluated by computing the neutrality test statistics F S 74 , R 2 75 and D 76 , by using DNaSP. The significance of these statistics was assessed using 1000 coalescent simulations.
Microsatellites data analysis. We tested for deviations from the expected Hardy-Weinberg and linkage equilibria using FSTAT 2.9.3. The software GENETIX was used to obtain estimates of genetic diversity based on the observed and the expected heterozygosity and the number of alleles per locus.
The extent of population genetic structure across the study area was investigated by the Bayesian clustering algorithm implemented in TESS 2.3 77 . This method is advantageous in that it takes the spatial distribution of sampled individuals into account, and has been shown to outperform similar methods, particularly in cases of shallow population structures 77 . The TESS analysis was carried out by modelling admixture using the conditional autoregressive (CAR) model. A preliminary analysis with 20 000 runs (of which the first 5 000 were discarded as burn-in) and 10 replicates for each K value (i.e. the number of clusters) between 2 and 20, was carried to assess model performance, the need to fine-tune the model parameters, and to reduce the range of plausible K values. The final analysis contained 100 replicates with K 5 2 to 7, each with 50 000 runs (of which the first 20 000 were discarded as burn-in). The spatial interaction parameter was set to 0.6 (i.e. the default value) and the option to update this parameter was activated. The model that best fit the data was selected using the deviance information criterion (DIC). DIC values were averaged over the 100 replicates for each K value, and the most likely K value was selected as the one at which the average DIC reached a plateau. The estimated admixture proportions of the 10 runs with the lowest DIC values were averaged using CLUMPP 1.1.2 78 . Finally, to obtain an easily interpretable representation of the estimated admixture proportions, we generated posterior predictive maps by the spatial interpolation procedure available in TESS, and with the maps being visualised by TESS Ad-Mixer 79 . Finally, we used the posterior estimate of the allele frequencies at each locus in each cluster to compute the mean expected heterozygosity of each cluster.