The wolverine (Gulo gulo) is the largest terrestrial member of the mustelid family and is globally distributed throughout the boreal zone of the northern hemisphere. Wolverine occurrences worldwide are strongly associated with cold areas offering persistent spring snow cover, and year-round habitat use of the species generally is restricted to this narrow bioclimatic niche (Aubry et al. 2007; Copeland et al. 2010; Inman et al. 2013). In the contiguous United States, historical wolverine populations were extirpated by the early 20th century, and contemporary populations in the lower 48 States were reestablished by immigrants from the north (McKelvey et al. 2014). The species still is absent from Utah and only single individuals have been confirmed in California (Moriarty et al. 2009), Colorado (Packila et al. 2017), and Oregon (Oregon Wild 2019), all of which were part of the historical distribution of the species (McKelvey et al. 2014). Currently, the species occupies insular, high-elevation habitat patches in Idaho, Montana, Washington, and Wyoming (Banci 1994; Aubry et al. 2007; McKelvey et al. 2014). In this southernmost part of their distribution, wolverines live at low densities (e.g., Inman et al. 2012a) and have an estimated effective population size of only 28–52 individuals (Schwartz et al. 2009). Compared to wolverines in Canada and Alaska, wolverines in the contiguous United States also show lower levels of genetic diversity, and higher levels of genetic differentiation (Kyle and Strobeck 2001, 2002; Cegelski et al. 2003, 2006; Schwartz et al. 2007; McKelvey et al. 2011). To ensure the long-term viability of wolverines in the lower 48 States, it is particularly important to understand and maintain population connectivity and gene flow across broad spatial scales (Carroll et al. 2001, 2015; Schwartz et al. 2009; Inman et al. 2013).

Landscape genetic approaches make it possible to study the effects of landscape heterogeneity on dispersal and population connectivity in highly mobile, low-density species such as the wolverine (Storfer et al. 2007; Waits et al. 2016). A typical way for evaluating such landscape-genetic relationships is through landscape resistance modeling (Spear et al 2010, 2016). In this approach, the study landscape is modeled as a raster layer where each cell receives a value reflecting the hypothesized resistance of the landscape to animal movement and gene flow, for example due to physical or behavioral constraints (Zeller et al. 2012). This model then is used to estimate spatial proximity among sampling locations via “effective distances” (e.g., through least-cost paths), which account for the hypothesized landscape resistances. Finally, various statistical methods can be used to test for relationships between measures of genetic connectivity (e.g., genetic distances) and the estimates of landscape connectivity (i.e., effective distances—Wagner and Fortin 2016; Fortune et al. 2018; Shirk et al. 2018). A significant relationship between genetic and effective distances indicates that the modeled landscape resistance predicts spatial-genetic structures, a pattern often referred to as “isolation-by-resistance” (IBR—Cushman et al. 2006; McRae 2006). IBR often predicts genetic structures better than classic “isolation-by-distance” (IBD—Wright 1943), where genetic distances increase with increasing geographic (i.e., straight-line) distances among sampling locations. While IBD assumes that the landscape is homogeneous from the perspective of the study species, IBR accounts for species-specific landscape heterogeneity.

Such a landscape genetic approach was used to show that the distribution of persistent spring snow cover predicts wolverine gene flow in the contiguous United States (Schwartz et al. 2009). Wolverines have a strong association with the area defined by persistent spring snow cover, and the species depends on deep, persistent snow for denning and food caching (Copeland et al. 2010; Inman et al. 2013). Areas with persistent spring snow cover also predict wolverine occurrences during summer months, highlighting the year-round association of the species with a narrow bioclimatic niche (Copeland et al. 2010).

However, in addition to its narrow bioclimatic niche, gene flow in wolverines likely can be predicted by other landscape variables; for example, those related to terrain or anthropogenic infrastructures (e.g., May et al. 2006; Webb et al. 2016; Heim et al. 2017; Kortello et al. 2019). To understand landscape effects on wolverine dispersal and gene flow, and to effectively maintain meta-population connectivity (e.g., through conservation corridors), it is therefore important to evaluate the influence of additional landscape variables on wolverine dispersal and gene flow.

Dispersal and gene flow in mammals can arise from multiple underlying processes that act across different spatial scales (e.g., Travis et al. 2013). Specifically, mammals usually disperse over small spatial scales, while long-distance dispersal in contrast is rare (Sutherland et al. 2000; Whitmee and Orme 2013). Occasional long-distance dispersal events often require individuals to venture into unfamiliar areas, where they encounter landscape elements and influencing factors not present in their typical habitat (e.g., Walton et al. 2018; Bartoń et al. 2019). On account of these factors, several recent studies have shown that it is challenging to predict gene flow from habitat suitability models (e.g., Mateo-Sánchez et al. 2015; Roffler et al. 2016), because these models capture the suitability of the landscape for movement and occurrence of the species, but not necessarily the resistance of the landscape to long-distance dispersal movements. Often, individuals still are able and willing to disperse through certain areas, even if these areas are unsuitable for settlement and breeding (Mateo-Sánchez et al. 2015; Keeley et al. 2016; Abrahms et al. 2017). Current landscape genetic approaches do not distinguish among these different processes leading to gene flow, and rarely test for scale effects on landscape-genetic relationships.

Here, we extend the analysis of Schwartz et al. (2009) by including additional landscape variables in a multivariate, multiscale landscape genetics framework that can detect environmental predictors of genetic structure resulting both from typical, small-scale genetic exchange within a habitat, as well as exceptional, long-distance dispersal among habitats. Specifically, our study addresses three interrelated questions: (1) What is the relative importance of pure space (IBD) versus landscape heterogeneity (IBR) for predicting wolverine genetic structure? (2) What is the relative importance of different landscape variables for predicting wolverine genetic structure? (3) Are scale effects detectable in the data set and, if so, how do they affect the answers to questions (1) and (2)?

Methods

Genetic data set.—The genetic data set used in this study previously was analyzed by Schwartz et al. (2009) and includes 210 wolverine tissue samples gathered across northern and central Idaho, western Montana, and northwestern Wyoming. The study area covers approximately 330,000 km², including the entire Rocky Mountain range of wolverines in the contiguous United States (Fig. 1). Most samples are from harvested animals (N = 142), with the remainder stemming from different research projects. Samples were genotyped at 16 microsatellite loci with an expected and observed heterozygosity of 0.504 and 0.473, respectively. Mean number of alleles per locus was 4.06; a total of 65 alleles were detected. Additional details about sampling and laboratory analyses are detailed by Schwartz et al. (2009). This data set was used to calculate Miller's genetic distance D_ij for all pairs of individuals i and j with the software Alleles in Space (Miller 2005). This distance is the individual-based equivalent to the commonly used measure D for estimating genetic differentiation among populations (Nei et al. 1983) and takes on values between 0 (identical genotypes) and 1 (no allele in common). Since landscape genetic conclusions can depend somewhat on the choice of the genetic distance metric employed (Shirk et al. 2017), we also estimated alternative measures and compared results obtained with each genetic distance. Specifically, we estimated Rousset's a (Rousset 2000) in software Spagedi (Hardy and Vekemans 2002), and Smouse and Peakall's r (Smouse and Peakall 1999) in software Genalex (Peakall and Smouse 2005).

Fig. 1.

The study area includes states of Idaho (ID), Montana (MT), and Wyoming (WY) and borders Canada (CAN) in the north. Darker colors depict higher wolverine-specific landscape resistances inferred from the multivariate, broad-scale analysis (see main text for details). Locations of genetic wolverine samples are shown as white circles. The inset shows the location of the study area in the western United States.

Landscape data and suitability for providing connectivity.—As a first step for landscape resistance modeling, we chose eight digital landscape layers of variables that we hypothesized to be predictors of wolverine functional connectivity (Table 1). These variables were among those used by Inman et al. (2013) to predict habitat suitability for wolverines across the western United States. We chose these variables because wolverines in our study area occupy rugged, high-elevation habitats that generally follow the timberline (Hornocker and Hash 1981; Copeland 1996; Edelmann and Copeland 1999; Squires et al. 2006; Aubry et al. 2007; Copeland et al. 2007; Inman et al. 2012a, 2013). These habitats include conifer forests as well as forest edge, e.g., areas where forested and nonforested areas are adjacent (Brock et al. 2007). Wolverines also appear sensitive to anthropogenic influences (Copeland 1996; Rowland et al. 2003; Heim et al. 2017), and require deep spring snow-pack for denning (Magoun and Copeland 1998; Aubry et al. 2007). Variables chosen for this study thus included a climatic variable (snow depth), two topographic variables (terrain ruggedness, and latitude-adjusted elevation), three anthropogenic variables (human population density, housing density, and road density), and two vegetative variables (forest edge, and conifer cover). Further information on the different data layers can be found in  Supplementary Data SD1.

Table 1.

GIS data layers and sources used .in the dispersal suitability models.

We did not know a priori whether all of these landscape variables actually predict wolverine genetic structure, nor could we translate the landscape data into exact resistance values. However, we were able to hypothesize whether the relationship between the variable and genetic structure should be positive (increasing genetic connectivity with increasing variable values), or negative (decreasing connectivity with increasing variable values). For example, we hypothesized that if housing density predicts genetic structure in wolverines, this association most likely would be negative, with higher housing densities predicting reduced functional connectivity for wolverines. Thus, we rescaled the housing density layer so that highest values became 0 (lowest suitability for providing connectivity) and lowest values became 1 (highest suitability for providing connectivity). We rescaled the other anthropogenic data layers in the same way, thus hypothesizing that a potential relationship between these variables and wolverine functional connectivity would be negative. For all other data layers, we generally hypothesized that a potential relationship between the variable and connectivity would be positive (i.e., data were rescaled so that highest values became 1 and lowest values became 0). The only exception was the latitude-adjusted elevation index, which we used in two different ways. First, we calculated the absolute elevational difference between any cell in the landscape and the elevation of the alpine timberline at the latitude of that cell and rescaled the data to range from 0 (highest elevational difference) to 1 (elevation of timberline, “absdiff” in Table 1). This rescaling assumes that wolverine dispersal and gene flow generally follow the alpine timberline, as the occurrence of the species is associated with high-elevation habitats along this line (Squires et al. 2006; Aubry et al. 2007; Copeland et al. 2007; Inman et al. 2012a). However, this rescaling also assumes that even small differences between actual elevation and the elevation of the timberline reduce the suitability for successful wolverine dispersal, which is a rather restrictive assumption. We thus also rescaled latitude-adjusted elevation in an alternative way, where all areas within 300-m elevational difference to the timberline received a value of 1 (highest suitability for providing connectivity, “timberline” in Table 1). This use of elevation data for landscape resistance modeling is similar to that of other landscape genetic studies (e.g., Cushman et al. 2006; Shirk et al. 2010; Wasserman et al. 2010).

The rescaling approach only makes very general assumptions about the potential suitability of an area for providing connectivity via dispersal and ensures that the different variables are measured on the same scale (0 for lowest suitability, 1 for highest suitability), which is important for assessing their relative importance for predicting genetic structures.

Converting landscape data into resistance models.—Next, we converted the rescaled data into actual landscape resistance values, using the concept of weighted distances (Singleton et al. 2002). In this approach, each cell of the landscape is assigned a weighted distance value that reflects the hypothesized resistance of that cell to movement and gene flow. When the landscape in a cell is not assumed to provide resistance, the weighted distance is simply the length of the shortest possible path through the entire cell. In our case, with a cell size of 1 km², this shortest path has a length of 1,000 m. In contrast, when the landscape is assumed to have high resistance for successful dispersal, the weighted distances is 100 times greater than the shortest possible path through the cell (i.e., 100,000 m in our case). We used this approach to convert our rescaled landscape data into resistance values, so that landscape cells with low suitability for providing connectivity (rescaled data value = 0) received highest resistance values. In contrast, landscape cells with highest potential for successful dispersal (rescaled data value = 1) received lowest resistance values. Following Mateo-Sánchez et al. (2015), we used three different translations for this step.

First, we used a linear conversion to derive weighted distances from our rescaled data (see Fig. 2, dashed graph):

where WD is the resistance assigned to a specific cell (i.e., the weighted distance), CS is the length of the shortest possible path through the cell (“cell size,” i.e., 1,000 m), and DS is the value of the rescaled data (“dispersal suitability,” range 0–1).

In addition to this linear translation, we also defined two exponential conversions which account for possible nonlinear relationships between landscape data values and actual resistance (Fig. 2, dotted graph):

This second conversion assumes that decreased suitability for connectivity via dispersal increases landscape resistance most strongly when suitability is low. This is the kind of translation that has been recommended for deriving landscape resistances from habitat suitability values (Keeley et al. 2016), because animals often are still able and willing to move through nonhabitat areas, even if these areas are not suitable for settlement or breeding (e.g., Mateo-Sánchez et al. 2015).

Fig. 2.

Calculation of landscape resistances based on hypothesized dispersal suitability using linear (dashed line), upper exponential (solid line), and lower exponential (dotted line) conversions.

Finally, we used an inverse-symmetrical version of the previous formula, which assumes that reduced suitability for dispersal quickly results in higher landscape resistances (Fig. 2, solid graph):

Estimating effective distances from landscape resistance models.—Based on the resulting 27 resistance surfaces (9 rescaled landscape layers × 3 weighted distance translations), we calculated effective distances using two approaches. First, we calculated pairwise least-cost distances among individual wolverines with the ArcView extension PATHMATRIX (Ray 2005). Least-cost distances measure effective separation among all pairs of samples based on a single, optimal (i.e., least-costly) path through the resistance landscape (Adriaensen et al. 2003). Thus, the least-cost approach implicitly assumes that effective dispersal and resulting genetic structures are influenced by optimal animal choices that are based on substantial knowledge of the landscape.

In addition to least-cost paths, we also calculated pairwise effective resistances in program CIRCUITSCAPE 3.4 (McRae and Shah 2013). Effective resistances are based on circuit theory and estimate effective separation distances in terms of current flows across the resistance landscape (McRae 2006; McRae and Beier 2007). This approach assumes that movement of individuals is not optimal and based on incomplete knowledge of the landscape, so that multiple paths contribute to effective dispersal. Since inferences derived from the two approaches can differ (e.g., Avon and Bergès 2016; McClure et al. 2016; Marrotte and Bowman 2017), we applied both approaches and used objective criteria to choose the best approach for our data (see next section).

Model optimization and exclusion of spurious variables.—We next conducted an exploratory analysis to assess how best to represent each variable (i.e., which resistance translation to use and whether to use least-cost distances or effective resistances), and to select variables for the final, multivariate analysis (i.e., excluding spurious variables). For this, we first used simple Mantel tests (Mantel 1967) to evaluate the correlation of each effective distance with the genetic distances. For each environmental variable, the representation giving highest correlation coefficients was chosen. Any variable that did not lead to significant correlations was excluded from further analyses. Next, we used partial Mantel tests (Smouse et al. 1986) to determine whether the best representation of each variable is significantly correlated with genetic distances after accounting for the effects of space (i.e., geographic distances). Thus, we controlled for the effects of geographic distance and then assessed whether effective distances based on environmental variables have an influence on genetic structure “beyond IBD.” Thus, we essentially used IBD as a null model, which assumes that environmental heterogeneity does not impact genetic structure (Jenkins et al. 2010). Finally, simple Mantel tests among all effective distances were calculated, and if two effective distances showed a correlation of r > 0.7, one of them was excluded from further analyses (i.e., the one that gave lower correlations with the genetic data—Dormann et al. 2013). Because each effective distance was used for multiple statistical tests, inferences were based on Bonferroni-corrected P-values. This approach ensures that only effective distances are retained for further analyses that are (1) significantly correlated with genetic distances “beyond IBD,” and (2) not significantly correlated with any other effective distances.

Assessing relative variable importance.—Simple and Mantel tests are appropriate for the exploratory data analysis described above, but they are not suitable to assess the relative importance of different landscape variables (Cushman et al. 2013a; Prunier et al. 2015). Hence, we used a combination of alternative methods to analyze the final, multivariate data set. All analyses were based on multiple regression on distance matrices (MRM), which is essentially a multivariate regression performed on the unfolded, lower-triangle distance matrices (Legendre 2000; Lichstein 2007). The method is one of the most effective ways for analyzing landscape genetic data (Balkenhol et al. 2009; Wang 2013), and ideally is combined with variance partitioning techniques to infer relative variable importance (Prunier et al. 2015, 2017). Here, we combined MRM with two types of variance partitioning. First, to distinguish IBD effects from actual landscape influences, we used the spatial variance partitioning approach of Borcard et al. (1992). In this approach, the variation explained by a regression model is partitioned into three different components. Spatial variation is the amount of variation that is solely due to endogenous factors, such as landscape-independent space-use or species-specific dispersal limitations, i.e., IBD, which is represented by straight-line geographic distances. Environmental variation is the component that can be attributed to purely exogenous factors, such as the heterogeneous landscape, i.e., IBR, which is represented by effective distances (i.e., least-cost paths or circuit-theoretic resistances). Finally, there is a spatial component of environmental factors, so that there also is a shared component, which cannot be separated into purely environmental versus purely spatial contributions (e.g., because of correlations between straight-line and effective distances). Separating these three components via spatial variance partitioning commonly is used in many other ecological disciplines (e.g., Wagner 2004; Wagner and Fortin 2005; Lichstein 2007). In landscape genetics, spatial variance partitioning can be used to answer the question: What is the relative importance of space (IBD) versus environment (IBR) for predicting spatial-genetic structure?

Second, to assess the relative importance of the different effective distances for predicting genetic distances, we used hierarchical partitioning (Chevan and Sutherland 1991; Mac Nally 2000, 2009), which is one the most effective approaches for assessing relative variable importance in ecological studies (Murray and Conner 2009). The approach assumes that a hierarchy of models exists, ranging from the simplest, univariate models to the most complex model involving all variables. The method considers all possible models that can be constructed with the set of candidate variables, and estimates the increase in model fit generated by a certain variable for each model in which the variable is included. Hierarchical partitioning can be seen as an alternative to model selection (e.g., based on AIC values), which cannot be applied to pairwise distance matrices (Franckowiak et al. 2017). In landscape genetics, hierarchical partitioning can be used to answer the question: What is the relative importance of different landscape variables for predicting spatial-genetic structure? Since we were particularly interested in assessing landscape-genetic relationships independently from any IBD effects, we conducted hierarchical partitioning with the original data, and with residuals from regressions of geographic distances on all other distances (i.e., partial regression—Wagner and Fortin 2016).

Testing for scale effects.—Finally, we tested for scale effects in the data set based on spatial autocorrelograms. We partitioned the data set into spatial lags, which included only those pairwise comparisons within a certain distance class or “bin.” We used Sturges' rule to objectively determine the optimal number of bins (Sturges 1926), and assessed the correlation between genetic and geographic distances separately for each distance bin (Smouse et al. 1986; Goslee and Urban 2007). If spatial-genetic structure is homogeneous across scales, all bins should show a similar pattern, i.e., the correlogram should show a relatively straight line and little variation in spatial-genetic correlation coefficients. Because we found two distinct spatial scales in the wolverine data set where spatial-genetic correlation coefficients differed significantly (see “Results”), we repeated the spatial and hierarchical variance partitioning for these two distinct scales, thus assessing the influence of scale effects on our inferences.

All statistical analyses were conducted in the R software, using packages ecodist (Goslee and Urban 2007) and hier. part (Walsh and Mac Nally 2009). When significance tests were needed, P-values were estimated via 9,999 permutations, and 95% confidence intervals were based on 1,000 bootstraps (Goslee and Urban 2007).

Results

Our optimization approach produced models that performed well in predicting wolverine genetic structure, with Mantel correlations and amounts of explained variation in the multiple-matrix regressions being at the high end of values reported for individual-based studies (see, e.g., Wang et al. 2008). Our results were consistent across the different genetic distances we used, with all of them being significantly correlated to each other (r > 0.93). We therefore only report results based on Miller's D, as it gave slightly higher Mantel r-values than the other genetic distances.

Model optimization and exclusion of spurious variables.—All simple Mantel tests were highly significant (P = 0.0001), and the highest Mantel correlations obtained for each variable are shown in Table 2. Effective resistances based on circuit theory generally led to higher correlation coefficients than least-cost distances, except for timberline elevation and conifer cover (Table 2;  Supplementary Data SD2). Resistance models derived from the higher exponential translation led to slightly higher correlation coefficients for most variables (Table 2), but strong effects of the different translations only were apparent in three variables (i.e., road density, timberline, conifer cover;  Supplementary Data SD2). In the partial Mantel tests, all variables were significant after Bonferroni correction, with the exception of conifer cover and absolute elevation difference (Table 2). These two variables were therefore dropped from further analyses.

Table 2.

Highest Mantel correlations obtained for each variable. *“Best” refers to the conversion and effective distance that led to the highest r-values. CT = effective resistances based on circuit theory; LCD = least-cost distances. Variable abbreviations as in Table 1, geo is the geographic (i.e., straight-line) distance. Note that in the partial Mantel tests, a Bonferroni-corrected P-value of 0.0056 is required for significance at alpha = 0.05, meaning that treecov and absdiff are insignificant.

Effective distances based on the remaining seven environmental variables showed no collinearity with each other (i.e., Mantel r < 0.7, P ≥ 0.05), except for population density and housing density (Mantel r = 0.74, P = 0.0001). We chose to delete population density and retain housing density for further analyses, because the latter resulted in a higher correlation coefficient with genetic distances. Thus, we retained six environmental variables for our final analyses (i.e., elevation—timberline, terrain ruggedness—tri, snow depth—snow, housing density—hden, road density—roaden, and forest edge—edge), plus straight-line geographic distance (geo).

Scale-dependent results.—When using all pairwise comparisons, the MRM model explained 13.5% of the variation, and was highly significant (P < 0.001). The spatial variance partitioning of the full MRM model showed that pure environmental factors had a strong impact on the explained variation ( Supplementary Data SD3). About 64.7% of the explained variation can be attributed to pure environmental components (IBR), 23.9% to shared spatial-environmental components, and only 11.4% to pure spatial components (IBD). In terms of individual variable importance, hierarchical partitioning suggested that snow depth, terrain ruggedness, and housing density, best predict genetic distances across the 127 possible models ( Supplementary Data SD4; maps of these three most influential variables are shown in  Supplementary Data SD5). However, strong scale effects were detected in the genetic structure, as the Mantel correlogram (Fig. 3) suggested positive autocorrelation in the first five distance classes (i.e., up to ∼230 km), while negative autocorrelation existed in the last five distance classes (i.e., beyond distances of ∼420 km). In between these distances, there is no significant spatial autocorrelation. Thus, wolverines were genetically more similar than expected under spatially random mating at small spatial scales (up to ∼230 km distance), while they were less similar than expected at large spatial scales (beyond ∼420 km distance). Analyzing these two spatial scales separately is particularly interesting, because different landscape-genetic relationships could exist within each scale.

Based on these results, we partitioned the data set into two groups: one including only pairwise comparisons belonging to the first five distance bins (“small scale”), and the other including only those distances belonging to the last five distance bins (“broad scale”). We then repeated the two variance partitioning analyses for small and broad scale separately.

The spatial variance partitioning analyses suggested substantial differences in the relative importance of spatial and environmental variables between the two scales (Fig. 4). While both MRM models explained similar percentages of variation (14.7% for small scale, and 13.5% for broad scale;  Supplementary Data SD3), the combined effect of space and environmental heterogeneity were much more important at larger scales. Pure environmental components (IBR) were responsible for 78.9% of the explained variation at small scales, and 64.6% at broad scales. Pure space (IBD) accounted for 17.1% and 11.4% of the variation at small and broad scale, respectively. Shared components accounted for only 4.1% at fine scales, but for 24% at broad scales.

The relative importance of the different environmental variables also varied with scale (Fig. 5). At the small scale, snow depth was the most important variable, with a 42.9% contribution to the variance explained. Terrain ruggedness contributed slightly more than straight-line geographic distances (i.e., IBD) at this scale (16.6% and 13.2%, respectively), while the contribution of road density was similar to that of geographic distance (13.6% and 13.2%, respectively). All other variables contributed less than geographic distance to the genetic structure at this small scale.

At the broad scale, the contribution of snow depth to explained variance became much weaker (7.4%), while the importance of terrain ruggedness increased (24.9%). The most pronounced difference between the two scales occurred for housing density, which contributed only little at the small scale (12.1%), but was the most important variable at the broad scale (54.6%). Final resistance values for long-distance gene flow among wolverine in our study area are shown in Fig. 1, as inferred from the broad-scale model results. To create this final resistance layer, we added the best resistance layers of the six environmental variables, after weighting them by their relative importance at the broad scale (in % contribution to R²; Fig. 5).

Fig. 3.

Spatial Mantel correlogram used to infer scale effects. The number and width of the equally sized spatial bins was determined using Sturges' rule (N = 15, width of bin = 47 km). Black dots indicate Mantel correlation coefficients (r) that are significantly different from zero after Bonferroni correction; white dots indicate Mantel r values not different from zero.

Fig. 4.

Spatial variance partitioning for small and broad scales. Shown as percentages of explained variation.

Discussion

Dispersal is a key process in ecology and evolution (Clobert et al. 2012; Travis et al. 2013) and various approaches have been developed to study dispersal in the wild (i.e., Cayuela et al. 2018). Landscape genetic approaches are particularly useful for analyzing environmental influences on effective dispersal, i.e., dispersal leading to gene flow (Broquet and Petit 2009; Robertson et al. 2018). Since dispersal can be impacted by multiple environmental drivers and their relative importance can vary across spatial scales (e.g., Morton et al. 2018), it is important to evaluate the effects of multiple environmental variables in landscape genetic studies and to quantitatively test and account for scale effects on landscape-genetic relationships.

Fig. 5.

Relative variable importance for partial models at two distinct scales based on hierarchical partitioning. Variable importance is given in % contribution of each variable to overall variation explained. The black dotted line shows the contribution of straight-line geographical distance in the full models at the small scale (i.e., 13.2%; see  Supplementary Data SD4). The gray dashed line shows the contribution of geographical distance in the full models at the broad scale (i.e., 12.7%; see  Supplementary Data SD4).

Here we used a multivariate, multiscale approach for a detailed investigation of landscape-genetic relationships in wolverines inhabiting the contiguous United States. Several previous studies have demonstrated that wolverines in the contiguous United States exhibit lower levels of genetic diversity, and higher degrees of genetic structure than their conspecifics in Canada and Alaska (Kyle and Strobeck 2001, 2002; Cegelski et al. 2003, 2006). While the homogeneous boreal and consistently snowy areas in the north are less fragmented and can support larger effective population sizes of the species (Schwartz et al. 2009), suitable wolverine habitat in the contiguous United States only can be found in several high-elevation mountain ranges, which are separated by low-elevation basins with limited amounts of snow and nonpersistent spring snow cover (Aubry et al. 2007; Inman et al. 2013). However, the comparably low levels of gene flow in the southernmost parts of the species' range sometimes also have been attributed to increased anthropogenic impacts, which potentially have led to additional habitat fragmentation and reduced wolverine densities in this area (Guillot et al. 2005; Cegelski et al. 2006).

Our results suggest that wolverine genetic structure in the region can indeed be predicted by a combination of bioclimatic, anthropogenic, and topographic, landscape characteristics, that these environmental factors are much more important for predicting wolverine gene flow than pure space (IBD), and that the relative influence of these variables is scale-dependent. Deriving this inferential detail only was possible through our multivariate and multiscale analyses, which greatly increased the ecological meaningfulness of our study.

Methodological considerations for inferring landscape effects on gene flow.—Using MRM and hierarchical partitioning allowed us to avoid the use of significance testing (e.g., based on Mantel tests) for our final inferences, and instead assess the relative importance of each individual variable in a multivariate context, as recommended by Prunier et al. (2015, 2017). From an ecological standpoint, this is arguably more informative than trying to find the single “best” landscape resistance model for predicting genetic structures, as often is done in landscape genetics. Furthermore, testing and accounting for scale effects proved vital for our analyses and led to a more detailed understanding of the potential processes underlying successful dispersal of wolverine within versus among habitat complexes. For the wolverines in our study area, our multiscale approach enabled us to differentiate between landscape influences on small-scale gene flow (likely reflecting typical dispersal and mating movements within habitat complexes) versus landscape influences on broad-scale gene flow (reflecting long-distance dispersal or multigenerational gene flow among habitat complexes). This kind of inferential detail would not have been derived from the analyses involving all scales, as the model involving all scales could not clearly disentangle the relative importance of individual variables for predicting genetic structure ( Supplementary Data SD4). Our analytical framework thus can help to disentangle environmental impacts on the different processes leading to gene flow and results in a biologically meaningful distinction between habitat suitability and landscape resistance. These findings allow us to discuss possible causes of detected landscape-genetic relationships at fine and broad spatial scales.

Fine-scale predictors of wolverine gene flow.—In our fine-scale analyses (i.e., up to ∼230 km), pure spatial and pure environmental components explained most of the variation in genetic distances. This scale is well within reported wolverine dispersal distances of, e.g., up to 170–178 km (Vangen et al. 2001; Inman et al. 2012a). Accordingly, the strong impact of pure space at this small scale likely is not due to an actual, physical dispersal limitation (classical IBD—Wright 1943), but probably is caused by wolverine space-use and mating patterns (Dalerum et al. 2007). Wolverines show the typical space-use and polygamous mating patterns of mustelids, with male home ranges overlapping those of several females (Banci 1994; Dalerum et al. 2007; Persson et al. 2010). Also, wolverines may establish their home ranges in proximity of their natal range (Hedmark et al. 2007; Bischof et al. 2016; Aronsson and Persson 2018). These behavioral characteristics can lead to strong spatial-genetic correlations at fine spatial scales, because related individuals tend to be closer geographically than unrelated individuals (see, e.g., Wang et al. 2017; Kristensen et al. 2018, for other species). Because these are endogenous (i.e., landscape-independent) processes, they explain the relatively high percentage of pure spatial components in our small-scale analysis, despite the high mobility of the species.

In terms of landscape predictors of spatial-genetic structure, our small-scale results are congruent with those of Schwartz et al. (2009), who used data on the presence or absence of persistent spring snow cover to show that areas lacking such persistent snow cover have a resistance to gene flow approximately 20 times higher than areas with persistent spring snow cover. Because our weighted resistance approach for modeling resistance values required us to work with continuous variables, we used snow depth data in our study instead of categorical persistent spring snow cover data. Despite this difference in modeling snow-related landscape resistances, our results support the conclusion that deep, persistent snow cover predicts successful wolverine dispersal and resulting genetic structures. Indeed, snow was the most important predictor for wolverine gene flow across all scales, and especially up to distances of ca. 230 km. Most likely, gene flow at this scale is primarily determined by mating movements and natal dispersal of wolverine offspring, and these processes largely occur within primary habitat complexes, i.e., within mountain ranges. These complexes are typical wolverine habitats and often characterized by low housing densities and high topographic heterogeneity (e.g., May et al. 2006; Inman et al. 2013). Thus, at this spatial scale, wolverine gene flow cannot be predicted very well by housing and topography, as they simply are not “limiting factors” (sensu Cushman et al. 2013b) of dispersal within typical wolverine habitats. Instead, successful wolverine dispersal mainly is determined by the availability of high-quality breeding habitat, with deep, persistent spring snow cover as a very important predictor (Magoun and Copeland 1998; Copeland et al. 2010).

Other environmental variables that we investigated (i.e., conifer cover, forest edge habitat, elevation, and road densities) showed statistically significant, but relatively weak, relationships with genetic distances. However, these results do not contradict previous studies that identified these variables as important for wolverine occurrence and movement (Aubry et al. 2007; Copeland et al. 2007; Inman et al. 2013; Webb et al. 2016; Kortello et al. 2019; Sawaya et al. 2019), because these studies were conducted more locally or in different regions.

Broad-scale predictors of wolverine gene flow.—At scales beyond typical wolverine dispersal distances (i.e., > 420 km), housing density and terrain ruggedness were the best predictors of spatial-genetic structures in wolverines. These spatial distances likely are only covered by multigenerational gene flow or occasional long-distance dispersal events (Gardner et al. 1986; Moriarty et al. 2009). For example, an exceptional straight-line dispersal distance of 826 km for an adult male was reported by Packila et al. (2017), and male dispersal up to 500 km was detected using genetic means by Flagstad et al. (2004). To maintain connectivity over such large distances, wolverines eventually have to disperse outside of their typical, mountainous habitat, and potentially cross areas that are characterized by higher housing densities and areas that are topographically more homogeneous (e.g., Packila et al. 2017). Our results suggest that these factors strongly influence population genetic structure in wolverines at broad spatial scales. Given that anthropogenic features can affect wolverine space-use behavior, including habitat selection, at large spatial scales (e.g., Rowland et al. 2003; May et al. 2006; Stewart et al. 2016; Heim et al. 2017), it seems plausible that effective wolverine dispersal is higher across areas with low human housing density and associated human presence. In contrast, it is unlikely that wolverine dispersal actually is limited by topographic homogeneity per se, because wolverines show high levels of genetic connectivity in areas without much topographic variation (i.e., taiga and tundra—Wilson et al. 2000; Kyle and Strobeck 2001, 2002). Instead, the broad-scale relationship between genetic connectivity and topographic ruggedness observed in our study might be caused by comparably little human influences and associated infrastructures in areas with high topographic heterogeneity. For example, wolverines in Canada preferred more rugged terrain with less human development (Fisher et al. 2013), just as wolverines in Norway avoided areas with infrastructural development (May et al. 2006).

The increasing availability of suitable sites for food caches in more rugged terrain is an alternative or additional explanation for the importance of terrain ruggedness in predicting spatial-genetic structure in wolverines. Wolverines often cache food in cold, structured sites protected from competitors, such as bacteria, insects, and other scavengers, for example, under downed logs or large boulders (Magoun and Copeland 1998; Inman et al. 2012b). Inman et al. (2013) showed that terrain ruggedness was a strong predictor of wolverine resource selection and hypothesized that areas with increased terrain ruggedness offer more cliffs, boulders, and talus areas, and generally more topographic habitat structure than areas with low ruggedness, and showed that terrain ruggedness was indeed a strong predictor of wolverine resource selection. Thus, it seems possible that long-distance gene flow across topographically more variable areas is higher than across flat areas, because more rugged areas offer better caching opportunities, hence are more suitable for the species.

Overall, the relationship between terrain ruggedness and genetic connectivity detected in our study warrants further investigation such that potential mechanisms behind this landscape-genetic relationship can be clarified. Nevertheless, based on our results we can conclude that terrain ruggedness is a strong predictor of wolverine genetic structure at broad spatial scales, although not as strong as housing density.

Conclusions for wolverine dispersal and population connectivity.—Our study provides insights into gene flow among wolverine at different spatial scales, and improves our understanding of dispersal ecology and functional connectivity in the species. Our methodological approach allowed us to analyze effective dispersal (i.e., gene flow) in the species at distinct spatial scales, thereby enabling us to detect the different environmental predictors of genetic structure at small and broad scales. While previous studies had already shown the strong association between gene flow and deep, persistent snow (Schwartz et al. 2009), or suggested anthropogenic and topographic effects on population genetic structure without explicitly quantifying them (Cegelski et al. 2003, 2006; Guillot et al. 2005), we demonstrated the importance of terrain ruggedness and housing density for predicting broad-scale wolverine gene flow. Had we only analyzed the all-scale model, we would have concluded that snow depth is the best variable for predicting genetic structure. While this finding is not incorrect, it provides an incomplete picture, because deep snow in the region only occurs within typical wolverine habitat, where most gene flow occurs across usual dispersal and mating distances. In contrast, low housing density and terrain ruggedness occur outside of wolverine habitat and individuals only have to move through these areas during occasional long-distance dispersal. While not common, long-distance dispersal among habitat complexes is vital for maintaining meta-population viability and species persistence (Trakhtenbrot et al. 2005; Heinz et al. 2006). Our finding that housing density best predicts gene flow across broad spatial scales thus is of particular relevance for maintaining and improving connectivity among wolverine habitats. Specifically, our results suggest that connectivity among wolverine will be particularly sensitive to housing developments and other human impacts in rugged areas located between typical wolverine habitat. Even if these areas are not usually inhabited by wolverines, they are potentially used during dispersal and can therefore offer crucial pathways for gene flow across broad spatial scales. Efforts to identify optimal corridors among wolverine habitats in the region should accordingly consider housing density and terrain ruggedness as predictive variables of successful wolverine dispersal.

Our multivariate, multiscale analyses led to a more detailed understanding of wolverine dispersal in a heterogeneous environment. We suggest that the approach is generally useful for landscape genetic research in mammals, particularly for large scale studies involving highly mobile species, as such studies likely will capture gene flow at different spatial scales, and both within and between typical habitat of the study species.