Genetic attributes of midwife toad (Alytes obstetricans) populations do not correlate with degree of species decline

Genetic diversity is crucial for long-term population persistence. Population loss and subsequent reduction in migration rate among the most important processes that are expected to lead to a reduction in genetic diversity and an increase in genetic differentiation. While the theory behind this is well-developed, empirical evidence from wild populations is inconsistent. Using microsatellite markers, we compared the genetic structure of populations of an amphibian species, the midwife toad (Alytes obstetricans), in four Swiss regions where the species has suffered variable levels of subpopulation extirpation. We also quantified the effects of several geographic factors on genetic structure and used a model selection approach to ascertain which of the variables were important for explaining genetic variation. Although subpopulation pairwise FST-values were highly significant even over small geographic scales, neither any of the geographic variables nor loss of subpopulations were important factors for predicting spatial genetic structure. The absence of a signature of subpopulation loss on genetic differentiation may suggest that midwife toad subpopulations function as relatively independent units.


Introduction
The maintenance of population genetic diversity is of crucial importance for long-term population persistence (Frankham 2005;Evans and Sheldon 2008) Reducing the number of migrants into a recipient population is a proven mechanism for reducing its effective population size and causing a reduction in measurable genetic diversity through genetic drift (Frankham 1995(Frankham , 2005. Immigration rate is expected to scale with the availability of source populations, so population loss should have a direct and measurable effect on the genetic variability and genetic similarity of the remaining populations. Given sufficient time (i.e. 20-50 generations; Anderson et al. 2010), in areas where population loss is significant, among-population differentiation should be greater and within-population genetic variability lower than in areas where local extirpation has been less frequent (Lande et al. 1998;Frankham 2005). While the theory is well-developed (Gilpin 1991), empirical evidence in support of theory is inconsistent (Keyghobadi 2007). Populations of taxa that were classified as threatened by the International Union for Conservation of Nature (IUCN) and have experienced significant population loss exhibit lower population genetic diversity than do unthreatened sister taxa (Spielman et al. 2004;Evans and Sheldon 2008). However, other studies have found no such relationship (Gibbs 2001;Keyghobadi 2007). Determining why there are inconsistencies is hampered by the sampling design of many studies of population genetics and population loss. Most are unreplicated, conducted at a single site or within a single geographic area (Matocq and Villablanca 2001) and hence may be unrepresentative of the general relationship between rapid decline and loss of genetic diversity. A better approach would be comparative, measuring the consistency of the relationships between population loss and genetic variability in replicated regions.
The global decline of amphibians has resulted in the loss of amphibian populations over a matter of decades, but declines are not homogeneously distributed (Houlahan et al. 2000;Stuart et al. 2004). The amphibian extinction crisis is certainly one of the most topical conservation issues of our generation, but it could also be viewed as an opportunity for investigating the links between decline and population genetic structure. Previous studies have illustrated a relationship between population loss and genetic diversity on a limited geographic scale Beebee 1997, 1998;Beebee 2005). We know of no study that has made an explicit attempt to relate population loss to its effects on gene flow and population genetic variability in a replicated fashion.
Adopting the classic definitions of Wright (1965) from hereon, a region across which amphibian breeding ponds are distributed is equivalent to a population and each pond within a population to a subpopulation. A region refers to a geographic area harboring a number of subpopulations, but not necessarily representing a closed population without genetic exchange with other such populations. A suitable amphibian species to study the link between decline and genetic structure would have the following traits: (1) discrete subpopulations; (2) data available on subpopulation persistence across comparable populations; and (3) variable rates of subpopulation extinction across comparable populations. Such a model exists in the common midwife toad, Alytes obstetricans, in Switzerland, where the distribution of the species is well known (Grossenbacher 1988;Borgula and Zumbach 2003). A. obstetricans exhibits a complex life history with an aquatic larval stage and a terrestrial adult stage, so breeding ponds can be treated as defined subpopulations.
Alytes obstetricans has suffered strong declines in Switzerland: The species was first red-listed in 1982 (Hotz and Broggi 1982) because population declines were already observed in Switzerland in the 1960s (Escher 1972). The most recent update of the Swiss amphibian red list showed that since the mid-1980s~50% of subpopulations have been extirpated (Schmidt and Zumbach 2005). Therefore, A. obstetricans is categorized as "endangered" in the most recent Swiss Red List (Schmidt and Zumbach 2005). Very few colonizations of unoccupied ponds have been observed and Swiss A. obstetricans subpopulations are small (Borgula and Zumbach 2003;Schmidt and Zumbach 2005;Tobler et al. 2012). Given that A. obstetricans has a generation time of 1-2 years (B€ oll et al. 2012), >20 generations have passed since the onset of declines, sufficient time for subpopulation loss to have affected genetic structure (Anderson et al. 2010). Rate of subpopulation loss has varied among populations and in most cases subpopulation extirpation cannot be attributed to habitat loss (Borgula and Zumbach 2003;Schmidt and Zumbach 2005).
In this study we take advantage of existing knowledge regarding the spatial distribution of recent extirpations of A. obstetricans (Fig. 1) subpopulations and ask whether population genetic structure and diversity measured in subpopulations and populations varies and if this variation can be attributed to among-population variation in rates of recent subpopulation loss. To do this, we sampled four regions where the common midwife toad is found and where quantitative evidence of variation in subpopulation loss is available. We used microsatellite polymorphisms to measure within-subpopulation and within-population genetic variability, as well as gene flow among subpopulations. We distinguished between the  of subpopulation loss and region-specific characteristics on genetic diversity and population size by modeling whether other factors such as geographic connectivity, elevation, or location close to a stream contributed to the observed genetic patterns. We expected that the loss of subpopulations that occurred over a time frame of roughly 20 generations would be paralleled by decreasing population sizes and have caused an increase in genetic differentiation among subpopulations and a reduction in genetic diversity due to increased drift (Lande 1988). We predicted similar consequences due to the effect of weaker geographic connectivity (Cushman 2006). Because suitable habitat becomes scarcer towards the distribution limits of a species, higher elevation subpopulations should also be more isolated than lowland ones so we predicted similar effects as those resulting from range periphery on high-elevation subpopulations (Giordano et al. 2007). Thus, we expected reduced genetic diversity and heterozygosity, and increased subpopulation differentiation at higher elevation. Finally, because midwife toads sometimes occur in stream or inhabit ponds near streams (Barandun 2007), we also assessed the effect of streams on the genetic structure. Streams may act as corridors for gene flow either downstream (tadpoles) or upstream (adults; Morrissey and de Kerckhove 2009;Grant et al. 2010;Mullen et al. 2010), hence we expected weaker genetic differentiation among subpopulations along identical catchments (Table 1).

Study populations
We collected tissue samples for genetic analyses from 45 subpopulations clustered within four populations in Switzerland during spring and summer 2007: Baselland (BL, 15 subpopulations), Bern (BE, 11 subpopulations), Lucerne (LU, 11 subpopulations), and St. Gallen (SG, 8 subpopulations; Fig. 2). We chose these populations based on available knowledge on midwife toad subpopulation trends (Table 2) and the fact management practices have aimed to increase and improve the quantity and quality of breeding sites in all four populations. Therefore, there was no net loss of available breeding sites for A. obstetricans across our study system. BL is located in the Jura Mountains and supports a comparatively dense network of A. obstetricans subpopulations. The mean Euclidian distance between all known  Ryser et al. 2003). In BE, A. obstetricans is found predominantly in man-made ponds. Lucerne is a pre-alpine region where the distances between subpopulations are relatively large (mean Euclid-ian distance of 17.7 km among subpopulations; A. Borgula, pers. comm.; 13.1 km among study subpopulations). Relatively more subpopulations of A. obstetricans in LU have been lost (51 subpopulations reported in 1980 down to 23 reported in 2002) than in BL and BE. There are only three reports of colonization of new subpopulations (Borgula and Zumbach 2003). Breeding habitats for A. obstetricans in LU are man-made ponds located in prealpine and alpine meadows, and one stream relatively unchanged by human activity. BE and LU are adjacent but are separated from each other by prealpine topography and are part of different drainage basins. The colonization of prealpine areas after the last glaciation followed the retreat of the ice starting from the lowlands (Hewitt 1996;Mart ınez-Solano et al. 2004), hence we consider LU and BE genetically distinct populations, and test this (see below). SG is the second region where A. obstetricans has experienced severe subpopulation loss and is also a prealpine region. Since 1980, A. obstetricans is absent from 68 of 118 recorded subpopulations and another 23 are potentially extirpated. No colonizations are reported (Barandun 2004). Breeding habitats in SG are streams and man-made ponds; mean distance among subpopulations in this region is 19.8 km (J. Barandun, pers. comm.; 11.1 km among study subpopulations). Both LU and SG are located on the periphery of the species distribution in Switzerland (Fig. 2).
In each population we sampled 8-15 ponds previously reported to harbor A. obstetricans subpopulations. In LU, we sampled two clusters of subpopulations of which the southernmost was located at the distribution border. Subsequent analyses confirmed that, although no known populations exist between the southern and northern cluster, these clusters could still be treated as one population (see below; Fig. 3). The measure of population differentiation we used is relatively insensitive to unsampled populations (Koen et al. 2013) and the populations that were not sampled were small and most likely do not contribute 4.3 km (0.6-9.0) 5.7 km (0.9-13.2) 13.1 km (1.2-25.9) 11.1 km (0.9-20.8) many migrants. Hence, we do not expect that the presence of unsampled populations caused bias in our estimates of population differentiation and isolation-bydistance. We caught tadpoles by dip-netting and collected tissue for DNA extraction by cutting off less than 3 mm from the tail tip. We sought to sample a minimum of 25 tadpoles per site. All tadpoles were released into their source ponds immediately thereafter. Standard hygiene protocols to avoid the spread of infectious diseases were followed during field work (Schmidt et al. 2009).

Microsatellite development
Microsatellite primers were commercially developed by ecogenics GmbH (Zurich, Switzerland) using A. obstetricans tissue samples from central Spain, France and Switzerland. An enriched library was made from the DNA of one Swiss A. obstetricans: size selected genomic DNA was ligated to SAULA/SAULB-linker (Armour et al. 1994) and enriched by magnetic bead selection with biotin-labeled (GT) 13 , (CT) 13 , (GATA) 7 , (GTAT) 7 , (ACAG) 7 , and (GCGT) 7 oligonucleotide repeats (Gautschi et al. 2000). Of 1893 recombinant colonies screened, 241 gave a positive signal after hybridization. Plasmids from 187 positive clones were sequenced and primers were designed for 29 microsatellite inserts, of which 21 were tested for polymorphism. We selected a set of 12 primers that exhibited clear and reliable amplification, polymorphism and no evidence of null alleles in preliminary tests for generating population genetics data (Table 3).

Microsatellite amplification
We extracted DNA from tail clips using the BioSprint 96 DNA Blood Kit (Qiagen, Hombrechtikon, Switzerland), following the protocol for tissue extraction. Polymerase    (Table 3). Each well contained 1.25 lL Multiplex PCR Kit (Qiagen), the respective primer volumes (Table 3) and 1 lL of template DNA. Forward primers were color-labeled to allow multiplexing. PCRs were carried out on a TC-412 Thermal Cycler (Barloworld Scientific, Staffordshire, UK) with polymerase activation at 95°C for 15 min, followed by 33 cycles (multiplex 1) or 30 cycles (multiplexes 2 and 3) of denaturing for 0.30 min at 94°C, annealing for 1.30 min at 52°C (multiplex 1) or 56°C (multiplexes 2 and 3) and extension for 1.00 min at 72°C, followed by a final extension for 30 min at 60°C. PCR products were run on an ABI 3730 Avant capillary sequencer (Applied Biosystems, Rotkreuz, Switzerland) with internal size standard Gene-Scan-500 LIZ; peaks were visually scored using GENEM-APPER 3.7 (Applied Biosystems 2004).

Statistical analysis
Microsatellite loci were tested for the presence of null alleles, stuttering and allelic dropout of larger alleles using MICROCHECKER (Van Oosterhout et al. 2004). Because of the potential presence of siblings in our dataset (Goldberg and Waits 2010), we tried to identify siblings using COLONY 3.1 (Wang 2004). We did not find within-population sibship structure, indicating there are no confounding effects of relatedness on our data (J. Wang, pers. comm.). We used ARLEQUIN 3.1 (Excoffier et al. 2005) to test for linkage disequilibrium and deviations from Hardy-Weinberg (HW) equilibrium and tested all microsatellites for evidence of locus-specific selection using the program FDIST (Beaumont and Nichols 1996).
Measures of genetic diversity per subpopulation and population (observed heterozygosity [H o ], expected heterozygosity [H e ], and the frequency of private alleles) were calculated using GENETIX 4.05.2 (Belkhir et al. 1996(Belkhir et al. -2004 and a sample-size corrected estimate of allelic richness (A) was obtained using FSTAT 2.9.3.2. To obtain an estimate of number of breeding individuals per population (N b ), we ran the programs ONeSAMP (Tallmon et al. 2007) with prior information between 2 and 100 for each population, and Colony (Wang 2004) assuming a polygamous breeding system for both males and females, and using the full likelihood model with medium precision and no prior information. As both analyses returned very similar values, we only report the results obtained using Colony.
To test whether the regions could be treated as distinct genetic units, we performed a Bayesian clustering analysis using the software STRUCTURE (Pritchard et al. 2000). To infer the number of clusters in our data, we used the DK method proposed by Evanno et al. (2005), modified using a stepwise approach according to Coulon et al. (2008). Specifically, we repeated the estimation of the number of clusters with the DK method on each of the K groups inferred in the previous step. We repeated this process until the number of clusters inferred was 1, or individuals from any single pond were split. Each simulation was run with a burn-in period of 50,000 iterations and a sampling period of 250,000 iterations. For every simulation step, we set the bounds of K from 1 to 15 and repeated the simulation for each K 10 times.
Because a proportion of young-of-the-year Alytes tadpoles hibernate as tadpoles and only metamorphose in the year after hatching (Thiesmeier 1992), we first calculated F ST values treating overwintered tadpoles from 2006 and young-of-the-year tadpoles from 2007 as separate units ("populations" in the terminology of the program) in all ponds where data on two cohorts was available. This analysis tested whether reproduction among cohorts was strongly skewed towards a few parental individuals (Savage et al. 2010). Because all F ST -values between cohorts were not significant (results not shown), we assumed that our measures of genetic diversity with all samples from a subpopulation pooled were representative. We next calculated global F ST -values per region and pairwise F ST -values of subpopulations within regions as a measure of genetic population structure among and within regions and tested the significance of pairwise F ST within regions with an exact test with 1000 permutations (Excoffier et al. 2005). We tested for isolation by distance (IBD) across and within populations using F ST /(1ÀF ST )transformed F ST values and log-transformed Euclidian distances (Excoffier et al. 2005). To determine at what degree genetic variation was distributed within and among subpopulations and populations, we used analysis of molecular variance (AMOVA) implemented in Arlequin 3.1 (Excoffier et al. 2005).
To test our hypotheses regarding the effects of severity of declines, geographic connectivity, elevation, or location near a stream on population size and differences in genetic structure we combined linear mixed effects models with an Akaike's information criterion (AICc) based inference framework (Burnham and Anderson 2002).
Severity of subpopulation loss was classified as "none," "moderate," or "strong" (Table 2). We extracted data on elevation and pairwise Euclidian distances among subpopulations within a region from the national map (1:25,000; Swisstopo 2003-2009). Geographic isolation within a population was defined as mean pairwise Euclidian distance of a subpopulation to all other subpopulations in the same region. We considered subpopulations as stream subpopulations if they were within 200 m of a stream or when the breeding site was a stream. We used A, H e , and the mean subpopulation F ST as response variables in separate linear mixed effects models; we did not use all variables of genetic diversity that are commonly used in population genetic studies because all measures of allelic diversity (number of private alleles, fixed loci, and A) were correlated (all r > 0.45), and the same was true for H e and H o (r = 0.89). We calculated a mean subpopulation F ST as the mean of all pairwise F ST values of a subpopulation to all other subpopulations within the same population.
Our independent variables, all defined above, were (1) severity of subpopulation loss, (2) elevation, (3) location close to a stream, and (4) geographic isolation within a population. We used an intercept-only model as the null model. We included one main effect per model only. Population was defined as a random effect to accommodate non-independence of the sub-populations in all models (Rhodes et al. 2009) We standardized elevation and isolation for the use in linear mixed effects models using a z-transformation. After fitting the models, we used spline correlograms to investigate auto-correlation in the data (Rhodes et al. 2009).
We ranked models based on AICc for small sample sizes (Burnham and Anderson 2002) and report parameter estimates for models within 2 DAICc units of the best model. Model assessment was also based on inspection of the parameter estimates and their confidence intervals; that is, only parameters of which the confidence interval did not include 0 were considered important (Burnham et al. 2011). All linear mixed models were fitted to the data using maximum likelihood in R 2.8.1. using the package lme4 (R Development Core Team 2008; Bates et al. 2011).

Genetic diversity
We did not detect null alleles, large allele dropout or stuttering at any microsatellite locus (all P > 0.05). Locus Al-yobs20 was monomorphic in all studied subpopulations and excluded from further analyses. We found no deviations from HW equilibrium after Bonferroni correction and very little and inconsistent linkage at few markers in few subpopulations. With the application of a conservative significance level of P = 0.001, 24 signatures of linkage were found in a total of 4730 comparisons (11 markers times 43 subpopulations); linked markers were found in 8 of the 43 subpopulations. This is further indication that the (possible) presence of siblings in the genotype data did not negatively affect our results (Rasmussen 1979). Four markers were identified as being potentially non-neutral. Alyobs23 was designated as a candidate for directional selection, while Alyobs3, Alyobs4, and Alyobs8 were identified as candidates for balancing selection (Beaumont and Nichols 1996). Excluding these markers from the calculation of our measures of genetic diversity estimates did not alter the model selection results (results not shown) hence we concluded that potential non-neutrality did not bias our results. We thus report the estimates of genetic diversity based on all 11 polymorphic markers. Average allelic richness (A) within subpopulations was low. Six loci, on average, were fixed in every subpopulation and we detected on average one private allele in every second subpopulation (22 private alleles in 45 populations; Table 4). Average H e was also low and the difference between H e and H o never exceeded 0.12 (Table 4). Mean breeding population size was low with 16.1 individuals, ranging from 8 to 31 breeders (Table 4).
With only one exception we could assign subpopulations to the correct population in the first step of the STRUCTURE analyses (Figure 3). This initial analysis only split the samples into three clusters, with two of the populations (BE and LU) grouped together. As well, subpopulation BL-Reig, located geographically in BL, clustered with SG. This wrong assignment is probably due to a lack of power with only~2 alleles per locus: Out of 24 different alleles at 11 markers in BL-Reig, seven are shared between BL-Reig and the other BL subpopulations, and only one is shared between BL-Reig and the other SG subpopulations, but not with the rest of BL. The other 16 are shared among BL, BL-Reig, and SG. Given the higher accordance of BL-Reig alleles with other BL subpopulations than with SG, it is unlikely that the assignment of BL-Reig to the SG region is due to the release of individuals from SG at the site. Subsequent analysis of the firstlevel clusters split BE and LU into two separate groups and split BL-Reig from the SG cluster. Further reanalyses generated a total of 32 clusters (Fig. 3).
The AMOVA showed that genetic diversity was predominantly distributed within subpopulations (62%) but still considerable amounts of variation around 20% each occurred among subpopulations within populations and among populations (Table 5).
We could not derive a single model as the most parsimonious explanation for all measures of genetic ª 2013 The Authors. Ecology and Evolution published by John Wiley & Sons Ltd variability and structure. Allelic richness was best explained by variation in elevation (Table 6), but the confidence interval of the estimate overlapped with zero, suggesting that elevation did not provide a robust explanation for the observed pattern. The intercept-only model was the AICc-best model for H e , suggesting that none of the explanatory variables explained variation in H e well (Table 6). Two other models were within~2 DAICc values (Table 6), but the confidence interval of the estimates of the coefficients of the explanatory variables in these models included zero (Table 7) indicating that these alternative models are not reliable explanations for the data. F ST among subpopulations within populations was best explained by decline (Table 6), but the confidence intervals for the coefficients of the decline effect included zero (Table 7), and, as above, we did not accept this top-ranking model as a good fit. The same was true for two other models that ranked within 2 DAICc values of the best model. Inspection of the spline correlograms showed that there was no spatial autocorrelation.

Discussion
We predicted that we would find differences in subpopulation and population genetic structure and diversity that correlated with population declines and connectivity. Specifically, we expected to see differences among subpopulations attributable to different population decline histories and habitat variables expected to influence amphibian dispersal. Although these factors varied among the four regions, we detected no evidence that variation in genetic structure and diversity was influenced by any geographical  or historical factor (see also Dudaniec et al. 2012). However, estimates of numbers of breeding adults (N b ) were small and levels of genetic variation were very low. Moreover, private alleles were common in subpopulations, supporting our estimates of small breeding population sizes (N b ). Additionally, we identified strong genetic differentiation among subpopulations even across small geographical scales. Distances of less than 500 m were associated with significant pairwise F ST -values. Relatively large pairwise F ST values indicate moderate rates of gene flow. A pairwise F ST value of 0.2 implies that subpopulations exchange about one migrant per generation. Nevertheless, given the very small breeding population sizes that we estimated (Table 4), drift may have a stronger effect than gene flow and may quickly lead to the loss of rare and immigrant alleles and insignificant IBD (Fig. 4). A high proportion of genetic variation occurring among subpopulations within regions (Table 5) suggests that gene flow is not strong enough to homogenize populations.
The fact that local extirpation did not lead to changes in genetic structure of declining populations can be attributed to different causes. We can rule out insufficient time since declines as an explanation because 20 or more generations have passed since the earliest recorded declines (B€ oll et al. 2012). One possibility is that extirpations affected predominantly sink subpopulations (Gill 1978;Semlitsch 2000;Hels 2002), and loss of sink subpopulations should not result in a reduction in allelic variation (Hanski and Ovaskainen 2000).

Conservation implications
Matocq and Villablanca (2001) and Short Bull et al. (2011) pointed out the importance of suitable reference groups when interpreting population genetic data. Our study emphasizes the need for a comparative approach when ascertaining the effect of any threatening process on population genetic parameters. Had we only analyzed the subpopulation genetic structure in populations where strong declines of A. obstetricans had been reported, we might have concluded that declines were associated with low levels of genetic diversity and strong subpopulation differentiation. On the basis of this, we might have recommended translocations from other populations to enhance genetic diversity. The inclusion of reference populations where the species has not declined or experienced reduced rates of decline allowed us to show that decline and genetic structure were not associated. Certainly, our results showed that researchers should not simply assume that there will be a genetic signature of species decline and (sub)population loss, even if declines have been ongoing for a sufficient amount of time for genetic signatures of decline to manifest.
There is still cause for concern for common midwife toads in Switzerland based on our, and other, results. Allelic richness and heterozygosity were both extremely low across all four regions and A. obstetricans subpopulations in Switzerland are generally composed of relatively few breeding adults (this study and Tobler et al. 2012). Small population size and isolation of subpopulations mean that demographic stochasticity is a potential threat to population survival (Lande 1993). Small and isolated subpopulation sizes coupled with low levels of genetic variability are not ideal situations for long-term subpopulation persistence. Indeed, the fact that in two of our study populations subpopulation loss rate exceeded 50% clearly shows that A. obstetricans is a species at risk in Switzerland (Schmidt and Zumbach 2005). Given our analyses of the genetic structure of A. obstetricans populations, we believe that translocations to allow for genetic rescue would not be a suitable conservation strategy for this species. Our results showed that most genetic diversity is found within rather than among subpopulations, and that this is the case for stable and declining populations alike. Hence, the translocation of individuals from one subpopulation to another seems unlikely to benefit the long-term survival of subpopulations.
Increasing size (i.e., viability) of subpopulations appears to be the most promising conservation strategy. Larger subpopulations are less prone to environmental and demographic stochasticity and have a higher chance of long-term persistence. Increasing the number of subpopulations, for example, by the creation of new breeding sites, would also benefit long-term survival of populations. However, given the low dispersal rate and the resulting low chance for natural colonizations to occur, conservation management of Alytes should not rely solely on this approach. Thus, in the short term, increasing subpopulation size seems the most promising strategy to reduce local extinction risk while increasing connectivity through the establishment of new subpopulations may be a longterm goal.