Comparative landscape genetics of two endemic torrent salamander species, Rhyacotriton kezeri and R. variegatus: implications for forest management and species conservation

Comparative landscape genetic studies provide insights into whether relationships between landscape features and patterns of spatial genetic structure differ among populations, species, habitat types, and regions. For species with fragmented distributions, especially when management practices contribute to fragmentation, tests of the factors structuring population connectivity are particularly important for understanding continued risks. We determined levels of genetic diversity and tested the relationships of landscape-scale vegetative, geographic, and climate variables with genetic distance in two congeneric, endemic salamander species with status of concern. Using microsatellite data for 326 Rhyacotriton kezeri and 557 Rhyacotriton variegatus individuals collected from 17 to 29 localities, respectively, we implemented a model of landscape resistance based on circuit theory. The northernmost portions of each species’ range is more fragmented than areas to the south, leading to the prediction that these areas would have relatively lower genetic diversity in response. Due to reliance of both species upon cold-water habitats, we predicted that landscape variables maintaining cool, moist microhabitats would be correlated with gene flow. Genetic structure was high overall and trended toward increasing with the proportion of the forested landscape. Based on maximum likelihood population effects models across genetic clusters and species, land cover and roads were the best predictors of genetic distance, even though the degree of fragmentation differed across each species’ geographic range. Our results suggest that forest cover is essential for dispersal in these salamanders, indicating negative effects of fragmentation resulting from timber harvest and other forest disturbances.


Introduction
Landscape genetic studies integrate population genetics, landscape ecology, and spatial statistics to develop and test predictions about the influence of landscape features on spatial population genetic structure and gene flow (Manel et al. 2003;Storfer et al. 2007). As an alternative to a null model of isolation by distance (IBD; Wright 1943), these analyses may identify certain landscape-scale variables as suitable macro-habitat conditions and also indicate features that may act as barriers to or facilitators of gene flow. Landscape genetics studies can, in turn, inform more targeted management approaches to promote species conservation than those that simply use existing information about habitat preferences (Segelbacher et al. 2010). For example, recent studies on a range of taxa, including the hazel dormouse (Muscardinus avellanarius; Bani et al. 2017), the alpine newt (Ichthyosaura alpestris ;Luqman 2018), and wild tigers (Panthera Electronic supplementary material The online version of this article (https ://doi.org/10.1007/s1059 2-019-01172 -6) contains supplementary material, which is available to authorized users. 1 3 tigris; Thatte et al. 2018) have leveraged landscape genetic results to prioritize conservation goals for these species. Thus, there is a growing precedent for the incorporation of landscape genetics into conservation planning.
The majority of landscape genetic studies have focused on a single species in one portion of its range (Segelbacher et al. 2010;Storfer et al. 2010;Keller et al. 2014;Richardson et al. 2016). However, comparative studies can provide insights into whether relationships between landscape features and genetic distance differ among taxa, habitat types, and geographic areas. Comparative studies examining either ecologically similar or distinct species in the same landscape have become more common in recent years, providing novel perspectives on different species' relationships with the landscape (Engler et al. 2014;Whiteley et al. 2014;Mims et al. 2015;Kierepka et al. 2016;Burkhart et al. 2017;Gutiérrez-Rodríguez et al. 2017). In one such study, Engler et al. (2014) determined that among three Hesperid butterfly species in the genus Thymelicus, recent human land conversion restricted gene flow in a poorly dispersing habitat specialist, while a generalist with moderate dispersal abilities was restricted by climate, and another generalist that was the strongest disperser exhibited IBD. Similarly, in three amphibian species (Hyla arenicolor, Anaxyrus punctatus, and Spea multiplicata), Mims et al. (2015) determined that more philopatric species with longer larval periods exhibit greater genetic differentiation. Whiteley et al. (2014) found differences in the level of genetic differentiation of two ecologically similar Ambystomatid salamanders (Ambystoma maculatum and A. opacum), pointing to differential responses to habitat fragmentation between species. Furthermore, Burkhart et al. (2017) identified three distinct landscape resistance layers as the top models explaining genetic distance in A. maculatum, and A. opacum, A. annulatum), yet all unexpectedly indicated movement through dry areas. Thus, it is not always predictable how landscapescale connectivity will relate to site-specific habitat studies, and closely-related organisms may respond differently to the same landscape. Hypothesis testing based on species' biology may lead to the emergence of generalities as these comparative landscape genetic studies become more common.
In contrast to assessing multiple species in the same study area, other studies have focused on a single species across different portions of its range (Short Bull et al. 2011;Trumbo et al. 2013;Castillo et al. 2016). This approach can be particularly revealing for species with complex life histories, such as those that have distinct habitats for breeding, rearing, foraging, and dispersal, especially as conditions for each may vary with landscape-scale factors (Steele et al. 2009). Stream-breeding amphibians, which emerge and disperse across upland terrestrial habitats, may exemplify this type of ecological complexity. In one case, Trumbo et al. (2013) contrasted landscape use by the Cope's giant salamander, D. copei, in three distinct regions across its range, discovering that gene flow and the relative importance of particular landscape features varied greatly with regional habitat conditions. Whereas streams were essential for gene flow in one study region, in another area forest cover and solar radiation took precedence, demonstrating that application of ecological insights or conservation strategies from a single area to the species' distribution would be inappropriate (Trumbo et al. 2013). This example demonstrates that comparative studies can provide key insights about ecological heterogeneity that would not be possible with study designs focusing on a single species or a limited study area.
Herein, we apply a comparative landscape genetic approach to two closely related, geographically proximate, and ecologically similar torrent salamander species (Rhyacotritonidae) across multiple portions of their heterogeneous geographic ranges. Torrent salamanders are endemic to the U.S. Pacific Northwest (Oregon, Washington, and northern California), occurring in small, cool, shaded stream habitats in moist coniferous forests (Sheridan and Olson 2003;Olson and Weaver 2007). They are highly sensitive to desiccation, warm temperatures, and ground disturbances that result in sedimentation of their instream breeding habitats, stressors which have been coincident with past timber harvest practices in the region (Adams and Bury 2002). Adverse effects of timber harvest on torrent salamanders have been reported in multiple studies (Diller and Wallace 1996;Welsh and Lind 1996;Sheridan and Olson 2003;Olson and Weaver 2007;Olson and Burton 2014).
We focus our comparison on the Columbia torrent salamander, Rhyacotriton kezeri, which is endemic to coastal southwest Washington and northwest Oregon and has a restricted geographic range compared to other amphibians (Whitton et al. 2012), and its congener the southern torrent salamander, R. variegatus, which has a somewhat larger yet still restricted range extending south to northern California from their shared range boundary in coastal Oregon (Fig. 1). Both species are of conservation concern (ODFW 2008; CDFW BDB 2011), with R. kezeri proposed for listing under the U.S. Endangered Species Act (WSDOT 2017). Forest disturbances are key concerns for both species, as their ranges are highly managed for wood production in the Pacific Northwest, where timber harvest and associated roads may degrade and fragment stream breeding and upland dispersal habitats (Bury and Corn 1988;Corn and Bury 1989). Other main land uses in these species' ranges include paved roads and rural development. Also, fire has periodically affected forest habitat in these species' ranges. For example, with the 1987 Silver Fire Complex (566 km 2 ) and 2002 Biscuit Fire (2000 km 2 ) in the R. variegatus range in Oregon, and four fires from 1933 to 1951 collectively known as the Tillamook Burn (1400 km 2 ) in the R. kezeri range in Oregon. Approximately 28% of the historical area of old-growth conifer forest remains in this region of the U.S. as of 2000, which is estimated to be 46,656 km 2 (Strittholt et al. 2006). More recently, climate change projections have raised concerns for R. variegatus, as its southern and interior distribution may be limited by historically warm, dry climates that may become warmer and dryer in the future (Bury 2015). A comparative landscape genetic study of R. kezeri and R. variegatus may help inform management efforts for these species of concern by revealing landscapescale associations with genetic structure, which can be used to identify forest landscape planning priorities in the area.
In this study, we expand sampling relative to previous work (Emel and Storfer 2015), allowing us to take a comparative approach across portions of each species' geographic range and across species, as well as examine the roles of additional landscape variables. This approach facilitates a  comparison of the potential effects of relative differences in forest fragmentation on the two species, which is a conservation concern. Specifically, the goals of this study are to: (1) compare levels of genetic diversity in relation to forest cover between the species and among different portions of each species' distribution; and (2) determine environmental predictors of genetic distance within each species by examining a mix of land cover type, geographic, and climate factors. We hypothesize that genetic diversity would be lower in R. kezeri due to smaller geographic range and restricted gene flow among populations, and that gene flow relies on dispersal through landscape variables minimizing desiccation due to both species' associations with cool aquatic habitats.

Population sampling
Sampling localities (putative populations) were selected based on the habitat specificity of Rhyacotriton and the availability of previous locality information obtained from the California Natural Diversity Database, the Oregon Department of Fish and Wildlife, the Washington Department of Fish and Wildlife, and Miller et al. (2006). For R. kezeri, sampling occurred in two regions: north and south of the Columbia River in Washington and Oregon. For R. variegatus, we also sampled in two regions: the northernmost portion of the species' range, and an area in the southern portion of its range in northern California, above the phylogeographic barrier of the Smith River identified by Miller et al. (2006; Fig. 1). We opportunistically sampled 17-20 individuals from each of 17 localities for R. kezeri and 29 localities for R. variegatus between 2010 and 2013 (general locality information available in Online Resource 1; due to conservation concern for these species, detailed information may be made available upon request). Adult and larval salamanders were captured by hand or net. Tissues were sampled in the form of tail clips (< 1 cm) taken with a sharp razor blade to minimize suffering and preserved in 95% ethanol. All salamanders were released after tissue collection. This study was carried out under Animal Subjects Protocol #04146-001 of the Washington State University Institutional Animal Care and Use Committee (IACUC) and under the guidance of the US Forest Service prior to its institution of IACUC for amphibians.

Population genetic analysis
DNA was extracted from all samples using Qiagen DNeasy© Blood & Tissue Kits (Qiagen, Inc., Germantown, MD) according to the manufacturer's protocol. Ten microsatellite loci were amplified using fluorescently-labeled primers in 326 individuals of R. kezeri and 557 individuals of R. variegatus, following a standard protocol (Emel and Storfer 2014). Fragment analysis was performed on an ABI 3730XL automated sequencer (Life Technologies, Inc., Carlsbad, CA, USA) at Washington State University. Genotyping was carried out using GeneMapper® Software v4.1 and Peak Scan-ner™ Software v1.0 (Life Technologies, Inc.). We used the software COLONY v2.0 (Wang 2004;Wang and Santure 2009) to identify full siblings using a maximum likelihood approach. We then removed individuals such that full sibling families had no greater than two members with a sibship probability greater than 0.9 represented for each sampling locality, as including siblings may artificially inflate inbreeding and pairwise genetic distance estimates (Goldberg and Waits 2010). Consequently, one R. kezeri and four R. variegatus individuals were removed from subsequent analyses, yielding final sample sizes of 325 and 558, respectively. We assessed locus quality using MICROCHECKER v2.2.3 (van Oosterhout et al. 2004) and FREENA (Chapuis and Estoup 2007;Chapuis et al. 2008). We tested for departures from Hardy-Weinberg equilibrium among loci and sites and linkage disequilibrium between pairs of loci, both using a Bonferroni correction for multiple comparisons in GENEPOP v4.2 (Raymond and Rousset 1995;Rousset 2008).
We performed a spatial genetic mixture analysis of localities (subsequently referred to as populations, due to their concordance with genetic structure) in BAPS v6.0, which incorporates geographic coordinates as priors in a Bayesian model, improving its ability to assign populations to clusters when the genetic data are not conclusive (Corander et al. 2007). We selected this model specifically due to its ability to incorporate spatial information, which is not possible in programs such as STRU CTU RE. We varied K (number of clusters) between 2 and 17 for R. kezeri, and between 2 and 29 for R. variegatus. Evanno et al. (2005) suggested that the K with the lowest maximum likelihood (logML) found by Bayesian clustering algorithms was rarely the correct number of clusters and proposed an alternative method for selecting K; unfortunately, this method is not possible with the output of BAPS, thus we identified the best supported K by the greatest change in log maximum likelihood (ΔlogML). Based on preliminary results, we performed an individual-based admixture BAPS analysis for R. kezeri populations 10 and 11. Because results of Bayesian clustering approaches must be considered alongside other genetic and ecological information for a species (Janes et al. 2017), we thus took prior phylogeographic findings (Miller et al. 2006) and extreme barriers (the Columbia river) into account when selecting the genetic clusters for downstream analysis.
We calculated summary genetic diversity statistics per population using GenAlEx v6 (Peakall and Smouse 2006), including sample size (N), average number of alleles per locus (A R ), average number of expected alleles per locus (A e ), observed heterozygosity (H o ), Nei's unbiased estimate of expected heterozygosity (H e ); and the inbreeding coefficient, (F IS ). To determine differences between the species for these summary measures, we performed one-tailed t-tests to include the direction of the difference in R v3.4.4 (R Core Team 2018).
We used two measures of pairwise genetic distance among all populations: F′ ST , which is F ST divided by its maximum possible value, calculated with GenoDive v2.0b27 (Meirmans and van Tienderen 2004) and D PS , which is one minus the proportion of shared alleles, calculated in Microsatellite Analyzer v4.05; (Dieringer and Schlötterer 2003). F ST is frequently used in population and landscape genetic studies (Emel and Storfer 2012), and the standardized measure F′ ST allows for direct comparison between species by correcting for within-population variation (Meirmans and Hedrick 2011). D PS is less frequently used, but it is not subject to the equilibrium assumptions inherent in using F ST and its analogs to estimate gene flow, and it is thus potentially more appropriate for this study due to likely non-equilibrium conditions caused by recent human disturbance throughout the study area.

Forest fragmentation analysis
To compare the habitat fragmentation with genetic diversity for each species, we used land cover data based on a 16-category scheme from the National Land Cover Database (Homer et al. 2007). We trimmed this layer to rectangles bounded by the maximum spatial extent of the localities for each genetic cluster identified per species, plus a 10-km buffer, and combined the three forest categories into one. Next, we calculated four class metrics based on all forest patches (contiguous forest areas) within each rectangle (landscape) using FRAGSTATS v4.2 (McGarigal et al. 2012). The first metric is the proportion of the landscape covered by forest, ranging from no forest (0) to entirely forested (1). Next, we examined patch density, the number of patches per 100 hectares, which ranges from no forest (0) to a value constrained by the cell size of the landscape data and is informative as a comparative measure because resolution is held constant across landscapes. Additionally, we used two metrics based on the cumulative distribution of patch sizes. The landscape division index, defined as the probability that two randomly selected pixels of the same patch type are in different patches, ranges from entirely forested (0) to almost entirely non-forest with only a single small forest patch as it approaches 1. Finally, we calculated the splitting index, which is the number of patches of equal size that could yield the same landscape division index, ranging from entirely forest (1) to the number of cells in the landscape squared. Thus, higher values reflect increased forest fragmentation for all of these metrics. The results were compared qualitatively with the genetic diversity measures for each genetic cluster.

Landscape genetic analysis
As in forest fragmentation analysis, we created rectangular resistance surfaces for each genetic cluster bounded by the maximum spatial extent of the localities plus a 10-km buffer for each of ten different landscape and climate variables that were important in previous amphibian landscape genetic studies (see Storfer 2008, 2010;Trumbo et al. 2013;Zancolli et al. 2014;Emel and Storfer 2015). Three of ten landscape variables analyzed were the discrete variables of land cover (described above; Homer et al. 2007), road cover (North American Detailed Streets, ESRI, Redlands, CA), and stream cover (Streamnet GIS Data 2003), which we downloaded as vector data and converted into raster data. The seven remaining variables examined were continuous. These variables included percent canopy cover based on tree canopy density (Homer et al. 2007), annual frost-free period (Rehfeldt 2006), annual growing season precipitation (Rehfeldt 2006), compound topographic index (Moore et al. 1993;Gessler et al. 1995), heat load index (McCune and Keon 2002), roughness (Blaszczynski 1997;Riley et al. 1999), and slope. Canopy cover, frostfree period, and growing season precipitation were downloaded in the form of raster files. Compound topographic index, heat load index, roughness, and slope were calculated from the National Elevation Dataset digital elevation model (Gesch et al. 2002;Gesch 2007) using the Geomorphometry and Gradient Metrics toolbox v2.0-0 for ArcGIS (Evans et al. 2014). Compound topographic index is an estimate of soil wetness based on the expected movement of water downhill, heat load index reflects the amount of solar radiation reaching a surface, and roughness quantifies the change in elevation per unit area (Evans et al. 2014). To test for nonlinear relationships between these continuous variables and genetic distance, after rescaling all variables to range from 0 to 1 we created two transformed resistance surfaces for each variable: T1 = 100 (originalcostsurface) , and T2 = 100 − 100 (1−originalcostsurface) . T1 represents a scenario in which the slope of the resistance function is shallow at low resistance and steep at high resistance, whereas T2 has the opposite effect on the resistance function (Trumbo et al. 2013).
All resistance layers created for the landscape variables were at 30-m resolution with the exception of frost-free period and growing season precipitation, which were at 0.008333 decimal degrees, or approximately 825-m resolution. We created resistance surfaces for each variable using ArcGIS 10.0 (ESRI, Redlands, CA). For canopy cover, frostfree period, growing season precipitation, and topographical roughness, each raster cell reflected the value measured for that variable. We predicted that high canopy cover, high compound topographic index, long frost-free period, and high growing season precipitation would promote gene flow, hence we subtracted the value of each cell from the maximum value across all regions for each variable, so that low values would result in high cost. For heat load index and slope, we used the raw value as the cost, because we predicted high values of these variables to be associated with low gene flow. Because the effect of land cover type on these salamanders beyond simply the degree of canopy cover is unknown, we examined categorical land cover in addition to percent canopy cover described above by assigning a resistance of 10 or 100 for developed or barren land, 5 or 50 for non-forested vegetated land (e.g. grassland/herbaceous), and 1 for forested habitat, with the expectation that genetic distance will increase for paths crossing non-forested or developed land compared to forest. For stream cover, we first converted the vector data, which represents streams as linear elements, to raster data at 30-m, where each cell denotes the presence or absence of a stream. We then developed three alternative cost ratios for the resistance associated with dispersal over land versus dispersal via stream corridors because it was unknown a priori how strongly land restricts movement in Rhyacotriton over stream, but likely that it would be important for shaping gene flow. In all cases, stream pixels were assigned a resistance value of 1, while land pixels were assigned a value of 2, 10, or 100, because movement on land would presumably be at least twice as costly as movement within streams, but potentially quite a bit more costly (see Storfer 2008, 2010). We followed a similar procedure for road cover, except that we assigned non-road pixels a resistance of 1, and road pixels a resistance of 2, 10, or 100 because presumably roads are barriers to gene flow for these species. Thus, the best-supported cost ratios for these linear variables were determined via model selection (see below). Finally, for all variables, pixels containing high order rivers (Gary et al. 2010) or ocean were assigned a high cost of 20,000 to prevent the creation of unrealistic paths that would assume movement through these water bodies.
For each landscape variable, we created models of movement through the landscape based on circuit theory in Circuitscape v4.0.5 (McRae et al. 2008), in which paths between populations are modeled similarly to parallel resistors in a circuit. This method allows for multiple paths to exist between each pair of populations, in stark contrast to the single idealized path identified by least-cost paths, another frequently used approach, or a limited number of paths as in least-cost corridors. Additionally, Circuitscape resistance distances are less sensitive than least-cost paths, to the number of pixels used to represent the landscape and to the Euclidean distance between localities (Marrotte and Bowman 2017). Within each genetic cluster, we determined the pairwise average Circuitscape resistance using each of the landscape resistance surfaces listed above, based on a network connecting cells to their eight immediate neighbors. We also created a null model of IBD using a surface with all cells set to a value of 1.
To assess the relationship between the landscape variables and genetic distance, we applied linear mixed-effects models with maximum likelihood-population effects (MLPE; Clarke et al. 2002). All statistical analyses were performed in R (R Core Team 2018). MLPE accounts for the dependency of pairwise genetic distances with a specific covariance structure. Following the procedure of van Strien et al. (2012), we centered each variable around zero and fitted MLPE models with restricted maximum likelihood (REML) using the R package lme4 (Bates et al. 2015). We created separate MLPE models using both D PS and F′ ST as the response variables and each individual Circuitscape resistance model as a fixed effect. Population was included as a random effect, because differences between the pairs of populations may explain some of the variance in the models. The top models for each genetic distance measure and within each cluster were selected using the lowest Akaike Information Criterion corrected for finite sample size (AICc), implemented in the R package MuMIn (Bartoń 2018) setting sample size equal to the number of localities as opposed to the number of pairwise comparisons, and the Bayesian Information Criterion (BIC), implemented in the R base stats package, both of which are appropriate for REML models (Gurka 2006). We considered all models within ΔAICc < 2 of the model with the lowest AICc to be the top models, because these all have similar support. Although Franckowiak et al. (2017) recently called into question the appropriateness of model selection via AIC and related information criterion approaches due to selection of overly complex models, this concern is not applicable to our univariate Circuitscape analysis. Additionally, we used the method developed by Nakagawa and Shielzeth (2012) implemented in the R package GLMM (Johnson 2014) to approximate R 2 as a measure of goodness-of-fit for fixed effects alone (marginal R 2 ) and both fixed and random effects combined (conditional R 2 ) from mixed-effects models.

Population genetic analysis
Analyses of locus quality by population in MICRO-CHECKER indicated no evidence of large allele dropout in either species, and potential scoring error due to stuttering in 1-3 loci for the majority of R. variegatus populations. It also indicated the presence of null alleles in nearly every population analyzed. However, there were insufficient samples for MICROCHECKER to run in all populations. Thus, we also assessed the presence of null alleles using FreeNA and 1 3 obtained pairwise F ST values using the ENA correction. We then determined that the ENA-corrected F ST is very highly correlated with uncorrected F ST (r > 0.99, p < 0.0001), and thus would have minimal impact on subsequent results.
For R. kezeri, there were 18 departures from Hardy-Weinberg equilibrium after Bonferroni correction, with no loci consistently in disequilibrium across different pairings of populations and loci. We removed one locus, Rhyvar_56006, from analyses of R. kezeri, due to its high number of populations out of Hardy-Weinberg Equilibrium. Additionally, we removed the locus Rhyvar_36876 from R. kezeri analyses, as it was monomorphic in 13 of 17 populations. Thus, all subsequent analyses were performed using 8 of 10 loci for R. kezeri. Based on Δlog ML from the BAPS analysis for R. kezeri, we found the most support for K = 3, however, when taking into consideration the biology of the species and the restriction of Δlog ML approach to K ≥ 3, it was determined that K = 2 is most appropriate (Online Resource 2). Two populations north of the Columbia River, 10 and 11, cluster with the populations to the south, despite the substantial barrier between them (Online Resource 3). Our individual-based analysis revealed that these populations are highly admixed, with fewer than 60% of individuals in each assigned to the northern cluster. Thus, landscape genetic analyses were performed among the populations dividing the sampling localities into two spatially distinct clusters separated by the Columbia River (Fig. 1). Pairwise D PS among localities were high overall while F′ ST was low on average (mean ± SE; R. kezeri south ( For R. variegatus, there were no departures from Hardy-Weinberg equilibrium after Bonferroni correction, yet there was evidence of linkage disequilibrium in 8 of 45 pairs of loci. However, no loci were consistently in disequilibrium and thus analyses of R. variegatus were performed with all 10 loci. All further analyses of R. variegatus were performed separately on the three clusters identified previously (Emel and Storfer 2015), which were maintained with the addition of data from new populations examined in this study ( Summary genetic diversity statistics indicated overall lower genetic diversity within the R. kezeri than in the R. variegatus populations sampled (Tables 1, 2). Mean A R for R. kezeri populations was lower than the mean for R. variegatus (p = 0.007). Similarly, mean H o was lower for R. kezeri than for R. variegatus (p = 0.005). Mean F IS was lower for R. kezeri than for R. variegatus (p < 0.001), indicating lower inbreeding in R. variegatus. In contrast, mean H e was higher in R. kezeri than in R. variegatus (p = 0.05). Additionally, mean A e did not differ significantly between species.
Additionally, genetic diversity varied among the genetic clusters within species. The northern cluster for R. kezeri (Cluster 2) trended toward higher mean A R , A e, H o , and H e , than the Cluster 1 (Table 1). There was no clear trend for these metrics among the R. variegatus clusters, although A R , A e, and H e were lowest in the northernmost cluster (Cluster 3; Table 2). There was no clear trend in F IS within species.

Forest fragmentation
Analyses of four class metrics of forest cover suggested a higher level of fragmentation in the northernmost portions of both species' ranges (R. kezeri Cluster 2 and R. variegatus Cluster 3) than in the southern portion of their ranges (Table 3). Specifically, the proportion of the landscape covered in forest was lower for these two clusters relative to the others, and concordantly, patch density, landscape division index, and splitting index were all higher in these clusters as well. The proportion of forested landscape trended toward increasing with decreasing genetic diversity, but there is no clear pattern (Tables 1,  2, 3). Overall, the class metrics do not suggest a trend toward higher degree of fragmentation in the R. kezeri study area overall, despite the lower average genetic diversity in this species.

Landscape genetic analysis
In all cases except for the northern cluster of R. variegatus (Cluster 2), Circuitscape models based on landscape variables outperformed the null model of IBD (Table 4). Across clusters and genetic distance measures, categorical land cover was the most frequently supported model, at various and often multiple cost ratios. Road cover also appeared multiple times among the top models, at 1:10 and 1:100 cost ratios. Additionally, models based on percent canopy cover, heat load index, growing season precipitation, and roughness were all supported. Compound topographic index, length of the frost-free period, and notably stream cover were variables not included in the top models.
There were also consistent differences in the most supported Circuitscape models across clusters as well as species. In southern R. variegatus (Cluster 1), categorical land cover was always the top model, whereas in northern Cluster 2, categorical land cover was equally weighted with the null model, suggesting that distance alone may drive genetic differentiation in that cluster. In contrast, road cover was the exclusive variable in the top models for the northernmost Cluster 3, at a cost ratio of 1:10 or 1:100, suggesting a strong effect of roads on gene flow. For R. kezeri, there was much more variation across clusters and genetic distance metrics than in R. variegatus. Heat load index with a higher emphasis placed on the high values (T1) was the sole top model for southern R. kezeri (Cluster 1) using D PS , yet several models were within ΔAIC < 2 for F′ ST , including categorical land cover, road cover, and both untransformed and T1 percent canopy cover. The northernmost cluster, R. kezeri Cluster 2, differed in that growing season precipitation and roughness (both T2) were the most supported models for D PS and F′ ST , respectively, indicating that these variables are particularly impactful at low values.

Discussion
Our comparison of two closely related, ecologically similar species suggests that landscape genetic patterns can be consistent across broad spatial scales, yet ecological differences in habitat characteristics and disturbance can differentially affect landscape genetic structure at small scales. Overall, we determined that land cover and roads are the strongest predictors of genetic distance in these two torrent salamanders, but within the genetic clusters in each species there is variation in the relative importance of these and other variables related to minimizing desiccation. Our results also suggest that decreased genetic connectivity and lower genetic diversity in R. kezeri may be associated with disturbances that reduce forest cover, however, similar patterns of forest fragmentation were found in the northern portion of the study areas for both species. Forest cover can be affected by timber harvest, other land uses, landslides, and wildfire, and in the area analyzed, timber harvest is the most pervasive disturbance to landscape-scale forest cover (Nickerson et al. 2011). Given the smaller geographic range of R. kezeri compared to R. variegatus, the effects of fragmentation in a large proportion of this area may have stronger negative effects on observed genetic diversity. Furthermore, the specific landscape context, variation in degree of human development, biotic interactions, or evolutionary history may play a role in shaping genetic diversity.

Comparative landscape genetics
By comparing two or more closely related species, there is potential to either uncover cryptic differences among populations or species, or support expected similarities related to shared evolutionary histories among lineages. In this study, Circuitscape models produced largely consistent results across two congeneric salamander species with similar habitat requirements. Of the 10 tested variables, two were frequently among the top models in both species: categorical land cover and road cover. All three resistance ratios tested for categorical land cover appeared among the top models, and in some cases all three were within ΔAICc < 2, indicating that they are equally supported. Interestingly, where it appears, road cover is always found at both a 1:10 and 1:100 cost ratio, but never 1:2, suggesting that the true resistance of roads is likely much higher than twice that of non-road areas. Nonetheless, identifying the exact resistance values of each variable is beyond the scope of this comparative study, and could be achieved within study regions in future work through the application of methods such as Resistance GA (Peterman 2017).
Overall, our results point to a strong link between land cover, particularly human-altered landscape features, and gene flow in these two torrent salamander species. However, our findings for variable importance for R. kezeri varied much more than those for R. variegatus, suggesting gene flow was more dramatically affected by variation in landscape and environmental disturbances as well. In addition to categorical land cover and road cover, percent canopy cover appears to drive patterns of genetic structure in the southern R. kezeri cluster. The prevalence of categorical land cover over percent canopy cover may highlight areas where land conversion from forest to rural communities acts as a major barrier to gene flow, because of the difference in how these two resistance layers are constructed. Categorical land cover involves high resistance values for non-forested vegetated land, and even higher resistance for human-developed land, despite both of these categories also having low canopy cover. Our results suggest that any forest cover, regardless of its density, is important for dispersal. Additionally, the topographic variables heat load index and roughness were only significant for R. kezeri, suggesting that either R. kezeri individuals are more vulnerable to solar radiation, dry substrates, and energetically expensive changes in elevation than those of R. variegatus, or that these environmental factors vary more across the R. kezeri distribution. Furthermore, these conditions become increasingly important in the absence of forest cover, as desiccation becomes more likely. The importance of growing season precipitation could also reflect a similar pattern. However, IBD was selected in the northern R. variegatus cluster in addition to categorical land cover, suggesting that Euclidean distance alone can be an important driver of genetic structure. Finally, F′ ST and D PS yielded similar results within all R. variegatus clusters with moderate-to-high marginal R 2 values for all models (marginal R 2 = 0.51-0.76; Table 4). These findings would suggest that genetic distances in this study are robust to departures from the assumptions of F ST and its analogs. Still, these two genetic distance measures produced varying results in both R. kezeri clusters, with a sharp contrast in marginal R 2 indicating that the landscape data fit D PS (marginal R 2 = 0.63-0.93) better than F′ ST (marginal R 2 = 0.25-0.43). This finding may relate to the lower genetic diversity in R. kezeri.

Importance of landscape variables
The landscape genetic results for R. variegatus presented here differ from previously published work that investigated fewer landscape variables over a smaller geographic area using both Circuitscape and least-cost paths (Emel and Storfer 2015). Previously, Emel and Storfer (2015) reported that low percent canopy cover, low stream cover, high heat load index, steep slope, and a long annual frost-free period predict high genetic distance within two R. variegatus genetic clusters including the same southern cluster used here and a subset of the populations in the northern cluster. Interestingly, IBD was selected as a top model both in the northern cluster here and in the subset of the populations analyzed in the previous study. However, in the present study with a larger sample size over a greater geographic extent, we also included additional landscape variables, including compound topographic index and roughness, as well as categorical land cover and road cover (each at 3 alternative cost-ratios). Although by limiting our categorical variables to 3 cost-ratios it is possible that we do not capture the true cost of each cover type, doing so would be beyond the scope of this study. Our goal is to identify the variables that best explain genetic distance within genetic clusters and thus may be of importance to these species' conservation in the future, using cost-ratios motivated by the species' biology and previous research Storfer 2008, 2010;Trumbo et al. 2013). Furthermore, three of four of these new variables were present in the top landscape genetic models across species, with categorical land cover being the most supported model overall, and a higher proportion of the variation in genetic distance explained. We find that categorical land cover is nearly always a better model than percent canopy cover, suggesting the overall importance of this variable. Thus, we present a refined set of predictor variables of R. variegatus gene flow compared to Emel and Storfer (2015), and expand to include a species of greater conservation concern, R. kezeri. Our results demonstrate the importance of testing multiple variables in landscape genetic studies, so long as they are chosen a priori based on the species' biology . If the variables that are truly driving spatial genetic patterns are omitted, researchers could make spurious conclusions based on the subset of variables that they examine, which may be correlated with the true variables.
We predicted that streams would facilitate gene flow due to their importance in reducing desiccation risk during dispersal. However, unlike Emel and Storfer (2015), we did not find streams to be an important variable in explaining genetic structure in this study. There are a number of potential explanations for this discrepancy. First, rain may facilitate overland movement to the point that there is a stronger signal for other variables than for streams. In a study of a similarly desiccation-sensitive salamander (Plethodon albagula), Peterman et al. (2014) found unexpectedly that genetic distance was negatively correlated with increased desiccation risk across the landscape, leading to conclusions that rapid nighttime dispersal (associated with higher humidity) through unsuitable habitat may strongly affect gene flow. Moreover, the role of stream influence on gene flow may depend on stream size, with headwater streams facilitating gene flow, but higher order-streams and rivers acting as barriers due to predation by larger salamanders and fish Olson 2001, 2003) and/or lack of suitable sites for egg deposition. The role of stream order in upstream versus downstream dispersal, and the potential compounding effects of road crossings, could be tested by stratifying sampling within individual streams and using capture-recapture methods to model unidirectional gene flow. Alternatively, it is possible that Circuitscape resistance does not adequately model linear stream elements in the same way that least-cost paths do, however, stream Circuitscape models were still highly ranked in Emel and Storfer (2015). Therefore, this counterintuitive result demonstrates the value of landscape genetic approaches to reveal new insights about species' biology and inspire topics for further investigation.

Genetic diversity and conservation implications
Torrent salamanders (Rhyacotriton spp.) have specific habitat requirements, and are typically associated with forests and the cool, highly oxygenated water of headwater streams (Diller and Wallace 1996;Welsh and Lind 1996;Sheridan and Olson 2003;Olson and Weaver 2007;Olson and Burton 2014). This specificity makes them particularly vulnerable to habitat alteration and the potential effects of climate change. Our results suggest that both R. kezeri and R. variegatus use forested habitat for dispersal, avoiding roads, and exhibit high levels of genetic structure. Thus, fragmentation of these habitats by disturbances such as timber harvest can further isolate populations from one another, as landscape resistance in harvested areas increases beyond current levels.
Rhyacotriton kezeri has a much more restricted geographic range than R. variegatus, and significantly lower average genetic diversity within populations (based on lower A R , H e, and F IS ; Tables 1, 2). Decreases in gene flow due to habitat fragmentation can promote inbreeding, which may be exacerbated by population declines due to reductions in the quality of habitat present in the remaining habitat patches. Class metrics revealed evidence of greater forest fragmentation within the northernmost study regions in both species' ranges (Table 3), despite the lower genetic diversity in R. kezeri overall, and at the cluster level, genetic diversity does not closely reflect levels of fragmentation. However, the overall pattern of decreased fragmentation in the other two R. variegatus genetic clusters indicates that fragmentation could pose a greater overall threat to R. kezeri. Still, we cannot discount the potential effect of ascertainment bias due to applying markers developed for one species to another, albeit closely related, species. The reduction of the marker set from 10 to 8 loci due to lack of diversity at one locus and a high proportion of populations out of Hardy-Weinberg equilibrium in the other indicates that genetic diversity in this species may be underestimated, especially in A R . Further analysis with a larger set of microsatellite or SNP markers could rule out this factor, as well as more clearly define the genetic clusters in this species. Although the presence of null alleles may inflate heterozygosity in most populations, lower H o in R. kezeri also points to a potential increase in homozygosity due to genetic isolation. Thus, a combination of variation in landscape genetic relationships between the two species and variation in the degree of habitat fragmentation and quality may explain the difference in genetic diversity, and further study is warranted.
Nonetheless, these findings suggest an uncertain future for both species as they occur within a highly managed forest landscape where timber harvest is the dominant land-use activity (Nickerson et al. 2011), particularly when the species' aquatic habitat and upland forest dispersal habitat are not addressed in landscape-scale planning. Importantly, this work can aid designs of connectivity pathways in managed forest landscapes. Burnett (2009, 2013) provided conceptual designs for landscape-scale dispersal of these forested headwater species, extending and connecting protected riparian corridors up and over ridgelines at strategic locations across watersheds. Specific environmental correlates identified here provide empirical support for identification of more efficient dispersal routes for these species. Considering the projected threat of climate change on these species (Bury 2015), and the need to design climate-smart forests in the region (Kim et al. 2017), maintenance of forested stream and streamside habitats and over-ridge dispersal habitats will likely be essential to the conservation of these unique salamanders. Insights from this study may inform the status reviews for these species, especially for R. kezeri, which is under consideration for listing under the U.S. Endangered Species Act (WSDOT 2017).