Functional connectivity in a continuously distributed, migratory species as revealed by landscape genomics

Maintaining functional connectivity is critical for the long-term conservation of wildlife populations. Landscape genomics provides an opportunity to assess long-term functional connectivity by relating environmental variables to spatial patterns of genomic variation resulting from generations of movement, dispersal and mating behaviors. Identifying landscape features associated with gene flow at large geographic scales for highly mobile species is becoming increasingly possible due to more accessible genomic approaches, improved analytical methods and enhanced computational power. We characterized the genetic structure and diversity of migratory mule deer Odocoileus hemionus using 4051 single nucleotide polymorphisms in 406 individuals sampled across multiple habitats throughout Wyoming, USA. We then identified environmental variables associated with genomic variation within genetic groups and statewide using a stepwise approach to first evaluate nonlinear relationships of landscape resistance with genetic distances and then use mixed-effects modeling to choose top landscape genomic models. We identified three admixed genetic groups of mule deer and found that environmental variables associated with gene flow varied among genetic groups, revealing scale-dependent and regional variation in functional connectivity. At the statewide scale, more gene flow occurred in areas with low elevation and mixed habitat. In the southern genetic group, more gene flow occurred in areas with low elevation. In the northern genetic group, more gene flow occurred in grassland and forest habitats, while highways and energy infrastructure reduced gene flow. In the western genetic group, the null model of isolation by distance best represented genetic patterns. Overall, our findings highlight the role of different seasonal ranges on mule deer genetic connectivity, and show that anthropogenic features hinder connectivity. This study demonstrates the value of combining a large, genome-wide marker set with recent advances in landscape genomics to evaluate functional connectivity in a wide-ranging migratory species.


Introduction
Maintaining natural levels of functional connectivity is essential to sustaining healthy wildlife populations (Cosgrove et al. 2018, Tucker et al. 2018. Functional connectivity describes the influence of the environment on the movement of individuals (Taylor et al. 1993). It can be difficult to manage functional connectivity of highly mobile terrestrial organisms, especially those that migrate seasonally, because they require extensive, contiguous habitat and often move across jurisdictional boundaries (Bolger et al. 2008, Middleton et al. 2020. Various approaches to estimate functional connectivity exist, each of which requires different data types and leads to different inferences. GPS tracking and habitat selection studies reveal the influence of environmental variables on daily and seasonal movements of migratory species (Lindström et al. 2014, Monteith et al. 2018. Landscape genetic studies aim to understand the influence of environmental variables on patterns of genetic variation which can result from generations of movement and mating behaviors (Manel et al. 2003), thereby providing longer-term evaluations of functional connectivity than direct tracking methods (Kool et al. 2013). As habitat loss and fragmentation increasingly reduce movement options for highly mobile species, it is important to understand how natural and human-impacted environments influence gene flow (Tucker et al. 2018).
Mule deer Odocoileus hemionus are migratory ungulates that occupy diverse habitats throughout their North American range (Mackie 1994). Mule deer populations have declined or fluctuated considerably over the past several decades due to factors including extreme weather, habitat loss and degradation, and predation (Bishop et al. 2009, Hurley et al. 2011, Monteith et al. 2014. Some mule deer migrate seasonally to capitalize on spatial variation in food resources (Aikens et al. 2017), but whether to migrate and how far to migrate can vary across and even within populations (Sawyer et al. 2016). Despite mixed migration strategies, individual mule deer exhibit higher fidelity to their seasonal ranges and migratory routes than many other migratory ungulates (Sawyer et al. 2019). High fidelity to seasonal ranges and migratory routes could limit dispersal propensity and therefore reduce levels of genetic connectivity among mule deer populations. The state of Wyoming, USA includes diverse habitats that support different mule deer migration strategies. Western Wyoming encompasses several mountain ranges and valleys, where mule deer tend to migrate seasonally from higher-elevation summer ranges to lower-elevation winter ranges (Heffelfinger et al. 2003). Eastern Wyoming transitions to the Great Plains where mule deer are less likely to migrate seasonally due to lower spatiotemporal heterogeneity in food resources (Heffelfinger et al. 2003). The biogeographic variation in Wyoming provides an opportunity to investigate the influence of environmental variables, and reliant migration strategies, on genetic connectivity.
Habitat selection studies have highlighted the influence of topography, vegetation type and other environmental variables on mule deer movement (Nobert et al. 2016, Coe et al. 2018, Monteith et al. 2018, as well as anthropogenic disturbances such as roadways, residential sprawl and energy development (Northrup et al. 2015, Wyckoff et al. 2018, Dwinnell et al. 2019, Sawyer et al. 2020. Despite considerable movement research throughout the species' range, landscape genetic studies are limited to California where both natural ecological variation and human developed infrastructure are associated with mule deer gene flow (Pease et al. 2009, Mitelberg and Vandergast 2016, Fraser et al. 2019. Landscape genetic studies in other regions of the mule deer range, where demographics and landscapes are much different, would further understanding of mule deer ecology and inform management.
Identifying landscape features associated with gene flow at large geographic scales for highly mobile species is becoming increasingly possible due to more accessible genomic approaches, improved analytical methods and enhanced computational power (Manel andHolderegger 2013, Spear et al. 2016). The ability to generate large, genome-wide datasets of single nucleotide polymorphisms (SNPs) for non-model organisms (Peterson et al. 2012) enhances detection of finescale genetic variation, which is particularly useful for continuously distributed, highly mobile animals that have subtle levels of genetic structure (e.g. woodpeckers, Aguillon et al. 2018, mountain lions, Trumbo et al. 2019and marine fish, D'Aloia et al. 2020. Analytical improvements include methods to optimize resistance values associated with environmental variables based on genetic data itself (Graves et al. 2013, Peterman 2018, rather than relying on expert opinion, which can be biased or lacking for a study system (Zeller et al. 2012). Optimization strategies have also been improved by the ability to calculate resistance distances along multiple paths between individuals (McRae et al. 2008) and by the inclusion of nonlinear resistance surfaces that likely reflect relationships more accurately (Shirk et al. 2010, Zeller et al. 2017. Increasing computational power permits researchers to perform these memory-intensive analyses for larger study areas and at finer spatial resolutions, though computational feasibility remains a challenge (Zeller et al. 2012, Peterman et al. 2019. We combined these advances to determine the influence of the environment on the genetic structure of mule deer across a large geographic area in the core of the range, Wyoming, USA. We evaluated the influence of natural and anthropogenic environmental variables on gene flow at multiple geographic scales (i.e. within genetic groups and statewide) using individual-based comparisons. We hypothesized that fidelity to seasonal migration strategy would reduce genetic connectivity across the state, and therefore predicted that we would detect genetic structure related to seasonal ranges. We also hypothesized that environmental variables found to influence daily and seasonal movements also would affect gene flow.

Sample collection
Our study area spanned the state of Wyoming, USA (Fig. 1). During the autumn hunting seasons in 2014-2016, the Wyoming Game and Fish Department (WGFD) collected retropharyngeal lymph nodes from harvested mule deer as part of their disease surveillance program and later donated these samples to our laboratory. The WGFD recorded collection date, hunt area and harvest location (e.g. GPS coordinates, descriptive location), and used ArcGIS to assign GPS coordinates to samples with sufficient descriptive locations to confidently identify harvest location. Mule deer mating and seasonal migration both occur during the sample collection period, so we cannot discern which seasonal range animals were harvested on, though the large geographic scale of the study should compensate for variation in harvest location.
We selected 600 samples ( Fig. 1) to extract DNA from based on sex ratio (aiming for equal representation of males and females, though most harvested deer are males), sample location, and location data specificity (preference for GPS locations over descriptive locations). We included samples lacking GPS locations in genetic structure and diversity analyses but excluded them from landscape genomic analyses because the latter requires specific location data.

Library preparation and sequencing
We used Qiagen DNeasy Blood and Tissue Kits (Qiagen, Valencia, California) to extract DNA from selected samples, then quantified DNA concentrations using a BioTek Synergy HTX Multi-Mode Microplate Reader (BioTek Instruments, Winooski, Vermont). From the extracted samples, we selected 480 samples (297 males, 183 females) with the highest DNA concentrations and normalized each sample to 50 ng µl −1 for library preparation.
We prepared samples for reduced-representation sequencing following a protocol modified from Parchman et al. (2012). Briefly, we digested the DNA using EcoRI and MseI enzymes, ligated Illumina adapters and unique barcodes to digested fragments, and amplified the fragments using PCR. We used AMPure XP beads (Beckman Coulter, Brea, California) to clean and concentrate the PCR product and then size selected for 350-450 base pair fragments with a BluePippin (Sage Science, Beverly, Massachusetts). The University of Oregon Genomics and Cell Characterization Core Facility sequenced the prepared libraries on five lanes of an Illumina HiSeq4000 (Illumina, San Diego, California) with 150 base pair single-end reads.

Assembly and data processing
We filtered out contaminants and Illumina adapters using bowtie2 (Langmead and Salzberg 2012), then de-multiplexed reads. We aligned reads to a mule deer reference genome (16 208 scaffolds, N50 = 838 758, GenBank Accession: GCA_004115125.1, Russell et al. 2019) with the Burrows-Wheeler Alignment maximal exact matches algorithm in BWA (Li and Durbin 2009). After alignment, we identified and filtered SNPs using Samtools ). To exclude sequencing errors, we discarded SNPs with a mapping quality score below 20, with a Phred quality score below 19 or with a minor allele frequency below 0.01. We required at least three reads per sample to call a genotype at a given locus to maximize the number of loci retained with sufficient genotype certainty to calculate inter-individual genetic distances based on all loci combined (Buerkle andGompert 2013, Fumagalli 2013). To address linkage disequilibrium, we only retained the first SNP per read. To remove potential paralogs, we discarded SNPs with more than two alleles, with excess heterozygosity or with a sequencing depth greater than 100 (Willis et al. 2017). Because our analyses assumed neutral processes, we identified and removed loci potentially under selection using BayeScan with a false discovery rate of 0.1 (Foll and Gaggiotti 2008). Finally, we required loci to have data for at least 60% of samples, and discarded samples with data missing for more than 60% of loci. We checked for close relatives using VCFtools (Danecek et al. 2011).

Genetic structure and diversity
We tested for isolation by distance (IBD) using the individual-based spatial autocorrelation test in GenAlEx (Peakall and Smouse 2012). We binned sample pairs into 12 distance classes at 50-km intervals, and then calculated genetic distances among sample pairs in each distance class (Smouse and Peakall 1999). We calculated 95% confidence intervals using 999 bootstrap iterations and generated a null model of no spatial autocorrelation with 999 data permutations. We performed this test for all samples together, for samples within each genetic group, and for each sex separately because mule deer exhibit male-biased dispersal which could lead to different IBD results for each sex (Powell et al. 2013).
We characterized the genetic structure of mule deer in Wyoming using STRUCTURE, TESS3 and a principal components analysis (PCA). We used StrAuto (Chhatre and Emerson 2017) to run STRUCTURE (Pritchard et al. 2000) on the University of Wyoming 'Teton Computing Environment, Intel x86_64 group' (Advanced Research Computing Center 2018). Our analysis used the admixture and correlated allele frequencies models, included a burnin of 100 000 iterations and 500 000 Markov chain Monte Carlo iterations, and tested K = 1-10 possible genetic groups with 20 replicates of each K. We determined the most likely number of genetic groups by evaluating mean likelihood scores and ΔK in Structure Harvester (Evanno et al. 2005, Earl andvonHoldt 2012). We also ran the spatially explicit clustering program TESS3 in the R package tess3r (Caye et al. 2018), which accounts for geographic distance among samples to estimate ancestry coefficients (Caye et al. 2016). We tested K = 1-10 possible genetic groups with 20 replicates of each K and used default parameter settings. We performed a PCA based on pairwise genetic covariance calculated from loci shared by each sample pair to avoid potential bias associated with replacing missing data with population mean allele frequencies.
Using the R package diveRsity, we calculated pairwise F ST among groups with 95% confidence intervals calculated using 1000 bootstrap iterations (Weir andCockerham 1984, Keenan et al. 2013). We calculated observed and expected heterozygosity for each genetic group and statewide with the R package hierfstat (Goudet 2005). We estimated effective population size (N e ) with 95% confidence intervals for each genetic group using the linkage disequilibrium method in NeEstimator v2.1 (Do et al. 2014), and applied a correction for the number of chromosomes (n = 35; Wurster and Benirschke 1967) as suggested for N e estimates based on genome-wide SNP data (Waples et al. 2016). We did not report N e at the statewide level because the presence of substructure can inaccurately deflate N e estimates (Wang 2016, Gagne et al. 2018.

Landscape genomics
We used individual-based landscape genomics to investigate the relationship between landscape resistance (i.e. how much an environmental variable promotes or inhibits gene flow; Peterman 2018) and pairwise genetic distances (i.e. genetic dissimilarity, a proxy for gene flow; Smouse and Peakall 1999) among samples both within genetic groups and at the statewide scale (Manel et al. 2003). We excluded from landscape genomic analyses samples without GPS locations, and for analyses within each genetic group, we excluded samples that genetically assigned to the group, but were harvested outside of the general geographic area of the group to avoid biasing results with these geographic outliers.
We assigned resistance values to environmental variables using seven functional transformations and competed them against each other to determine which set of resistance values best matched the genetic data (Zeller et al. 2017). The seven transformations tested included positive linear, negative linear, positive monomolecular convex, negative monomolecular convex, positive monomolecular concave, negative monomolecular concave and inverse Ricker (Supporting information; Zeller et al. 2017). Positive and negative transformations represent increasing or decreasing resistance with increasing values of the environmental variable, respectively. The monomolecular functions allow for nonlinear relationships that may reflect landscape resistance more accurately. The inverse Ricker transformation represents low resistance at intermediate values of the environmental variable. For line features (e.g. highways), we only tested positive and negative linear transformations because these variables were binary (i.e. presence/absence of the feature). We used the transformed resistance rasters to calculate resistance distances among all pairs of samples using the commuteDistance function in the R package gdistance ( van Etten 2017), which applies circuit theory to average multiple paths between sample locations (McRae et al. 2008). We represented our null model of IBD with a pairwise geographic distance matrix calculated using the Haversine formula, which accounts for the earth's curvature, using the R package fossil (Vavrek 2011).
We selected ten environmental variables for our landscape genomic analyses based on previous mule deer studies focusing on landscape genetics (Pease et al. 2009, Fraser et al. 2019 or habitat selection (Sawyer et al. 2006, Lendrum et al. 2012, Northrup et al. 2015, Nobert et al. 2016, Coe et al. 2018, Monteith et al. 2018. We tested topographic features (elevation, slope, topographic roughness; Jarvis et al. 2008), land cover (proportion shrubland, proportion forest, proportion grassland; Homer et al. 2020), natural features (rivers; U.S. Geological Survey 2010) and anthropogenic features (highways, interstate highways, oil and gas wells; Wyoming Oil and Gas Conservation Commission 2012, U.S. Geological Survey 2014) (Fig. 2, Supporting information). We cropped all rasters to our study area plus a 50-km buffer around the state to reduce edge effects for samples located near study area boundaries. To standardize resistance values across rasters, we resampled rasters to a resolution of 270 m 2 (the finest computationally feasible resolution) and scaled resistance values to 1-100 (resampling methods in Supporting information; Hijmans 2020).
We used maximum-likelihood population-effects (MLPE) mixed models in the R package lme4 (Bates et al. 2015) to determine which functional transformation performed best for each environmental variable (Clarke et al. 2002, van Strien et al. 2012). The MLPE approach accounts for non-independence of pairwise comparisons by including individual as a random effect and outperforms other modeling approaches for individual-based landscape genetics (Shirk et al. 2018). We evaluated model performance using the Bayesian information criterion (BIC), which outperforms R 2 for ranking landscape genetic MLPE models (Row et al. 2017). To control for geographic distance inherent in all resistance-distance calculations and to reduce the likelihood of type I error, we included IBD as a covariate in all models (Row et al. 2017). Once we found the best performing functional transformation for each variable for each genetic group and across the state, we tested the resistance-distance matrices for collinearity. We retained variables with Pearson's correlation coefficients below 0.7 and variable inflation factors below five because collinearity can lead to incorrect interpretation in multivariate models (Dormann et al. 2013). Because we only included variables in our study that we thought would influence mule deer gene flow, we ran MLPE models for all variable subsets. All models with ΔBIC below five were considered top models. Top models were then refit using restricted maximum likelihood (REML) to estimate unbiased beta coefficients with 95% confidence intervals (Clarke et al. 2002, van Strien et al. 2012. Environmental variables with positive beta coefficients and confidence intervals not overlapping zero were considered significant. Row et al. (2017) found that variables with significant negative beta estimates typically indicated a non-true relationship in the context of landscape genetic modeling because if there was truly a negative relationship between the response and explanatory variables, then a negative functional transformation should have been selected as the top model and resulted in a positive beta estimate. Because we tested both positive and negative transformations of resistance values for every variable, we considered variables with negative beta estimates to be nonsignificant. We calculated marginal R 2 (mR 2 ) as a metric of model fit (Nakagawa and Schielzeth 2013) using the R package MuMIn (Bartoń 2020).

Assembly and data processing
Of the 480 samples sequenced, we discarded 19 samples with fewer than 12 000 reads after de-multiplexing; the next lowest sample had over 500 000 reads. After removing contaminants, the remaining 461 samples had 95.0-99.2% of reads align to the mule deer reference genome. Initial quality filtering resulted in two million SNPs, which was halved after filtering for maximum read depth and number of alleles. After filtering for minimum read depth and minor allele frequency, we had 340 000 SNPs, and then had 4069 SNPs after removing loci with data missing in more than 60% of samples. We removed two outlier loci and 16 loci with excess heterozygosity. Finally, we discarded 55 samples that had missing data at more than 60% of loci. We did not detect any close relatives to remove. Our final dataset included 406 individuals (250 males, 156 females) and 4051 loci.

Genetic structure and diversity
We detected low levels of positive spatial autocorrelation (r < 0.007) at distances up to 200-km, indicating the presence of IBD at shorter distances among sample pairs (Supporting information). We did not detect significantly different spatial autocorrelation among genetic groups or between sexes, so we report results for all samples together. Our STRUCTURE analysis showed the highest likelihood of 1-3 genetic groups (Supporting information) and the ΔK method supported three genetic groups (Fig. 3a, Supporting information). When K = 3, STRUCTURE admixture proportions revealed geographic clustering of genetic groups, roughly representing western, northern and southern regions of Wyoming, with admixture towards the center of the state and several geographic outliers that assigned highly to one genetic group but were harvested outside of the general geographic area of the group (Fig. 3c). When K = 2, the western and southern genetic groups collapsed into one group (Supporting information), suggesting the presence of some hierarchical structure. The PCA showed a single genetic cluster, but the three genetic groups identified by STRUCTURE were distinguishable with overlap among them (Fig. 3b). TESS3 could not clearly resolve the number of genetic groups (Supporting information), however, when we selected K = 3, the TESS3 genetic map (Supporting information) looked similar to the STRUCTURE assignment map (Fig. 3c). We proceeded with the three genetic groups designated in STRUCTURE for comparative landscape genomics, hereafter referred to as the western group, northern group and southern group.
F ST values were low, but significantly differed from zero among all pairs of genetic groups. F ST between the western and northern groups was 0.0051 (95% CI = 0.0037-0.0064), between the western and southern groups was 0.0042 (95% CI = 0.0028-0.0057), and between the northern and southern groups was 0.0049 (95% CI = 0.0036-0.0063). Observed and expected heterozygosity were similar among groups, and the western group had the highest estimated N e while the southern group had the lowest (Supporting information).

Landscape genomics
Our statewide modeling included 370 individuals with sufficient harvest location information. After removing geographic outliers for the analyses within genetic groups (Supporting information), modeling included 107 individuals in the western group, 136 individuals in the northern group and 100 individuals in the southern group. The dominant land cover type across the state was shrubland, however,  highest density of highways while the western group had the lowest (Supporting information). The northern group had the highest density of oil and gas wells (Supporting information). We report the best-performing functional transformation for each environmental variable in each analysis in Supporting information, and we report the global model for each analysis (after removing collinear variables) in Table 1.
For the statewide analysis, we identified three top models, which included proportion grassland, elevation and a combination of proportion grassland and elevation (Table 2, Supporting information). All the covariates in top models had positive beta estimates with confidence intervals not overlapping zero. Elevation had a positive linear function, meaning that more gene flow occurred at lower elevation. The selection of the inverse Ricker function for grassland indicated that more gene flow occurred in areas of mixed land cover or transitions from grassland to either shrubland or forest.
For the western group, the top model was the null model of IBD, and no other models had ΔBIC below five (Table 2, Supporting information). This result indicates that geographic distance among samples explained patterns of gene flow better than resistance from environmental variables in this part of the state. The model with the next lowest BIC (ΔBIC = 7.7) included proportion grassland, but the grassland variable was not significant.
For the northern group, we identified five top models that included proportion grassland, proportion shrubland, oil and gas wells, proportion forest and highways as significant covariates; the IBD covariate in each model was not significant (Table 2, Supporting information). Proportion grassland had a negative monomolecular convex function, indicating that more gene flow occurred in areas with high proportions of grassland. Proportion shrubland had a positive monomolecular convex function, indicating that more gene flow occurred in areas with low proportions of shrubland. Oil and gas wells had a positive monomolecular concave function, indicating that more gene flow occurred in areas with low to moderate densities of wells. Proportion forest had a negative monomolecular convex function, indicating that more gene flow occurred in areas with high proportions of forest. And lastly, highways had a positive linear function, indicating that more gene flow occurred in areas with fewer highways.
For the southern group, the one well-supported model included elevation as a significant covariate, and the IBD covariate was not significant (Table 2, Supporting information). Elevation had a positive linear function, indicating that more gene flow occurred in lower elevation areas.

Discussion
Detecting environmental variables associated with genetic patterns can be difficult for highly mobile, continuously distributed species (e.g. lobster, Benestan et al. 2015, black bears, Puckett andEggert 2016), but here we identified natural and anthropogenic landscape features associated with gene flow in mule deer. These features varied statewide and within each genetic group, demonstrating the presence of scale-dependent and regional variation in functional connectivity. Our top landscape genomic models included a variety of functional transformations of resistance values, bolstering support for the need to consider non-linear relationships between environmental variables and gene flow (Shirk et al. 2010). Overall mR 2 values were somewhat low (mR 2 = 0.003-0.06), but landscape genetics studies of other highly mobile mammals report similar values (Thatte et al. 2019, Trumbo et al. 2019) and the environment is one of many influences on spatial patterns of genetic variation (Manel et al. 2003).
We detected three genetic groups of mule deer with admixture among groups. The presence of geographically clumped genetic groups aligns with the high fidelity of mule deer to seasonal ranges and migration routes observed in GPS tracking studies (Sawyer et al. 2019), but admixture among the three groups, especially at the center of the state, shows that interbreeding or movement occurred between neighboring groups. We found several individuals that assigned highly to one group but were harvested in the geographic area of another group, which indicates some longer-range dispersal movements. We reported observed and expected heterozygosity and N e , which serve as benchmarks for comparison with future genome-wide SNP datasets.
Our modeling revealed that low elevation (statewide and in the southern group), grasslands (in the northern group) and transitional habitats (statewide) promoted gene flow, suggesting that winter range habitat is important for genetic connectivity. Most mule deer breed on their winter range, which is typically low elevation, relatively open habitat such as shrublands and grasslands (Heffelfinger et al. 2003). Table 1. Global landscape genomic models relating inter-individual genetic distances (GD) to landscape resistance for Wyoming mule deer collected in 2014-2016, using 4051 SNPs. Analyses included 370 deer statewide, 107 deer for the western genetic group, 136 deer for the northern genetic group and 100 deer for the southern genetic group. Environmental variables included elevation (elev), slope, proportion shrubland (shrub), proportion forest (forest), proportion grassland (grass), highways, oil and gas wells (wells) and geographic distance (IBD). We included individual (ind) as a random effect in every model to account for the nonindependence of pairwise comparisons. We tested all possible subsets of the global model for each analysis. We report ΔBIC as a metric of performance for each global model as compared to top models reported in Table 2.

Genetic group
Global model ΔBIC Statewide GD ~ elev + shrub + forest + grass + highways + wells + IBD + ind 37.7 Western group GD ~ slope + shrub + forest + grass + rivers + highways + wells + IBD + ind 47.2 Northern group GD ~ scrub + forest + grass + highways + wells + IBD + ind 32.9 Southern group GD ~ elev + shrub + forest + grass + highways + wells + IBD + ind 28.6 We also found, however, that high proportions of forest promoted gene flow in the northern region, which could suggest that summer range habitat is important for genetic connectivity because forested areas are more common in mule deer summer range (i.e. mountain ranges). During mild autumns, deer may breed while intermixed with other populations on summer range. Mule deer generally show high fidelity to their seasonal ranges and migratory routes (Sawyer et al. 2019), but the mixing of deer on summer range makes dispersal (i.e. young animal leaves natal region; Jakopak et al. 2019) or switching (i.e. adult animal changes seasonal migration strategy) events more plausible. In the western group, where more mule deer migrate than the rest of the state, the only top model was the null model of IBD. The western region includes several mountain ranges representing shared summer ranges for deer that spend winters in different foothills and valleys (Fig. 3). It is possible that deer mixing on shared summer range outweighs any influence of the landscape on gene flow in this area. Investigation of genetic structure in GPS-tracked deer could further elucidate the influence of different seasonal ranges on genetic connectivity.
We found negative effects of anthropogenic disturbance in the form of highways and oil and gas wells on gene flow in the northern group. Though highways are a relatively new addition to the landscape in evolutionary terms and were not detected as a genetic barrier in another Wyoming ungulate, pronghorn (LaCava et al. 2020), most wildlife-vehicle collisions in Wyoming and across the western United States involve mule deer (Huijser et al. 2009, Coe et al. 2015. Further, some roadways in Wyoming act as movement barriers to mule deer (Kauffman et al. 2018) and highways in southern California have also reduced mule deer gene flow (Mitelberg andVandergast 2016, Fraser et al. 2019). Our findings suggest that through time, such barrier and mortality effects may reduce genetic connectivity for mule deer.
The negative relationship between gene flow and well density in the northern group aligns with studies showing that mule deer consistently avoid energy infrastructure (Sawyer et al. 2017, Dwinnell et al. 2019. The northern group had the highest density of oil and gas wells by a large margin (Fig. 2i, Supporting information), which likely accounts for the inclusion of well density in a top model for Table 2. Top landscape genomic MLPE mixed models (ΔBIC < 5) relating inter-individual genetic distances to landscape resistance statewide and within genetic groups of Wyoming mule deer collected in 2014-2016, using 4051 SNPs. Functional transformations of environmental variables include positive linear (+ linear), positive monomolecular convex (+ convex), positive monomolecular concave (+ concave), negative monomolecular convex (− convex) and inverse Ricker (Supporting information). Marginal R 2 (mR 2 ) reported as a metric of model fit. Beta (β) estimates and 95% confidence intervals (95% CI) were calculated by refitting the model with restricted maximum likelihood. Positive 95% CIs not overlapping zero were considered significant. IBD = isolation by distance. this group specifically. Interestingly, the positive relationship between low elevation and genetic connectivity seen in other parts of the state was not present in top models for the northern group. Development is often concentrated in low elevation areas (Copeland et al. 2009), so it is possible that the negative development effect we detected in the northern group counteracted the benefit of low elevation habitat. As development continues to expand in Wyoming and throughout North America (Leu et al. 2008, Trainor et al. 2016, functional connectivity for mule deer and other migratory ungulates may be more difficult to sustain. Our study highlights environmental variables that promote and inhibit mule deer gene flow, providing crucial information for land-use planning and conservation efforts. Maintaining functional connectivity is central to wildlife management (Fletcher et al. 2016) and requires consideration of numerous biological patterns and processes (Kool et al. 2013). Combining genetic data with other lines of evidence (e.g. GPS tracking data) can improve decision making related to management such as habitat preservation (Zeller et al. 2017) or mitigating barrier impacts of roads (Balkenhol and Waits 2009).
We employed a stepwise approach to set resistance values empirically and consider non-linear functional relationships (Zeller et al. 2017), however, we were unable to perform full optimization of the resistance values (Peterman 2018) in our large study area at the desired resolution due to computational limitations, despite having access to a large computing cluster (Advanced Research Computing Center 2018). In addition, some of the best performing functional transformations for covariates seemed biologically unintuitive (e.g. in the statewide analysis, highways had a positive relationship with genetic distance while interstate highways had a negative relationship; Supporting information), however, the improbable transformations did not end up in top multivariate models (Table 2). If a variable does not influence gene flow, then the process of selecting the best functional transformation essentially involves selecting among functions that do not represent a significant relationship, and therefore it seems reasonable that an unintuitive function could be selected. In the future, researchers could consider either excluding improbable functional transformations completely or removing variables before multivariate modeling if an improbable functional transformation is selected as performing best. Despite these limitations, the methods used here can be applied to study functional connectivity of other highly mobile, continuously distributed species across large geographic areas. The variation in environmental variables important to gene flow among our analyses adds to a growing body of literature that demonstrates the need to evaluate functional connectivity at different geographic scales and in regions with different environmental characteristics (Short Bull et al. 2011, Robertson et al. 2018, Kozakiewicz et al. 2019).