Wildlife managers are challenged to manage spatially structured populations efficiently and effectively, therefore dispersal and gene flow are vital to understand and manage, particularly for a harvested species. We used a genetic approach to describe the metapopulation structure of Rocky Mountain elk (Cervus elaphus) in Idaho to assess past patterns of population distribution and influences of harvest. We used elk tissue and DNA samples (n = 216) to examine genetic dissimilarity between 7 regions and 9 elk management zones throughout Idaho using microsatellite loci (n = 11). Using 5 approaches, including pairwise FST-values, assignment tests, and a Bayesian model–based clustering of genotypes, we examined the distribution of genetic variation. The distribution of genetic variation between elk populations indicated low levels of genetic differentiation among regions (expected heterozygosity [HE] = 0.55–0.61, overall FST = 0.011) and elk management zones (HE = 0.54–0.60, overall FST = 0.017). Assignment tests and migration rates indicated directional gene flow between elk populations. A patchy metapopulation best describes the distribution of genetic variation among Idaho elk populations because likely enough individual interchange occurs between geographically separated populations. The elk populations we sampled could be part of a geographically larger patchy metapopulation potentially stretching from Yellowstone National Park through Idaho into western Canada. Because of historical translocations of elk from Yellowstone National Park, insufficient time may have passed to detect differences in genetic variation. Subtle differences in the distribution of genetic variation were observed in 2 of the 9 elk management zones within 2 different regions of the state. Our findings indicate management of Idaho elk populations and dispersal are maintaining sufficient gene flow. Metapopulation structure of a harvested species based on the distribution of genetic variation is an indicator of potential genetic consequences of harvesting and sustainable harvest levels.
The term “metapopulation,” first used by Levins (1970, p. 105), was defined as “a population of populations.” Levins' (1970) classical metapopulation (i.e., a set of populations that go extinct locally, but persist because of recolonization at the metapopulation level—Harrison and Taylor 1997) has been criticized as too limiting (Harrison 1991; Doak and Mills 1994) and expanded based on empirical evidence to suggest that the structure of metapopulations falls on a continuum of patch size and patch isolation to fully capture the diversity of spatially structured populations (Fig. 1; Stith et al. 1996; Harrison and Taylor 1997; Pannell and Charlesworth 2000). Because patch size and isolation or connectivity, or both, are dynamic processes varying temporally and spatially, a single model may not describe a species' metapopulation structure (Harrison and Taylor 1997; Hanski and Gaggiotti 2004).
Fig. 1.
Empirical models of metapopulation structure expanded from the classical metapopulation structure. Models fall along continua of patch size (i.e., size of area[s] occupied by a population) and patch connectivity or isolation, which is related to dispersal strategy. The classical metapopulation structure includes both occupied and unoccupied patches of equal size with patches sufficiently connected for recolonization of unoccupied patches to occur and the metapopulation to persist (Harrison 1991). The nonequilibrium metapopulation structure has all small patches that are so highly isolated no dispersal occurs. Each population is highly extinction prone and each acts as its own independent population and metapopulation (Harrison 1991; Hanski 1997). The patchy metapopulation structure has both small and large patches that are highly connected by dispersers and patches are united as a single population and metapopulation with a low probability of populations within patches going extinct (Harrison 1991). The mainland–island metapopulation structure is characterized by a large patch combined with smaller patches. The extent of patch isolation and dispersal strategy determines how small patches are combined into a population with the large patch or if each patch acts as an independent population (Harrison 1991). Figure is modified from Harrison and Taylor (1997) and Stith et al. (1996).
Use of the term metapopulation has evolved to include research and management that views populations as relatively discrete units (i.e., patches) on the landscape, frequently distributed over a large area, among which dispersal and gene flow occurs (Harrison 1994; Wells and Richmond 1995; Pannell and Charlesworth 2000; Hanski and Gaggiotti 2004; Hanski 2012). Four types of metapopulation structure have been described, including classical (i.e., moderate genetic differentiation between small patches because of moderate levels of connectivity that leads to recolonization of unoccupied patches), mainland–island (i.e., low to moderate genetic differentiation between patches), patchy (i.e., low genetic differentiation because of highly connected small to large patches), and nonequilibrium (i.e., high genetic differentiation because of small isolated patches [Fig. 1; Table 1]). Understanding the metapopulation structure of a species, particularly a harvested species, sheds light on the amount of connectivity (i.e., dispersal and gene flow) between populations. Relative similarity in genetic variation among populations can indicate which populations are interacting (i.e., exchanging dispersers) and detect key components of a species' metapopulation structure (Table 1; Cosentino et al. 2012; Roderick et al. 2012).
Table 1.
Genetic characteristics of metapopulation structure models and the associated genetic methods we used to empirically evaluate each model. Examples of each model are listed with citations. See Fig. 1 for further details on metapopulation structure models and see text for further explanation of genetic methods. Range of FST-values for low, moderate, and high genetic differentiation is based on Hartl and Clark (1997). HE = expected heterozygosity; HO = observed heterozygosity.
A species' metapopulation structure can be revealed by integrating information on levels of genetic similarity or dissimilarity among populations; patterns of movement and geographic barriers to movement; the spatial distribution of individuals, populations, or habitat; and statistical correlation in demographic rates between populations (Garton 2002; Garton et al. 2012). All of these factors merit research on how each individually and in combination can be used to describe a species' metapopulation structure. In the absence of habitat characteristics (e.g., forage quality or habitat fragmentation) or geographic barriers (e.g., mountain ranges or rivers), current genetic similarity captures past patterns of population distributions and successful immigration, which may or may not be reflected in current demographic rates (Bossart and Prowell 1998; Mills 2007; Lowe and Allendorf 2010). For example, spatially disjunct populations can have demographic rates fluctuating independently, but these populations could exchange dispersers, which would be evident in the distribution of genetic variation (Garton 2002; Lowe and Allendorf 2010). Furthermore, metapopulation structures, except for the nonequilibrium, could exhibit source–sink population dynamics, which could be revealed by the distribution of genetic variation among populations (Fig. 1; Pulliam 1988; Kawecki 2004; Runge et al. 2006; Andreasen et al. 2012; Mills 2013).
Large mammals such as the Rocky Mountain elk (Cervus elaphus) exemplify the challenge of landscape-level management because their seasonal and annual ranges often cover areas the size of watersheds or subbasins (Wisdom and Cook 2000). Therefore, the metapopulation structure of elk using a genetic approach could inform landscape-level management because spatially disjunct populations could arise as elk populations are highly influenced by harvest and variation in harvest pressure and hunting regulations (e.g., adult males–only harvest or either sex harvest [Law 2001; Compton 2007]). The Rocky Mountain elk is a highly sought-after game species, which accentuates the importance of understanding its metapopulation structure because of its economic benefits and recreational opportunities (Bryant and Maser 1982; Edge et al. 1986).
Previous genetic studies to determine interactions among elk populations have found minimal genetic differentiation between geographically separated populations (New Mexico [Hicks et al. 2007]; Ontario, Canada [Polziehn et al. 2000]; and South Dakota [Williams et al. 2002]). These patterns could indicate that elk populations over large areas (i.e., multiple landscapes) are interacting. Niedziałkowska et al. (2012) found low to moderate genetic differentiation among harvested red deer (Cervus elaphus) populations inhabiting forests in central Europe, but they also found less interaction among populations as distance increased between forests. These studies suggest that connectivity and metapopulation structure can be discerned through measures of relative similarity in genetic variation among populations.
Wildlife managers and conservation biologists are increasingly challenged to conserve and manage wildlife populations over large spatial extents consisting of multiple landscapes. Developing landscape-scale management recommendations requires understanding the spatial structure of populations (i.e., metapopulation structure) across the landscape, which can indicate levels of dispersal (i.e., gene flow) between populations (Kauffman et al. 2004; Davies et al. 2005; Epps et al. 2005; Garton et al. 2012; Hanski 2012). Without this understanding, a management plan could fail to notice management problems and their causes (Cooper and Mangel 1999).
Reed and Frankham (2003) showed maintaining genetic variation is important for a species' long-term survivability and indicates a species' potential to adapt to future environmental change. Because they are spread broadly over large landscapes, elk populations undergo variable harvest pressure that removes potential breeders from the population and may limit the number of successful dispersers, which in turn could limit gene flow and influence the distribution of genetic variation among populations (Law 2001; Peek et al. 2002; Stalling et al. 2002; Hard et al. 2006). Harvest levels of elk populations differ across management units and vary widely across their range (Mohler and Toweill 1982; Compton 2007; Aycrigg 2009; Rachael 2010). However, elk are capable of dispersing long distances through a variety of habitats, which may keep populations connected, but large areas of unsuitable habitat could limit their dispersal distance or reduce their survival during dispersal, thereby resulting in unoccupied habitat patches (Shoesmith 1980; Lyon and Christensen 2002; O'Gara 2002; Stalling et al. 2002). Therefore, we used a genetic approach to examine the distribution of genetic variation in elk populations across Idaho to describe metapopulation structure of a harvested species at 2 management extents (regions and elk management zones), and to evaluate whether population management through variation in harvest regulations influenced the metapopulation structure of a harvested species. Based on our knowledge of elk populations in Idaho in relation to variable harvest pressure, dispersal capabilities, and potential limitations to dispersal distance, we hypothesized a classical metapopulation would best describe the metapopulation structure exhibited by elk in Idaho except that patch sizes would not all be small as described by Levins (1970), but instead vary in size (Fig. 1; Compton 2007; Zager et al. 2007; Aycrigg 2009). We expected to observe a moderate level of genetic differentiation between regions and elk management zones, particularly those widely separated geographically (Table 1).
Materials and Methods
Historical populations and harvest levels
Historically, elk occurred throughout Idaho; however, they were most common in and around the mountainous areas of northern, central, and eastern Idaho and infrequently observed in the arid plains of southern Idaho (O'Gara and Dundas 2002). In the mid-1800s an influx of gold miners led to unregulated year-round hunting and by the early 1900s elk populations were few and far between (O'Gara and Dundas 2002). In response, a translocation program ran from 1915 to 1946 moving 675–800 elk from Yellowstone National Park (YNP) and Jackson Hole, Wyoming, to numerous locations throughout Idaho (Humbird 1975; Robbins et al. 1982; O'Gara and Dundas 2002). After these translocations, Idaho elk populations increased until about 1960. After 1960, they began to decrease. Subsequently, antlerless hunting (i.e., removal of females and males with no antlers present) was suspended in 1975 and populations increased again until the late 1980s (Zager et al. 2007). The general trend over the last 50 years has been toward increasing elk populations; however, pressure on elk populations in the form of human development and increased human access to elk habitat has increased at a greater rate (Compton 2007).
Harvest of elk populations increased steadily from about 2,000 elk annually in 1935 to approximately 16,000 in 1960 (Compton 2007). After 1960, harvest declined to a low of 4,000 in about 1976. After 1976, harvests levels climbed to approximately 28,000 over the next 20 years (Compton 2007). Over the last 10–15 years, harvest levels have ranged between 16,000 and 25,000 annually (Compton 2007; Rachael 2010).
Idaho Department of Fish and Game (IDFG) divided the state into 7 regions (approximately 15,400–43,500 km2) and 29 elk management zones (approximately 1,550–35,600 km2; Fig. 2). The 7 IDFG regions, all of which were included in our analysis, were Panhandle, Clearwater, Southwest, Salmon, Magic Valley, Southeast, and Upper Snake River. We chose these 2 management extents for our analysis because IDFG regions were the units within which our elk tissue samples were collected by IDFG wildlife managers, but elk management zones were the units within which elk population management objectives were applied (i.e., harvest levels).
Fig. 2.
Rocky Mountain elk (Cervus elaphus) tissue (muscle and ear punch) and DNA samples collected from harvested and live-captured individuals obtained during 1996–1998 and 2000–2004 from 7 Idaho Department of Fish and Game (IDFG) regions (left) and 9 elk management zones (right) in Idaho. Each circle represents a single sample. Regions are Clearwater (n = 52), Magic Valley (n = 15), Panhandle (n = 31), Salmon (n = 22), Southeast (n = 18), Southwest (n = 38), and Upper Snake River (n = 33). Elk management zones included in our analysis are Boise River (n = 15), Diamond Creek (n = 17), Lolo (n = 26), Palisades (n = 19), Palouse (n = 18), Panhandle (n = 36), Salmon (n = 19), Smoky Mountains (n = 15), and Weiser River (n = 16). Sample sizes included here differ from total sample sizes used in analysis (n = 216) because not all samples had specific location information and some locations were duplicated.
DNA sample collection and amplification
Elk muscle tissue samples (n = 216), via harvested and livetrapped individuals, were obtained during 1996–1998 (n = 33) and 2000–2004 (n = 183) from the 7 regions across Idaho (Fig. 2). IDFG wildlife managers and biologists collected muscle tissue samples from harvested elk brought to hunter check stations within each of the 7 regions in the state (n = 144). Furthermore, IDFG personnel and volunteers collected tissue samples from live-captured elk (n = 21 in 2002 and 2004). IDFG's Wildlife Health Lab also provided DNA samples from elk muscle tissue (n = 51). Each sample represented a unique individual. All samples were preserved by freezing or placing in ethanol until DNA extraction. Location information recorded with each tissue sample allowed us to assign 181 of the 216 tissue samples to an elk management zone. Samples were collected from each region and assigned to 9 of the 29 elk management zones, which covered most of the spatial extent occupied by elk in Idaho. We used both regions and elk management zones as our units of analysis for assessing the distribution of genetic variation among elk populations and describing the metapopulation structure of elk.
The DNA was extracted using DNAzol (Invitrogen Life Technologies, Grand Island, New York). Each extraction yielded approximately 25–150 μg of DNA suspended in 50 μl of deionized and distilled water (ddH2O). Fourteen polymorphic microsatellite loci (BL42, BM203, BM415, BM848, BM888, BM1009, BM4107, BM4208, BM5004, BM6506, BOVIRBP, FCB193, MAF109, and RM006 [Borst et al. 1989; Swarbrick and Crawford 1992; Buchanan and Crawford 1993; Kossarek et al. 1993; Bishop et al. 1994; Talbot et al. 1996) were amplified using polymerase chain reaction. These microsatellite primers, which were either bovine- or ovine-derived primers, were selected based on polymorphism (2–14 alleles/locus) and consistency of scoring. Amplification reactions, using polymerase chain reactions, totaled 10 μl and contained 2 μl of DNA, 0.2–0.8 μl of locus-specific primer (4 pmol/μl except BM6506 and RM006, which were 2 pmol/μl; MWG-Biotech, Hunstville, Alabama), 1 μl of 10X buffer (1.5 mM MgCl2; Promega, Madison, Wisconsin), 0.2–0.3 μl of 10 mM deoxynucleotide triphosphates (Promega), 0.2–0.4 μl Promega Taq (1.0 U Taq/reaction; Promega), and ddH2O.
Each amplification was denatured at 95°C for 5 min, then cycled 30–50 times, depending on the primer, at 94°C for 30 s, annealed at 50–60°C for 30 s (annealing temperature was optimized for each locus), and held at 72°C for 30 s. After cycling each sample was held at 72°C for 5 min. This amplification procedure was the same for all loci except BOVIRBP, in which a touchdown sequence was used. The 1st cycle started at 60°C and decreased by 0.5°C each cycle until 50°C. Multiplex polymerase chain reaction amplifications were BM1009, BM4107, and BM5004; RM006 and MAF109; and BM4208 and BM888. The remaining primers were amplified in singleplexes. Negative controls (i.e., ddH2O) were used to detect for polymerase chain reaction contamination and leakage between gel lanes during allele scoring.
No polymerase chain reaction products were diluted for electrophoresis except for those products with BM203 (diluted 1:3), BL42 (diluted 1:3), and BOVIRBP (diluted 1:2) primers. We combined 0.8 μl of polymerase chain reaction product with 1.6 μl of loading buffer, which consisted of 1.5 μl of formamide, 0.35 μl of ABI loading dye (blue dextran, 50 mg/ml; Applied Biosystems, Grand Island, New York), and 0.35 μl of GENESCAN-500 TAMRA (Applied Biosystems). Only 0.7 μl of this solution was loaded onto a single tooth of a Quick-Comb-96 (0.2 mm thick with 96 teeth; Sigma-Aldrich, Munich, Germany) for electrophoresis on a 6% acrylamide gel. To test for scoring consistency, 5–10 sample reactions were repeated on each 96-sample gel.
Polymerase chain reaction products were detected and visualized using an ABI 377 DNA automated sequencer (Applied Biosystems). GENESCAN 3.2.1 and GENOTYPER 2.5 software (Applied Biosystems) were used to size alleles. On each electrophoresis gel, 5–10 sample reactions were repeated and used to determine genotyping error rates. These repeated samples were done for all loci and made up about 5–10% of the total reactions. Genotyping error rates (i.e., errors per allele) were calculated by counting the number of errors per allele for each locus. A mean was calculated over all alleles for each locus and over all loci for an overall genotyping error rate. Errors were determined over multiple polymerase chain reaction runs (2–7 runs/locus) and multiple electrophoresis gels (2–7 gels/locus) when inconsistent scoring was apparent at the same allele and locus. The “true” allele(s) at each locus was defined as the allele most frequently or clearly observed over multiple polymerase chain reaction runs, was consistently scored, and its electropherogram had a scaled height ≥ 75. All loci were genotyped multiple times, but only 56 (86%) of 65 alleles were amplified multiple times. Loci success rates (i.e., number of successfully scored reactions over all reactions) also were calculated. Three microsatellite loci (BM848, MAF109, and BOVIRBP) with low success rates (< 80%) were dropped from any further analysis. Any individuals with at least 8 of the remaining 11 loci genotyped were included in our analyses (percent missing genotypes per locus ranged from 0.9% to 11.6%). In addition to calculating overall genotyping error rates, we randomly selected 27 samples from each locus to calculate allelic dropout and false allele error rates. Allelic dropout was calculated for each locus as the number of dropouts per total number of heterozygous repeats and total allelic dropout was averaged across all loci. For false allele error rates for each locus, we divided the number of false alleles detected by the total number of polymerase chain reactions and total false allele error rate was averaged across all loci. Both calculations provided unbiased estimates of genotyping error frequencies (Broquet and Petit 2004).
Genetic analysis
All genetic analyses were conducted over 216 individual samples for regions and 181 samples for elk management zones. However, our isolation by distance analysis among regions only used 209 samples, because some samples had duplicate locations and all duplicates were removed for this analysis.
Tests for Hardy–Weinberg equilibrium per locus and pairwise locus tests for linkage disequilibrium were conducted for all regions and elk management zones using GENEPOP 4.0.7 (Rousset 2008). A sequential Bonferroni correction for multiple tests was applied (Rice 1989). A test of the Wahlund effect over all samples was performed using GENEPOP 4.0.7 (Rousset 2008). We used an exact G-test to test for equal allelic distributions between the 2 sampling periods (i.e., 1996–1998 and 2000–2004). To evaluate distribution of genetic variation among regions and elk management zones we used 5 a priori approaches, which included both direct and indirect measures. First, an exact G-test was calculated using GENEPOP 4.0.7 (Rousset 2008). We tested the null hypothesis of equal allelic distributions among the regions and elk management zones. All region and elk management zone pairs and all loci were included in the exact G-tests at alpha of 0.05.
Second, average expected heterozygosity (HE) and observed heterozygosity (HO), number of private alleles, and mean allelic richness between all regions and elk management zones were used as indirect measures of the distribution of genetic variation. The program GIMLET 1.3.3 (Valiere 2002) was used to calculate HE and HO. HE was reported in addition to HO because it is less sensitive to sample size than HO (Frankham et al. 2002). We also calculated standard errors for HE-values and standard deviations for HO-values. Number of private alleles and allelic richness were calculated using FSTAT 2.9.3.2 (Goudet 2002; updated from Goudet 1995), based on minimum sample sizes for each analysis. We also calculated mean number of alleles per locus and allelic richness. Allelic richness is a standardized value, not influenced by sample size, and can be compared across different sample sizes (Allendorf and Luikart 2007).
Third, pairwise FST-values (i.e., a measure of the total genetic variation partitioned among regions or elk management zones), calculated as theta values based on Weir and Cockerham (1984), were used as a direct measure to examine levels of gene flow and genetic differentiation between regions and elk management zones (GENEPOP 4.0.7—Rousset 2008). According to Hartl and Clark (1997), guidelines for the interpretation of FST-values for gauging genetic differentiation are: 0–0.05 little; 0.05–0.15 moderate; 0.15–0.25 great; and > 0.25 very great. High genetic differentiation indicates restricted or low gene flow between regions or elk management zones. Using FSTAT 2.9.3.2, we tested whether pairwise FST-values among regions (n = 216) and elk management zones (n = 181) were significantly different from zero (Goudet 2002; updated from Goudet 1995). For each pairwise FST-value, 95% confidence intervals were calculated based on bootstrapping. Overall FST-values and their 95% confidence intervals were calculated using jackknifing and bootstrapping.
Fourth, individuals were assigned to the most likely population, based on their genotypes, using the direct measure of assignment tests developed by Paetkau et al. (1995). This test determined from which region or elk management zone each individual most likely originated (Paetkau et al. 1995, 1997). Population assignments were made using the Doh assignment test calculator (Brzustowski 2002), based on Paetkau et al. (1995, 1997). Furthermore, a resampling algorithm was used to compute a probability estimate that an individual originated from its region or elk management zone of origin (Paetkau et al. 2004). Lower probability estimates indicated fewer migrants within a region or elk management zone, whereas higher probability estimates indicated more than a few migrants.
Last, we used a Bayesian model–based clustering of genotypes to infer population structure (Pritchard et al. 2000). Individuals were assigned to a region or elk management zone based on their genotypes and population allele frequencies were calculated. The number of populations (K) was chosen a priori to range from 1 to 7 regions or 1 to 9 elk management zones with the software program structure 2.2 (Pritchard et al. 2000) attempting to cluster the genotypes over 10 runs for each value of K. We averaged results over the 10 runs to obtain mean values. For each of the 10 runs, we used an admixture model with correlated allele frequencies, along with a burn-in of 10,000 iterations and a run length of 100,000 iterations. Posterior probabilities of K closer to zero indicated a better fit of the data.
We tested for isolation by distance with a simple Mantel test with 10,000 randomizations using zt 1.1 (Bonnet and Van de Peer 2002) between regions (n = 209) and elk management zones (n = 181) for all individuals. Pairwise FST-values calculated as theta values based on Weir and Cockerham (1984) in GENEPOP 4.0.7 (Rousset 2008) were used as genetic distances. Minimum and mean distances in meters between samples from region pairs and elk management zone pairs were calculated using ArcGIS 9.2 (Environmental Systems Research Institute 2006) and used as geographic distances.
In addition to evaluating the distribution of genetic variation of elk, we estimated recent migration rates (i.e., within the last 1–3 generations) between regions (n = 216) and elk management zones (n = 181) using a Bayesian multilocus genotyping method (Wilson and Rannala 2003). We used BayesAss 3.0.3 (Wilson and Rannala 2003; Rannala 2007) to estimate the mean posterior distribution of migration rates for all region and elk management zone pairs. Both the proportion of residents and the proportion of immigrants into each region or elk management zone were estimated. Using our data on elk genotypes, we initially ran the Markov chain Monte Carlo analysis 5–7 times to optimize the input values for mixing parameters of inbreeding coefficients and allele frequencies that produce acceptance rates between 20% and 60% (Rannala 2007). When acceptance rates fell into this range, we ran the Markov chain Monte Carlo analysis for 3 × 107 iterations, sampling every 2,000 iterations after 1 × 107 burn-in iterations to look for convergence of the log-posterior probabilities over the last 2 × 107 iterations. After convergence was observed in the Markov chain Monte Carlo analysis, we completed 6 independent runs of BayesAss and averaged results for migration rates and standard deviations.
Results
DNA sample collection and amplification
A total of 216 samples of elk muscle tissue (females = 107, males = 101, unknown = 8) were used in our analyses. All 216 samples were assigned to an IDFG region, whereas 181 samples (females = 88, males = 88, unknown = 5) were assigned to 1 of 9 elk management zones. Based on simulations to detect genetic differentiation by Ryman et al. (2006), we believe our sample sizes and level of polymorphism at each locus were sufficient to detect genetic differentiation between IDFG regions and elk management zones. Our genotyping error rate was 0.1 errors/allele (range = 0.01–0.22 errors/allele) and our loci success rate was 93.3% after 3 loci with low success rates (< 80%) were removed from our analysis. Our overall allelic dropout rate for the remaining 11 loci after removing the loci with low success rates was 0.13 (range = 0.02–0.3) and our false allele error rate was 0.07 (range = 0–0.25).
Genetic analysis
No deviations from Hardy–Weinberg equilibrium were found among regions (n = 216) or elk management zones (n = 181) for all individuals sampled after applying a sequential Bonferroni correction. No deviations from linkage equilibria were detected at any locus in any region or elk management zone for all individuals after a Bonferroni correction was applied. Because no test deviated from linkage equilibrium, we assumed data from different loci were independent. Based on test results of the Wahlund effect, 2 loci suggested subpopulation structure over all individuals (n = 216). The exact G-test between the 2 sampling regions indicated no significant difference in allele distributions (P < 0.05).
Our 5 approaches for evaluating directly and indirectly the distribution of genetic variation indicated genetic variation was more or less evenly distributed between populations across our 2 management extents (regions and elk management zones). However, across regions and elk management zones a few patterns were evident. Results for the exact G-test (i.e., our 1st approach to evaluate distribution of genetic variation) for allele distributions showed 10 of 21 region pairs (P-value range: 0–0.043) and 18 of 36 elk management zone pairs (P-value range: 0–0.048) had significantly different allele distributions (P < 0.05). Among the regions, Clearwater and Southwest occurred most frequently among the 10 pairs that had significantly different allele distributions. Among the elk management zones, Boise River and Lolo occurred most frequently among the 18 pairs that had significantly different allele distributions.
For our 2nd approach, regions and elk management zones had a range of 15–54 individuals in our analysis with 3.5–4.8 mean number of alleles per locus, 0.50–0.61 HO-values, 0.54–0.61 HE-values, 0–3 private alleles, and 3.3–3.9 mean allelic richness values (Table 2). Among regions and elk management zones for all individuals our values were comparable, suggesting minimal differences in genetic diversity and therefore genetic variation between regions and elk management zones.
Table 2.
Total number of alleles, mean observed and expected heterozygosity (HO and HE, respectively), mean number of alleles per locus, mean allelic richness, number of private alleles, and number of individual Rocky Mountain elk (Cervus elaphus) in 7 regions and 9 elk management zones in Idaho. Results are from tissue and DNA samples collected during 1996–1998 and 2000–2004 and are based on 8–11 polymorphic microsatellite loci. HO- and HE-values were calculated using GIMLET 1.3.3 (Valiere 2002) and mean allelic richness was calculated using FSTAT 2.9.3.2 (Goudet 2002; updated from Goudet 1995).
Region pairs with significantly different allele distributions from the exact G-test had pairwise FST-values < 0.05 (i.e., our 3rd approach), which suggested little genetic differentiation; however, many of them were statistically significant from zero (Table 3). The pairwise FST-value between Boise River and Lolo elk management zones was 0.0872, indicating a moderate amount of genetic differentiation and was statistically significant from zero. Within regions, most pairwise FST-values that were significantly different from 0 included the Clearwater or Southwest region and within elk management zones included Lolo or Boise River, which were located within the Clearwater and Southwest regions (Table 3; Fig. 2). After applying a sequential Bonferroni correction, 4 and 1 pairwise FST-values were significantly different at P < 0.05 for regions and elk management zones, respectively (Table 3). Panhandle and Salmon elk management zones had the single pairwise FST-value that was significant, whereas Panhandle and Salmon regions were 1 of the 4 regions with significantly different pairwise FST-values. Most of the Salmon region tissue samples for genetic analysis were obtained within the Salmon elk management zone (Fig. 2).
Table 3.
Pairwise FST-values for all individual Rocky Mountain elk (Cervus elaphus) within 7 regions and 9 elk management zones in Idaho. FST-values were calculated using GENEPOP 4.0.7 (Rousset 2008). Overall FST-values are based on jackknifing over 11 loci. Values in parentheses are 95% confidence intervals based on bootstrapping. Both bootstrapping and jackknifing were done using FSTAT 2.9.3.2 (Goudet 2002; updated from Goudet 1995). Sequential Bonferroni corrections were calculated at P < 0.05. Results are from tissue and DNA samples collected during 1996–1998 and 2000–2004 and based on 8–11 polymorphic microsatellite loci.
Table 3.
Continued.
Our 4th approach based on assignment test results further supported minimal differences in genetic variation between regions and elk management zones. Overall only 29.6% and 23.8% of all individuals were assigned to the region or elk management zone of origin (Table 4). Assignment test results showed 3 of 7 regions (Clearwater, Southwest, and Upper Snake River) and 5 of 9 elk management zones (Boise River, Diamond Creek, Lolo, Palisades, and Panhandle) having most of their individuals assigned back to the region or elk management zone of origin (Table 4). Again, the Clearwater and Southwest regions and the Lolo and Boise River elk management zones had most of their individuals assigned back to them. The probability estimates based on randomizations supported the pattern observed within our assignment test results (Table 4). Fewer migrants (i.e., lower probability estimates) were indicated within the Clearwater and Southwest regions and the Lolo and Boise River elk management zones.
Table 4.
Assignment of individual Rocky Mountain elk (Cervus elaphus) to 7 regions and 9 elk management zones within Idaho. Values are number of all individuals from each region or elk management zone in rows assigned to each region or elk management zone in columns. Individual assignments were based on genotypes using 8–11 polymorphic loci. The value in parentheses is a probability estimate of observing at least as many cross assignments as would be seen if the null hypothesis of elk in regions or elk management zones being one well-mixed population in Hardy–Weinberg equilibrium. Lower probability estimates indicate fewer migrants between that pair of regions or elk management zones. Results are from tissue and DNA samples collected during 1996–1998 and 2000–2004. Assignments and probability estimates were calculated using the Doh assignment test calculator (Brzustowski 2002).
The mean posterior probability for the number of populations using the Bayesian model-based clustering of genotypes (i.e., our 5th approach) showed the strongest support for K = 1 (Pritchard et al. 2000; Table 5). Our results showed the proportion of individuals assigned to each value of K to be approximately equal (∼ 1/K), which indicated genetic variation was evenly distributed among our sampled individuals (Pritchard et al. 2007). This indicates low genetic differentiation among all groups of individuals across the 7 regions and 9 elk management zones.
Table 5.
Posterior probability and associated variance of number of populations (K) from Bayesian model–based clustering of Rocky Mountain elk (Cervus elaphus) genotypes from 7 regions and 9 elk management zones in Idaho using structure (Pritchard et al. 2000). Values of K were chosen a priori and the posterior probability of K closest to 0 indicates a best fit of the data. Results are from tissue and DNA samples collected during 1996–1998 and 2000–2004 and based on 8–11 polymorphic microsatellite loci.
Our isolation-by-distance tests showed no significant correlation between genetic (pairwise FST-values) and minimum and mean geographic distances among samples within regions (P = 0.36 for minimum distances and P = 0.27 for mean distances) and elk management zones (P = 0.41 for minimum distances and P = 0.42 for mean distances).
Our mean posterior distribution of migration rates using the Bayesian multilocus genotyping method showed that most individuals within each region and elk management zone were residents (Table 6). However, the Panhandle region and the Panhandle and Salmon elk management zones had mean migration rates ≥ 0.10. Additionally, the Clearwater region had mean migration rates > 0.05. These values suggest directional migration from these regions and elk management zones.
Table 6.
Means of the posterior distribution of migration rates and associated standard deviations for Rocky Mountain elk (Cervus elaphus) within 7 regions and 9 elk management zones in Idaho. Migration rates were calculated using BayesAss 3.0.3 (Wilson and Rannala 2003; Rannala 2007). Region or elk management zone from which each individual was sampled is shown in rows, whereas the region or elk management zone from which an individual migrated is shown in columns. Values along the diagonal are the proportion of individuals from the source region or elk management zone. Migration rates ≥ 0.10 are in boldface type and standard deviations ≥ 0.05 are italicized. All values are means over 6 independent runs of BayesAss using a Markov chain Monte Carlo (MCMC) analysis. Results are based on tissue and DNA samples collected from 1996 to 1998 and 2000 to 2004 and are based on 8–11 polymorphic microsatellite loci.
Table 6.
Continued.
Discussion
Our goal was to describe the metapopulation structure of elk populations at 2 management extents (regions and elk management zones) in Idaho using a genetic approach. We hypothesized that variable harvest pressure and the potential for unsuitable habitat to limit dispersal would support a classical metapopulation structure (Fig. 1) with variable patch sizes as the metapopulation structure of elk in Idaho. However, our results support a patchy metapopulation structure (Fig. 1; Table 1) because sufficient interchange of individuals is occurring that only slight differences in genetic variation among regions and elk management zones were detected. The high level of patch connectivity suggested by the patchy metapopulation structure and our results indicates local elk populations have a low risk of losing genetic variation (Harrison 1991; Stith et al. 1996; Harrison and Taylor 1997; Mills 2007). Therefore, elk populations in Idaho, from a genetic perspective, should be resilient to varying harvest levels over time because individual movement maintains genetic variation across the entire metapopulation.
Our number of individuals, HO, HE, mean number of alleles per locus, and mean allelic richness fall in the range of values obtained in other studies of elk populations (number of individuals = 10–56, HO = 0.41–0.61, HE = 0.26–0.60, mean number of alleles per locus = 1.9–5.0, mean allelic richness = 2.9–4.0 [compare to Table 2; Polziehn et al. 2000; Williams et al. 2002; Hicks et al. 2007; Hundertmark and Van Daele 2010]). Elk have similar genetic variation between different populations that are widely separated geographically from South Dakota (Williams et al. 2002) to New Mexico (Hicks et al. 2007) to Ontario, Canada (Polziehn et al. 2000) to Washington and Alaska (Hundertmark and Van Daele 2010). These results suggest that our sampled individuals are representative of genetic variation in other elk populations.
The patchy metapopulation structure supported by our results may not be restricted to Idaho. The elk we sampled could be part of a geographically larger metapopulation because our results are consistent with those of Polziehn et al. (2000), who found minimal differences in genetic variation among 6 elk populations in the Canadian Rocky Mountains and elk in YNP. Hicks et al. (2007) also found minimal genetic variation among reintroduced elk populations (Arizona, Oklahoma, Oregon, New Mexico, and North Dakota) and their founder population in YNP. Comparisons between elk populations in YNP and Custer State Park, South Dakota, which were reintroduced from YNP, indicated little difference in genetic variation (pairwise FST = 0.030), even though the reintroduction was > 80 years ago (Williams et al. 2002). Polziehn et al. (2000) proposed that elk in the Rocky Mountains are part of a continuum with individuals moving through corridors connecting parks from YNP into Canada. This continuum could represent a patchy metapopulation structure stretching at least from YNP through Idaho into western Canada. Even though not all, or even most, elk populations have been sampled over this extensive area, our results combined with those of Hicks et al. (2007), Polziehn et al. (2000), and Williams et al. (2002) do provide support for this idea.
Historical reintroductions of elk from 1915 to 1918 when 200 elk were moved from YNP to many areas in Idaho and during 1935–1939 when 416 elk were brought to Idaho may have influenced the distribution of genetic variation among regions and elk management zones (Robbins et al. 1982). Potentially, the influx of several hundred elk from YNP may have diluted genetic variation among elk populations present at that time (Leberg and Ellsworth 1999). Alternatively, the translocated elk may have increased the genetic variation within the existing elk populations, but the genotypes of the translocated elk were relatively uniform because YNP source populations were genetically similar (Hicks et al. 2007). Both these hypotheses could explain our results.
Even though > 60 years have passed since the largest elk translocation and the collection time of our genetic samples, the influence of genetic drift, which can lead to differences in the distribution of genetic variation among populations, has not been observed among the sampled regions or elk management zones. Historical elk populations in Idaho were described as small and scattered (O'Gara and Dundas 2002), but our results do not indicate founder effects (i.e., loss of genetic variability because of small founding group size) occurred, which could be attributed to high gene flow posttranslocation or high population growth rates soon after elk were translocated. We do not have the population data needed to tease apart these 2 hypotheses. From our findings and those of Hicks et al. (2007), gene flow seems to have played a major role in maintaining genetic variation among elk populations of different sizes moved into a wide variety of places, thereby reducing the influence of founder effect. Also, Hundertmark and Van Daele (2010) found that high population growth rates in an introduced population of elk (Cervus elaphus roosevelti) minimized heterozygosity loss and thereby minimized the founder effect.
DeYoung et al. (2003) and Leberg and Ellsworth (1999) found that the genetic contribution of translocated white-tailed deer (Odocoileus virginianus) in the southeastern United States was present in recipient populations for many years and explained, in part, the observed genetic similarity among different white-tailed deer populations. Additionally, Hicks et al. (2007) hypothesized insufficient time had passed to see differences in the distribution of genetic variation among 5 elk populations reintroduced from YNP 21–93 years prior to analyzing their genetic variation. Our results support this hypothesis because differences in the distribution of genetic variation among our regions and elk management zones were not distinct. Even the 2 main source populations of elk from YNP showed no genetic differentiation (Hicks et al. 2007), further influencing the distribution of genetic variation among elk in Idaho. We do not know the distribution of genetic variation nor the population demographics among elk populations before any reintroductions, but our findings suggest that elk populations in the early 1900s already exhibited minimal differences in genetic variation. This suggests a patchy metapopulation structure may have been maintained by elk throughout Idaho over many years.
For a harvested species, differences in harvest regulations across large landscapes can influence the distribution of genetic variation; for example, the Southwest and Clearwater regions and the Boise River and Lolo elk management zones exhibited subtle differences in the distribution of genetic variation based on assignment test results and pairwise FST-values (Tables 3 and 4). Interchange of individuals probably still occurred within these regions and elk management zones, but at a comparably lower rate. The elk population in the Boise River elk management zone within the Southwest region is known to have been small for many years. In 2008, there was a 1:1 ratio of hunters to elk, which was relatively high since the population was low in the 1990s, in part, because male harvest exceeded male-calf recruitment (Compton 2007; Rachael 2010). The Lolo elk management zone within the Clearwater region has elk populations below its stated management objectives (50–77% fewer females and 25–54% fewer males—Compton 2007; Rachael 2010). This is attributed to low calf productivity and recruitment since the late 1980s (Compton 2007; Zager et al. 2007). Genetic characteristics of a population can be altered by harvest when harvest regimes bias adult sex ratios heavily toward females and reduce adult life span (Hartl et al. 1991; Harris et al. 2002; Hard et al. 2006). Both effects are likely observed in these 2 elk management zones. In the Boise River elk management zone because adult male harvest exceeded male-calf recruitment, the population has been low in the past, and the ratio of hunters to elk is high suggests reduced adult life span and a female-skewed adult population (Rachael 2010). In the Lolo elk management zone, harvest has been skewed toward males since 1976 and low calf productivity and recruitment suggest a female-skewed adult population and reduced adult life span (Zager et al. 2007). However, gene flow within a patchy metapopulation can counteract these effects to some extent, particularly when mature adult males are contributing to gene flow (Hard et al. 2006). To further counteract these effects, managers could reduce harvest levels through permit harvests (i.e., reduce harvest of adult males—Hard et al. 2006; Allendorf et al. 2008; Coltman 2008).
Both metapopulation structure and source–sink population dynamics focus on the role of dispersal among populations (Kawecki 2004). Source–sink population dynamics might occur within a patchy metapopulation structure whereby some populations are net “contributors” to the metapopulation, whereas other populations act as sinks because they do not contribute to the metapopulation (Runge et al. 2006). Our assignment test results indicated directional movement of elk between some regions and elk management zones (e.g., from the Clearwater region to Magic Valley, Salmon, and Southwest regions, and from the Panhandle elk management zone to Lolo and Palouse elk management zones [Table 4]). Furthermore, the migration rates among regions suggest that the Panhandle region may act as a source to other regions (Table 6). Among elk management zones, both Panhandle and Salmon could be considered sources. Even though source–sink population dynamics are suggested based on our genetic analysis, Aycrigg (2009) estimated the intrinsic population growth rate to be stable over 20 years within elk management zones, thereby suggesting, at least within elk management zones, there are no source or sink populations. Without knowing the recruitment and emigration (i.e., successful dispersal) contribution of each population to the patchy metapopulation, variability in emigration could exist between populations that we observed in our genetic analysis (i.e., assignment tests and migration rates [Tables 4 and 6]). Andreasen et al. (2012) found source–sink population dynamics occurring among harvested cougar (Puma concolor) populations across Nevada and California. Therefore, because elk are harvested throughout Idaho under varying harvest regulations, variability in emigration between populations is likely. Our results regarding source–sink population dynamics also could be influenced by not having sampled all populations exchanging migrants, an assumption of the Bayesian multilocus genotyping method we used (Wilson and Rannala 2003).
Most of our tissue samples, collected over several years, were from the fall harvest of elk populations. Our samples may be inherently biased toward describing the distribution of genetic structure of elk populations in the fall (Latch and Rhodes 2006). Hicks et al. (2007) collected genetic samples in a similar ad hoc approach and also found minimal differences in distribution of genetic variation between 5 elk populations. Ideally, samples should be collected randomly with regard to elk population ecology and life history; however, this rarely occurs with harvested species because more samples can be collected opportunistically during the hunting season (Waples 1998; Latch and Rhodes 2006; Schwartz and McKelvey 2009). Schwartz and McKelvey (2009) found sampling design can influence how well clustering algorithms, such as structure (Pritchard et al. 2000), perform. Furthermore, pairwise FST-values and assignment test results are influenced by sample sizes and number of loci; however, the patterns of genetic differentiation and proportion of individuals assigned to each population provide useful insights into the distributional patterns of genetic structure (Manel et al. 2005). To increase the power to detect genetic differentiation in elk populations, samples collected during different time periods, larger samples sizes, and more loci should be taken into account (Manel et al. 2005; see also Ryman et al. 2006; Schwartz and McKelvey 2009).
We attempted to minimize genotyping errors (i.e., the differences observed between ≥ 2 genotypes obtained from the same tissue sample) by using blind samples, repeating amplifications for a portion of our samples, and using negative controls to monitor contamination (Bonin et al. 2004; Kelly et al. 2011). Despite these efforts, our genotyping error rate was high (0.1 errors/reaction), but we observed no deviations from Hardy–Weinberg equilibrium, which detects homozygous excess from allelic dropouts or null alleles (Hoffman and Amos 2005). Kelly et al. (2011) found that genotyping errors led to populations appearing to be more genetically distinct than they actually were. Our results from both direct and indirect methods indicate lows levels of genetic differentiation (Tables 2–4). However, the subtle genetic differentiation we detected in the Southwest and Clearwater regions and the Boise River and Lolo elk management zones may not be substantiated with our high genotyping error rate. We acknowledge our genotyping error could have been reduced by more rigorous sampling and analytic efforts, but we believe our results regarding genetic variation of elk in Idaho are applicable to their conservation and management.
Many species, including those with more continuous spatial distributions and high mobility, have a metapopulation structure (Hanski and Gaggiotti 2004). Here we focused on the level of genetic similarity among elk populations in Idaho to describe their metapopulation structure. Our results suggest a network (i.e., patchy metapopulation structure) of elk populations linked by extensive gene flow and exhibiting little to no differences in genetic variation, which is common among cervids (Coltman 2008). The lack of geographic barriers and the dispersal capabilities of elk support our results based on genetics (Shoesmith 1980; O'Gara 2002). Knowledge of a harvested species' genetic structure is a good indicator of potential genetic consequences of harvest and whether the harvest is sustainable through time. Our findings indicated that current management with variable harvest regulations of elk populations in Idaho is maintaining sufficient gene flow between populations. Using our genetic analysis results as a baseline, future elk population management could include periodic genetic sampling and analysis to ensure a sustainable harvest. Without information on the amount of genetic connectivity within a species' metapopulation, wildlife managers could unknowingly manage populations as isolated from each other and potentially lower the genetic connectivity between populations (Cooper and Mangel 1999). However, connectivity observed through genetic approaches does not always directly translate to demographic connectivity because genetic connectivity relies on dispersing individuals to successfully contribute to gene flow (Palsbøll et al. 2006; Lowe and Allendorf 2010). If demographic connectivity is low (i.e., few dispersers), there is a risk of local extinction within a metapopulation even with sufficient gene flow. Therefore, combining genetic results with movement patterns and demographic rates of a harvested species will strengthen the understanding of its distribution of genetic variation and improve its population management and conservation (see Lande 1988; Tallmon et al. 2002; Aycrigg 2009; Lowe and Allendorf 2010).
Acknowledgments
Financial support was provided by IDFG, University of Idaho, Rocky Mountain Elk Foundation, Boone and Crockett Club, Idaho Cooperative Fish and Wildlife Research Unit, Charles DeVlieg Foundation, and the United States Geological Survey Gap Analysis Program. We owe our gratitude to numerous people for their assistance, effort, advice, expertise, and support, including the Regional Wildlife Managers for IDFG: J. Hayden, J. Rachael, J. Rohlman, J. Crenshaw, C. Anderson, R. Smith, B. Compton, and T. Keegan. We also thank M. Lucia, D. Spicer, B. Palmer, G. Vecellio, J. Jusa, K. Warner, K. Laudon, D. Parrish, J. Naderman, J. Rydlach, W. Seybold, L. S. Mills, P. Zager, J. M. Scott, C. Anderson, L. Waits, J. Rachlow, J. Peek, A. Metge, C. Miller, J. Horne, M. DeBarba, J. Adams, M. Murphy, J. Hicks, B. Rankin, J. Keehner, D. Bates, M. Moses, R. Dobson, D. Roon, M. Scott, E. Strand, B. Gilbert, B. Dennis, K. Steinhorst, C. Williams, K. Strickler, L. Vierling, L. Svancara, and the Waits Lab Genetics group. This manuscript has been improved by several anonymous reviewers.
