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

Abstract 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.


| INTRODUC TI ON
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(UNDESA, , 2015 has been accompanied by increases in the percentage of people living in urban areas from 29. 4% in 19504% in to 52.1% in 20114% in and up to 70% anticipated by 20504% in (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.
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 poolbreeding 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.
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 m 2 , Ousterhout and Burkhart, 2017;wood frogs: up to 32,165 m 2 , 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.
1. 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 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., TA B L E 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  (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.

| 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-km 2 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) (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. When fewer than 25 individuals were collected at a site in one year, we returned to sample in the subsequent year.

| 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 (A O ) 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 × 10 7 iterations, a burn-in period of 1 × 10 6 iterations, and a thinning interval of 1 × 10 3 (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 logtransformed 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 logtransformed 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;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 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 abovedescribed 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 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.
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 (AIC C ; Burnham & Anderson, 2002). AIC C compares relative support of candidate models including a penalty for the number of variables incorporated, thereby encouraging parsimony. Models with ΔAIC C value <2 were considered equally supported. We used all possible combinations of included variables as candidate models and calculated AIC C values and their relative weights using R package glmulti (Calcagno & Mazancourt, 2010).

| 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).
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 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).
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 humanpopulated area around Portland, Maine; however, a similar pattern was not observed for wood frogs ( Figure 5).

| Regression and multivariate analyses
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 F I G U R E 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  (Table 4).
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).

Distance class (km) Mantel correlation coefficient
Wood frogs Spotted salamanders -0.05 0.00 0.05 0.10 0.15 0.20 F I G U R E 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 -10 0.0 10 E ective migration rates among demes (m) 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).

| 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 TA B L E 4 Results of multiple regression models assessing effects of three road types on allelic richness, expected heterozygosity, and site isolation  Total length of light roads within 1km (m) Population isolation (Average IBD residual) y = 5.78*10 -7 x -2.76*10 -3 Adjusted r 2 = 0.106 p = .001 y = 6.30*10 -7 x -2.75*10 -3 Adjusted r 2 = 0.310 p < .001 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.
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 AIC C 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.

| D ISCUSS I ON
Natural and anthropogenic landscape features contribute to interpopulation genetic structuring for both spotted salamanders and F I G U R E 7 Linear regression analyses illustrating negative relationships between sampling site isolation and allelic richness and expected heterozygosity for analyzed spotted salamanders and wood frogs

| 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). 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.  (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 TA B L E 7 Results of additive landscape resistance models ranked based on the parsimony-weighted AIC C 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

| 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 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 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). 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 ( 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 S1-S4) for comparison of different studies or species at overlapping IBD inference scales.

| CON CLUS IONS
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. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

CO N FLI C T O F I NTE R E S T S
None declared.

AUTH O R CO NTR I B UTI O N S
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.

DATA AVA I L A B I L I T Y S TAT E M E N T
All R code and data associated with this manuscript are available individually and as an R package at https ://github.com/jared homol a/