Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Genetic Differentiation, Isolation-by-Distance, and Metapopulation Dynamics of the Arizona Treefrog (Hyla wrightorum) in an Isolated Portion of Its Range

  • Meryl C. Mims ,

    mmims@usgs.gov

    Current address: United States Geological Survey, Forest and Rangeland Ecosystem Science Center, Corvallis, Oregon, United States of America

    Affiliation School of Aquatic and Fishery Sciences, University of Washington, Seattle, Washington, United States of America

  • Lorenz Hauser,

    Affiliation School of Aquatic and Fishery Sciences, University of Washington, Seattle, Washington, United States of America

  • Caren S. Goldberg,

    Affiliation School of the Environment, Washington State University, Pullman, Washington, United States of America

  • Julian D. Olden

    Affiliation School of Aquatic and Fishery Sciences, University of Washington, Seattle, Washington, United States of America

Abstract

Population attributes such as diversity, connectivity, and structure are important components of understanding species persistence and vulnerability to extinction. Hyla wrightorum, the Arizona treefrog, is native to the southwestern United States and Mexico, and an isolated group of populations exists in the Huachuca Mountains and Canelo Hills (HMCH) of southeastern Arizona, USA. Due to concerns about declining observations of the species within the isolated HMCH portion of its range, the HMCH group is currently a candidate for federal protection under the U.S. Endangered Species Act. We present results of a genetic study examining population diversity, structure, and connectivity within the HMCH region. We sampled DNA from H. wrightorum larvae and adults from ten distinct locations, 8 of which were breeding sites and 4 of which were previously undescribed localities for the species. We developed and genotyped 17 polymorphic microsatellite loci and quantified genetic diversity, population differentiation, and landscape influences on population genetic structure. We found evidence of larger than expected effective population sizes, significant genetic differentiation between populations, and evidence of distance being the primary driver of genetic structure of populations with some influence of slope and canopy cover. We found little evidence of recent genetic bottlenecks, and individual-based analyses indicate admixture between populations despite significant genetic differentiation. These patterns may indicate that the breeding sites within the Huachuca Mountains constitute a metapopulation. We suggest that the HMCH region may contain larger and more connected breeding populations than previously understood, but the dynamics of this system and the limited geographic extent of the HMCH group justify current concern for the persistence of the species in this region. Efforts to ensure availability of high-quality breeding habitats and control for local threats such as effects of invasive predators may be critical to the persistence of these unique populations of H. wrightorum.

Introduction

Effective conservation of species vulnerable to extinction throughout some or all of their range requires knowledge of population attributes such as diversity, connectivity, and structure [1]. For decades, population genetic approaches have helped successfully identify scenarios that may compromise the health, resilience, or persistence of a species and its populations [2]. Population genetic approaches have proven particularly useful in evaluating the status of amphibian species [3], many of which are declining globally due to threats including habitat loss, invasive species, disease, and climate change [45].

Hyla wrightorum (Taylor, 1938), the Arizona treefrog, is an anuran native to the southwestern United States and Mexico and is currently a candidate species for federal protection by the U.S. Fish and Wildlife Service (USFWS) within a subset of its range [6]. Hyla wrightorum is most commonly associated with streams, cienegas (wetlands), and manmade ponds, particularly during the summer monsoon season when intermittent pools are used for breeding [78]. The species’ range includes three disjunct regions, with the majority of the range occurring in two disparate regions in the United States and in Mexico. Spatially intermediate to these two regions is an isolated portion of the species’ range in the Huachuca Mountains and Canelo Hills (HMCH) region of southeastern Arizona, USA (Fig 1) [910]. This distinct portion of the species’ range is notable for its geographic isolation, and genetic evidence suggests that the HMCH region has been isolated from the other portions of the species’ range since the late Pleistocene [11].

thumbnail
Fig 1. Hyla wrightorum range map and sampling locations.

The known range extent of H. wrightorum [6, 9] (A, B), and H. wrightorum sampling locations in the Huachuca Mountains and Canelo Hills region (C). Symbols indicate 2014 sampled locations (blue) or unsampled locations at which H. wrightorum were observed in 2015 (white). Symbols for sampled locations are proportional to sample size, and color indicates life stage of sampled individuals (dark blue = adults; light blue = larvae). Population numbers are shown in black font, and information for each sampling location included in Table 1.

https://doi.org/10.1371/journal.pone.0160655.g001

Hyla wrightorum populations within the HMCH region are of conservation concern. In the past two decades, confirmed observations have been made at only half (8 of 16) of previously reported locations [6, 8]. The geographic range occupied by H. wrightorum in the HMCH is small compared to the larger two portions of its distribution (Fig 1), with known breeding sites occurring within an area no larger than 85 km2 [6]. Based on previous sampling efforts, HMCH population sizes were presumed to be small, with breeding choruses typically including 30 or fewer adults [6, 7]. Potential threats to local persistence of this species include disease [12] and predation by invasive American bullfrogs (Lithobates catesbeianus [13]). In addition, available breeding habitat is generally scarce in the region due to topographic constraints limiting suitable shallow breeding habitat [7]. Contemporary H. wrightorum breeding habitat in the HMCH region largely consists of manmade ponds originally built to provide water for livestock. Much of the region’s wetland habitat has been dewatered, and these “stock ponds” now serve as surrogate aquatic habitat for many native aquatic species. However, stock ponds also host non-native species, such as bullfrogs, which compete with and prey upon a host of native aquatic and terrestrial species [14]. Taken together, these factors contribute to the concern regarding the vulnerability of H. wrightorum populations in the HMCM region to regional processes such as dewatering, catastrophic fires, and other manmade disturbances likely to increase in frequency and intensity with climate change and increasing human development [1517].

These concerns, coupled with the phenotypic, genetic, and geographic uniqueness of this group, have led to candidacy for federal protection of H. wrightorum in the HMCH region [6]. However, the status of populations in this isolated range is generally unknown. It is difficult to determine whether apparent absences over the last two decades at up to half of historically reported sites are due to true declines, natural population fluctuations, or failure to detect the species. Because H. wrightorum may form breeding choruses only one or a few nights in a given year, traditional census and survey efforts may fail to provide reliable information, particularly if populations are in fact in decline [18]. Genetic approaches provide promising alternatives to traditional census and survey efforts for evaluating population structure, size, and connectivity for many amphibians [3] and may offer valuable insight into the status of H. wrightorum in the Huachuca Mountains region.

We present results of a population genetic study for H. wrightorum in the HMCH region of southeastern Arizona. The objectives of our study were to address current knowledge gaps regarding population genetic diversity, population structure, and effective population sizes (Ne) of this species in the region. To achieve this, we used a combination of population-level and individual-based approaches. Based on current knowledge of the species, we hypothesized small (< 30) effective population sizes for H. wrightorum in the HMCH region. This was based on both small observed breeding choruses for this species [6, 11] as well as recent reported estimates of Ne for a sympatric congener (H. arenicolor, average Ne = 30.7 [19]). We also searched for evidence of any recent genetic bottlenecks. Finally, we tested five hypotheses regarding the effects of landscape attributes on gene flow. First, if populations are sufficiently small, genetic drift may result in complete decoupling of genetic structure and spatial and landscape factors, producing isolated populations. This has been observed for the pond breeding salamander Ambystoma mavortium stebbinsi, sympatric to H. wrightorum in the Huachuca Mountains region [20]. Alternatively, an isolation-by-distance pattern may emerge in which populations in closer geographic proximity are genetically more similar than those farther apart [21]. We also tested three additional hypotheses of landscape effects: isolation-by-slope, in which steep topography inhibits gene flow; connectivity-by-canopy, in which gene flow is higher along more shaded areas that may provide more cover from predation or desiccation; and connectivity-by-stream in which gene flow is higher along riparian corridors. We conclude by discussing how the findings of this research can help inform management and conservation efforts for this species.

Methods

Sampling

Adult and larval H. wrightorum individuals were sampled during the summer monsoon season of 2014. We selected sampling locations using a combination of historical records of H. wrightorum occurrences as well as opportunistic visits to potential habitat (ponds and wetlands) in the known range of the species (Fig 1, Table 1). Sites were visited in the evenings, typically following thunderstorms, to maximize the chance of hearing breeding choruses. Sites were also surveyed during daylight hours using a combination of visual search and dip-netting. DNA was collected from each sampled individual via tail clip (larvae) or buccal swab (adult, following Goldberg et al. [22]). Tail clips were immediately stored with a desiccant (Drierite), and buccal swabs were immediately placed in a vial with buffer ATL from a DNeasy Blood & Tissue DNA extraction kit (QIAGEN). Tail clips were then kept at room temperature and buccal swabs at -10°C until DNA was extracted. Sampling protocol was approved by the University of Washington’s Institutional Animal Care and Use Committee (Protocol # 4172–03). All efforts were made to minimize animal handling time and suffering. All sampling was conducted on public lands, and sampling permits were obtained from the U.S. Forest Service (Special Use Permit SIE0150), Arizona Game and Fish Department (SP685479), and from the U.S. Army Garrison Fort Huachuca Environmental and Natural Resources Division (IMHU-PWB 200–1).

DNA extraction, genotyping, and marker screening

Whole genomic DNA was extracted using DNeasy 96 Blood & Tissue Kit (QIAGEN) at the Molecular Ecology Research Lab at the University of Washington’s School of Aquatic and Fishery Sciences. Hyla wrightorum microsatellite markers were developed by the Evolutionary Genetics Core Facility at Cornell University and are described in S1 Table, including GenBank accession numbers. Polymerase chain reaction (PCR) was used to amplify DNA for multiplexed loci using Multiplex PCR kits (QIAGEN). Reactions consisted of 0.2 μM of each primer, 1X Qiagen Multiplex PCR Master Mix, and 1μl DNA in a 10 μl reaction. PCR conditions followed QIAGEN guidelines and included an initial activation step of 15 minutes at 95°C, followed by 25 cycles through three steps: denaturation (30 seconds at 94°C), annealing (90 seconds at 60°C), and extension (90 seconds at 72°C). All PCR reactions were performed on C1000 Touch or S1000 thermal cyclers (Bio-RAD). PCR products were genotyped using 3730xl 96-Capillary Genetic Analyzer (Applied Biosystems) at Yale University’s DNA Analysis Facility (New Haven, CT). Genotypes were analyzed using the software GENEMAPPER 4.1 (Applied Biosystems). Individuals with poor-quality genotype data were processed again through PCR and genotyping, and in a few cases DNA extractions were performed again provided sufficient tissue remained. Individuals with > 25% missing data after three genotyping attempts were discarded from downstream analyses.

Alleles were binned using the program TANDEM [23]. Any alleles occurring in < 3 individuals and > 1 repeats from other alleles at the locus were confirmed by re-amplifying and genotyping the samples. Loci were screened for the presence of linkage disequilibrium using the log-likelihood ratio statistic for each pair of loci in each population (GenePop 1.2 [24]). Loci were also screened for deviations from Hardy-Weinberg equilibrium (HWE) using Fisher’s exact test as implemented in GenePop, and the presence of null alleles was evaluated with Micro-Checker [25]. All screening procedures were performed first using adult samples only (three populations) and then again with all population samples (adults and larvae) following full sibling removal (description follows).

Larval samples can bias population genetics findings by artificially inflating genetic differentiation due to family structure [26]; therefore, we screened all larval samples for full siblings using the program COLONY 2.0 [27]. All genotype data (with and without siblings removed) are available via figshare (doi: 10.6084/m9.figshare.3365293). Microsatellite loci were first screened for linkage disequilibrium or deviations from HWE using adult genotypes to identify any loci that may affect detection of siblings. When full siblings were detected, one sibling was retained from each family in the final dataset.

Genetic diversity and population structure

We calculated genetic diversity estimates of expected heterozygosity (HE), observed heterozygosity (HO), allelic richness (AR) rarified to the smallest number of sampled individuals per population, and FIS [28]. We also calculated two global measures of genetic differentiation, global FST [28] and G'ST [29]. A hierarchical analysis of molecular variance (AMOVA) among and within sampling locations was performed using ARLEQUIN 3.1 [30] with 10,000 permutations to assess significance. Pairwise genetic distance (between each pair of sample sites) was examined in three ways: Dps, a method of measuring genetic differentiation based on proportion of shared alleles [31], FST [28] bootstrapped with 10,000 replicates to determine significant difference from zero, and Slatkin’s linearized FST, calculated as FST /(1- FST) [32]. Here we report Slatkin’s linearized FST as linFST. Genetic diversity, global genetic differentiation, and pairwise genetic differentiation metrics were calculated using MSA 4.05 [33].

We estimated effective population size (Ne) for each sampling location and for all individuals combined. Ne for each sampling location was estimated using the linkage disequilibrium method (LDNe) of Waples and Do [34] as implemented in NeEstimator V2 [35]. LDNe has been shown to perform well in simulations, even in the presence of migration [36]. However, migration can inflate population-level Ne estimates, and overall Ne estimates can be downwardly biased if sub-structure exists among populations [37]. In order to account for these potential biases, we used a heirarchical estimation of total Ne that accounts for population structure. We estimated total Ne following Wright [21] as implemented for hierarchical estimation of total Ne in Holleley et al. [38]: Total Ne was calculated with mean local Ne estimates and upper and lower 95% jackknifed confidence interval values. A value of 10,000 was used in the case of local Ne estimates with infinite upper confidence values [19].

Evidence of recent bottlenecks was evaluated using the program BOTTLENECK 1.2.02 [39]. We performed a Wilcoxon signed rank test (1,000 iterations), which is appropriate for tests evaluating fewer than 20 loci [39], to determine if observed heterozygosity was higher than expected in a population at mutation drift equilibrium, indicating that a genetic bottleneck had occurred. We followed recommended parameters for microsatellite markers (100% of mutations following the infinite allele model and a variance of 0.36 for the geometric distribution of the model) [39], and a Bonferroni correction was applied to adjust the critical p-value for multiple comparisons. We also performed a mode-shift test to determine whether distortion of the expected L-shape frequency distribution was detected and attributable to a recent bottleneck in any populations [40].

To evaluate whether individuals sampled at a location likely originated from the local population, we used an individual assignment analysis as implemented in GeneClass2 [41]. We used the Bayesian approach of Rannala and Mountain [42] with 1,000 Monte Carlo-simulated individuals per sample. The designation of individuals as migrants was based on the likelihood of the individual genotype within the population where the individual was sampled (L_home), which is appropriate in cases where source populations may be missing from the samples (i.e., some populations may not have been sampled) [41]. Individual-based hierarchical population structure was evaluated using the Bayesian clustering program STRUCTURE 2.3.4 [43]. Each sampling site was treated as an independent putative population with a total of n putative populations. Ten iterations of each K from 1 to n + 1 were run for 500,000 cycles with a burn-in of 50,000 cycles. Given the close geographic proximity of these populations and the likelihood that they are relatively closely related, we allowed admixture and correlated allele frequencies. Incorporating sampling information using the LOCPRIOR model has been shown to perform better than models without LOCPRIOR in cases of weak but significant genetic structure [43]. Because we anticipated weak population structure for H. wrightorum based on previous studies of sympatric pond-breeding species [19, 44], we conducted STRUCTURE analyses with LOCPRIOR. We then compared our results to those obtained without using LOCPRIOR. We determined the most likely K using the delta-K method [45] in which the most likely value of K is assessed by the second-order rate of change in the log-likelihood. A delta-K value cannot be calculated for K = 1. However, K = 1 is assumed most likely for runs in which K = 1 has the greatest log-likelihood [46]; α, the Dirichlet parameter for degree of admixture, varies throughout a run rather than converging [47]; and assignment to genetic clusters when K > 1 tends to be highly admixed within individuals [47]. We identified terminal clusters (K = 1) by first examining log-likelihood and secondarily by visually inspecting α and individual admixture results. This analysis was repeated for genetic clusters in which both K > 1 and n > 1 to identify hierarchical population structure until terminal clusters were described [45]. Individuals from a given sampling location were not subdivided in hierarchical analyses. Rather, sampling location groups were kept intact to evaluate hierarchical structure across rather than among sampling locations. STRUCTURE output visualizations were constructed using the program DISTRUCT 1.1 [48].

Landscape genetic analysis

To examine our hypotheses of landscape effects on population genetic structure, we first constructed four landscape connectivity surfaces using CIRCUITSCAPE [49]. CIRCUITSCAPE uses circuit theory to simulate gene flow (i.e., “current”) through a resistance surface in which landscape features hypothesized to promote gene flow are assigned low resistances, and landscape features hypothesized to inhibit gene flow are assigned high resistances. CIRCUITSCAPE allows gene flow across multiple pathways and reports pairwise summations of resistance between sampling locations. Modeling multiple pathways is appropriate for dryland anurans with potentially high dispersal ability [19, 44]. To generate pairwise resistance data, we built raster maps of resistance (low to high) using spatial data describing stream networks, canopy cover, and topography (slope). Four resistance surfaces were created to represent four simple models of landscape effects on gene flow: isolation-by-distance (uniform resistance across the landscape); isolation-by-slope (high resistance for steep regions, and low resistance across flat land), connectivity-by-canopy (low resistance with high canopy cover, high resistance in areas with low or no canopy), and connectivity-by-stream (low resistance along riparian corridors, defined as a 100 m buffer surrounding streams). All spatial data were obtained from publicly available sources and are described in detail in the supporting information, S8 Table. A geographic information system (ArcGIS 10.1, Environmental Systems Research Institute) was used to catalog and manipulate landscape data and generate resistance raster maps. Each resistance map was scaled for hypothesized landscape resistance to gene flow from 1–100 where 1 indicates low resistance and 100 indicates high resistance. Values of slope and percent canopy cover were scaled linearly from 1–100 to create resistance surfaces, and stream resistance surface consisted of values of 1 (low resistance) within stream buffers and 100 (high resistance) outside stream buffers. The scale of resistance values was arbitrary and was designed to reflect hypothesized relationships between landscape features and genetic connectivity. All resistance surfaces were scaled to 30 m resolution—the minimum resolution available for canopy cover (see S8 Table). The spatial extent of resistance rasters ensured a buffer of at least 7 km from the edge of the raster to a given sample location. We selected a 7 km buffer distance because 7 km is the farthest nearest-neighbor distance between any sampling locations. With this grain and extent, we were able to perform all CIRCUITSCAPE analyses in the pairwise source/ground modeling mode and using a cell connection scheme of eight neighbors, allowing maximum freedom of current flow. Finally, we evaluated correlations between pairwise landscape resistances, reported these correlations in the supporting information (S7 Table), and did not include resistances with r > 0.7 in the same models.

We evaluated relationships between landscape resistance and population genetic differentiation (as measured by Dps and linFST) in two ways. First, we first employed a causal modeling framework as implemented in Cushman et al. [50]. By calculating a series of Mantel tests [51] and partial Mantel tests, we assessed the relationship between each landscape resistance matrix and the genetic distance matrix while sequentially partialling out the influence of the other landscape resistance matrices. This approach is effective in identifying spurious correlations between genetic distance and landscape processes, particularly when landscape factors may be correlated [52].

Additionally, we evaluated relationships between pairwise genetic distance and pairwise landscape resistances using a mixed-effects modeling approach [53]. Through mixed-effects modeling, explanatory variables (pairwise landscape resistances) are treated as fixed effects, and sampling locations are treated as random effects in a pairwise covariance structure to account for non-independent values in distance matrices. All resistance variables (fixed effects) were mean-centered, and we evaluated model fit with three assessment criteria. We calculated Akaike’s Information Criterion corrected for small sample sizes (AICc [5455]) and Bayesian Information Criterion (BIC [56]). We also evaluated models using the R2β statistic [57]. R2β compares a model with fixed and random effects (pairwise landscape distance or resistance and sampling location) to a null model with only the random effect (sampling location) and an intercept. All analyses were performed in R version 2.14.0 (R Development Core Team, 2012), using a modified version of lme4 [58] and PBKRTEST [59] for R2β calculation as described in van Strien et al. [53].

Results

We sampled H. wrightorum individuals at 10 locations for a total of N = 299 individuals, 82 of which were adults (Figs 1 and 2, Table 1). Four of the 10 sampling locations had no previous record of H. wrightorum presence (sites 2, 4, 5, and 10). Two locations had insufficient sample sizes for population-level analyses (populations 2 and 5); among the remaining 8 populations, the mean sample size was 36.6 individuals, with a minimum of 19 and a maximum of 50. We retained only 1 full sibling from each family identified, reducing the final sample size to a total of N = 221 (Table 2). Among the 8 populations with sufficient sample sizes for population level analyses, mean sample size with full siblings excluded was 27.6, with a minimum of 19 and maximum of 40 (Table 1).

thumbnail
Fig 2. Hyla wrightorum individuals and breeding sites.

(A) H. wrightorum breeding pair at Site 1; (B) Site 10; (C) Site 2; (D) Site 4; and (E) H. wrightorum larvae sampled at Site 4. Sites in B, C, and D are all new localities for this species.

https://doi.org/10.1371/journal.pone.0160655.g002

thumbnail
Table 2. Results of Hyla wrightorum sibling identification among larval samples.

https://doi.org/10.1371/journal.pone.0160655.t002

All individuals were genotyped for 17 novel polymorphic microsatellite loci (Table 3; microsatellite loci described in S1 Table, including GenBank accession numbers). Genotyping failed for nine larval individuals (not included in total sample size) in which individuals were missing > 25% genotype data despite multiple re-run attempts. In some cases, a second DNA extraction attempt was performed. Some of these larvae were either very small when sampled or were found dead among live larvae and salvaged as voucher specimens. We suspect that species identification was incorrect for these individuals (in the case of very small larvae that may have been H. arenicolor) or DNA was too degraded for successful genotyping (in the case of salvaged individuals). We found deviations from Hardy-Weinberg expectations at 21 of the 136 sample/locus combinations (15%) prior to application of a Bonferroni correction. However, with a Bonferroni correction applied, we found only 5 instances of significant deviation from HWE (S2 Table). Because these deviations occurred across five different loci and three different populations, we elected to retain all loci in our analyses. Significant LD for two markers pairs were found in only one of 8 populations each, and no evidence of consistent LD was observed for any marker pair (S3 Table). We also found no evidence of null alleles or large allele dropout from Micro-Checker.

thumbnail
Table 3. Polymorphic microsatellite loci for Hyla wrightorum.

https://doi.org/10.1371/journal.pone.0160655.t003

Estimated total Ne was 972.1 (mean), with a jackknifed 95% confidence interval of 401.7 to 30402.4. Population-level Ne averaged 116.7 individuals per population, with a minimum estimate of 32.5 and a maximum estimate of 277.9. Confidence intervals for Ne ranged from 23.1 to infinite values (Table 1). Observed heterozygosity averaged across all loci ranged from 0.647 (Site 4) to 0.695 (Site 7), and expected heterozygosity averaged across all loci ranged from 0.671 (Site 4) to 0.736 (Site 7). Allelic richness ranged from 5.47 (Site 4) to 6.24 (Site 8). Information for population and locus pairs, including observed heterozygosity, expected heterozygosity, allelic attributes, and FIS are included in the supporting information, S2 Table. We observed a negative relationship between Ne and FIS (Fig 3, R2 = 0.72), which may indicate some inbreeding in populations with smaller effective population sizes. Wilcoxon tests for bottlenecks found evidence of a recent bottleneck among individuals at only one location (Site 7), and mode-shift tests found no evidence of bottlenecks (Table 4).

thumbnail
Fig 3. Population-level Ne and FIS for Hyla wrightorum.

Relationship between effective population sizes (Ne) derived from LDNe (Ne, x-axis) and FIS (y-axis) for 8 H. wrightorum populations sampled in the Huachuca Mountains and Canelo Hills, Arizona. The black line and R2 reflect a logarithmic regression with the equation y = 0.03ln (x) + 0.19. All genetic diversity metrics are derived from 17 microsatellite loci described in the text.

https://doi.org/10.1371/journal.pone.0160655.g003

Global genetic differentiation measures provided evidence for genetic structure among populations with values of FST = 0.04 and G’ST = 0.19. A hierarchical AMOVA revealed modest but significant genetic structure with 4.2% variation explained among populations and 3.4% explained among individuals (supporting information, S4 Table). Pairwise FST revealed significant differentiation among populations; 25 of 28 pairwise comparisons were significantly different from zero (FST, Bonferroni correction applied, Table 5). The three non-significant pairwise comparisons all involved Site 6, one of the most central sampling sites (Fig 1). Site 6 was not significantly different from Sites 1, 7 or 8.

thumbnail
Table 5. FST, linFST, Dps values for Hyla wrightorum populations.

https://doi.org/10.1371/journal.pone.0160655.t005

Despite significant genetic differentiation between populations, 50.2% of individuals were assigned to populations other than their sampling location (Table 6, S5 Table). In addition, 5 individuals had low probabilities (<5%) of assignment to all reference populations, indicating these individuals may be recent migrants from unsampled populations. Admixture between populations was also evident in the Structure analyses, which provided support for modest hierarchical clustering among H. wrightorum individuals with evidence of isolation-by-distance (Fig 4). Support for K = 2 and K = 3 clusters across all individuals was similar (S6 Table). For both K = 2 and K = 3, we found support for a general shift from one group (blue, the dominant group for Sites 1, 2, 3, and 4) to a second group (orange, dominant for Sites 9 and 10). This pattern followed a north to south trend, with individuals from centrally located sites assigned to both clusters (Fig 4). In the K = 3 scenario, we identified a third genetic group within Site 10. Fig 4 shows both K = 2 and K = 3 for all individuals, and hierarchical results are shown for the K = 3 scenario. Site division for hierarchical analysis included sites 1, 2, 3, and 4, which had few individuals assigned to the second group (orange), and sites 6, 7, 8, and 9, which had a modest to dominant presence of the second group. Hierarchical results for the alternate K = 2 scenario are included in the supporting information, with the only difference from the K = 3 scenario being the inclusion of site 10 in the second group (S6 Table). At the secondary level (within major clusters identified across all individuals), we found support for separate genetic clusters in Sites 3 and 9. Sites 5, 6, 7 and 8 were ultimately grouped together with strongest support for K = 1. We found support for K = 3 clusters for the group encompassing Sites 1, 2, and 4; however, visual inspection of structure within this group does not support strong differentiation between populations but rather a possible isolation-by-distance pattern. STRUCTURE analyses without LOCPRIOR largely agree with the results of the analyses performed with LOCPRIOR (supporting information, S7 Table). There was only one discrepancy between analyses with and without LOCPRIOR: analyses with LOCPRIOR found modest support for genetic clusters among sites 1, 2, 3, and 4, whereas analyses without LOCPRIOR suggest that those sites may be panmictic. Delta-K tables for all species and genetic clusters are included in the supporting information (S6 and S7 Tables).

thumbnail
Fig 4. Structure results for Hyla wrightorum.

Numbers correspond to sampling locations (Fig 1, Table 1). Each vertical bar represents one individual. Colors indicate the most likely genetic cluster assignments. Black vertical bars denote individuals from the same sampling locations. Each cluster was hierarchically analyzed for nested structure; nested structure results are shown directly below the original cluster. Hierarchical analyses were repeated until terminal clusters (K = 1) were reached. Note that results for K = 2 (upper-most Structure plot) and K = 3 (just below K = 2) clusters are shown for all individuals, because both had similar support (S6 Table). Hierarchical analyses are shown following K = 3 for all individuals. Nested, colored outlines on the map correspond to population clusters.

https://doi.org/10.1371/journal.pone.0160655.g004

We found multiple lines of evidence that, of the spatial attributes examined, distance is a primary driver of population structure. We also found evidence that slope and canopy may affect population structure. Of the four Mantel tests examining hypothesized relationships between spatial and landscape variables (Distance, Slope, Canopy, and Stream) and genetic distance (Dps and linFST), only Distance was significantly univariately related to Dps (Mantel r = 0.50, p-value = 0.015). LinFST was not significantly univariately related to any landscape variables (Table 7). Distance remained significantly related to Dps in three partial Mantel tests controlling for Slope, Canopy, and Stream, and Distance was significantly related to linFST in two of three partial Mantel tests (controlling for Slope and Canopy). These findings support the expected linear relationship between resistance (of a uniform landscape) and genetic distance in an isolation-by-distance scenario [49] (S1 Fig). For both measure of genetic distance, isolation-by-slope had support as a driver of genetic structure after accounting for Canopy, and there was also a significant effect of canopy after accounting for slope (Table 7; but note the high correlation between canopy and distance, S7 Table). There was no support for stream as a driver of genetic distance in this dataset.

thumbnail
Table 7. Mantel and partial Mantel results for Hyla wrightorum.

https://doi.org/10.1371/journal.pone.0160655.t007

Mixed-effects models also support distance as the primary driver of genetic structure for this dataset (Table 8). The information criteria used (AICc and BIC) were consistent in indicating that the best model was the model with Distance only for both Dps and linFST (Table 8). The standard threshold for identifying the best models in a suite of candidate models is a difference of 2 [51], and the ΔAICc and ΔBIC values for the next best model of isolation-by-slope were well above this threshold for both measures of genetic distance. Distance had the highest R2β value among univariate models for Dps, and Distance was tied with Slope for the highest R2β value for linFST. Of the more complex models, those that included Canopy outperformed those with Distance for mixed-effects models according to the R2β metric (Table 8), but not according to AICc and BIC.

thumbnail
Table 8. Mixed-effects models and landscape genetic results for Hyla wrightorum.

https://doi.org/10.1371/journal.pone.0160655.t008

Discussion

We found a higher number of breeding sites than previously described and evidence of larger than expected population sizes for this isolated group of populations of H. wrightorum within the HMCH region. We found support for significant genetic differentiation between many populations despite evidence of admixture. The influence of spatial and landscape variables on population structure was evident, with significant support for isolation-by-distance and some support for the effects of slope and canopy on genetic structure. Taken together, these results suggest that populations of H. wrightorum in the Huachuca Mountains may constitute a metapopulation in which the spatial and temporal variability of available breeding habitat may contribute to local extinctions, colonizations, and exchange of individuals. Similar evidence exists for many pond breeding amphibians [60].

The discovery of previously unknown breeding sites for this species as part of this research effort supports two alternative hypotheses for why this species has not been observed at some historical sites over the last few decades. First, previous survey efforts may simply have missed active breeding sites due to spatial or temporal mismatches between survey efforts and breeding activity. Alternatively, natural population fluctuations or metapopulation dynamics may account for the absence of this species at historical breeding sites and the discovery of the species at new sites over the last two decades. Over the course of an intensive 3-week field season in 2014, we found H. wrightorum in 4 new locations (Sites 2, 4, 5, and 10). Site 4 is a confirmed breeding site (larvae present), and Site 10 is a presumed breeding site (> 20 adults sampled in a breeding chorus). Site 10 extends the known range of contiguous breeding sites in the Huachuca Mountains southeast by roughly 7 km. Additionally, limited observations of H. wrightorum have been reported at a small number of wetland habitats at Rancho Los Fresnos, Sonora, Mexico, south of the Huachuca Mountains [6]. The Los Fresnos population is presumed to be small and is typically described as an “outlier” from the Huachuca and Canelo Hills populations due to its isolation and distance from other known breeding sites in the Huachuca Mountains. However, the newly discovered Site 10 is approximately intermediate in distance between Los Fresnos and previously described breeding sites, representing a significant reduction in the spatial gap between these sites. Survey efforts for this species in 2015 (after the conclusion of sampling efforts for this study) confirmed H. wrightorum at three additional breeding locations, all previously undescribed (Fig 1; CSG, K. Strickler unpublished data). Two 2015 sites are also intermediate to previously known sites and the Los Fresnos location. In general, the discovery of new sites during this study supports the hypotheses that these populations exist within a metapopulation in the region, and that more intensive survey efforts are required for detection of this species. Results from the GeneClass analysis revealed that some individuals may have originated from one or more unsampled populations in the HMCH region.

In consideration of these observations, continued survey efforts exploring additional plausible breeding sites within the region are recommended. For example, within a 7 km buffered area encompassing all sites in this study, over 90 ponds with intermediate hydroperiods were identified using satellite imagery (publically available imagery, accessible via Google Earth). Many of these ponds were not visited in 2014 due to logistical and time constraints, and many have not been surveyed for this species. Temporal dynamics such as the frequency, timing, and duration of the period during which these ponds are wet versus dry is not known. Future research might aim to quantify the hydrologic dynamics of these ponds, and how they may change in the future, in order to help managers better understand the spatiotemporal variability of available breeding habitat. Additionally, H. wrightorum are thought to prefer or require emergent vegetation in breeding habitat [7], and they likely avoid ponds with predatory or competitive species present [61]. These two criteria alone may exclude many of the intermittent ponds within 7 km of known breeding sites; however, specific habitat associations and requirements of H. wrightorum remain uncertain. Additional survey efforts aimed at detecting this species and refining knowledge of breeding habitat will be critical in identifying other potential breeding sites for this species within the HMCH region and for optimal management of known breeding sites [62].

Significant pairwise FST comparisons, AMOVA results, and evidence of hierarchical structure among populations indicate considerable genetic differentiation between populations, with a significant influence of spatial and/or landscape variables. Individual-based cluster analyses of H. wrightorum most strongly supported a predominant pattern of isolation-by-distance. Additional hierarchical analyses reveal some genetic clusters strongly associated with a certain locality (Sites 3, 9, and 10). However, the presence of heterogeneity in assignment of individuals to clusters (Sites 9 and 10) or the lack of support from analyses that did not incorporate sampling location information (Site 3) indicates weak clustering at these locations. Admixture between sites was also supported with the assignment test results. In cases of more complete genetic isolation, we might expect to see near-perfect correlation of cluster assignments with sampling location, as observed in the one of the region’s aquatic macroinvertebrates, Abedus herberti, which requires perennial aquatic habitat [63].

Our findings suggest significant genetic structure of H. wrightorum populations in the HMCH region, which is in contrast to some sympatric pond-breeding anurans, including Anaxyrus cognatus, Scaphiopus couchii, and Spea multiplicata, found to have panmictic population structure across similar spatial scales [19, 44]. However, the pond-breeding anurans in those studies are more desiccation tolerant, more mobile, have higher fecundity and larger population sizes, and are able to better utilize more ephemeral breeding ponds due to short larval periods (i.e., one- to three-week larval periods for species with the lowest larval requirements versus one to a few months for hylids in the region) [6465]. These life history differences may explain the greater genetic differentiation observed for H. wrightorum in this study. By contrast, Storfer et al. [20] report much greater genetic differentiation among populations of an endangered pond breeding salamander, Ambystoma mavortium stebbinsi, in the San Rafael Valley and foothills of the HMCH region. They found no support for isolation-by-distance among populations and instead found evidence that genetic drift and small population sizes are responsible for high levels of genetic differentiation among populations. We found that H. wrightorum population structure was most similar to H. arenicolor in the Huachuca Mountains [19], a predominantly stream-dwelling amphibian and H. wrightorum’s only congener in the region.

Our analyses comparing the relationship between genetic structure and landscape variables of distance, slope, canopy cover, and streams provide additional support for the role of spatial and landscape factors in the genetic structure we observed. Taken together, we found strong support for physical distance as a driver of genetic distance. We also found some support for the role of slope and canopy cover, indicating that populations that are closer together are more highly related, and that flatter, more forested areas may play a role in facilitating gene flow. Correlation between canopy cover and distance highlights the difficulty of interpreting collinear landscape distances [52]. However, our results provide strong support for distance as the most likely factor affecting genetic structure. Although we did not detect an influence of streams or a strong influence of canopy cover in this study, the study region receives considerable precipitation during most summers and has high canopy cover throughout this portion of the species range. If favorable landscapes (i.e., those that promote gene flow) occur at a high density on the landscape, detecting their influence can be difficult using genetic approaches [66]. Additionally, we used simple linear or binary resistance values for the landscape variables in this study, and it is possible that genetic distance is related to landscape variables in non-linear ways. Finally, we hypothesized and tested three landscape variables we expected might influence genetic structure of H. wrightorum populations, but additional landscape variables not considered in this study may also explain some portion of H. wrightorum genetic structure. Thus, these findings do not rule out the possibly important role of landscape variables other than distance in connecting populations.

Total and population-level Ne estimates had a high degree of uncertainty as evidenced by wide confidence intervals. This uncertainty may be due in part to migration and other demographic complexities that complicate the estimation of Ne, particularly for metapopulations. Although LDNe is considered one of the most reliable methods for estimating Ne with a single sample in time [36], migration and population structure can bias estimates. Migration can upwardly bias Ne due to immigrants creating a larger total pool of parents than a sample with only local breeders [37]. Given the high degree of gene flow observed between populations, migration may have inflated LDNe estimates within populations. By contrast, samples with sub-structure may downwardly bias estimates of Ne calculated with LDNe due to linkage disequilibrium created by the mixture of multiple gene pools [37]. We attempted to avoid this bias by implementing the hierarchical calculation of Wright [21] which accounts for population structure and has been shown to outperform non-hierarchical total Ne estimates [38]. However, migration and population dynamics likely make characterizing “true” Ne values difficult regardless of the method employed, and the Ne estimates in this study should be interpreted with caution.

Despite uncertainty in our Ne estimates, our results suggest that previous assumptions of consistently small population sizes, with Nc estimates near or below 30 individuals [6], likely underestimate population sizes of this species. For highly fecund species or r-selected species, such as amphibians, differential recruitment may result in a Ne / Nc ratio that is quite low [67]. Thus, it is likely that Nc of H. wrightorum populations in the study region are larger than the estimated Ne values reported here. Still, concern regarding the size of H. wrightorum populations in the region is warranted. We found small mean family sizes among larval samples, indicating that variance in reproductive success may be modest in this species. Thus, the ratio of Ne / Nc may be closer to 1 than for other amphibians. Mean Ne was less than 100 for the overall estimate and for half the populations we sampled, and the 95% confidence interval of Ne included values < 100 for all populations, as is typical for many pond breeding amphibians [68]. Ne < 100 may indicate these populations are at risk of genetic depletion due to genetic drift and/or inbreeding [3], particularly if gene flow is sufficiently low. In this case, the relationship between Ne and FIS indicates that inbreeding may occur in some populations.

Despite genetic connectivity between populations and evidence of higher than expected Ne, the dynamics of this system and small number of known populations indicate that these populations are of conservation concern. Ne values reported here should be considered in the context of temporal uncertainty, particularly because environmental instability, demography, and other factors that vary through time can influence population dynamics [69]. Although the breeding habitat of H. wrightorum in the HMCH region occurs in a largely remote area with little human activity or development, many of the breeding ponds are frequented by cattle. The effects of grazing and other disturbances associated with livestock are currently not known for this amphibian. Predation by the invasive American bullfrog (Lithobates catesbeianus) could significantly reduce the number of adults at breeding ponds—particularly if choruses are small [13]. Site 7, the only H. wrightorum population with any evidence of a bottleneck, may have faced decades of predation from bullfrogs that established at a perennial pond just meters away. An intensive bullfrog eradication program began in the area in the early 2000’s [14], and it is possible that historical bullfrog predation at Site 7 may have contributed to or driven the bottleneck we observed (T. Jones, personal communication). Batrachochytrium dendrobatidis (Bd), a disease associated with the decline and/or extinction of many amphibians worldwide [70], is documented in the region [71] and is known to infect H. wrightorum [72]. Hyla wrightorum from northern Arizona showed no increase in mortality risk when exposed to Bd in a controlled, experimental setting [73], but the effects of Bd populations in the HMCH region are not currently known.

On a regional scale, climate change will likely bring drier conditions to the southwestern United States [16] with the possibility of spatial and temporal reduction in aquatic breeding habitat as well as changes in disease dynamics [74]. Shifts in climate along with other human modifications to the landscape will also likely result in larger and more intense wildfires in the region [15]. Such fires are often accompanied by major flooding and erosion events that can scour or fill breeding sites with sediment. We suggest continued survey efforts to determine occupancy of breeding habitats by H. wrightorum, and we suggest additional research aimed at quantifying both the habitat requirements of and threats to this species to characterize suitable breeding habitat. Such efforts will help inform efficient and effective management of intermittent ponds to promote the persistence and continued connectivity of H. wrightorum populations within the HMCH region.

Supporting Information

S1 Fig. Relationship between genetic distance and uniform landscape resistance for H. wrightorum.

https://doi.org/10.1371/journal.pone.0160655.s001

(DOCX)

S1 Table. Primer sequences and GenBank Accession numbers for 17 H. wrightorum microsatellite loci.

https://doi.org/10.1371/journal.pone.0160655.s002

(DOCX)

S2 Table. Microsatellite loci information for each population of H. wrightorum sampled.

https://doi.org/10.1371/journal.pone.0160655.s003

(DOCX)

S3 Table. Results for test for linkage disequilibrium for H. wrightorum microsatellite loci.

https://doi.org/10.1371/journal.pone.0160655.s004

(DOCX)

S4 Table. Results for an analysis of molecular variance (AMOVA) for H. wrightorum.

https://doi.org/10.1371/journal.pone.0160655.s005

(DOCX)

S5 Table. GeneClass2 individual assignment probabilities for H. wrightorum.

https://doi.org/10.1371/journal.pone.0160655.s006

(DOCX)

S6 Table. Delta K calculations and log likelihoods for Structure output for H. wrightorum with LOCPRIOR.

https://doi.org/10.1371/journal.pone.0160655.s007

(DOCX)

S7 Table. Delta K calculations and log likelihoods for Structure output for H. wrightorum without LOCPRIOR.

https://doi.org/10.1371/journal.pone.0160655.s008

(DOCX)

S8 Table. Sources, descriptions, and details for spatial data used in landscape genetic analyses for H. wrightorum.

https://doi.org/10.1371/journal.pone.0160655.s009

(DOCX)

S9 Table. Pearson correlations coefficiens between all distance matrices (genetic and spatial).

https://doi.org/10.1371/journal.pone.0160655.s010

(DOCX)

S10 Table. Landscape resistance matrices generated using CIRCUITSCAPE.

https://doi.org/10.1371/journal.pone.0160655.s011

(DOCX)

Acknowledgments

We would like to extend our sincerest gratitude to the following funding sources and people for making this work possible. Funding for this research was provided by the Hall Conservation Genetics Graduate Research Award at the University of Washington’s College of the Environment and by the U.S. Army Award Number W9214A-14-P-0048. Debbie Brewer (US Army, Fort Huachuca), Sheridan Stone (US Army, Fort Huachuca), John Kraft (US Forest Service), and Tom Jones (Arizona Game and Fish Department) all provided logistical and planning support for this work. Katherine Strickler, Kody Cochrell, Trevor Eakes, and Jessica Hale assisted in fieldwork and sampling.

Author Contributions

  1. Conceived and designed the experiments: MCM LH CSG JDO.
  2. Performed the experiments: MCM.
  3. Analyzed the data: MCM.
  4. Contributed reagents/materials/analysis tools: MCM LH CSG JDO.
  5. Wrote the paper: MCM LH CSG JDO.

References

  1. 1. Williams SE, Shoo LP, Isaac JL, Hoffmann AA, Langham G (2008) Towards an integrated framework for assessing the vulnerability of species to climate change. PloS Biology 6:e325.
  2. 2. Allendorf F, Luikart G (2007) Conservation and the genetics of populations. Blackwell Publishing, Malden MA.
  3. 3. Beebee TJC (2005) Conservation genetics of amphibians. Heredity 95:423–427. pmid:16106261
  4. 4. Sodhi NS, Bickford D, Diesmos AC, Lee TM, Koh LP, Brook BW, et al. (2009) Measuring the meltdown: drivers of global amphibian extinction and decline. PLoS one 3:e1636.
  5. 5. Stuart SN, Chanson JS, Cox NA, Young BE, Rodrigues ASL, Fischman DL, et al. (2004) Status and trends of amphibian declines and extinctions worldwide. Science 306:1783–1736. pmid:15486254
  6. 6. Fish U.S. and Wildlife (2013) U.S. Fish and Wildlife Service species assessment and listing assignment form: Hyla wrightorum. 36 pp.
  7. 7. Gergus EWA (2000) Demography, population genetics, and systematics of Huachuca tree frogs (Hyla wrightorum): taxonomic and conservation implications. Arizona Game and Fish Department Heritage Grant Project #I96047, IIPAM Project Report. Tempe, AZ.
  8. 8. Sredl MJ, Wallace JE, Miera V (2000) Aquatic herpetofauna inventory of Fort Huachuca and vicinity. Phoenix, Arizona.
  9. 9. Duellman WE (2001) The Hylid frogs of Middle America. Society for the Study of Amphibians and Reptiles, St. Louis, MO.
  10. 10. Brennan TC, Holycross AT (2006) A field guide to amphibians and reptiles in Arizona. Arizona Game and Fish Department, Phoenix, Arizona.
  11. 11. Gergus EWA, Reeder TW, Sullivan BK (2004) Geographic variation in Hyla wrightorum: advertisement calls, allozymes, mtDNA, and morphology. Copeia 2004:758–769.
  12. 12. Bradley GA, Rosen PC, Sredl MJ, Jones TR, Longcore JE (2002) Chytridiomycosis in native Arizona frogs. J Wildl Dis 38:206–212. pmid:11838218
  13. 13. Jones TR, Timmons RJ (2010) Hyla wrightorum (Arizona treefrog): predation. Herpetol Rev 41:473–474.
  14. 14. Hale JR, Mims MC, Bogan MT, Olden JD (2015) Links between two interacting factors, novel habitats and non-native predators, and aquatic invertebrate communities in a dryland environment. Hydrobiologia 746:313–326.
  15. 15. Westerling AL, Hidalgo HG, Cayan DR, Swetnam TW (2006) Warming and earlier spring increase western U.S. forest wildfire activity. Science 313:940–943. pmid:16825536
  16. 16. Seager R, Ting M, Held I, Kushnir Y, Lu J, Vecchi G, et al. (2007) Model projections of an imminent transition to a more arid climate in southwestern North America. Science 316:1181–1184. pmid:17412920
  17. 17. Seager R, Ting M, Li C, Naik N, Cook B, Nakamura J, et al. (2012) Projections of declining surface-water availability for the southwestern United States. Nat Clim Chang 3:482–486.
  18. 18. Chadès I, McDonald-Madden E, McCarthy MA, Wintle B, Linkie M, Possingham HP. (2008) When to stop managing or surveying cryptic threatened species. Proc Natl Acad Sci 105:13936–13940. pmid:18779594
  19. 19. Mims MC, Phillipsen IC, Lytle DA, Hartfield Kirk EE, Olden JD (2015) Ecological strategies predict associations between aquatic and genetic connectivity for dryland amphibians. 96:1371–1382.
  20. 20. Storfer A, Mech SG, Reudink MW, Lew K (2014) Inbreeding and strong population subdivision in an endangered salamander. Conserv Genet 15:137–151.
  21. 21. Wright S (1943) Isolation by distance. Genetics 28:114–138. pmid:17247074
  22. 22. Goldberg CS, Kaplan ME, Schwalbe CR (2003) From the frog’s mouth: buccal swabs for collection of DNA from amphibians. Herpetol Rev 34:220–221.
  23. 23. Matschiner M, Salzburger W (2009) TANDEM: integrating automated allele binning into genetics and genomics workflows. Bioinformatics 25:1982–3. pmid:19420055
  24. 24. Raymond M, Rousset F (1995) GENEPOP (version 1.2): population genetics software for exact tests and ecumenicism. J Hered 86:248–249.
  25. 25. Van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P (2004) MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Mol Ecol Notes 4:535–538.
  26. 26. Goldberg CS, Waits LP (2010) Quantification and reduction of bias from sampling larvae to infer population and landscape genetic structure. Mol Ecol Resour 10:304–313. pmid:21565025
  27. 27. Jones O, Wang J (2009) COLONY: a program for parentage and sibling inference from multilocus genotype data. Mol Ecol Resour 10:551–555. pmid:21565056
  28. 28. Weir BS, Cockerham CC (1984) Estimating F-statistics for the analysis of population structure. Evolution 38:1358–1370.
  29. 29. Hedrick PW (2005) A standardized genetic differentiation measure. Evolution 59:1633–1638. pmid:16329237
  30. 30. Excoffier L, Laval G, Schneider S (2005) Arlequin ver. 3.1: An integrated software package for population genetics data analysis. Evol Bioinf Online 1:47–50.
  31. 31. Bowcock AM, Ruiz-Linares A, Tomfohrde J, Minch E, Kidd JR, Cavalli-Sforza LL (1994) High resolution human evolutionary trees with polymorphic microsatellites. Nature 368:455–457. pmid:7510853
  32. 32. Slatkin M (1995) A measure of population subdivision based on microsatellite allele frequencies. Genetics 139:457–462. pmid:7705646
  33. 33. Dieringer D, Schlötterer C (2003) Microsatellite analyser (MSA): a platform independent analysis tool for large microsatellite data sets. Mol Ecol Notes 3:167–169.
  34. 34. Waples RS, Do C (2008) LDNe: a program for estimating effective population size from data on linkage disequilibrium. Molecular Ecology Resoruces 8:753–756.
  35. 35. Do C, Waples RS, Peel D, Macbeth GM, Tillett BJ, Ovenden JR (2014) NeEstimator v2: re-implementation of software for the estimation of contemporary effective population size (Ne) from genetic data. Mol Ecol Resour 14:209–214. pmid:23992227
  36. 36. Gilbert KJ, Whitlock MC (2015) Evaluating methods for estimating local effective population size with and without migration. Evolution 59: 2154–2166.
  37. 37. Waples RS, England PR (2011) Estimating contemporary effective population size on the basis of linkage disequilibrium in the face of migration. Genetics 189:633–644. pmid:21840864
  38. 38. Holleley CE, Nichols RA, Whitehead MR, Adamack AT, Gunn MR, Sherwin WB (2014) Testing single-sample estimators of effective population size in genetically structured populations. Conserv Genet 15:23–35.
  39. 39. Piry S, Luikart G, Cornuet JM (1999) BOTTLENECK: a computer program for detecting recent reductions in the effective population size using allele frequency data. J Hered 90:502–503.
  40. 40. Luikart G, Allendorf FW, Cornuet JM, Sherwin WB (1998) Distortion of allele frequency distributions provides a test for recent population bottlenecks. J Hered 89:238–247. pmid:9656466
  41. 41. Piry S, Alapetite A, Corneut JM, Paetkau D, Baudouin L, Estoup A (2004) GENECLASS2: A software for genetic assignment and first-generation migrant detection. J Hered 95:536–539. pmid:15475402
  42. 42. Rannala B, Mountain JL (1997) Detecting immigration by using multilocus genotypes. Proc Natl Acad Sci USA 94:9197–9201. pmid:9256459
  43. 43. Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics 155:945–59. pmid:10835412
  44. 44. Chan LM, Zamudio KR (2009) Population differentiation of temperate amphibians in unpredictable environments. Mol Ecol 18:3185–3200. pmid:19573030
  45. 45. Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol 14:2611–20. pmid:15969739
  46. 46. Spear SF, Crisafulli CM, Storfer A (2012) Genetic structure among coastal tailed frog populations at Mount St. Helens is moderated by post-disturbance management. Ecol Appl 22:856–869. pmid:22645816
  47. 47. Pritchard JK, Wen X, Falush D (2010) Documentation for structure software: version 2.3. Available: http://pritchardlab.stanford.edu/structure.html.
  48. 48. Rosenberg NA (2004) Distruct: a program for the graphical display of population structure. Mol Ecol Notes 4:137–138.
  49. 49. McRae BH (2006) Isolation by resistance. Evolution 60:1551–1561. pmid:17017056
  50. 50. Cushman SA, McKelvey KS, Hayden J, Schwartz MK (2006) Gene flow in complex landscapes: testing multiple hypotheses with causal modeling. Am Nat 168:486–499. pmid:17004220
  51. 51. Mantel N (1967) The detection of disease clustering and a generalized regression approach. Cancer Research 27:209–220. pmid:6018555
  52. 52. Cushman SA, Landguth EL (2010) Spurious correlations and inference in landscape genetics. Mol Ecol 19:3592–3602. pmid:20618896
  53. 53. Van Strien MJ, Keller D, Holderegger R (2012) A new analytical approach to landscape genetic modelling: least-cost transect analysis and linear mixed models. Mol Ecol 21:4010–4023. pmid:22738667
  54. 54. Akaike H (1973) Information theory and an extension of the maximum likelihood principle. Proc 2nd Inter Symposium on Information Theory. Budapest. 267–281.
  55. 55. Burnham KP, Anderson DR (2002) Model selection and multimodel inference: a practical information-theoretic approach. 2nd ed. Springer, New York, NY.
  56. 56. Schwartz GE (1978) Estimating the dimension of a model. Annals of Statistics 6:461–464.
  57. 57. Edwards LJ, Muller KE, Wolfinger RD, Qaqish BF, Schabenberger O (2008) An R2 statistic for fixed effects in the linear mixed model. Stat Med 27:6137–6157. pmid:18816511
  58. 58. Bates D, Maechler M, Bolker B, Walker S (2014) lme4: Linear mixed-effects models using Eigen and S4. R package version 1.1–7. Available: http://CRAN.R-project.org/package=lme4.
  59. 59. Halekoh U, Højsgaard S (2012) PBKRTEST: parametric bootstrap and Kenward Roger based methods for mixed model comparison. R package version 0.3–2. Available: http://CRAN.R-project.org/package=pbkrtest.
  60. 60. Marsh DM, Trenham PC (2001) Metapopulation dynamics and amphibian conservation. Cons Bio 15:40–49.
  61. 61. Collins JP (1994) Final report: a status of three species of endangered/sensitive amphibians in Arizona. Report to the Arizona Game and Fish Department, Heritage Fund IIPAM #I92014. Phoenix, Arizona.
  62. 62. Rosen PC, Schwalbe CR (1998) Using managed waters for conservation of threatened frogs in Environmental, economic, and legal issues related to rangeland water developments. Arizona State University, College of Law, Tempe, AZ, pp 180–202.
  63. 63. Phillipsen IC, Lytle DA (2013) Aquatic insects in a sea of desert: population genetic structure is shaped by limited dispersal in a naturally fragmented landscape. Ecography 36:731–743.
  64. 64. Zweifel R (1961) Larval development of the treefrogs Hyla arenicolor and Hyla wrightorum. American Museum Novitates 2056:1–20.
  65. 65. Stebbins R (2003) Western reptiles and amphibians. Houghton Mifflin. Singapore.
  66. 66. Cushman SA, Shirk AJ, Landguth EL (2013) Landscape genetics and limiting factors. Conserv Genet 14:263–274.
  67. 67. Hedrick PW (2005) Large variance in reproductive success and the Ne/N ratio. Evolution 59 1596–1599. pmid:16153045
  68. 68. Beebee TJC, Griffiths RA (2005) The amphibian decline crisis: a watershed for conservation biology? Biol Conserv 125:271–285.
  69. 69. Palstra FP, O’Connell MF, Ruzzante DE (2007) Population structure and gene flow reversals in Atlantic salmon (Salmo salar) over contemporary and long-term temporal scales: effects of population size and life history. 16:4504–4522.
  70. 70. Berger L, Speare R, Daszak P, Green DE, Cunningham AA, Goggin CL, et al. (1998) Chytridiomycosis causes amphibian mortality associated with population declines in the rain forests of Australia and Central America. Proc Natl Acad Sci 95:9031–9036. pmid:9671799
  71. 71. Game Arizona and Department Fish (2006) Rana subaquavocalis. Unpublished abstract compiled and edited by the Heritage Data Management System. Phoenix, AZ.
  72. 72. Sredl MJ, Field KJ, Peterson AM (2003) Understanding and mitigating effects of Chytrid fungus to amphibian populations in Arizona. Phoenix, AZ.
  73. 73. Gervasi SS (2013) All hosts are not created equal: variation in amphibian responses to an emerging fungal pathogen and why it matters. PhD Dissertation, Oregon State University.
  74. 74. Pounds JA, Bustamante MR, Coloma LA, Consuegra JA, Fogden MPL, Foster PN, et al. (2006) Widespread amphibian extinctions from epidemic disease driven by global warming. Nature 439:161–167. pmid:16407945