The population structure of organisms across a landscape is a fundamental part of the natural history of any species that wildlife managers will want to know. Because the amount of genetic divergence among subdivided populations is directly related to the movement of individuals (Hartl & Clark 2007), estimating population structure provides an indirect quantification of connectivity across the landscape. Previous investigations of genetic divergence in mule deer Odocoileus hemionus have found that population structure corresponds to subspecies boundaries (Latch et al. 2009, Pease et al. 2009).

Mule deer are distributed across western North America, ranging from northern Canada to central Mexico, and from the Pacific Ocean to east of the Rocky Mountains (Anderson & Wallmo 1984, Mackie et al. 2003). By the end of the 19th century, mule deer populations had experienced declines in both size and extent across much of their range (Mackie et al. 2003). Within the state of Montana, these population declines led to the establishment of a translocation program for mule deer (Picton & Lonner 2008), which occurred concurrently with improvements in habitat quality and reductions in both predator densities and human harvest of deer (Mackie et al. 1998, Picton & Lonner 2008). Populations of mule deer peaked between the 1950s and early 1970s and now appear to be regulated by environmental conditions, hunting and predation pressure (Mackie et al. 1998, Picton & Lonner 2008).

Chronic wasting disease could greatly impact the future management of mule deer. Chronic wasting disease is expected to reduce mule deer population sizes (Dulberger et al. 2010), but might also reduce hunting pressure due to perceived health hazards (Williams et al. 2002). As of 2012, chronic wasting disease is present in mule deer populations in the states of South Dakota and Wyoming, USA, and the Canadian provinces of Saskatchewan and Alberta, all of which border Montana. Control strategies for reducing chronic wasting disease infections focus on disrupting the routes through which the disease is spread (Miller et al. 2006). The goal of this investigation was to describe the genetic population structure of mule deer sampled across Montana, with the ultimate goal of better understanding potential chronic wasting disease transmission across the landscape, and thus helping to inform potential chronic wasting disease management actions. Additionally, our study will contribute to the understanding of mule deer genetic structure across large geographic regions.

Material and methods

We collected samples of lymph node tissue, on a volunteer basis, at hunter check stations from 370 mule deer harvested during fall 2007 and 2008 across a study area aimed to encompass the entire state of Montana, with additional samples collected in Wyoming, Colorado and Utah, USA (Fig. 1). We also collected 19 samples from white-tailed deer Odocoileus virginianus harvested in northwestern Montana. At the time of collection, we recorded GPS coordinates for the approximate location of harvest, along with the sex of the deer. ArcMap version 9.3.1 (Environmental Systems Research Institute, Inc., Redlands, California, USA) was used to reproject all coordinates into NAD83 UTM zone 12. For an indication of the size of our study area, the largest distance between two deer in our analysis was 1,432 km.

We isolated DNA from lymph node tissue using Qiagen DNeasy 96 Tissue Kits and suspended DNA extracts in Buffer AE (Qiagen, Massachusetts, USA). We performed polymerase chain reactions (PCR) on 16 microsatellite loci in a multiplex format. Multiplex PCRs amplified four unique loci, each of which was labeled with a forward primer at the 5′ end for a specific dye (6FAM, VIC, NED and PET; Life Technologies, California, USA). Loci included in the four multiplexes were: 1) M, P, K and N, 2) D, Q, O and R, 3) E, BM4107, Rt30 and Rt7 and 4) G, Ovir, Rt24 and Cervid3 (Bishop et al. 1994, DeWoody et al. 1995, Wilson et al. 1997, Jones et al. 2000). A typical multiplex PCR reaction consisted of 1 μM of each primer, 5 μl 2X Qiagen Multiplex PCR Kit, approximately 50 ng DNA and enough water for a final volume of 10 μl (Qiagen). The thermoprofile consisted of one activation step at 95°C for 15 minutes followed by 40 cycles of 94°C for 30 seconds, 60°C for 90 seconds and 72°C for 60 seconds, and a final extension step at 72°C for 30 minutes. In addition, we sequenced 489 base pairs (bp) of the D-loop, or control region, of the mitochondria using primers developed by Latch et al. (2008). PCR chemistry consisted of 1 μM of each primer, 2 μl 5X MyTaq Reaction Buffer, 0.375 Units MyTaq HS DNA polymerase, approximately 50 ng DNA and enough water for a final volume of 10 μl (Bioline, Massachusetts, USA). The thermoprofile consisted of one activation step at 95°C for four minutes followed by 40 cycles of 94°C for 30 seconds, 55°C for 15 seconds and 72°C for 30 seconds and a final extension step at 72°C for one minute. We used ExoSAP-IT for PCR cleanup (Affymetrix). Sequencing reaction and subsequent cleanup were done using Big Dye Terminator v3.1 Cycle Sequencing Kit and ethanol/EDTA/sodium acetate precipitation as suggested by the manufacturer (Life Technologies). Both microsatellite fragments and D-loop DNA sequences were visualized using a 3100-Avant Genetic Analyzer (Life Technologies). Scoring of genotypes was performed using Genemapper v. 3.7 (Life Technologies), and sequences were manually aligned with the aid of Sequencher v. 4.1 (Gene Codes Corporation, Michigan, USA). Mitochondrial D-loop sequences were deposited in GenBank (accession numbers: JN040634-JN040649, JN040651-JN040684, JN040686-JN040687, JN040691-JN040693 and JN040697-JN040716).

We screened for the presence of hybridization between mule deer and white-tailed deer using the admixture model in the program STRUCTURE version 2.3.3 (Pritchard et al. 2000, Falush et al. 2003, 2007). We tested for the number of genetic clusters (K) using five independent runs of STRUCTURE for K between one and five, each run for 100,000 Markov Chain Monte Carlo (MCMC) iterations after a burn-in of 10,000 iterations. A recessive alleles model was used with missing genotypes entered as null allele homozygotes. We did not set prior membership for the 19 white-tailed deer collected in northwestern Montana. The ΔK method of Evanno et al. (2005) was used to select the appropriate K, after which we ran 10 independent runs of STRUCTURE using the same model parameters, but with a fixed K. Results from these final runs were compiled into a single consensus cluster membership using the FullSearch method in CLUMPP version 1.1.2 (Jakobsson & Rosenberg 2007), with results visualized using DISTRUCT version 1.1 (Rosenberg 2004). We classified hybrids as those individuals who had a consensus estimated genetic contribution from white-tailed deer of > 10%. Following the removal of hybridized and misclassified individuals, we estimated the number of genetic clusters in our study area using the same STRUCTURE parameterization as the initial analysis of the number of genetic clusters.

We compared observed genotype frequencies with Hardy-Weinberg expectations using the exact test of Guo & Thompson (1992). Population structure can affect genotype frequencies, so we analyzed three broad-scale sampling regions (Montana/Wyoming, Utah and Colorado) independently. We used the Markov chain exact test in GENEPOP on the Web (Raymond & Rousset 1995, Rousset 2008) to perform these calculations. Allele frequencies and unbiased expected heterozygosity for each sampling region were estimated using GENALEX version 6.41 (Peakall & Smouse 2006). We tested the homogeneity of mean expected heterozygosities using linear combinations. A total of three tests were run to compare the mean expected heterozygosity of each broad-scale sampling region to the common mean of the other regions. Allelic richness was calculated for a constant sample size of 20 genes in each sampling location using rarefaction with the program HP-RARE (Kalinowski 2005). Global F_ST was estimated, using the methodology of Weir & Cockerham (1984), among the three sampling regions with the program F-STAT version 2.9.3.2 and a 95% confidence interval was calculated using bootstrap resampling across loci (Goudet 2002).

To test for the presence of isolation by distance, we first calculated a count of dissimilar alleles between two deer at a given locus using equation 1 of Eding & Meuwissen (2001), except that we defined a success as an event where two alleles were different. Allele sharing distance (ASD) was then the estimated by summing these counts of allelic dissimilarity across loci and dividing by the total number of comparisons. Binomial logistic regression was used to estimate the linear relationship between Euclidian distance and the log odds of having dissimilar alleles. A Mantel test (Mantel 1967) with 5,000 random permutations was used to calculate a P-value for the test of no relationship between geographic distance and the log odds of having dissimilar alleles. A 95% confidence interval for the estimated change in the log odds of having dissimilar alleles as a function of Euclidean distance was calculated by inverting the Mantel hypothesis test as described in Manly (2007), using 10,000 random permutations.

We further tested for a relationship between genetic and geographic distance using the multivariate spatial autocorrelation methods of Smouse & Peakall (1999), as implemented in GENALEX version 6.41 (Peakall & Smouse 2006). We tested the null hypothesis of no spatial structure within equal distance classes of 10 km using 10,000 random permutations, and corrected P-values for false discovery (Benjamini & Hochberg 1995). We ran a principal coordinates analysis for the three sampling locations using the covariance-standardized method in GENALEX version 6.41 (Peakall & Smouse 2006).

We used a resampling test to explore whether translocations may have reduced the genetic distance between donor and recipient locations relative to historic levels. Translocation zones were set to include all individuals within 20, 40, 60, 80 and 100 km of the centroid of Montana counties that received translocated individuals (Picton & Lonner 2008). These removal distances were selected to account for the expected increase in the geographic area affected by a translocation event due to dispersal over time. We first calculated the observed slope of the logistic regression line between the log odds of having different alleles and Euclidean distance for the data after removing all individuals located within a specified distance of a translocation release site. We then randomly removed the same number of individuals and reestimated the slope of the regression line. The P-values for these tests were calculated as the proportion of 5,000 random removals with estimated regression slopes greater than or equal to the slope that was observed after removing individuals within the translocation zone.

We used a combination of mitochondrial DNA (mtDNA) and microsatellite markers to test for sex-biased dispersal in mule deer. We tested for historical sex-biased dispersal by comparing the genetic divergence observed between sampling groups at a sex-linked marker (mtDNA) to that observed at biparentally inherited markers (Prugnolle & de Meeus 2002). We analyzed a total of 489 bp of the mitochondrial control region in 76 samples from Montana that were pooled into five distinct sampling groups (see Fig. 1). Samples were pooled into these five groups because of their geographic separation within the larger sampling region of Montana. We used ARLEQUIN version 3.5 (Excoffier & Lischer 2010) to estimate F_ST for microsatellite genotypes and for the mitochondrial sequences following a multiple hits correction (Tamura 1992). Estimates of female-specific dispersal are confounded by the fact that uniparental inheritance reduces the effective population size of this marker to one quarter the size of biparentally inherited markers. This increases the genetic divergence expected at mitochondrial loci. We calculated the expected increase in genetic divergence based on both an island model of migration and a model where populations evolve in complete isolation, based on equations provided in Zink & Barrowclough (2008). Isolation by distance was estimated for these samples using Mantel tests with 5,000 random permutations (Manly 2007) comparing the ratio of genetic distance F_ST/(1-F_ST ) to the natural logarithm of geographic distance as proposed by Rousset (1997). We tested for contemporary sex bias in dispersal using pairwise ASD from samples located in Montana and Wyoming. Samples were split into female (N = 106) and male (N = 213) groups, and the difference in mean ASD between the sexes was assessed using a permutation test, the null distribution of which was developed by randomly permuting sex and recalculating the difference in mean ASD. This test is analogous to the F_ST-based test for sex-biased dispersal of Goudet et al. (2002). Unless otherwise stated, all analyses were performed in the statistical computing package R (R Development Core Team 2009).

Results

STRUCTURE analysis provided strong evidence that the tissue samples of putative mule deer collected in the field included white-tailed deer and hybrids between mule and white-tailed deer. When we included genotypes of deer known to be white-tailed in a STRUCTURE analysis, and varied the number of genetic clusters (K), the best supported number of clusters was two (one for each species as shown by the mode of ΔK; Evanno et al. 2005). This analysis suggested that seven samples labeled as mule deer in the field were actually white-tailed deer and that four additional deer appeared to be post first-generation hybrids between white-tailed deer and mule deer. The hybridized individuals in the data set had admixtures of 12.5%, 14.5%, 29.8% and 40.9% white-tailed deer genetic contribution. We removed all misclassified and hybridized samples from subsequent analyses, reducing the data set to 359 individual deer, 320 of which were located in the Montana/Wyoming sampling region, 29 in the Utah sampling region and 10 in the Colorado sampling region. A subsequent STRUCTURE analysis indicated a single genetic cluster existed across this region. For this analysis, the number of genetic clusters was selected based on the mean log posterior probabilities of the data given K.

Two of the loci we genotyped (O and Q) had an excess of homozygotes and were removed from subsequent analyses. Previous investigations have also observed an excess of homozygotes at these loci (Pease et al. 2009), which suggests there may be null alleles at these loci. One locus (P) was identified as out of Hardy-Weinberg equilibrium in the Montana/Wyoming sampling region from a Bonferroni corrected P-value ≈ 0.00. This locus was not removed from analysis because the F_IS at this locus (0.03) was within the range (0.007-0.078) of positive F_IS-values that were not identified as out of Hardy-Weinberg equilibrium, and this locus was in Hardy-Weinberg equilibrium in both other sampling locations.

Estimates of heterozygosity were similar across all sampling locations (Table 1; all P-values ≥ 0.47), with an average expected heterozygosity of 0.701 across our entire study region (95% confidence interval: 0.662-0.740). Average allelic richness was also similar among the broad-scale sampling regions (5.49 in Montana/Wyoming, 5.12 in Utah and 5.64 in Colorado). These broad-scale sampling regions were characterized by low levels of genetic divergence, with an estimated global F_ST of 0.012 (95% confidence interval: 0.008-0.016).

Data are consistent with weak isolation by distance in mule deer across this region (P-value = 0.0334; Fig. 2). We estimated that a 100 km increase in the distance between two mule deer was associated with a 1.0057-fold increase in the odds that randomly chosen alleles at the same locus in the two deer were different (95% confidence interval: 0.9995-1.0118).

We found no evidence in the spatial autocorrelation analysis to support spatial structuring of mule deer genotypes within our smallest distance class of 10 km (P-value = 0.033, adjusted α-level = 0.004). There was also no apparent clustering among the broad-scale regions in the principle coordinate analysis (Fig. 3); the first two axes of which explained 38.1% of the variation in genotypic distance.

We did not detect any evidence that historical translocations of mule deer in Montana have homogenized the population structure of this species. That is, a directed removal of samples located within a given proximity to translocation release sites did not increase the observed level of isolation by distance (smallest P-value was 0.057 for removing all individuals within 80 km of the translocation site, adjusted α = 0.01).

A total of 39 mitochondrial haplotypes were observed in the sample of 76 mitochondrial sequences (Table 2). Nucleotide diversity appeared to be reduced in samples collected from central and southwestern Montana relative to other sampling groups across the state (Table 3). We observed higher pairwise F_ST-values for mtDNA than for nuclear DNA (Table 4), four comparisons of which indicated recent geographic structuring of groups (Zink & Barrowclough 2008) or historic sex-biased dispersal (Prugnolle & de Meeus 2002) based on the fact that only the mtDNA F_ST indicated considerable genetic structure with estimates > 0.2 (Zink & Barrowclough 2008). We did not detect isolation by distance for either the microsatellite or mitochondrial genotypes in these smaller sampling groups (P-values = 0.100 and 0.197, respectively). However, a higher mean ASD at microsatellite loci was observed for female mule deer (ASD_FEMALE - ASD_MALE = 0.013, P-value = 0.04 from 5,000 permutations) across Montana and Wyoming.

Discussion

The weak isolation by distance we observed in Montana is consistent with previous dispersal estimates (Wright 1943, 1946), which indicate that both male and female mule deer often move long distances (Anderson & Wallmo 1984, Mackie et al. 2003). The fact that there are few, if any, barriers to dispersal, whether complete or permeable, is further indicated by the lack of observed autocorrelation of genotypes across our study area. Finally, the suggestive evidence of potential sexually dimorphic patterns of genetic divergence suggests that males are probably more important than females for the long-distance spread of chronic wasting disease into new regions. This possibility was suggested in a previous study of infected mule deer populations in Colorado (Miller & Conner 2005).

Levels of genetic structure similar to our results have previously been found in studies of mule deer population structure conducted at broad spatial scales. For example, Latch et al. (2009) found no phylogenetic structure below the subspecies level when examining mtDNA. Similarly, while Pease et al. (2009) found five genetic clusters of mule deer within the state of California, USA, these clusters corresponded with previous subspecies classifications. Cullingham et al. (2011) detected spatial autocorrelation of female mule deer extending up to 2 km in Alberta and Saskatchewan, Canada, but their estimated global F_ST (0.008) indicated similar levels of genetic structure to our study area. However, these authors did find evidence of two genetic clusters within their study area, and individuals assigned to clusters logically based on patterns of isolation by distance (Cullingham et al. 2011). Therefore, our study appears to fit well with previous work on mule deer, and in general indicates that populations are characterized by low levels of genetic structure below the subspecies level.

Cross et al. (2005) showed that disease transmission across a spatially-structured population is a function of migration rate and the infectious period of the pathogen. Duration of infectiousness is related to host immune responses, lifespan and environmental persistence of the pathogen. The observed high levels of connectivity among mule deer likely represent a conservative estimate of actual dispersal, because molecular methods only record successful reproduction following dispersal events (Cushman et al. 2006). Therefore, our findings of limited genetic isolation suggest high connectivity across the sampling area, which combined with the potentially long-term survival of prions in the environment (Brown & Gajdusek 1991, Seidel et al. 2007), suggest few, if any, barriers to chronic wasting disease spread.