Landscape genetics reveals unique and shared effects of urbanization for two sympatric pool-breeding amphibians.

Metapopulation-structured species can be negatively affected when landscape fragmentation impairs connectivity. We investigated the effects of urbanization on genetic diversity and gene flow for two sympatric amphibian species, spotted salamanders (Ambystoma maculatum) and wood frogs (Lithobates sylvaticus), across a large (>35,000 km2) landscape in Maine, USA, containing numerous natural and anthropogenic gradients. Isolation-by-distance (IBD) patterns differed between the species. Spotted salamanders showed a linear and relatively high variance relationship between genetic and geographic distances (r = .057, p < .001), whereas wood frogs exhibited a strongly nonlinear and lower variance relationship (r = 0.429, p < .001). Scale dependence analysis of IBD found gene flow has its most predictable influence (strongest IBD correlations) at distances up to 9 km for spotted salamanders and up to 6 km for wood frogs. Estimated effective migration surfaces revealed contrasting patterns of high and low genetic diversity and gene flow between the two species. Population isolation, quantified as the mean IBD residuals for each population, was associated with local urbanization and less genetic diversity in both species. The influence of geographic proximity and urbanization on population connectivity was further supported by distance-based redundancy analysis and multiple matrix regression with randomization. Resistance surface modeling found interpopulation connectivity to be influenced by developed land cover, light roads, interstates, and topography for both species, plus secondary roads and rivers for wood frogs. Our results highlight the influence of anthropogenic landscape features within the context of natural features and broad spatial genetic patterns, in turn supporting the premise that while urbanization significantly restricts interpopulation connectivity for wood frogs and spotted salamanders, specific landscape elements have unique effects on these two sympatric species.

INTRODUCTION

Landscape alterations that accompany increases in human population density (i.e., urbanization) commonly influence ecological and evolutionary processes. For instance, landscape fragmentation and habitat loss can lead to reduced connectivity among wildlife populations, consequently disrupting demographic support from metapopulation dynamics that would otherwise improve population stability via immigration (Andrén, 1994; Crosby, Licht, & Fu, 2008). Reductions in connectivity are potentially associated with reduced levels of gene flow and genetic diversity (Crawford, Peterman, Kuhns, & Eggert, 2016; Ortego, Aguirre, Noguerales, & Cordero, 2015) and an increased risk of extirpation from lost demographic support (Bascompte & Sole, 1996; Haddad et al., 2015; Wilcox & Murphy, 1985). Populations that do persist may be subject to decreased fitness associated with inbreeding depression (Andersen, Fog, & Damgaard, 2004; Lopez, Rousset, Shaw, Shaw, & Ronce, 2009). These issues are of immediate concern as urban areas are becoming larger and more prevalent worldwide. For instance, growth in global human population size (UNDESA, 2012, 2015) has been accompanied by increases in the percentage of people living in urban areas from 29.4% in 1950 to 52.1% in 2011 and up to 70% anticipated by 2050 (UNDESA, 2012). During the past century, cities have also become increasingly diffuse, leading to a greater proportion of the landscape being affected by their growth (Seto, Sánchez‐Rodríguez, & Fragkias, 2010; Theobald, 2010). Given these trends, the ecological and evolutionary effects of urbanization on wildlife are likely to intensify in the coming decades.

Metapopulation‐structured species are especially vulnerable to negative effects of habitat fragmentation associated with urbanization (Graham, Haines‐Young, & Field, 2017). These species are spatially arranged in discrete subpopulations that are spread across a heterogeneous landscape and are dependent to a degree on dispersal among constituent subpopulations (Hanski, 1998). In classical metapopulation theory, subpopulations experience bouts of extinction and recolonization while maintaining overall metapopulation stability (Hanski, 1991; Levins, 1969). A loss of connectivity among subpopulations has been predicted to result in local extirpation and abundance declines for metapopulations (Grilli, Barabás, & Allesina, 2015; Reigada, Schreiber, Altermatt, & Holyoak, 2015; Schnell, Harris, Pimm, & Russell, 2013), a result that has been observed in several natural systems (Crooks et al., 2017). In some cases, habitat fragmentation has allowed evolutionary forces to generate changes in phenotypic traits, influencing characteristics such as dispersal propensities (Cheptou, Hargreaves, Bonte, & Jacquemyn, 2017) or life history traits (De Roissart, Wybouw, Renault, Leeuwen, & Bonte, 2016), as well as broader eco‐evolutionary processes (Fronhofer & Altermatt, 2017). The metapopulation concept has proven broadly applicable in urbanization‐associated habitat fragmentation scenarios, facilitating an improved understanding of how landscape alterations can affect regional population processes for many taxa, such as birds (Millsap, 2018; Padilla & Rodewald, 2015) and amphibians (Cox, Maes, Calster, & Mergeay, 2017; Hale et al., 2013; Heard, McCarthy, Scroggie, Baumgartner, & Parris, 2013).

Landscape genetics provides a framework for evaluating support for potential environmental correlates of observed interpopulation structure, allowing for the generation of inferences of spatially explicit drivers of gene flow. This approach has revealed allele frequency changes associated with diminished gene flow due to habitat fragmentation (Epps & Keyghobadi, 2015; Zellmer & Knowles, 2009). Studies focused on disentangling effects of natural versus anthropogenic landscape elements, and using multiple species to obtain a more comprehensive understanding of connectivity, have been highlighted as critical areas for advancing landscape genetics research (Manel & Holderegger, 2013; Richardson, Brady, Wang, & Spear, 2016). However, detecting landscape genetic effects of anthropogenic fragmentation is not without challenges. Unlike many natural landscape features present since deglaciation, the widespread appearance of anthropogenic features has occurred relatively recently, limiting the opportunity for drift or gene flow to affect gene frequencies, and placing constraints on statistical power. This necessitates sampling adequate numbers of populations at geographic scales where disruptions to gene flow are most apt to affect background patterns of drift–migration equilibrium.

Habitat losses, such as the landscape changes coincident with urbanization, have been singled out as a leading threat to amphibian species, with empirical evidence mounting for population declines associated with changes in land cover (Price, Dorcas, Gallant, Klaver, & Willson, 2006), canopy cover (Clark, Reed, Tavernia, Windmiller, & Regosin, 2008), and roadways (Andrews, Gibbons, & Jochimsen, 2008). Ecological studies have consistently suggested that pool‐breeding amphibians may be particularly susceptible to negative effects of landscape fragmentation (Baldwin & de Maynadier, 2009; Semlitsch, 2003). For instance, because they require access to both wetland and upland environments to complete their semiaquatic life cycle (Semlitsch, 2008), any barriers between those two environments could impair a population (Homan, Windmiller, & Reed, 2004). Additionally, the ability to occasionally disperse among populations is important for these metapopulation‐structured amphibian species due to highly variable interannual recruitment success at the local scales of pools (Baldwin, Calhoun, & de Maynadier, 2006b; Green, Hooten, Grant, & Bailey, 2013). Tracking studies suggest altered habitats between pools can reduce overall dispersal propensities (Cline & Hunter, 2014, 2016). Despite these observed ecological effects, signals of urbanization‐related influences on the genetic structure of pool‐breeding amphibians have not been consistently detected (e.g., Coster, Babbitt, Cooper, & Kovach, 2015; Peterman et al., 2015; Richardson, 2012; Table 1). In some of these cases, a lack of effect could be associated with the surveyed spatial extents that were relatively small and that may not be commensurate with the scales over which drift–migration equilibrium is most disrupted and detectable.

Table 1

Studies of the effects of urban landscape elements (e.g., roads and developed lands) on connectivity of pool‐breeding amphibians based on microsatellite genotyping. Negative (↘) and negligible effects (↔) are noted

Species	Approx. study area (km2)	No. populations	Effect	Citation
Blanchard's cricket frogs (Acris blanchardi)	2,320	28	↘	Youngquist, Inoue,Berg, and Boone (2017)
Columbia spotted frogs (Rana luteiventris)	213	8	↘	Goldberg and Waits (2010)
Eastern tiger salamander (Ambystoma tigrinum)	220	26	↔	Titus, Bell, Becker, and Zamudio (2014)
Long‐toed salamanders (Ambystoma macrodactylum)	213	4	↘	Goldberg and Waits (2010)
Natterjack toad (Epidalea calamita)	40	23	↘	Cox et al. (2017)
Ringed salamander (Ambystoma annulatum)	35	20	↔	Peterman et al. (2015)
Spotted salamander (Ambystoma maculatum)	35	23	↔	Peterman et al. (2015)
Spotted salamander (Ambystoma maculatum)	2,080	23	↘	Coster, Babbitt, Cooper, et al. (2015)
Spotted salamander (Ambystoma maculatum)	21,000	22	↘	Richardson (2012)
Wood frog (Lithobates sylvaticus)	200	9	↔	Peterman et al. (2013)
Wood frog (Lithobates sylvaticus)	375	65	↔	Gabrielsen, Kovach, Babbitt, and McDowell (2013)
Wood frog (Lithobates sylvaticus)	2,080	20	↔	Coster, Babbitt, Cooper, et al. (2015)
Wood frog (Lithobates sylvaticus)	21,000	22	↘	Richardson (2012)

We investigated the effects of urbanization on metapopulation processes by examining the landscape genetics of two metapopulation‐structured pool‐breeding amphibians, spotted salamanders (Ambystoma maculatum) and wood frogs (Lithobates sylvaticus), in an area of overlap between their native ranges. Several characteristics make these species excellent subjects for evaluating the effects of urbanization with a landscape genetics approach. For instance, the species have small annual home ranges (spotted salamanders: up to 301 m2, Ousterhout and Burkhart, 2017; wood frogs: up to 32,165 m2, Blomquist & Hunter, 2010; Groff, Calhoun & Loftin, 2017), relatively short generation times (spotted salamanders: maturity in 2–7 years, Flageole & Leclair, 1992; wood frogs: maturity in 2–3 years, Sagor, Ouellet, Barten, & Green, 1998), and high rates of philopatry (Vasconcelos & Calhoun, 2004). Both species are vulnerable to degradation of the breeding sites they share within their overlapping ranges in the northeastern United States (Harper, Rittenhouse, & Semlitsch, 2008). These two amphibians have several differentiating life history and behavioral attributes that likely affect how urbanization influences their interpopulation dynamics. For instance, wood frogs tend to be shorter lived, have larger home ranges, and are more vagile than spotted salamanders (Berven & Grudzien, 1990; Madison, 1997; Semlitsch, 1998). Due to these characteristics, we expected the magnitude and dynamics of urbanization‐related effects to differ between the species. We tested three sets of hypotheses to examine effects of urbanization on individual populations, interpopulation dynamics, and genetic structure more broadly.

Broadscale patterns of genetic structure: We hypothesized broadscale patterns of genetic structure will be much stronger for wood frogs than spotted salamanders. Previous work to characterize the isolation‐by‐distance (IBD) relationships for these two species has supported this hypothesis through observations of clear positive correlations between geographic and genetic distances for wood frogs (Crosby et al., 2008; Peterman, Feist, Semlitsch, & Eggert, 2013; Richardson, 2012; Squire & Newman, 2002) and either high variance positive correlations (Burkhart et al., 2017; Peterman et al., 2015; Richardson, 2012; Zamudio & Wieczorek, 2007) or nonsignificant relationships (Purrenhage, Niewiarowski, & Moore, 2009; Whiteley, McGarigal, & Schwartz, 2014) for spotted salamanders. Additionally, we hypothesize that IBD relationships within species are not absolute, but are instead scale‐dependent, such that the strongest correlations between genetic isolation and distance will occur at some intermediate geographic scales of analysis. Below these scales, the strength of IBD is expected to be weaker due to the highly variable patterns of local dispersal and gene flow overwhelming drift, and limited sample size. Above these scales, sample sizes are greater, but drift–migration equilibrium may take longer to occur, and stochastic processes weaken correlations. Understanding the scale dependence of IBD may be important to understanding why different studies on a specific species obtain different IBD inferences and help to inform whether an IBD study has sufficient power to identify effects of landscape features on population connectivity.

Urbanization's influence on isolation and genetic diversity: We hypothesized that greater urbanization will affect metapopulation dynamics by increasing among‐population isolation and reducing within‐population genetic diversity for both species, relative to patterns in less‐urbanized areas (Frankham, 2015; Kenney, Allendorf, Mcdougal, & Smith, 2014; Pavlova et al., 2017). We further hypothesize that the substantive differences in life history and vagility of our focal species will contribute to different effects of the same landscape features (Moyle, 2006; Phillipsen et al., 2015). Wood frogs are often more mobile than spotted salamanders and may be more capable of traversing urbanization features than salamanders. However, this increased mobility may increase wood frog exposure to deleterious features of the broader urbanized landscape relative to spotted salamanders.

Landscape features influencing gene flow: We predict signals of diminished gene flow associated with anthropogenic landscape features will be present; however, they will be weaker than signals from natural features (e.g., rivers, elevation gradients). Landguth et al. (2010) used simulations to illustrate a lag time of approximately 200 generations when using a landscape genetics approach to detect signals of barriers, although they note the magnitude of this effect likely varies with each population's effective size. The influence of lag time is also supported by empirical research that considers effects of both natural and anthropogenic features on amphibian gene flow (Garcia, Ivy, & Fu, 2017; Peterman et al., 2015; Richardson, 2012). Due to this lag, most natural features will likely have stronger signals, such as rivers appearing more important than interstate highways. However, given the magnitude of expected effects for heavily urbanized areas, we anticipate that the landscape features that coincide with the most intense urbanization (e.g., road density) will be the strongest predictors of isolation.

MATERIALS AND METHODS

Sample collection

Larval and embryonic wood frogs and spotted salamanders were collected from vernal pools throughout Maine, USA (Figure 1). We sampled across a 35,000‐km2 region that contains several natural (e.g., topographic) and anthropogenic (e.g., urbanization) landscape gradients. Sampling effort was concentrated in areas of urbanization to ensure adequate power to detect effects of relatively recent urban‐associated landscape features (Balkenhol, Cushman, Waits, & Storfer, 2016) with additional sampling conducted away in rural areas to facilitate appropriate comparisons to natural landscape features. When a site (treated throughout the study as a population) was visited prior to egg hatching, we sampled up to 40 egg masses, collecting one embryonic individual from each mass to reduce the likelihood of sampling siblings because the inclusion of closely related individuals has been found to bias some genetic analyses (Goldberg & Waits, 2010, Rodríguez‐Ramilo & Wang, 2012, Peterman, Brocato, Semlitsch, & Eggert, 2016, Wang, 2018, but see Waples and Anderson (2017)). If larvae were free‐swimming upon sampling, a small dipnet was used to collect individuals from throughout the pool and full siblings were later removed based on sibship analyses (see below). Sampling occurred during April and May 2014, 2015, and 2016. When fewer than 25 individuals were collected at a site in one year, we returned to sample in the subsequent year.

Figure 1

Locations of 90 wood frog and 87 spotted salamander vernal pool sampling sites in Maine, USA. Inset maps illustrate densely sampled regions around Portland (thin outline) and Bangor, Maine (heavy outline). Solid black line in largest extent maps represents interstate highways. Red background coloration indicates high levels of nighttime light intensity based on NASA Visible Infrared Imaging Radiometer Suite (VIIRS) data and is provided as a proxy for human population density

Genetic data collection and quality control

Genomic DNA was extracted from whole embryonic or larval individuals using Qiagen DNeasy kits following the manufacturer's instructions. We analyzed variability at 10 microsatellite loci to evaluate spatial genetic structure for each species. PCR components, thermal cycler profiles, and citations for loci primer sequences are described in Appendix S1. Negative controls were included in each 96‐well PCR to allow for detection of reagent contamination. Microsatellite fragment analysis was conducted using an ABI 3730 automated genetic analyzer (Applied Biosystems, Inc.) at the University of Maine DNA Sequencing Facility. Genotyping was performed using Geneious v7.1.9 with fragment sizes based on GeneScan 500 LIZ Size Standard (Applied Biosystems) and all allele calls confirmed visually. A random 10% of individuals were genotyped a second time to evaluate genotype error rates.

A series of data filtering steps was performed to reduce the potential influence of sampling bias and to ensure conformance to assumptions of population genetic analyses. First, individuals with fewer than five successfully amplified loci were removed. Peterman et al. (2016) found five microsatellite loci to be as informative as both 10 and 15 loci for estimating heterozygosity and allelic richness in other spotted salamander populations. Next, to reduce the likelihood of mischaracterizing allele frequencies due to small sample sizes, we eliminated sites with fewer than ten individuals successfully genotyped. Finally, we performed sibship reconstruction for all individuals sampled at each site using COLONY (v2.0.5.9; Jones & Wang, 2010; Wang, 2004) and subsequently haphazardly removed all but one individual from any apparent full‐sibling family. COLONY analyses assumed polygamy in both sexes, no inbreeding, and were performed using a long run with the full likelihood method. In addition to minimizing the degree of family structure present in our sample set, this post hoc removal of siblings improves congruence in sampling design between populations sampled at the egg stage and those sampled as free‐swimming larvae where inadvertently collecting siblings is more likely.

We estimated the frequency of null alleles for each locus and tested for Hardy–Weinberg equilibrium for each locus–sampling site combination using PopGenReport (Adamack & Gruber, 2014) in R v3.4.1 (R Core Team, 2016). Independent sorting of genotypes (i.e., linkage disequilibrium) was evaluated using exact testing in Arlequin v3.5.1.2 (Excoffier & Lischer, 2010). Alpha levels to determine statistical significance for tests of Hardy–Weinberg proportions and independent sorting of genotypes were adjusted using the false discovery rate (FDR) approach of Benjamini and Hochberg (1995) based upon a 0.05 alpha level.

Genetic diversity and differentiation

We quantified genetic diversity within each site and genetic differentiation among sites using multiple measures. Average number of alleles per site (AO) was estimated using PopGenReport, and allelic richness (i.e., allelic counts rarefied based on smallest sample size per species; spotted salamander: 10, wood frog: 12; AR), expected heterozygosity (H E), and Wright's inbreeding coefficient (F IS) were estimated using the R package hierfstat (Goudet & Jombart, 2015). Genetic differentiation was calculated using G ST (i.e., Nei, 1973; Nei & Chesser, 1983) and G″ ST (Meirmans & Hedrick, 2011). G ST (commonly reported as F ST) summarizes the amount of diversity contained among populations relative to the diversity of all populations combined (Nei, 1973), whereas G″ ST provides a scaled maximum value of G ST based on the genetic diversity within a measured population (Meirmans & Hedrick, 2011). Both G ST and G″ ST were estimated using the R package mmod (Winter, 2012). Statistical significance of pairwise population differentiation was evaluated with an exact G test implemented using the genetic differentiation option in Genepop v4.2 (Raymond & Rousset, 1995; Rousset, 2008) with a FDR correction for type I error rates.

The spatial arrangement of effective genetic diversity was visualized using estimated effective migration surfaces (EEMS; Petkova, Novembre, & Stephens, 2016). Effective genetic diversity reflects the expected genetic dissimilarity of two individuals sampled within each deme assuming a generally IBD‐driven system and a stepping‐stone dispersal pattern (Petkova et al., 2016). EEMS constructs a dense, regular grid across the study range and assigns sampling sites to the nearest grid intersection (node), often resulting in a set of fewer demes than the actual number of sampling sites. Diversity values are then interpolated among the demes to create a continuous surface for visualizing spatial patterns. Our starting grid provided 500 potential nodes for deme assignment, of which 462 were incorporated into the analysis due to the irregular landscape boundaries. Previous work has demonstrated EEMS results to be qualitatively robust when various numbers of nodes were used in analyses (Petkova et al., 2016). EEMS analysis parameters were adjusted to achieve the recommended 20%–40% acceptance rates before running the analysis using 1 × 107 iterations, a burn‐in period of 1 × 106 iterations, and a thinning interval of 1 × 103 (Combs, Puckett, Richardson, Mims, & Munshi‐South, 2017; Petkova et al., 2016). All EEMS plotting was performed using rEEMSplots R package (Petkova et al., 2016).

We also evaluated all populations for the presence of bottlenecks that may be associated with urbanization using the program Bottleneck (Cornuet & Luikart, 1996; Piry, Luikart, & Cornuet, 1999). We used the two‐phase model of microsatellite mutation (TPM; Di Rienzo et al., 1994) with variance set to 12 and the probability of single‐step mutations set to 95% as recommended by Piry et al. (1999). Significance was evaluated using a one‐tail Wilcoxon test with an FDR‐adjusted alpha level.

Isolation by distance

To examine IBD relationships, we compared each pairwise measure of genetic differentiation with between‐site Euclidean geographic distance. Genetic differentiation measures were linearized (G ST/(1 − G ST)) as suggested by Slatkin (1995), and geographic distances were measured as straight‐line Euclidean distances using the “distm” function of the R package geosphere (Hijmans, 2017). We examined relationships between linearized G ST and both log‐transformed and nontransformed geographic distances. The slope of IBD relationships based on log transformation of geographic distance is useful for understanding the dispersal kernel relationships in scenarios of two‐dimensional movements (Rousset, 2000), whereas nontransformed distances are helpful for understanding broadscale patterns of IBD (sensu Hutchison & Templeton, 1999). While we provide the slope of the IBD relationship based on log‐transformed geographic distances, it is important to note that our study design does not meet the assumptions for actually estimating dispersal kernel size per se (e.g., sampling extent greater than 0.56σ/√2μ, where σ is the parent–offspring axial distance and μ is the mutation rate of the loci; Rousset, 2004), and key parameters are unknown (i.e., D, the effective density). Therefore, the regression slopes we report should be considered a broad approximation of observed increases in genetic differentiation with geographic distance (Dσ 2) and useful only for comparisons between species within this specific study. Associations between the distance matrices were tested using regression and Mantel tests (Mantel, 1967) that were implemented in the R package vegan (Oksanen et al., 2017). We evaluated Mantel tests for significance based on 9,999 permutations.

We examined the relationship between genetic and geographic distances as a function of the spatial scale of analysis using two methods. First, we constructed a Mantel correlogram (Borcard & Legendre, 2012; Legendre & Legendre, 2012; Oden & Sokal, 1986) to quantify the strength of the relationship between genetic and geographic distances within various distance classes using the “mantel.correlog” function in vegan. Distance class breakpoints were placed every 20 km, and larger distance classes that did not contain every sampling site were omitted to avoid bias (Wagner et al., 2005). Statistical significance of correlations was assessed using 10,000 permutations and a FDR correction based on a 0.05 alpha level. Next, we estimated the slope (β) of variable intercept IBD regressions that were performed repeatedly using expanding datasets based on distance between sampling sites, which generated IBD scaling profiles for each species. For example, the first iteration of the analysis was conducted using the 20 shortest pairwise geographic distances, the next iteration with the 21 shortest, and so on, until all pairwise comparisons were included. We performed 1,000 bootstrapped replicates of each regression using the “Boot” function in the R package car (Fox & Weisberg, 2011) to estimate each beta coefficient and its 95% confidence intervals. The slope confidence intervals for each IBD regression were then plotted using the maximum analyzed distance as a response variable to generate the IBD scaling profiles.

We also used EEMS to visualize spatial patterns of connectivity among sampling sites. When assessing connectivity, EEMS identifies areas with greater differentiation than expected between neighboring demes assuming a generally IBD‐driven system and a stepping‐stone dispersal pattern (Kimura & Weiss, 1964; Petkova et al., 2016). The number of effective migrants among sites is then interpolated to construct a graphic depiction of connectivity across the landscape.

Regression and multivariate analyses

We assessed the influence of urbanization on a sampling site's degree of isolation and genetic diversity. The intensity of urbanization near a site was quantified using six environmental characteristics measured in ArcGIS v10.2 (ESRI): distance to nearest roadway, percent impervious surface within one km, length of roads within 1 km for light, secondary, and primary road types, and percent canopy cover within 1 km. Road type and classification was determined using the State of Maine's NG911 Roads dataset (http://www.maine.gov/megis/catalog/, accessed Feb 18, 2016); impervious surface extent was based on the Maine Department of Inland Fisheries and Wildlife's July 2016 impervious surface dataset (J. Czapiga, Maine Department of Inland Fisheries and Wildlife, unpublished data); and percent canopy cover data were drawn from the 2011 National Land Cover Database (Homer et al., 2015). Collinearity between the urbanization‐related explanatory variables was evaluated, and one variable was selected at random to be retained from each set with a correlation coefficient exceeding 0.7.

We quantified the degree of isolation experienced at each sampling site by averaging the residuals of the IBD data points that include that site. This approach is similar to the decomposed pairwise regression analysis to detect outlier populations described by Koizumi, Yamamoto, and Maekawa (2006) and essentially provides an index of genetic differentiation corrected for geographic distance. Sites with the largest average residual values were presumed to be more isolated than those with smaller average residuals. Due to strong correlation between G ST and G″ ST for both species (spotted salamander r = .993, wood frog r = .979), we quantified isolation using only the G ST‐based IBD relationships. Genetic diversity relationships were assessed using H E and AR. We conducted three statistical analyses to test hypotheses concerning the relationship between these three factors: multiple regression between each measure of genetic diversity and all retained urbanization variables, multiple regression between urbanization variables and degree of site isolation, and simple linear regression between each measure of genetic diversity and degree of isolation.

We examined the influence of urbanization and spatial proximity to observed interpopulation genetic differentiation with two approaches complementary to the above IBD and regression analyses: distance‐based redundancy analysis (dbRDA) and multiple matrix regression with randomization (MMRR). We conducted dbRDA using the “capscale” function and examined the significance of individual model terms using 10,000 permutations with the “anova.cca” function in vegan. Our global dbRDA model included pairwise G ST values in the response matrix, the above‐described urbanization‐related metrics in an explanatory matrix, and a conditional matrix containing the latitude and longitude of each site in decimal degrees. Model terms were eliminated using a backward optimization procedure where nonsignificant terms were removed, and a simplified model was tested until all remaining terms were significant. MMRR provides a multivariate method for examining the relationships between a response matrix (e.g., interpopulation genetic divergence) and multiple explanatory matrices (e.g., environmental characteristics) while accounting for interpopulation geographic distances (Wang, 2013). With data included in the dbRDA global model, we implemented MMRR using the “lgrMMRR” function in PopGenReport, which involved 10,000 permutations to allow statistical significance to be evaluated based on the pseudo‐t statistic of Legendre, Lapointe, and Casgrain (1994).

Landscape resistance modeling

We tested support for a series of resistance surface models to determine the relative influence of ten landscape features on the genetic structuring of each species. The modeling was a two‐step process. First, we optimized resistance values for each feature, and then, we conducted generalized additive modeling to determine which features were most influential for each species. Features to be analyzed were generally selected based on data availability and previous resistance modeling for these species (Richardson, 2012). Land cover was based on the 2011 National Land Cover Database (NLCD; Homer et al., 2015) and merged into three class that generally describe forests (land cover A), open areas and agriculture (land cover B), and developed areas and open waterbodies (land cover C; Table 2; Richardson, 2012). Road data were derived from the State of Maine NG911 Roads dataset and sorted into three classes describing limited‐access interstate freeways, secondary roads (e.g., state highways), and light roads. We subset river data from the National Hydrography Dataset (USGS, https://nhd.usgs.gov/, access Feb. 18, 2016) into two classes based on the Strahler numbering system. Medium rivers included order 4 and 5 streams; large rivers included order 6 and 7 streams; and lower order waterways were not considered. Railway data were based on the Maine Department of Transportation's RailRouteSys dataset (Johnson et al., 2011). We calculated a terrain ruggedness index (TRI; Riley, DeGloria, & Elliot, 1999) using the “tri” function in the R package spatialEco (Evans, 2017) to characterize topographic heterogeneity. Raster processing was performed using ArcGIS v10.2. Processing included buffering all linear features to ensure their continuity following conversion to a raster and the resampling of all rasters to a 90‐m resolution, which was necessary given computation constraints owing to the extent of the landscape being processed.

Table 2

Reclassification of National Land Cover Database fields into three categories for use in resistance surface modeling

NLCD category	NLCD descriptions	Assigned category
41	Deciduous forest	A
42	Evergreen forest	A
43	Mixed forest	A
90	Woody wetlands	A
95	Emergent herbaceous wetlands	A
21	Developed, open space	B
52	Shrub/scrub	B
71	Grassland/herbaceous	B
81	Pasture/hay	B
82	Cultivated crops	B
11	Open water	C
22	Developed, low intensity	C
23	Developed, medium intensity	C
24	Developed, high intensity	C
31	Barren land	C

Pairwise effective resistance between each sampling site was measured based on a circuit theory approach in GFlow (Leonard et al., 2017). We conducted partial Mantel tests with 10,000 permutations using the vegan R package to evaluate correlations between pairwise effective resistance values and genetic differentiation (G ST) while controlling for effects of geographic distance between sites. The candidate resistance cost values that explained the most variation (largest R 2 value) were selected as optimal. Four to seven resistance values were tested for each landscape feature. These values were selected based upon the results of Richardson (2012) and always included a value of 1 to allow comparisons between the candidate resistance values and a simple IBD relationship. All nonfeature raster cells were assigned a value of 1 during the optimization procedure, and the terrain ruggedness index was optimized by adding various resistance values to the actual index values.

Optimized cost surfaces were used to inform a series of generalized linear additive models to assess the relative contribution of each landscape feature to overall patterns of genetic differentiation among sites for each species. We only considered landscape variables that explained genetic diversity patterns better than IBD alone. Models were compared using the small sample size‐corrected Akaike's information criterion (AICC; Burnham & Anderson, 2002). AICC compares relative support of candidate models including a penalty for the number of variables incorporated, thereby encouraging parsimony. Models with ΔAICC value <2 were considered equally supported. We used all possible combinations of included variables as candidate models and calculated AICC values and their relative weights using R package glmulti (Calcagno & Mazancourt, 2010).

RESULTS

Sampling and quality control

Due to a longer duration prior to hatching, all spotted salamanders were collected as embryos, reducing the likelihood of siblings being sampled, whereas wood frogs were occasionally collected as larvae. Therefore, sibship analyses and subsequent elimination of all but one member from each family group were performed only for the wood frogs. We sampled across multiple years for 31 spotted salamander and 27 wood frog populations. Allelic richness and expected heterozygosity did not differ depending on the number of years a site was sampled (p > .05). Resampled sites were never found to be unoccupied in any particular year. In total, we collected and genotyped 2,862 spotted salamander eggs and 2,935 wood frog eggs and larvae. Following removals of individuals based on genotype completeness, sample size, and sibship, 2,413 spotted salamanders from 90 sites and 2,439 wood frogs from 87 sites were included in our analyses (Table 3). Pairwise distances between sites ranged from 0.12 km to 320.55 km for spotted salamanders (mean = 120.32 km) and 0.07 to 337.88 km for wood frogs (mean = 120.47 km).

Table 3

Site location information, sample size, quality control, and genetic diversity information for 90 spotted salamander and 87 wood frog populations

Site name	Latitude	Longitude	Spotted salamanders							Wood frogs
			N	HW	LD	A O	AR	H E	F IS	N	HW	LD	A O	AR	H E	F IS
ALF‐1	43.50788	−70.76709	17	0	0	6.63	5.73	0.723	0.045	23	0	0	11.30	5.25	0.855	0.067
AMH‐1	44.856017	−68.415817	15	0	0	6.75	6.05	0.755	−0.049	23	0	0	11.30	5.04	0.798	0.004
AMH‐2	44.831733	−68.41081	44	0	0	7.88	5.97	0.744	0.049	40	0	2	13.00	5.06	0.805	0.058
AMH‐3	44.8625	−68.396077	33	0	0	7.63	5.67	0.704	−0.053	17	0	0	10.40	5.18	0.831	−0.004
AMH‐4	44.855217	−68.411587	29	0	0	7.50	5.95	0.757	0.014	17	0	0	10.10	5.14	0.811	0.079
BAN‐1	44.863786	−68.736084	13	0	0	6.00	5.63	0.756	0.011	–	–	–	–	–	–	–
BAN‐2	44.79044	−68.83422	18	1	1	6.63	5.66	0.726	0.047	27	0	0	11.40	4.91	0.789	−0.033
BAN‐3	44.79794	−68.83723	–	–	–	–	–	–	–	19	0	2	10.50	5.11	0.816	0.193
BAN‐5	44.8632	−68.75721	30	0	0	7.50	5.78	0.740	0.070	–	–	–	–	–	–	–
BAN‐6	44.86372	−68.75248	30	0	0	7.00	5.55	0.723	−0.007	–	–	–	–	–	–	–
BAN‐7	44.80195	−68.78882	–	–	–	–	–	–	–	26	0	1	8.10	4.25	0.745	0.041
BEN‐1	44.58715	−69.51235	23	0	0	7.13	5.72	0.693	−0.046	43	0	3	13.60	5.29	0.837	−0.009
BRU‐1	43.89274	−69.98714	68	0	0	8.00	5.71	0.716	0.058	58	1	6	13.90	5.41	0.863	0.035
BRU‐2	43.935533	−70.008217	35	0	1	6.63	5.23	0.700	0.042	28	0	0	12.60	5.24	0.851	0.023
BRU‐3	43.916783	−69.98465	26	0	0	7.63	5.79	0.733	−0.009	74	1	1	15.40	5.42	0.861	0.048
BUX‐1	43.59865	−70.469383	23	0	0	7.38	5.85	0.742	0.048	33	0	2	13.40	5.27	0.847	0.009
CUM‐1	43.819133	−70.252717	32	0	1	7.38	5.75	0.733	0.062	62	2	1	14.00	5.46	0.876	0.056
DIX‐1	44.6857	−69.134633	28	0	0	7.63	5.85	0.714	−0.011	–	–	–	–	–	–	–
EDG‐1	43.97081	−69.58153	61	0	0	8.00	5.57	0.732	0.002	27	1	3	10.90	5.02	0.834	0.056
EDM‐1	44.8853	−67.2789	–	–	–	–	–	–	–	24	0	0	10.70	4.86	0.768	0.060
ELL‐1	44.60868	−68.36858	26	0	0	6.88	5.66	0.727	0.011	26	0	2	9.70	4.73	0.777	0.016
FAL‐1	43.80073	−69.75086	23	0	0	5.88	4.98	0.698	0.134	19	0	0	8.00	4.61	0.794	−0.028
FRA‐1	44.627483	−68.904683	30	0	0	7.25	5.60	0.672	0.032	23	0	1	10.80	5.02	0.823	0.054
FRE‐1	43.856667	−70.080133	13	0	0	5.88	5.41	0.685	0.078	15	0	0	9.80	5.43	0.854	0.051
FRE‐2	43.87105	−70.110917	43	0	0	7.50	5.73	0.744	0.031	61	1	1	14.90	5.27	0.853	0.055
FRE‐3	43.8005	−70.1331	22	0	0	6.88	5.51	0.693	0.040	18	0	0	9.60	4.92	0.831	0.081
GLS‐1	45.17524	−67.82938	25	0	0	6.75	5.50	0.690	0.009	27	0	0	10.80	4.89	0.774	0.014
GLS‐2	45.17594	−67.86995	29	0	0	6.75	5.33	0.688	0.088	–	–	–	–	–	–	–
GLS‐3	45.18357	−67.965	–	–	–	–	–	–	–	21	2	0	10.90	4.90	0.784	0.092
GRE‐1	45.440053	−69.392208	31	0	0	7.38	5.52	0.712	−0.006	–	–	–	–	–	–	–
GRE‐2	45.538783	−69.543654	23	0	0	7.25	5.93	0.722	−0.013	–	–	–	–	–	–	–
GRE‐3	45.695279	−69.46204	30	0	0	7.38	5.81	0.744	0.074	58	0	0	13.80	5.20	0.831	0.037
GRE‐4	45.708404	−69.669413	20	0	0	6.25	5.33	0.708	0.001	39	0	0	12.90	5.14	0.832	−0.003
HAM‐1	44.78044	−68.79021	39	0	0	7.50	5.88	0.747	0.058	50	0	4	12.00	5.02	0.797	−0.009
HAM‐2	44.87251	−68.70609	39	0	2	5.88	4.99	0.680	0.010	17	0	2	9.90	4.92	0.805	0.022
HAM‐3	44.725234	−68.839386	12	0	0	5.88	5.59	0.737	−0.028	19	0	0	11.00	5.18	0.826	0.013
HAM‐4	44.72599	−68.83836	23	1	0	6.50	5.60	0.743	−0.022	16	0	0	10.40	5.11	0.795	0.030
HAM‐5	44.72611	−68.83833	17	0	0	6.13	5.51	0.764	0.108	–	–	–	–	–	–	–
HAM‐6	44.76887	−68.81382	–	–	–	–	–	–	–	16	0	0	10.40	5.20	0.808	0.123
JEF‐1	44.13455	−69.5714	23	0	0	7.38	5.78	0.717	0.044	–	–	–	–	–	–	–
JEF‐2	44.13133	−69.5752	60	0	1	8.63	5.73	0.714	−0.043	31	0	0	12.80	5.28	0.841	0.035
JEF‐3	44.22778	−69.46122	64	0	0	8.88	5.98	0.714	−0.014	37	0	2	12.90	5.24	0.839	−0.003
JEF‐4	44.1574	−69.57211	23	0	0	7.13	5.89	0.740	0.060	30	0	1	12.70	5.30	0.837	0.034
JON‐1	44.70088	−67.5263	23	0	0	7.88	6.37	0.764	0.005	29	0	0	11.30	4.89	0.769	0.028
LIN‐1	44.33146	−68.063	17	0	0	6.25	5.62	0.702	0.146	–	–	–	–	–	–	–
LOV‐1	44.2119	−70.87366	29	0	0	7.63	5.85	0.705	0.018	23	1	0	13.00	5.41	0.865	0.042
LOV‐2	44.18099	−70.9356	15	0	0	6.13	5.64	0.719	−0.122	25	0	0	12.00	5.30	0.861	0.024
MAC‐1	44.72039	−67.5263	–	–	–	–	–	–	–	19	0	0	8.70	4.56	0.748	0.012
MNT‐1	44.44843	−69.30656	31	0	0	7.00	5.48	0.676	0.016	16	0	0	10.50	5.24	0.835	0.010
NEW‐1	43.9913	−69.5885	–	–	–	–	–	–	–	26	0	0	11.30	5.14	0.835	0.015
OGU‐1	43.24324	−70.61281	22	0	0	7.50	5.77	0.717	0.013	28	1	0	13.40	5.34	0.860	0.063
ORO‐1	44.888583	−68.78225	20	0	0	6.50	5.55	0.743	0.083	26	0	0	11.40	5.11	0.816	0.007
ORO‐3	44.878717	−68.757317	19	1	0	6.75	5.92	0.721	−0.024	–	–	–	–	–	–	–
ORO‐4	44.895533	−68.722867	29	0	3	7.50	5.79	0.757	0.064	32	0	1	10.70	4.91	0.801	−0.031
ORO‐5	44.8933	−68.723917	15	0	0	6.25	6.03	0.762	0.034	–	–	–	–	–	–	–
ORO‐6	44.889383	−68.760517	18	0	0	6.63	5.58	0.754	0.056	–	–	–	–	–	–	–
ORO‐7	44.8981	−68.68752	12	0	0	5.75	5.78	0.728	−0.088	12	0	0	8.80	4.99	0.791	−0.032
ORO‐8	44.87286	−68.70531	37	0	0	8.25	5.43	0.723	0.055	23	0	0	11.20	5.09	0.808	0.000
ORO‐10	44.89615	−68.72778	23	0	0	7.13	5.97	0.746	−0.014	–	–	–	–	–	–	–
OT‐2	44.93884	−68.68915	–	–	–	–	–	–	–	26	0	0	11.60	5.05	0.805	0.032
OT‐3	44.93901	−68.6713	19	0	0	6.88	5.82	0.706	−0.025	–	–	–	–	–	–	–
OT‐5	44.93901	−68.6713	–	–	–	–	–	–	–	21	0	0	11.20	5.15	0.807	0.015
PHI‐1	44.82591	−70.40891	20	0	0	6.75	5.86	0.718	0.064	22	0	0	12.10	5.39	0.853	−0.019
PIT‐1	44.79525	−69.3681	19	0	0	6.63	5.72	0.733	0.036	44	3	1	13.60	5.49	0.863	0.091
PRO‐1	44.536833	−68.880817	16	0	0	6.88	6.00	0.748	0.032	–	–	–	–	–	–	–
POW‐1	43.94315	−70.203017	–	–	–	–	–	–	–	28	0	0	12.50	5.32	0.856	−0.014
SBR‐1	43.87555	−69.56431	53	0	0	7.25	5.48	0.683	0.020	33	0	1	10.70	4.85	0.798	−0.008
SCA‐1	43.55003	−70.36169	32	1	2	6.50	4.97	0.702	0.051	30	0	4	10.30	4.95	0.832	0.031
SCA‐2	43.60163	−70.38015	24	0	0	6.13	5.16	0.699	0.026	31	1	0	12.70	5.34	0.851	0.060
SEA‐1	44.49271	−68.93039	–	–	–	–	–	–	–	26	0	2	10.10	4.93	0.805	−0.003
SEB‐1	43.90471	−70.67239	49	0	0	7.63	5.53	0.731	0.013	27	0	0	13.40	5.31	0.849	−0.017
SHA‐1	43.50199	−70.79466	20	0	0	7.00	5.86	0.745	0.139	29	0	0	12.00	5.18	0.834	0.049
SID‐1	44.42603	−69.70481	–	–	–	–	–	–	–	26	1	0	12.20	5.21	0.845	0.064
SKO‐1	44.7757	−69.74332	27	0	0	7.75	5.91	0.749	0.000	28	0	0	12.30	5.22	0.838	0.008
SKO‐2	44.76549	−69.59269	–	–	–	–	–	–	–	20	0	2	10.90	5.18	0.834	0.014
STA‐1	44.77255	−69.91174	16	0	0	6.25	5.56	0.716	0.000	38	0	0	13.60	5.31	0.842	0.050
STA‐2	44.77752	−69.92757	20	0	0	6.38	5.42	0.717	−0.015	19	0	0	10.80	5.30	0.858	0.058
SUL‐1	44.52482	−68.16853	22	0	2	6.88	5.67	0.721	0.001	22	0	1	10.20	4.96	0.792	0.033
TAT‐1	43.28084	−70.69041	51	0	0	7.63	5.71	0.732	0.028	18	0	0	11.00	5.13	0.830	−0.002
TOP‐1	43.954483	−69.976483	23	0	0	6.25	5.31	0.710	−0.027	30	0	0	11.60	5.00	0.821	0.018
TRE‐1	44.80536	−67.15285	–	–	–	–	–	–	–	12	0	0	7.50	4.41	0.736	0.021
WAL‐1	44.180917	−70.024217	15	0	0	7.00	6.14	0.752	0.008	45	0	0	13.80	5.36	0.852	0.043
WAY‐1	44.366033	−70.035417	33	0	0	7.25	5.58	0.702	0.055	58	0	4	14.50	5.21	0.838	0.011
WB‐1	43.74574	−70.36393	–	–	–	–	–	–	–	21	1	0	11.60	5.40	0.863	0.042
WEL‐1	43.34105	−70.55378	18	0	0	6.75	5.75	0.734	−0.013	40	0	0	13.40	5.22	0.845	0.018
WEL‐2	43.340982	−70.55058	31	0	3	7.13	5.63	0.747	−0.015	28	1	1	12.50	5.42	0.868	0.045
WEL‐4	43.31988	−70.59468	32	0	0	6.50	5.20	0.717	0.059	19	0	0	11.40	5.10	0.817	−0.030
WEY‐1	45.39568	−70.0009	–	–	–	–	–	–	–	15	0	0	10.40	5.28	0.842	0.047
WEY‐2	45.044343	−69.990037	11	0	0	5.25	5.11	0.657	−0.004	25	0	4	12.30	5.38	0.853	−0.080
WEY‐3	45.09518	−69.78388	28	0	0	7.50	5.48	0.728	−0.025	26	0	0	12.40	5.23	0.841	0.031
WG‐1	44.22255	−69.893033	11	0	0	6.00	5.84	0.722	−0.036	34	0	0	13.90	5.49	0.868	−0.023
WG‐2	44.221783	−69.901083	12	0	0	5.13	4.89	0.680	−0.026	–	–	–	–	–	–	–
WHI‐1	44.78317	−67.54623	19	0	0	7.38	6.04	0.745	0.061	30	0	1	11.40	4.91	0.768	−0.023
WIL‐1	44.61431	−70.26837	17	0	0	6.00	5.27	0.719	−0.107	23	0	6	10.20	5.21	0.855	0.033
WIS‐1	43.960033	−69.693833	–	–	–	–	–	–	–	29	0	1	11.40	5.22	0.851	0.015
WLM‐1	45.3036	−69.34457	24	0	0	6.25	5.41	0.724	−0.067	–	–	–	–	–	–	–
WM‐1	45.023405	−70.454206	–	–	–	–	–	–	–	13	0	4	9.00	5.29	0.847	0.143
WM‐2	44.92064	−70.50822	19	0	0	7.38	5.99	0.710	−0.097	12	0	0	9.50	5.28	0.851	0.056
WM‐3	44.767724	−70.49833	10	0	0	5.75	5.75	0.732	−0.158	–	–	–	–	–	–	–
WM‐4	44.73882	−70.45364	17	0	0	5.88	5.29	0.716	−0.023	30	0	0	13.60	5.45	0.867	0.012
WOL‐1	43.9494	−69.807	29	0	0	7.88	6.23	0.776	−0.087	28	0	0	12.10	5.17	0.828	0.033
WOO‐1	44.4036	−70.548683	21	0	0	7.00	5.71	0.733	0.022	29	0	0	13.20	5.45	0.869	0.014
YAR‐1	43.79115	−70.205483	14	0	0	6.75	6.00	0.706	−0.040	20	0	0	11.80	5.34	0.856	0.026
YAR‐2	43.821117	−70.209033	47	0	0	6.75	5.19	0.693	0.124	–	–	–	–	–	–	–
YAR‐3	43.98565	−69.944283	–	–	–	–	–	–	–	21	0	1	11.10	5.24	0.850	−0.033
YAR‐4	43.812283	−70.176717	–	–	–	–	–	–	–	16	0	4	9.20	5.02	0.835	−0.068
YOR‐1	43.122728	−70.64477	30	0	0	5.75	4.92	0.681	−0.044	–	–	–	–	–	–	–
YOR‐2	43.175521	−70.71255	32	0	0	7.13	5.29	0.692	0.002	24	0	0	12.30	5.34	0.852	0.013
YOR‐3	43.177481	−70.640711	25	0	0	6.75	5.71	0.759	0.015	–	–	–	–	–	–	–
YOR‐4	43.230644	−70.631586	57	1	0	7.88	5.43	0.711	0.022	–	–	–	–	–	–	–
Average			26.81	0.06	0.18	6.91	5.64	0.72	0.01	28.03	0.20	0.83	11.58	5.15	0.83	0.03

Abbreviations: A O, average number of alleles across loci; AR, allelic richness; F IS, Wright's inbreeding coefficient; H E, expected heterozygosity; HW, number of Hardy–Weinberg equilibrium nonconforming samples; LD, number of samples with loci in linkage disequilibrium; N, sample size.

Two spotted salamander loci had high null allele frequencies (AmaD328: 0.272 and AmaD315: 0.283) and were therefore excluded from further analyses. Null allele frequencies for the remaining eight spotted salamander loci ranged from 0.005 to 0.027. Tests of nonrandom assortment of genotypes indicated 16 of 2,578 tests (0.6%; Table 3) were significant. Significant violations of Hardy–Weinberg proportions were observed in 5 of 728 tests (0.7%; Table 3). Wood frog null allele frequencies ranged from 0.012 to 0.055 among loci. Tests of nonrandom assortment of genotypes indicated 72 of 3,897 tests (1.8%; Table 3) for wood frogs were statistically significant. Significant violations of Hardy–Weinberg proportions were observed in 17 of 870 tests (2.0%; Table 3) for wood frogs. No clear patterns of significance were detected within loci or sampling sites for either nonrandom assortment of genotypes or Hardy–Weinberg testing for either species; therefore, no loci or sites were excluded on the basis of these tests. F IS averaged 0.01 (±0.0052 SE) for spotted salamanders and 0.03 (±0.0041 SE) for wood frogs (Table 3). Missing allele calls occurred for 1.3% of locus–sample combinations for spotted salamander and 2.1% for wood frogs. Genotyping error rates were observed in 0.8% of instances for wood frogs and 0.9% for spotted salamanders.

Genetic diversity, differentiation, and isolation by distance

Measures of genetic diversity, including H E, AR, and F IS, varied between the species but were of similar magnitudes. Across loci and among sites, spotted salamander AR averaged 5.64 (±0.298 SE) and H E averaged 0.72 (±0.024 SE), whereas wood frog AR averaged 5.15 (±0.233 SE) and H E averaged 0.83 (±0.032 SE; Table 3). The greater difference in values between species for AR relative to H E is unsurprising given the relative insensitivity of H E to the number of alleles observed (Maruyama & Fuerst, 1985). Following an FDR correction, no evidence of genetic bottlenecks was detected for either species at any population. EEMS analyses generated 54 spotted salamander and 53 wood frog demes. Strongly contrasting geographic patterns of genetic diversity were observed between the two species. For instance, across the range of study sites, spotted salamanders had several interspersed regions of high and low diversity, whereas wood frogs showed a clear gradient of high diversity to the west transitioning to lower diversity in the east (Figure 2). Because the analysis occasionally groups separate sampling sites into a single deme, and both species have strong spatial structuring, some locally high levels of diversity identified by the analysis may be due to the grouping of dissimilar populations into a single deme.

Figure 2

Spatially heterogeneous effective rates of genetic diversity among 54 spotted salamander and 53 wood frog demes. Black points indicate the location and relative sample size of each deme

Average genetic differentiation among sites was relatively low for both species but varied widely. For spotted salamanders, global G ST was 0.024 and G″ ST was 0.087. Pairwise G ST values ranged from −0.006 to 0.068 and G″ ST ranged from −0.048 to 0.384. Following FDR correction of alpha levels, 3,474 of 4,005 tests (86.7%) were significant. For wood frogs, global G ST was 0.032 and G″ ST was 0.189. Pairwise G ST values ranged from −0.002 to 0.068, and G″ ST ranged from −0.029 to 0.583. Following FDR correction of alpha levels, 3,639 of 3,741 tests (97.3%) were significant. All pairwise G′ST and G″ ST values are provided in Appendix S2.

Isolation‐by‐distance patterns differed between the two species. Despite relatively weak correlation using each genetic distance, IBD relationships were statistically significant for spotted salamanders based on nontransformed geographic distances (G ST: r = .196, p < .001; and G″ ST: r = .203, p < .001), as well as following log transformation (G ST: r = .18, p < .001; and G″ ST r = .181, p < .001). IBD patterns without the geographic distance log transformation were stronger for wood frogs for G ST (r = .628, p < .001) and G″ ST (r = .593, p < .001) and were marginally weakened following transformation for both G ST (r = .461, p < .001) and G″ ST (r = .433, p < .001). Because the IBD relationship for the wood frogs appeared nonlinear, we also fit a quadratic rather than linear model to the data. Due to similar patterns between the genetic distance measures, only plots based on G ST are shown (Figure 3 inset panels). The regression slopes for genetic distance versus log‐transformed geographic distances were significantly different from zero for both species (p < .001), and the slope estimate for spotted salamanders (β = 0.0018) was less than that for wood frogs (0.0064).

Figure 3

Associations between the slope (β) of the regressed isolation‐by‐distance (IBD) relationship and the maximum pairwise distance of the sample set considered for spotted salamanders and wood frogs. Shaded areas indicate the 95% confidence intervals of β coefficients for each iteration of the analysis. Inset figures depict pairwise relationships between geographic (km) and linearized genetic distances (G ST/1−G ST) indicated by ordinary least squares (OLS) regression (solid line) and 95th and 5th quantile regressions (dashed lines) for 90 spotted salamander and 87 wood frog populations. Wood frog OLS regression: y = 1.028 × 10−2 + −1.831 × 10−5 x + 4.581 × 10−7 x 2, R 2 = 0.445, Mantel's r = 0.628, p < .001. Spotted salamander: y = 9.416 × 10−3 + 2.535 × 10−5 x, R 2 = 0.038, Mantel's r = .196, p < .001

Our IBD scaling profiles indicated that β values ranged widely for each species depending on the maximum pairwise distance included in the analysis and the responses of β to maximum pairwise distances were strongly nonlinear for both species (Figure 3). The Mantel correlogram indicated that IBD relationships were strongest at shorter distance classes for both species, with spotted salamander associations becoming nonsignificant at distances greater than 60 km (Figure 4). For wood frogs, the greatest distance class had significant negative spatial autocorrelation, which aligns well with the particularly high levels of population differentiation at large scales that were detected with IBD regressions. EEMS identified several geographic regions with more and less gene flow than expected under an IBD scenario. For instance, the north‐central portion of the study area consistently had relatively low connectivity, whereas multiple coastal regions were more connected. An area of low connectivity was also noted for spotted salamanders in the most densely human‐populated area around Portland, Maine; however, a similar pattern was not observed for wood frogs (Figure 5).

Figure 4

Mantel correlograms indicating associations between genetic (G ST) and geographic distance among spotted salamander and wood frog site pairs. Filled symbols indicate statistical significance based on 10,000 bootstrap replicates and a false discovery rate correction for multiple testing based on an alpha level of 0.05

Figure 5

Spatially heterogeneous effective rates of migration among 54 spotted salamander and 53 wood frog demes. Black points indicate the location and relative sample size of each deme. Migration rates (m) are relative to background rates. For instance, a value of 10 is reflective of 10× greater migration than the background rate

We detected significant relationships among genetic diversity, urbanization, and isolation. Residuals were measured using the relationship among linearized G ST and nontransformed geographic distance because a strong correlation with residuals of the log‐transformed relationship was observed for both species (spotted salamanders: 0.994; wood frogs: 0.934). Sampled sited represented a broad range of the measured environmental attributes (Appendix S3). Nearby canopy cover, distance to roadway, and nearby amount of impervious surface were all highly correlated with the total distance of light roads within 1 km of each study site, allowing us to retain only the three road classes. Multiple regression models that sought to explain variation in allelic richness and expected heterozygosity based on the three road classes were generally nonsignificant. Statistical significance was detected for the wood frog allelic richness model, and the secondary roads term had a significant positive relationship with expected heterozygosity for spotted salamanders, although each of these relationships explained very little variation, suggesting limited biological relevance (Table 4).

Table 4

Results of multiple regression models assessing effects of three road types on allelic richness, expected heterozygosity, and site isolation

Parameter	Allelic richness				Expected heterozygosity				Site isolation
	Spotted salamanders		Wood frogs		Spotted salamanders		Wood frogs		Spotted salamanders		Wood frogs
	β	p	β	p	β	p	β	p	β	p	β	p
Light roads	−1.78 × 10−5	.099	−9.82 × 10−6	.103	−1.20 × 10−6	.173	1.33 × 10−7	.871	5.67 × 10−7	.003	5.99 × 10−7	<.001
Secondary roads	8.95 × 10−6	.768	−1.91 × 10−5	.399	5.22 × 10−6	.037	−5.24 × 10−6	.093	1.44 × 10−7	.781	1.75 × 10−8	.968
Primary roads	−2.97 × 10−5	.202	−1.79 × 10−5	.319	−2.30 × 10−7	.903	−2.86 × 10−6	.245	−5.31 × 10−8	.893	2.75 × 10−7	.431
Adjusted R 2	0.026		0.061		0.024		0.026		0.086		0.297
Full model p	.155		.041		.170		.162		.013		<.001

Variables and models with statistically significant values (p < .05) are indicated in bold font.

The degree of isolation (mean IBD residual) experienced by a site was significantly greater for locations with more nearby light roads for both species, with this effect being stronger for wood frogs (β = 5.99 × 10–7, p < .001) than spotted salamanders (β = 5.67 × 10–7, p = .003; Table 4). For the wood frog dataset involving nearby light road length, a single outlier site was removed due to having over twice the distance of nearby light roads than the next closest site. The influence of light roads on site isolation was also analyzed using a simple linear regression (Figure 6), which further enforced the positive relationship. Finally, expected heterozygosity and allelic richness declined as a site's degree of isolation increased for both species; however, these declines were stronger for spotted salamanders than for wood frogs (Figure 7).

Figure 6

Linear regression analyses indicating a positive relationship between site isolation and the length of light roads within one km of a breeding site for spotted salamanders and wood frogs

Figure 7

Linear regression analyses illustrating negative relationships between sampling site isolation and allelic richness and expected heterozygosity for analyzed spotted salamanders and wood frogs

Distance‐based RDAs for each species identified relationships between interpopulation genetic differentiation and measures of urbanization. Following backward optimization, each model contained the light roads variable while controlling for latitude and longitude. Density of light roads was significantly associated with G ST values for wood frogs (F = 3.95, p = .007), but statistical support for the spotted salamander model was marginal (F = 1.59, p = .097). Variance partitioning based on the RDAs' adjusted R 2 values revealed a much better overall model fit for wood frogs than spotted salamanders, mostly attributable to a stronger IBD signal for the wood frogs. For spotted salamanders, the light roads explained 2.4%, geography explained 5.4%, and the terms collectively explained 12.7% of the variation. For wood frogs, the amount of variation explained was 3.6% for light roads, 47.6% for geography, and 53.7% for the combined terms. MMRR further supported the existence of a significant relationship between interpopulation genetic divergence and both geographic distance and density of nearby light roads for both species (Table 5).

Table 5

Results of multiple matrix regression with randomization assessing effects of three road types on interpopulation genetic differentiation values

Parameter	Spotted salamanders			Wood frogs
	β	t	p	β	t	p
Light roads	6.12 × 10−7	11.93	.006	5.31 × 10−7	18.39	<.001
Secondary roads	4.53 × 10− 8	0.36	.921	2.22 × 10–7	1.69	.549
Primary roads	2.08 × 10− 7	2.66	.616	5.76 × 10–7	7.24	.123
Geographic distance	2.06 × 10−5	10.22	<.001	1.07 × 10–4	52.64	<.001

Variables with statistically significant values (p < .05) are indicated in bold font.

Landscape resistance

The resistance surface models that we constructed provided insight into the relative influence of numerous natural and anthropogenic landscape features on interpopulation connectivity for spotted salamanders and wood frogs. Optimization of resistance surfaces indicated several differences in which landscape features most influence connectivity for each species (Table 6; Appendix S4). For instance, light roads and interstates were suggested as important for spotted salamanders but not for wood frogs. Both river classes had less influence on genetic structure than distance alone for spotted salamanders, whereas rivers were strongly influential for wood frogs. Generally, anthropogenic features such as roads (particularly interstate highways) and developed land cover trended toward negative influences on connectivity in both species.

Table 6

Landscape features included in resistance surface models and optimized resistance values based on partial Mantel testing

Landscape feature	Spotted salamanders	Wood frogs
Land cover A	1	1
Land cover B	1	1
Land cover C	5	15
Interstates	1,000	500
Secondary roads	1	25
Light roads	25	10
Medium rivers	1	500
Large rivers	1	4,500
Railroads	1	1
Terrain ruggedness index	TRI + 500	TRI + 500

Values of 1 represent an isolation‐by‐distance model.

Generalized linear modeling based on effective resistance distances between populations provided strong support for the influence of multiple landscape features on connectivity of each species. By limiting inclusions of landscape features to only those with a stronger influence than distance alone, we assessed four variables for spotted salamanders and seven variables for wood frogs, in addition to an IBD‐only model. This resulted in 17 candidate models for spotted salamanders and 129 for wood frogs. Spotted salamander resistance values were best explained with a single top model that included land cover class C (development and open water), interstates, light roads, and terrain ruggedness (Table 7). For wood frogs, there were four models within two AICC points of each other, indicating they were equally strongly supported. These four models each included land cover class C, medium rivers, large rivers, light roads, and terrain ruggedness with interstates and secondary roads each occurring in two of these top four models (Table 7). The top eight models all included land cover class C, and terrain ruggedness was present in the each of the top 12 models, suggesting a major role for these features in determining genetic structure for wood frogs. The IBD‐only model was one of the least supported for both species.

Table 7

Results of additive landscape resistance models ranked based on the parsimony‐weighted AICC

Landscape resistance model	AICC	ΔAICC	Weight
Spotted salamanders
LcC + Interstates + LtRoads + TRI	−26,456.06	0.00	0.94
Interstates + LtRoads + TRI	−26,450.33	5.73	0.05
Interstates + LtRoads	−26,443.78	12.28	0.00
LcC + Interstates + LtRoads	−26,443.74	12.32	0.00
LcC + LtRoads	−26,383.67	72.39	0.00
LtRoads + TRI	−26,381.98	74.07	0.00
LcC + LtRoads + TRI	−26,381.77	74.29	0.00
LtRoads	−26,378.66	77.40	0.00
LcC + Interstates + TRI	−26,254.38	201.67	0.00
LcC + Interstates	−26,215.35	240.71	0.00
Interstates	−26,206.53	249.53	0.00
Interstates + TRI	−26,204.55	251.51	0.00
Wood frogs
LcC + LtRoads + MedRivers + LgRivers + TRI	−23,885.32	0	0.32
LcC + Interstates + LtRoads + MedRivers + LgRivers + TRI	−23,884.73	0.59	0.24
LcC + SecRoads + LtRoads + MedRivers + LgRivers + TRI	−23,884.19	1.13	0.18
LcC + Interstates + SecRoads + LtRoads + MedRivers + LgRivers + TRI	−23,883.61	1.71	0.14
LcC + MedRivers + LgRivers + TRI	−23,881.90	3.41	0.06
LcC + SecRoads + MedRivers + LgRivers + TRI	−23,880.89	4.43	0.03
LcC + Interstates + MedRivers + LgRivers + TRI	−23,879.95	5.37	0.02
LcC + Interstates + SecRoads + MedRivers + LgRivers + TRI	−23,878.94	6.37	0.01
SecRoads + LtRoads + MedRivers + LgRivers + TRI	−23,784.06	101.25	0.00
Interstates + SecRoads + LtRoads + MedRivers + LgRivers + TRI	−23,782.06	103.26	0.00
LtRoads + MedRivers + LgRivers + TRI	−23,780.59	104.72	0.00
Interstates + LtRoads + MedRivers + LgRivers + TRI	−23,778.59	106.73	0.00

The top 12 ranked models are shown. Models only included variables with resistance values that explained genetic distances better than the null isolation‐by‐distance scenario.

Abbreviations: Interstates, interstate highways; LcC, land cover class C; LgRivers, sixth‐ and seventh‐order rivers; LtRoads, light roads; MedRivers, fourth‐ and fifth‐order rivers; SecRoads, secondary roads; TRI, terrain ruggedness index.

DISCUSSION

Natural and anthropogenic landscape features contribute to interpopulation genetic structuring for both spotted salamanders and wood frogs; however, the strength of effects for specific landscape features differs markedly between the species. Both species experience increased population isolation in urban areas and decreased genetic diversity as population isolation increases; however, no direct connection between roads and genetic diversity is evident. Moreover, specific natural and anthropogenic landscape features generally affect each species' gene flow differently, which is likely a result of how differing life history and behavioral tendencies influence interactions of each species with the landscape. These local effects also occur in the context of contrasting broadscale patterns of the distribution of genetic variation for each species. Collectively, these results suggest that urban landscape elements are reshaping metapopulation‐level dynamics for spotted salamanders and wood frogs, although the effects are not necessarily consistent among these two sympatric species.

Effects of urbanization

Our work has identified elements of urban landscapes that are capable of influencing connectivity among spotted salamander and wood frog populations. Density of light roadways was identified in multiple analyses as an important factor in restricting connectivity among populations. However, it is important to put this result in the context of a high level of correlation among light roads, canopy cover, distance to nearest road, and amount of nearby impervious surface. As such, the effects of light roadways are likely an indicator of urbanization as a whole, rather than light roads exclusively. However, the effects of roadways themselves on gene flow should not be understated, as they have consistently been recognized as hazardous for migratory amphibian species (reviewed in Schmidt & Zumbach, 2008) and previously observed to diminish interpopulation connectivity for spotted salamanders (Coster, Babbitt, Cooper, et al., 2015; Richardson, 2012) and wood frog (Richardson, 2012).

Estimated effective migration surfaces and resistance surface modeling revealed distinctive effects of several landscape features on gene flow for spotted salamanders and wood frogs. For instance, our resistance surface analyses indicated very strong effects of rivers on wood frog connectivity, but no detectable effect for spotted salamanders. Interstate highways, another modeled linear landscape feature, were found to have very strong effects on both spotted salamander and wood frog connectivity, which was unexpected given the interstates in Maine have only been in place for 60 years or less (Ferris, 1979). Relatively rapid responses in genetic structure to the presence of roadways have been observed in other species, but the effect is inconsistent (Holderegger and Di Giulio, 2010). Wood frog gene flow generally conformed to our hypothesis of stronger effects of natural (rivers, terrain ruggedness) versus anthropogenic features (roadways, developed landscapes), whereas salamander gene flow did not, with terrain ruggedness being the only natural landscape feature with an influence appreciably greater than distance. EEMS gene flow models (Figure 5) indicated several clear differences between species. For instance, spotted salamanders had relatively little gene flow where several peninsulas were sampled in the south‐central region of the coast, an area that also coincides with some of the most dense urban development (Figure 2). Wood frogs did not share this pattern, instead having a large area of restricted gene flow in the north‐central region of the study area.

Comparing our study with previous work highlights the context dependency of landscape genetic inferences. For instance, our results often contrast those of Richardson (2012), who used similar landscape genetics approaches with the same two species in the Connecticut River Valley, a region approximately 250 km to the southwest of our study range. That study found a lower IBD slope for wood frogs relative to spotted salamanders; however, the maximum geographic extent of our study was greater (approximately 350 km vs. 225 km), and we found the strongest differentiation for wood frogs occurred at geographic distances greater than those examined by Richardson (2012). Because Richardson (2012) found wood frog slopes were less than slopes for spotted salamanders, he suggests that gene flow likely occurs more frequently for wood frogs than spotted salamanders, whereas our results suggest that such a conclusion likely reflects spatial context. Our resistance surface modeling results also often contrasted with those presented by Richardson (2012), who identified a strong influence of medium and large rivers on spotted salamander structuring and strong effects of railways on both species, whereas we found no detectable effect. Regional differences in the correlation of various landscape features could lead to different variables having stronger or weaker relationships to resistance, and it is likely that real‐world resistance results from compounding factors that are not easily isolated. For instance, Richardson (2012) occurred in a topographically diverse environment, where railways likely follow specific elevation contours, creating collinearity between slope and railway variables that would not be present in most of our study region. Similar disparities among the influence of specific landscape features have been observed across the range of other species. For instance, Cope's giant salamander (Dicamptodon copei) had varying responses to waterways and forest cover in different areas of its range (Trumbo, Spear, Baumsteiger, & Storfer, 2013). Short Bull et al. (2011) found variable effects of the same landscape features across 12 black bear study areas in Montana and Idaho. That study also found that the most variable features within a study area were more likely to receive support in their model. Future studies or efforts to manage landscapes for resistance that focus on features identified in other regions may overlook locally important features and contribute to incomplete or ineffective management. This shortcoming may be addressed through replicated study designs, whereby separate areas of a species' range are the units of replication, or through separate studies as demonstrated here via our comparisons with Richardson (2012). A standardized procedure for quantifying landscape heterogeneity would also improve researchers' ability to make equitable comparisons among species and landscapes.

Inferences from resistance surface modeling are sensitive to the spatial scale (i.e., grain size and study extent) of analyses. Previous work has demonstrated a relationship between the grain size and estimated resistance values of environmental variables (Zeller et al., 2014); however, very high sensitivity was uncommon and even when observed, it remains difficult to determine which grain sizes are most ecologically meaningful. Our grain size was 90 m, which is reasonable given the home range of these species. Moreover, the grain size was determined by the resolution of our most coarse raster dataset and represented a trade‐off with the large study extent to maintain computational tractability. The effects of study extent on resistance estimates are less well resolved (Zeller, McGarigal, & Whiteley, 2012). Given that our study extent was many orders of magnitude greater than the study species' range size, we likely captured some effects unrelated to short‐term dispersal. In that case, the effects of discrete landscape elements (e.g., roadways or rivers) are likely underestimated, as their effect would be diluted across the greatest analyzed geographic distances. The inclusion of multiple species in our analyses should buffer any spurious outcomes associated with our chosen spatial grain and extent (Richardson et al., 2016).

We used the IBD residuals to inform an index for population isolation, which provides a pairwise genetic distance measure standardized using geographic distance and a means of determining the factors that contribute to population‐wise departures from an IBD pattern. Using this metric, we identified a significant relationship between isolation and the urbanization indicator of nearby density of light roads for both species, which was further supported by the elevated resistance values assigned to light roads (Table 6) and the significant influence of light roads identified in our dbRDA. We also used this isolation metric to detect relationships between increasing levels of isolation and declining levels of genetic diversity for both species. A similar relationship was observed by Crawford et al. (2016) who quantified the connectivity among populations of the pool‐breeding Jefferson salamander (Ambystoma jeffersonianum) and found that less connectivity among populations resulted in significantly lower heterozygosities and allelic richness. Similarly, Cosentino, Phillips, Schooley, Lowe, and Douglas (2012) found populations of tiger salamander (A. tigrinum) that were smaller and more isolated had less genetic diversity than populations that were more connected to one another. This lack of a direct effect between genetic diversity and urbanization is not theoretically surprising considering some gene flow likely still occurs, and a lag is likely present between the time for urbanization to generate isolation and the subsequent effects of that isolation on loss of diversity through drift. Along those lines, it is reasonable to hypothesize that much of the reduced genetic variation with isolation in our dataset comes from natural isolation patterns on the landscape. However, genetic diversity losses associated with urbanization have been observed in other species (e.g., Munshi‐South, Zolnik, & Harris, 2016), although the effect is not ubiquitous. For instance, dwarf salamanders (Eurycea quadridigitata) experienced reduced allelic richness as nearby road density increased, although the effect was absent for southern leopard frogs (Lithobates sphenocephalus) in the same landscape (McKee, Maerz, Smith, & Glenn, 2017).

Broadscale genetic structure

The spotted salamander IBD relationship has high variance and a low y‐intercept (Figure 3 inset), making it unusual among commonly observed patterns (Phillipsen et al., 2015; Hutchison & Templeton, 1999). Similar patterns have been observed for spotted salamanders throughout their range, including in central Massachusetts (Whiteley et al., 2014), central Missouri (Burkhart et al., 2017; Peterman et al., 2015), and northeastern Ohio (Purrenhage et al., 2009). Other studies that did not report an IBD intercept did find high variance relationships (Coster, Babbitt, & Kovach, 2015; Zamudio & Wieczorek, 2007). One potential explanation for high variance IBD patterns is an increase in the influence of genetic drift due to consistently depressed effective population sizes caused by limited recolonization capacity associated with relatively small dispersal distances for spotted salamanders. However, the small intercept value suggests an appreciable role of gene flow, and our effective population size estimates had confidence intervals consistently including infinity, likely due to an insufficient number of individuals or loci being sampled per population (analyses not shown). An unidentified factor in the species' biology or ecology such as exceptionally high microsatellite mutation rates or undocumented dispersal processes may also be contributing to the observed pattern.

The strongly nonlinear IBD pattern observed for wood frogs was unexpected and is not typically observed. Generally, nonlinear IBD relationships have been suggested to indicate departures from dispersal–drift equilibrium, secondary contact, or a colonization event (Bradbury & Bentzen, 2007; Hutchison & Templeton, 1999). However, most observed and simulated nonlinear IBD relationships follow a pattern of decreasing slope as geographic scale increases, rather than the increasing slope that we observed (Bradbury & Bentzen, 2007). The observed IBD pattern may result from a combination of contemporary and past processes. First, the absolute range boundaries along one edge of our sampling area created by the Atlantic Ocean may be inflating the degree of genetic differentiation experienced for populations near the boundary. Those populations can be expected to have fewer potential sources for inbound dispersers, possibly limiting the overall contribution of gene flow to homogenizing allele frequencies among populations and amplifying the effects of genetic drift or selection (Eckert, Samis, & Lougheed, 2008). In our case, such boundary effects may be particularly strong with several of the most distant wood frog sampling sites located along the coastline (Figure 1), meaning both populations in the most spatially distant paired comparison may be experiencing this effect. However, if genetic drift alone were driving differentiation, we should see an increase in G ST variance for these sites as some drift to allele frequencies more similar to geographically distant sites (i.e., Hutchison and Templeton's (1999) case IV pattern). The observed pattern may also be influenced by postglaciation or postdeforestation recolonization patterns. For instance, recolonization of the region from multiple refugia may cause our most distant contemporary sampling sites to exhibit relatively high levels of differentiation, while secondary contact has eroded these differences for more centrally located sites (Durand, Jay, Gaggiotti, & François, 2009). Another consideration is that an increasing IBD slope may be more common than currently recognized and was only revealed here by the large number of sampled populations and the broad sampled extent relative to the species' dispersal distance (Table 1; Jenkins et al., 2010; Koen, Bowman, Garroway, & Wilson, 2013). Without concentrating sampling on the periphery of a study region, IBD data are inherently scant near study margins due to fewer possible pairwise population combinations at the largest distance classes, which causes anomalously high variance at the upper end of the distance distribution and potentially obscures ecologically relevant patterns.

We identified strong relationships between the spatial extent of analyses and the strength of IBD relationships using our IBD scaling profiles. By quantifying the IBD slope across the range of distance values in our dataset, we could assess the relative importance of gene flow to genetic structuring across scales for each species. Slopes were greatest for wood frogs up to about 6 km and to about 9 km for spotted salamanders, suggesting that opportunities for instances of substantive pairwise isolation and divergence quickly increase with distance as the strongest locally constraining effects of gene flow become less universal to all population pairs. Beyond these distances of inclusion, IBD slopes decrease as the incidence of isolated pairs that are substantively divergent becomes increasingly marginal. An analogous pattern of a strong influence of gene flow at relatively short distances was identified by van Strien, Holderegger, and Heck (2015), who performed a similar analysis to identify the distance of maximum correlation using simulated data. This result also emphasizes the importance of considering spatial scaling of inferences in landscape genetic studies, an issue that has recently been emphasized by other authors (Balkenhol et al., 2016). Cushman and Landguth (2010) conducted a series of Mantel tests to examine relationships between genetic and geographic distances using simulated data while varying the window size (i.e., extent) of their analysis. In their study, Mantel r values declined as window size increased, although the transition was relatively gradual. In our study, the slope of the regressed IBD relationship experienced significant nonlinear dynamics depending on the spatial extent that was evaluated. If such scaling patterns are prevalent for IBD, which seems likely, this would suggest that there is little validity in directly comparing overall IBD slopes across studies conducted at very different geographic scales. Although more in‐depth analyses are possible (e.g., Galpern & Manseau, 2013), a workable alternative is to use the provided R code to generate IBD scaling profiles (see Appendices [Link], [Link], [Link], [Link]) for comparison of different studies or species at overlapping IBD inference scales.

CONCLUSIONS

Our study identified critical differences and similarities in how two sympatric species with similar habitat requirements are affected by landscape context. At the scale of single populations, both species responded negatively to the effects of nearby urbanization, whereas interpopulation dynamics differed between the species depending on landscape features. These results can inform conservation of pool‐breeding amphibians, as well as metapopulation‐structured species more broadly. Species with small home ranges capable of satisfying all their life history requirements are sometimes protected using a core–habitat conservation approach (Baldwin, Calhoun, & de Maynadier, 2006a; Semlitsch & Jensen, 2001) that aims to protect species through preserving the structures and functions of requisite habitats. This approach is often applied to pool‐breeding amphibians, where breeding habitats and adjacent upland environments are targeted for protections (Baldwin et al., 2006b). While this conservation strategy protects the majority of individuals that are faithfully philopatric to natal sites (Vasconcelos & Calhoun, 2004), the dispersers that demographically and evolutionarily connect subpopulations are left vulnerable if the intervening landscape is unprotected. The importance of landscape protections aimed at preserving connectivity among subpopulations has been recognized and implemented for some large‐bodied species (e.g., establishment of wildlife corridors; Sharma et al., 2013), and a landscape genetic approach as applied here is well positioned to provide insight into dispersal routes for more cryptic species. Preservation or restoration of those landscape types that are highly permeable to gene flow could be particular effective for increasing metapopulation‐level stability in amphibian species, as a loss of connectivity among populations has been identified as a leading cause of biodiversity loss (Pittman, Osbourn, & Semlitsch, 2014), which is occurring worldwide (Dudaniec, Spear, Richardson, & Storfer, 2012; McCallum, 2007). Assessing the habitat and corridor requirements of sympatric species could also improve the efficacy of management actions by identifying features that similarly affect multiple species. In our case, similarly strong effects of interstate highways indicate that mitigation efforts targeted toward large roadways (e.g., wildlife underpasses or overpasses; Hamer, Langton, & Lesbarreres, 2015) may provide the strongest return on investment for management actions. However, given contrasting results between our study and those of Richardson (2012), the generalizability of these inferences in the context of regionally dependent correlates that may drive observed relationships is warranted.

This research also underscores the importance of scale dependency when considering spatially explicit interpopulation relationships, as highlighted by our IBD scaling profiles that revealed strong spatial variation in isolation‐by‐distance relationships. Research geared toward quantifying variation in gene flow across other studies (e.g., using IBD β values) would benefit from considering the extent over which the data were collected to generate equitable comparisons. Quantifying such variation can help provide an understanding of the heterogeneity in a species' dispersal propensity, consequently producing a more realistic understanding of how species interact with their surrounding environments. Overall, this research provides a strong example of the capacity of urbanization to shape species' interpopulation dynamics. A deeper understanding of the causes and consequences of these effects will provide a robust foundation for identifying and mitigating current and future risks to biodiversity.

CONFLICT OF INTERESTS

None declared.

AUTHOR CONTRIBUTIONS

JJH, CSL, and MTK designed the study. JJH collected samples, performed laboratory work, and conducted data analyses. JJH led the writing of the manuscript. All authors contributed critically to the drafts and approved of the final manuscript.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.