Sex-Biased Gene Flow Among Elk in the Greater Yellowstone Ecosystem Free

The Greater Yellowstone Ecosystem (GYE) supports world-renowned populations of elk (Cervus elaphus) that provide significant visitor enjoyment and benefits to local economies through guiding, hunting, and ecotourism. Elk are the most numerous large mammal in the GYE (N ∼ 50,000) and have strong effects on other species including predators and scavengers. Elk influence ecosystem characteristics and processes such as soil fertility, and vegetation production and diversity (Toweill et al. 2002). Elk maternal gene flow (where gene flow is defined as the exchange of alleles between populations) is important to understand because female gene flow and physical movement influence genetic and demographic rescue, colonization rates, demographic vital rates, and the spread of certain diseases (Thorne et al. 1979; Martin et al. 2000; Tallmon et al. 2004).

The GYE is one of a few areas where elk were not extirpated in North America by the early 1900s due to overharvest, competition with livestock, and perhaps disease (Houston 1974). Elk have not been translocated into or within the GYE. These observations make the GYE among the best (and few remaining) locations to study natural population genetic structure and patterns of gene flow in elk (Houston 1974; Boyce and Hayden-Wing 1979; Polziehn and Strobeck 1998). The GYE stretches approximately 400 km north–south and 300 km east–west, spanning portions of Idaho, Montana, and Wyoming with elevation ranges from 1,200 to 4,200 m for our study area (Figure 1).

Mitochondrial (mt) DNA is a useful marker for resolving maternal gene flow because it is a maternally inherited haploid marker with a single chromosome coming only from the mother whereas nuclear microsatellite DNA markers are biparentally inherited (with a single allele from each parent; Allendorf et al. 2013). Recent work from Hedrick et al. (2013) presents an equation to estimate the ratio of male to female gene flow from genetic distance estimates (F_ST calculated from mtDNA or microsatellites) using Wright's (1951) island model of gene flow. The equation is useful for studies in which markers for global gene flow for both sexes (e.g., microsatellites from nuclear DNA) and maternal gene flow (e.g., mitochondrial DNA) are available. Genetic distance estimates also allow for tests for patterns of genetic structure such as isolation by distance (e.g., genetic distance and geographic distance between populations are related such that genetic distance increases monotonically with geographic distance; Wright 1943) that can help managers understand patterns of gene flow.

We assessed maternal and overall gene flow patterns in elk by estimating the ratio of male to female gene flow using F_ST values calculated using mtDNA from eight populations and microsatellite markers from a subset (three of the eight) of populations in the GYE. We then assessed the generality of our findings by comparing our ratios of gene flow in elk to other large mammals using published studies reporting the F_ST for both maternal and overall genetic differentiation. Finally, we tested the hypothesis of isolation by distance in elk maternal genetic structure. From our results we suggest new and promising directions for future research in landscape genetics and sex-specific gene flow in elk and other species.

We collected samples of blood, tissue, or fecal pellets from 380 elk (223 females, 19 males, 138 unknown) in eight populations in the GYE (Table 2; Figure 1). We sampled more females than males to precisely quantify maternal gene flow to help predict spread of maternally transmitted diseases, and because we expected males to show little structure and thus be less informative for gene flow studies (Hicks et al. 2007). In addition, we preferentially captured females captured for radio-collar studies. We collected samples (n ≈ 20) during multiple years (2005–2008) from four populations (Paradise Valley, Madison Valley, Northern Range, and Muddy Creek) to test for temporal stability of haplotype frequencies (Table 1).

We defined populations as large groups or collections of individuals from a location where elk congregate, such as winter ranges with hundreds to thousands of elk. All collection of samples occurred within a 4-yr time frame (one generation) to reduce potential intergenerational effects (e.g., on spatial F_ST estimates). We collected blood or tissue from captured (Northern Range, Paradise Valley, and Madison–Firehole) or hunter-killed elk (Madison Valley, Shoshone River, Muddy Creek). We collected fecal pellets within 1–2 h after defecation in Sand Creek and Bench Corral and from 5 of 62 individuals from Muddy Creek (Figure 1). To prevent repeated sampling of the same individual, we only collected feces from individuals observed defecating, small groups (5–10) of individuals that were at least 0.5–1 km apart, or from individuals with distinctive natural marking or unique ear tags, or radio collars.

We isolated DNA from tissue and blood using the Qiagen QIAamp isolation kit (Chatsworth, California) and from feces using the QIAamp blood kit as described in Maudet et al. (2004). We conducted polymerase chain reaction (PCR) amplification on a 596–base pair fragment of the mtDNA control region using primers 275-294F (5′-CTCGTAGTACATAAAATCAA-3′) and 990-968R (5′-ATAAGGGGGAAAAATAAGAA-3′) and reaction conditions given in Polziehn and Strobeck (1998). The PCR and sequencing were conducted by the University of Washington High-Throughput Genomics Center, Department of Genome Sciences (Seattle, Washington). Each 100-mL PCR was performed on a 9600 Perkin-Elmer Cetus Thermocycler using the following conditions: a 3-min denaturing step at 94°C; 30 cycles at 94°C for 15 s, 56°C for 30 s, and 72°C for 30 s.

We conducted forward and reverse strand sequencing on all samples to ensure data quality. We performed each sequencing reaction using approximately 8 µL of PCR product, as described in the Perkin-Elmer Dye Terminator Cycle Sequencing Ready Reaction kit. Cycle sequencing reaction parameters on a 9600 Perkin Elmer Cetus Thermocycler were as follows: denaturation at 96°C for 15 s, annealing at 50°C to for 1 s, and extension at 60°C for 4 min. We separated sequencing reactions by electrophoresis on an ABI Prism 377 Perkin Elmer automated sequencer. The sequin file for the 380 mtDNA sequences can be found in GenBank and Dryad (accession numbers: JX125702 - JX126108; see Archived Material).

We corrected and aligned sequences of mtDNA using DNAstar 5.0 software package (DNASTAR Inc., Madison, Wisconsin). We double-checked sequences visually for quality and correctness, including every polymorphic site. We randomly re-extracted and resequenced 5% of all samples to monitor for potential errors; we found none. Amplification and sequencing success was relatively high even for fecal samples, where ∼85% of samples yielded useable sequences, considering the relatively long mtDNA fragment we amplified (596 base pairs).

We genotyped eight microsatellite loci for three populations: Muddy Creek, Northern Range, and Paradise Valley (Figure 1; see Archived Material), using the same individuals for both microsatellites and mtDNA sequencing from the Muddy Creek and Northern Range populations. The Paradise Valley samples for mtDNA were from hunter-killed elk, whereas microsatellite DNA samples were from live captured elk in the same geographical location. The captured elk provided better quality DNA from fresh blood, which typically yields more reliable microsatellite genotypes than hunter-killed samples and which we obtained after the mtDNA sequences. We used the following microsatellite DNA loci: BM5004, BM888, BM1009, BM4208, FCB193, OarkP6, RM006, BM415 (Buchanan and Crawford 1992; Kossarek et al. 1993; Bishop et al. 1994; Paterson and Crawford 2000).

All microsatellite PCRs consisted of an initial denaturation at 96°C for 15 s, annealing at 50 to 56°C (Table 1) for 1 s, and extension at 60°C for 4 min. We amplified the BM5004, BM888, and BM1009 loci together in one PCR; the BM4208 and FCB193 loci together in another (separate) PCR; and OarkP6 and RM006 in a third PCR. We amplified the BM415 locus alone at 50°C (Table 1). We extracted whole genomic DNA from elk tissue and blood samples using the QIAGEN Dneasy Tissue Kit (Qiagen) according to manufacturer's instructions. The reaction volume (10 µl) contained 1.0 µL DNA, 1× reaction buffer (Applied Biosystems, Foster City, California), 2.0 mM MgCl₂, 200 µM of each dNTP, 1 µM reverse primer, 1 µM dye-labeled forward primer, 1.5 mg/ml bovine serum albumin, and 1 U Taq polymerase (Applied Biosystems). We visualized the resultant products (PCR profiles) on a LI-COR DNA analyzer (LI-COR Biotechnology, Lincoln, Nebraska) and called genotypes visually (manually). We independently visualized and double-checked the profiles and genotypes for quality and correctness. We randomly re-extracted and regenotyped 5% of all samples to test for potential errors and found none.

We computed diversity indices (haplotype diversity and nucleotide diversity for each population), number of variable nucleotide sites, average number of haplotype differences, population pairwise F_ST values, and a global F_ST value using Arlequin 3.5 (Excoffier and Lischer 2010) and confirmed these values by GenAlEx 6.5 (Peakall and Smouse 2012). All F_ST values considered only the differences in haplotype frequencies using 10,000 permutations to test for statistical significance (i.e., to test if F_ST > 0.0). We used the Mantel test (Mantel 1967) to evaluate correlations between genetic distances (population pairwise F_ST values) and pairwise geographic distances (Euclidean distance) between sample populations. We tested for isolation by distance in maternal genetic structure using the R package “vegan” using 10,000 permutations to test for significance (Oksanen et al. 2013). We also tested for isolation by distance expressed as F_ST /(1− F_ST) against the natural log-transformed Euclidean distance (Rousset 1997).

For microsatellite data, we used Arlequin 3.5 to calculate pairwise and total F_ST weight averaged over all loci (using analysis of molecular variance tests and 10,000 permutations for significance and 10,000 replicates for confidence intervals bootstrapped over loci). We calculated per-population–based measures of allelic richness, observed and expected heterozygosity, deviation from Hardy–Weinberg proportions (F_IS), and linkage disequilibrium. We double-checked all values for microsatellites and confirmed them in GenAlEx 6.5.

We calculated the expected male to female gene flow ratio using a recently derived equation 7(c) in Hedrick et al. (2013):

In equation (1), m_m is male gene flow and m_f is female gene flow, F_ST is measured overall genetic differentiation for a population (considering both male and female gene flow), and F_ST(f) is the measured genetic differentiation for females in a population using maternally inherited mtDNA markers. Using equation (1), we calculated the ratio of male to female gene flow using global F_ST values calculated for elk populations in the GYE. These estimates of gene flow assume an island model of gene flow, with assumptions that may be violated in some populations. Thus, only relative levels of gene flow can be roughly approximated (Whitlock and McCauley 1999).

Hedrick (2005) and Jost (2008) noted the tendency of F_ST to be lower than expected for populations with high gene diversity (or high heterozygosity). Therefore, as an alternative analysis of male vs. female gene flow, we also computed the standardized genetic differentiation measure G″_ST to remove potential bias due to relatively high haplotype (gene) diversity for mtDNA (Meirmans and Hedrick 2011). We used the SMOGD program to calculate the value of global maternal gene flow (G″_ST(f); Crawford 2010). Confidence intervals are often calculated from bootstrapping over loci in programs such as GenAlEx and Arlequin. However, because mtDNA is treated as a single locus, bootstrapping for mtDNA (single locus) is not available in most genetic programs (Fstat, Arlequin, Genepop, GenAlex, etc.). Instead, we calculated confidence intervals for mtDNA by bootstrapping over individuals. This method of bootstrapping is available in the SMOGD program (which does not provide an F_ST calculation). For microsatellites, we calculated the value for global G″_ST in GenAlEx 6.5 (we arrived at confidence intervals at by bootstrapping over loci). Both programs SMOGD and GenAlEx 6.5 provided nearly identical estimates of both global values of G″_ST. This analysis assumed that G″_ST can replace F_ST in equation (1) (since G″_ST is an analog of the original F_ST), which Hedrick et al. (2013) point out might be more useful. However, Hedrick et al. (2013) also recommend that G″_ST should be investigated theoretically (as they did for F_ST) to ensure it is appropriate to use G″_ST in equation (1).

For mtDNA, haplotype diversity (gene diversity) ranged between 0.75 and 0.91 with an average of 0.85, and 13.0 haplotypes per population. The total number of variable mtDNA nucleotide sites was 27, which defined 30 haplotypes (Table 2). Most substitutions were transitions except for one transversion, which is typical for mammals, including ungulates (e.g., Luikart et al. 2001). The average number of differences between haplotypes was 3.0 nucleotide sites and the mean number of variable sites within populations was 14.6.

For mtDNA, the global F_ST(f) was 0.161 (P  =  0.001). Pairwise F_ST between populations ranged between 0.103 and 0.213 (P < 0.001 for all pairwise comparisons, Table 3). Temporal F_ST values were zero. We found no significant correlation between geographic distance and mtDNA genetic population pairwise distances (r  =  0.274, P  =  0.168). Similarly, a second Mantel test using the relationship of F_ST /(1 − F_ST) against the natural log-transformed Euclidean distance was not significant (r  =  0.202, P  =  0.228). Thus, there was no evidence of geographic isolation by distance among populations.

For microsatellite loci, mean heterozygosity ranged between 0.56 and 0.62 for the study populations, with an average allelic richness of 3.88 (Table 4). Populations and loci were all in or near Hardy–Weinberg proportions, with no significant gametic disequilibrium. The global F_ST from microsatellites was 0.002 (95% CI: 0.000–0.011; P  =  0.332) when averaged across loci (Table 5). The estimated ratio of male to female gene flow was m_m/m_f  =  46 using our global F_ST values derived from mtDNA (F_ST(f)  =  0.161) and from microsatellites (F_ST  =  0.002) in equation 1. Global maternal G″_ST(f) for mtDNA was equal to 0.277 (95% CI: 0.153–0.430). For microsatellites the global G″_ST was 0.005 (95% CI: 0.000–0.030; P  =  0.331). We again calculated the ratio of male to female gene flow and found it to be m_m/m_f  =  37, which was still high but nonetheless lower than the m_m/m_f  =  46 produced from using F_ST values.

Elk in the GYE had a high ratio of male to female gene flow compared to findings for other large mammals in the literature, including some species in the same taxonomic family and genus (Table 6). The global F_ST(f) we computed for mtDNA was 81 times larger than the global F_ST value for microsatellites, which yielded an estimated rate of gene flow that was 46 times higher for male elk than for female elk in the GYE. Similarly, the results from G″_ST (an analog of F_ST) also suggested male-biased gene flow with 37 times higher male-mediated gene flow and with nonoverlapping confidence intervals for G″_ST calculated from mtDNA and microsatellites. In comparison, red deer (Cervus elaphus) in the Scottish Highlands have gene flow for males that is 13 times higher than that for females (Table 6; Pérez-Espona et al. 2009, 2010). Also, Yellowstone bison (Bison bison), present in the same geographical region as our elk, have a moderately high ratio of male to female gene flow of 5 (Halbert et al. 2012).

These comparisons highlight the rather high male-biased gene flow in GYE elk compared to other ungulate species known or suspected to exhibit male-biased dispersal and female natal philopatry (when females return to their birthplace to breed). In this context, dispersal is defined as the movement of individuals from their place of birth to a spatially discrete or distant population, with permanent or long-term settlement for the purpose of breeding (Lowe and Allendorf 2010). Similar to our microsatellite data, Hicks et al. (2007) also report low F_ST values for elk microsatellites (F_ST  =  0.004; P  =  0.281) genotyped from tissue samples from the northern and southern portions of the GYE. Their study used five of the same microsatellite markers (BM5004, BM888, BM4208, FCB193, and BM415) and samples from a similar geographic area (∼260 km between the Hicks et al. [2007] study populations as compared to ∼300 km between our Paradise Valley and Muddy Creek populations from the northern and southern GYE, respectively; Figure 1).

Gene flow estimates from equation 1 assume equal effective population sizes (N_e; or equal variance in reproductive success) for males and females. However, high variance in male reproductive success can reduce local N_e for microsatellite loci (without reducing N_e for mtDNA) and result in higher local genetic drift and thus higher F_ST for microsatellites compared to mtDNA (Hedrick et al. 2013). For example, if only 20% of males reproduce (e.g., due to only a few males dominating reproduction), then the local N_e (for microsatellites) is reduced by approximately 50% (Allendorf et al. 2013:figure 7.3), and the global F_ST for microsatellites roughly doubles from the value expected with random male reproductive success. In other words, if there is high variance in male reproductive success for elk, which likely exists given male dominance and harems, then the expected difference between F_ST values for mtDNA vs. microsatellites would be smaller (all else being equal) because of the relatively reduced local N_e (and higher expected F_ST) for microsatellites. The likely reduced male N_e (compared to females) suggests that an even greater magnitude of male-biased gene flow (than reported here) is needed to explain the much lower F_ST observed for microsatellites compared to mtDNA.

Female gene flow (or genetic structure) was not significantly correlated (r  =  0.274, P  =  0.168) with the straight (Euclidian) geographic distance between populations. The lack of correlation between geographic and genetic distance was best illustrated by populations that have high pairwise F_ST values, but are located close together geographically. For example, Muddy Creek and Bench Corral were separated by a relatively small geographic distance (∼65 km), but have one of the largest pairwise genetic distances (F_ST  =  0.203; Figure 1). The lack of isolation by geographic distance raised the hypothesis that behavioral factors (e.g., seasonal patterns of migration) or landscape features might be more important in explaining female gene flow patterns than geographic (Euclidian) distance.

In summary, maternal gene flow among elk populations in the GYE was low compared to male gene flow, which resulted in high sex-biased gene flow estimates compared to other large mammals. The low female gene flow over distances of 50 to 325 km was an intriguing result for such a (potentially) highly mobile species. Future studies should use many additional populations and apply a landscape genetics approach to test for effects of landscape features on female gene flow because simple geographic (Euclidean) distance did not explain maternal genetic differentiation. The growing availability of genetic data to compare relative male to female gene flow provides many exciting opportunities to explore the magnitude, causes, and implications of sex-biased gene flow likely to occur in many species.

Please note: The Journal of Fish and Wildlife Management is not responsible for the content or functionality of any archived material. Queries should be directed to the corresponding author for the article.

To cite this archived material, please cite both the journal article (formatting found in the Abstract section of this article) and the following recommended format for the archived material.

Data A1. 2013-12-16_mtDNA_8pops. Arlequin (v3.5) formatted file for elk (Cervus elaphus) in Greater Yellowstone National Park. File includes raw mtDNA sequences, sampled between the years of 2005 and 2008, with 596 base pairs for 380 elk from eight populations.

Found at DOI: http://dx.doi.org/10.5061/dryad.8g118

Data A2. 2013-12-16_Microsats_3pops. Genepop formatted file for elk in the Greater Yellowstone Ecosystem. File contains nuclear microsatellite DNA for eight loci genotyped in 59 elk in three populations sampled between the years of 2005 and 2008.

Found at DOI: http://dx.doi.org/10.5061/dryad.8g118

The University of Washington High Throughput Sequencing Center (HTSeq) provided mtDNA sequencing, and K. Pilgrim helped isolate DNA and genotype microsatellites. Fred W. Allendorf provided helpful comments on an earlier version of this paper. Peter Smouse provided helpful advice on interpreting measures of genetic differentiation and computing confidence intervals. Two anonymous reviewers and the Subject Editor provided helpful reviews and comments. National Science Foundation grants supported GL and PC (DEB-0723928, DEB-1258203, DEB-1067613), BKH and MK (DGE-0504628, and RAG DEB-0716188). GL was also supported by grants from the US National Park Service (USA) and a grant (PTDC/BIA-BDE/65625/2006) from the Portuguese Science Foundation (FCT). ADM and MJK were supported by grants from the Wyoming Game and Fish Department, the Rocky Mountain Elk Foundation, the Wyoming Animal Damage Management Board, and the Wyoming Governor's Big Game License Coalition.

Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government.