https://orcid.org/0000-0002-8197-5138
Abstract
Yellowstone National Park is home to the only plains bison population that has continually existed as wildlife, on the same landscape, through the population bottleneck of the late 19th century. Nevertheless, by the early 1900s, only 23 wild bison were known to have survived poaching. Salvation efforts included the addition of 18 females from Montana and 3 bulls from Texas to augment this population. A century later, nuclear microsatellite-based population-level assessment revealed two genetically distinct bison subpopulations. However, in 2016, an analysis of mitochondrial haplotypes showed the two founding lineages were distributed throughout the park. This study is designed to delineate any current substructure in the Yellowstone bison population by strategically sampling the two major summer breeding herds and the two major winter ranges. Population-level metrics were derived using the same microsatellite loci as the original study along with a newly developed set of highly informative bison-specific single nucleotide polymorphisms. Our analyses reveal that the modern bison in Yellowstone National Park currently consists of one interbreeding population, composed of two subunits.
Introduction
Yellowstone National Park is home to arguably the most iconic bison population. Since its establishment in 1872 as the first National Park in the United States and the world, significant and unwavering efforts have been dedicated to the preservation of the bison population residing there. Not only has this population played a part in bison becoming a national icon but they are also a true testament to the resilience, hope, and strength bison symbolize. Yellowstone bison are the largest free-roaming bison population in the United States and the only plains bison to have continually existed in the wild and on same the landscape since prehistoric times.
The history of this population is a tale of near tragedy turned conservation success story. During the population bottleneck of the late 19th century, the Greater Yellowstone Area bison population narrowly avoided extinction, with an estimated 23 individuals remaining in 1902 (American Bison Society 1908). To preserve this population, additional bison were brought in from private ranches. Eighteen females from the Pablo-Allard herd in Montana, three males from the Goodnight herd in Texas, and four calves from the indigenous herd were used to establish a secondary “introduced” population (Meagher, 1973; Coder 1975). The introduced herd was moved to the Lamar Valley in 1907 and closely day herded or corralled through at least 1915 (Meagher 1973). Meanwhile the indigenous herd was isolated, wintering in the Pelican Valley and summering in the high-elevation grasslands of the Upper Lamar River.
The subsequent recovery of Yellowstone bison is well documented (Meagher 1973; White et al. 2015). Intermingling of the native and introduced herds increased after 1921, when managers fenced the introduced herd out of the Lamar Valley and moved them to areas used by the indigenous herd during summer. Managers also preferentially removed adult males from the introduced herd to increase the indigenous lineage. The annual report of Superintendent Toll in 1929 indicated “[t]here seems to be a gradual intermingling of the wild [Pelican] and tame [Lamar] herd. It has reached a point where it is difficult to distinguish the buffalo of the wild herd from those of the Lamar Valley herd” (Toll 1929). Therefore, it is likely the bison formed a single “Northern Herd” by the 1930s summering together and separating into two wintering units called the “Lamar” and “Pelican” bison (Meagher 1973). In 1936, managers relocated 71 bison from the Lamar bison to the Firehole and Hayden Valleys (Skinner and Alcorn 1942). The animals formed the “Mary Mountain Bison” or “Central Herd.” Population reductions and subsequent recovery likely kept the Northern and Central Herds separated through the 1970s. Movements between the herds were believed to increase through the 1980s when Northern Herd animals wintering in the Pelican Valley began moving to the Hayden Valley and integrating into the Central Herd (Meagher 1993, 1998). By the 1990s, the animals from the Central Herd began moving to wintering areas of the Northern Herd (Meagher 1989, 1993, 1998). Today, all Yellowstone bison roam relatively freely within Yellowstone National Park, with limitations on their distribution into areas outside the park boundaries. During the past decade, the population has fluctuated between 4,000 and 6,000 animals. Modern GPS technology confirms that bison congregate into two main herds (Northern, Central) during the breeding season, some animals from each herd share wintering areas, and some individuals switch breeding areas over time.
In a study by Halbert et al. (2012a), evaluating 661 Yellowstone bison sampled between 1997 and 2003, Yellowstone bison appeared to split into two genetically distinct subpopulations defined by microsatellite genotype diversity and allelic distributions. At that time, a concern raised by Halbert et al. (2012a) involved unequal culling across these two bison populations and the impact this could have on long-term genetic diversity. Subsequently, White and Wallen (2012) suggested “managers should be promoting the conservation of wildness and natural selection to retain adaptive capabilities, rather than preconceived notions of ‘natural’ genetic or population structures that were likely created or exacerbated by human actions.” As this discussion continued, (Halbert et al. 2012a) pointed out that because all modern bison populations are due to anthropogenic activities, it may be best to error on the side of caution and focus Yellowstone bison management to preserve genetic diversity in its, at the time, current state (Halbert et al. 2012a) . Although no definitive conclusion was reached at the time, it was agreed that there should be continued monitoring of genetic variation in Yellowstone bison and use the best available scientific data as a foundation for future bison management.
Following the study by Halbert et al. (2012a), a mitochondrial genome study was conducted to analyze haplotype diversity among Yellowstone bison (Forgacs et al. 2016). In this study, 25 mitochondrial genomes of Yellowstone bison were evaluated and compared to an additional 20 bison from diverse populations. This study revealed two distinct mitochondrial haplotype clades found within Yellowstone and the overall population of bison. This evidence was consistent with historical documentation of a second maternal lineage brought in from the Pablo-Allard herd in the early 1900s (Meagher 1973; Coder 1975). Unlike the previous study (Halbert et al. 2012b), there was no evidence, based on mtDNA haplotypes, to support population subdivision among the Yellowstone bison. In fact, this study was able to identify that bison of both mitochondrial haplotype clades, and therefore maternal lineages, were present in both present-day breeding areas (Forgacs et al. 2016).
However, both studies had limitations to their scope. Namely, neither study included Yellowstone bison sampled from both breeding populations during the same breeding season. In the study conducted by Halbert et al. (2012a), although a portion of the samples were collected within the park boundaries, over 90% of the samples were collected during the winter migration when bison were leaving the park. Considering the northern and central groups have an overlapping winter range in the Gardiner Basin, it was not possible to confidently assign home ranges without associated individual movement data such as radio telemetry (Meagher 1989; Halbert et al. 2012b). Therefore, no definitive conclusions could be drawn as to whether the identified genetic population substructure extended to the breeding groups. In the work of Forgacs et al. (2016), the mitochondrial DNA only portrays a portion of the genetic story. Although able to assign bison to maternal lineages and compare life histories, it failed to characterize levels of admixture within these individuals.
The current study was designed to evaluate population substructure in Yellowstone bison by strategically sampling the two major summer breeding populations, where one would expect to see genetic differentiation in ongoing population subdivision, and the two major winter ranges, where migration patterns of both breeding herds can overlap. Metrics used to evaluate population dynamics were determined using a set of 24 microsatellite loci (Halbert 2003), previously characterized in Yellowstone bison (Halbert et al. 2012), as well as a newly developed set of highly informative bison-specific single nucleotide polymorphisms (SNPs) (Stroupe and Derr 2024). Identifying current structure of the Yellowstone bison population is important to understanding the history of this iconic population, as well as developing strategies for its future conservation.
Methods
Biological material
Tissue biopsy samples were collected from 282 Yellowstone National Park bison. These samples were collected during the summer of 2019 and winter of 2021 in two major ranges within the park, Central and Northern (Fig. 1). The summer collection included 154 samples, 62 from the Central breeding area and 92 from the Northern breeding area. The winter collection included 128 samples, 44 from the Central wintering area and 84 from the Northern wintering area. Within each geographic area, Northern or Central, bison groups were opportunistically located during sample collection periods with up to 30% of the bison within a group randomly sampled. Tissue samples were collected using a Type P Pneu-Dart Biopsy RDD. Tissue samples were extracted from the dart using forceps cleaned in ethanol. Tissue samples were suspended in an ethanol solution and refrigerated prior to lab analysis. Biopsy tissue samples were collected according to the NPS IACUC IMR_YELL_White_Ungulates_2022.A3.
Map of sample collection sites.
Tissue samples were cut into smaller portions in an isolation hood and divided for microsatellite or SNP genotyping. Of the 282 samples, eight were identified as duplicates with SNP data using Sequoia v2.5.3 (Huisman 2017) and confirmed with microsatellite data. Duplicates were removed from final data analysis.
Microsatellites
DNA was extracted from a portion of each biopsy tissue sample using the Gentra Puregene kit (Qiagen) according to manufacture protocols. A set of 24 autosomal microsatellite markers were used in this study including the core set of loci used for parentage determination (Supplementary Table 1) (Schnabel et al. 2000; Halbert 2003). Microsatellite samples were genotyped according to established lab protocol (Halbert 2003) using an ABI 3730 Genetic Analyzer and STRand software (Toonen and Hughes 2001; https://vgl.ucdavis.edu/STRand). Three samples had a call rate below 90% and were removed from further analysis resulting in a final microsatellite data set of 271 samples and 24 microsatellite markers. Of those, 148 samples were collected during the summer breeding season.
Cervus v3.0.7 (Kalinowski et al. 2007) was used for initial overview of population genetic diversity and formatting data set for downstream analysis. A series of principal component analysis (PCA) were used to establish patterns in the data using collection location and time of year to group samples. Ade4 v1.7-22 (Dray and Dufour 2007) was used to calculate eigenvalues and eigenvectors on the final data set and plotted with 95% confidence level ellipses using ggplot (Wickham 2016) in R v4.1.2 (R Core Team 2022). Observed heterozygosity (HO), expected heterozygosity within populations (HS), allelic richness (AR), FST, and FIS with 95% confidence intervals were calculated with the Hierfstat R package (Goudet 2005). FST was calculated using the Weir and Cockerham (1984) equation. Values calculated per individual or per loci were averaged across each population. These analyses were carried out in R v4.1.2 (R Core Team 2022).
Population parameters were evaluated according to the following methods. Evidence of population structure was assessed using STRUCTURE (Pritchard et al. 2000). The data set was evaluated across 20 iterations at each K from 1 to 6. Notable changes in default program setting include a burn-in period of 40,000 replicates and 80,000 Markov chain Monte Carlo replicates. Structure Harvester was used to compile and format the resulting data (Earl and vonHoldt 2012). Results were visualized with all iterations at each K aligned and merged using the R package pophelper v2.3.1 (Francis 2017).
Single nucleotide polymorphisms
Portions of the biopsy tissue samples were sent directly to NeoGen Canada for processing, DNA extraction, and genotyping. Samples were genotyped according to standard Illumnia Infinium HD Ultra Assay protocol guideline on the GGP Equine-Bison chip (NeoGen). Initial data quality control was done by NeoGen using Illumina’s GenomeStudio with a call cutoff of 0.95. Details regarding the development of the SNP panel used in this study are described in Stroupe and Derr (2024). Briefly, a larger set of SNPs was filtered to a panel of 798 autosomal and 13 mitochondrial SNPs that was then evaluated in 995 bison across ten populations. This panel was found to be highly informative for individual identification and population distinction.
The Illumnia SNP Chip final report data was converted to plink lgen format using the script “illumina_to_lgen.R,” originally written by Ryan Schubert (github@RyanSchu), with adjustments made to include metadata specific to our data set. In addition, a plink fam file was created from the sample list. Data were then converted into VCF format using Plink v1.9 (Purcell et al. 2007). Nuclear and mitochondrial SNP genotypes were analyzed separately in downstream analysis.
The nuclear SNP data set was further thinned to remove monomorphic SNPs and SNPs with a genotyping call rate below 90% among the Yellowstone samples, using VCFtools v0.1.16 (Danecek et al. 2011). Two samples had a genotyping call rate below 90% and were therefore removed resulting in a final SNP data set of 272 samples and 725 SNPs. Of those, 151 samples were collected during the summer breeding season.
A series of PCA were used to establish patterns in the data. Plink v1.9 (Purcell et al. 2007; Chang et al. 2015) was used to calculate eigenvalues and eigenvectors on the final data set and plotted with 95% confidence level ellipses using ggplot (Wickham 2016) in R v4.1.2 (R Core Team 2022). Observed heterozygosity (HO), expected heterozygosity within populations (HS), allelic richness (AR), FST, and FIS with 95% confidence intervals were calculated with the Hierfstat R package (Goudet 2005). FST was calculated using the Weir and Cockerham (1984) equation. Values calculated per individual or per loci were averaged across each population. These analyses were carried out in R v4.1.2 (R Core Team 2022).
fastStructure 1.0 (Raj et al. 2014) was then used to determine population structure among the samples with values of K = 1 to K = 6 and visualized using the R package pophelper v2.3.1 (Francis 2017). fastStructure chooseK.py was used to test for the most likely number of subpopulations.
Mitochondrial haplotype assignments were made according to Stroupe and Derr (2024). Briefly, mitochondrial SNPs originally identified by Forgacs et al. (2016) were used to assign interspecific haplotypes (bison or domestic cattle) and intraspecific clade haplotypes (Clade I or Clade II).
Results
Biopsy tissue samples were collected from 282 bison from Yellowstone National Park during the summer breeding season of 2019 and the winter of 2021. Collection location was used to group bison into either the Central or Northern Herd according to previous studies, migration patterns, population history, and observations (Meagher 1973, 1989; Halbert et al. 2012; Geremia et al. 2014) (Fig. 1). Collection time of year was used to distinguish samples during summer breeding season and winter migration. Of these samples, eight were identified as duplicates and removed from final analyses. All duplicate bison were sampled on the same range; however, two bison were sampled during both the summer and winter on the same range. In those cases, the sample collected during the summer breeding season was kept and the winter sample was removed from the data set. Additionally, two samples from the SNP data and three samples from the microsatellite data were removed due to a genotyping call rate below 90%. The summer breeding season samples included 151 bison in the SNP data set and 148 bison in microsatellite data set.
Samples were genotyped at 24 microsatellites (Schnabel et al. 2000; Halbert 2003), 798 autosomal SNPs, and 13 mitochondrial SNPs (Stroupe and Derr 2024). Microsatellite marker frequencies for individual loci can be found in Supplementary Table 1. The autosomal SNP data were further filtered to remove monomorphic variants and variants with a genotyping call rate below 90% resulting in a final data set of 725 autosomal SNPs. All microsatellites and mitochondrial SNPs passed the 90% call rate threshold. Genotyping call rates for the final data sets were 99.94% for microsatellites and 99.67% for SNPs.
Evidence of population structure and differentiation was evaluated using PCA, admixture analysis, and calculations of FST (Figs 2 and 3, Table 1). In the PCA of only samples collected during the summer breeding season, there is not a definitive separation between the breeding herds (Fig. 2). In the PCA based on 725 SNPs, the central herd seems to only contain a subset of the genetic diversity found within Yellowstone, whereas the northern herd has a wider distribution (Fig. 2a). In the PCA based on 24 microsatellites, there is more overlap between the two breeding herds sample distribution (Fig. 2b). The inclusion of samples collected during the winter migration did not substantially change the distribution of samples in the PCA for either herd (Supplementary Fig. 1). In fact, there are minimal differences between the PCAs of summer samples only and all samples besides the density of sample points. Additionally, admixture analyses, fastStructure, with 725 SNPs, and STRUCTURE, with 24 microsatellites, did not separate the geographically defined herds into genetically distinct subpopulations (Fig. 3, Supplementary Fig. 2). In the fastStructure analysis ran from K = 1 to K = 6 for only the samples collected during the breeding season, both the model complexity that maximizes marginal likelihood and model components used to explain the data structure were revealed to be equal, meaning the best fit for this data is one population (Fig. 3a). This was the same when including winter samples (Supplementary Fig. 2a). In the microsatellite-based STRUCTURE analysis, there was again no distinction between the two summer breeding herds (Fig. 3b). Evaluation of multiple runs revealed the best fit as K = 3 according to the delta K (Evanno et al. 2005), meaning there is evidence for three genetically defined clusters. However, this did not extend to geographically defined breeding herds. In fact, all samples seem to have an even representation of all clusters at all instances of K. Inclusion of samples collected during the winter migration was consistent with the summer only analysis in admixture of individuals; however, the best fit was K = 2 instead (Supplementary Fig. 2b).
Measures of population differentiation (FST) between each grouping of Yellowstone bison defined by sample collection location and time of year. Estimates of FST are based on 725 SNPs or 24 microsatellite markers.
| SNPs | ||||
|---|---|---|---|---|
| FST | Central Herd—Summer | Central Herd—Winter | Northern Herd—Summer | Northern Herd—Winter |
| Central Herd—Summer | −0.0007 | 0.0049 | 0.0018 | |
| Central Herd—Winter | −0.0007 | 0.0050 | 0.0011 | |
| Northern Herd—Summer | 0.0049 | 0.0050 | 0.0013 | |
| Northern Herd—Winter | 0.0018 | 0.0011 | 0.0013 | |
| FST 95% CI | Central Herd—Summer | Central Herd—Winter | Northern Herd—Summer | Northern Herd—Winter |
| Central Herd—Summer | [−0.0019,0.0004] | [0.0039, 0.0058] | [0.0004, 0.0028] | |
| Central Herd—Winter | [−0.0019,0.0004] | [0.0037, 0.0067] | [−0.0001, 0.0020] | |
| Northern Herd—Summer | [0.0039, 0.0058] | [0.0037, 0.0067] | [0.0005, 0.0020] | |
| Northern Herd—Winter | [0.0004, 0.0028] | [−0.0001, 0.0020] | [0.0005, 0.0020] | |
| Microsatellites | ||||
| FST | Central Herd—Summer | Central Herd—Winter | Northern Herd—Summer | Northern Herd—Winter |
| Central Herd—Summer | −0.0043 | 0.0044 | 0.0031 | |
| Central Herd—Winter | −0.0043 | 0.0008 | 0.0000 | |
| Northern Herd—Summer | 0.0044 | 0.0008 | 0.0007 | |
| Northern Herd—Winter | 0.0031 | 0.0000 | 0.0007 | |
| FST 95% CI | Central Herd—Summer | Central Herd—Winter | Northern Herd—Summer | Northern Herd—Winter |
| Central Herd—Summer | [−0.0059, −0.0015] | [0.0017, 0.0077] | [−0.003, 0.0075] | |
| Central Herd—Winter | [−0.0059, −0.0015] | [−0.0014, 0.0031] | [−0.0036, 0.0043] | |
| Northern Herd—Summer | [0.0017, 0.0077] | [−0.0014, 0.0031] | [−0.0005, 0.0022] | |
| Northern Herd—Winter | [−0.003, 0.0075] | [−0.0036, 0.0043] | [−0.0005, 0.0022] | |
| SNPs | ||||
|---|---|---|---|---|
| FST | Central Herd—Summer | Central Herd—Winter | Northern Herd—Summer | Northern Herd—Winter |
| Central Herd—Summer | −0.0007 | 0.0049 | 0.0018 | |
| Central Herd—Winter | −0.0007 | 0.0050 | 0.0011 | |
| Northern Herd—Summer | 0.0049 | 0.0050 | 0.0013 | |
| Northern Herd—Winter | 0.0018 | 0.0011 | 0.0013 | |
| FST 95% CI | Central Herd—Summer | Central Herd—Winter | Northern Herd—Summer | Northern Herd—Winter |
| Central Herd—Summer | [−0.0019,0.0004] | [0.0039, 0.0058] | [0.0004, 0.0028] | |
| Central Herd—Winter | [−0.0019,0.0004] | [0.0037, 0.0067] | [−0.0001, 0.0020] | |
| Northern Herd—Summer | [0.0039, 0.0058] | [0.0037, 0.0067] | [0.0005, 0.0020] | |
| Northern Herd—Winter | [0.0004, 0.0028] | [−0.0001, 0.0020] | [0.0005, 0.0020] | |
| Microsatellites | ||||
| FST | Central Herd—Summer | Central Herd—Winter | Northern Herd—Summer | Northern Herd—Winter |
| Central Herd—Summer | −0.0043 | 0.0044 | 0.0031 | |
| Central Herd—Winter | −0.0043 | 0.0008 | 0.0000 | |
| Northern Herd—Summer | 0.0044 | 0.0008 | 0.0007 | |
| Northern Herd—Winter | 0.0031 | 0.0000 | 0.0007 | |
| FST 95% CI | Central Herd—Summer | Central Herd—Winter | Northern Herd—Summer | Northern Herd—Winter |
| Central Herd—Summer | [−0.0059, −0.0015] | [0.0017, 0.0077] | [−0.003, 0.0075] | |
| Central Herd—Winter | [−0.0059, −0.0015] | [−0.0014, 0.0031] | [−0.0036, 0.0043] | |
| Northern Herd—Summer | [0.0017, 0.0077] | [−0.0014, 0.0031] | [−0.0005, 0.0022] | |
| Northern Herd—Winter | [−0.003, 0.0075] | [−0.0036, 0.0043] | [−0.0005, 0.0022] | |
Measures of population differentiation (FST) between each grouping of Yellowstone bison defined by sample collection location and time of year. Estimates of FST are based on 725 SNPs or 24 microsatellite markers.
| SNPs | ||||
|---|---|---|---|---|
| FST | Central Herd—Summer | Central Herd—Winter | Northern Herd—Summer | Northern Herd—Winter |
| Central Herd—Summer | −0.0007 | 0.0049 | 0.0018 | |
| Central Herd—Winter | −0.0007 | 0.0050 | 0.0011 | |
| Northern Herd—Summer | 0.0049 | 0.0050 | 0.0013 | |
| Northern Herd—Winter | 0.0018 | 0.0011 | 0.0013 | |
| FST 95% CI | Central Herd—Summer | Central Herd—Winter | Northern Herd—Summer | Northern Herd—Winter |
| Central Herd—Summer | [−0.0019,0.0004] | [0.0039, 0.0058] | [0.0004, 0.0028] | |
| Central Herd—Winter | [−0.0019,0.0004] | [0.0037, 0.0067] | [−0.0001, 0.0020] | |
| Northern Herd—Summer | [0.0039, 0.0058] | [0.0037, 0.0067] | [0.0005, 0.0020] | |
| Northern Herd—Winter | [0.0004, 0.0028] | [−0.0001, 0.0020] | [0.0005, 0.0020] | |
| Microsatellites | ||||
| FST | Central Herd—Summer | Central Herd—Winter | Northern Herd—Summer | Northern Herd—Winter |
| Central Herd—Summer | −0.0043 | 0.0044 | 0.0031 | |
| Central Herd—Winter | −0.0043 | 0.0008 | 0.0000 | |
| Northern Herd—Summer | 0.0044 | 0.0008 | 0.0007 | |
| Northern Herd—Winter | 0.0031 | 0.0000 | 0.0007 | |
| FST 95% CI | Central Herd—Summer | Central Herd—Winter | Northern Herd—Summer | Northern Herd—Winter |
| Central Herd—Summer | [−0.0059, −0.0015] | [0.0017, 0.0077] | [−0.003, 0.0075] | |
| Central Herd—Winter | [−0.0059, −0.0015] | [−0.0014, 0.0031] | [−0.0036, 0.0043] | |
| Northern Herd—Summer | [0.0017, 0.0077] | [−0.0014, 0.0031] | [−0.0005, 0.0022] | |
| Northern Herd—Winter | [−0.003, 0.0075] | [−0.0036, 0.0043] | [−0.0005, 0.0022] | |
| SNPs | ||||
|---|---|---|---|---|
| FST | Central Herd—Summer | Central Herd—Winter | Northern Herd—Summer | Northern Herd—Winter |
| Central Herd—Summer | −0.0007 | 0.0049 | 0.0018 | |
| Central Herd—Winter | −0.0007 | 0.0050 | 0.0011 | |
| Northern Herd—Summer | 0.0049 | 0.0050 | 0.0013 | |
| Northern Herd—Winter | 0.0018 | 0.0011 | 0.0013 | |
| FST 95% CI | Central Herd—Summer | Central Herd—Winter | Northern Herd—Summer | Northern Herd—Winter |
| Central Herd—Summer | [−0.0019,0.0004] | [0.0039, 0.0058] | [0.0004, 0.0028] | |
| Central Herd—Winter | [−0.0019,0.0004] | [0.0037, 0.0067] | [−0.0001, 0.0020] | |
| Northern Herd—Summer | [0.0039, 0.0058] | [0.0037, 0.0067] | [0.0005, 0.0020] | |
| Northern Herd—Winter | [0.0004, 0.0028] | [−0.0001, 0.0020] | [0.0005, 0.0020] | |
| Microsatellites | ||||
| FST | Central Herd—Summer | Central Herd—Winter | Northern Herd—Summer | Northern Herd—Winter |
| Central Herd—Summer | −0.0043 | 0.0044 | 0.0031 | |
| Central Herd—Winter | −0.0043 | 0.0008 | 0.0000 | |
| Northern Herd—Summer | 0.0044 | 0.0008 | 0.0007 | |
| Northern Herd—Winter | 0.0031 | 0.0000 | 0.0007 | |
| FST 95% CI | Central Herd—Summer | Central Herd—Winter | Northern Herd—Summer | Northern Herd—Winter |
| Central Herd—Summer | [−0.0059, −0.0015] | [0.0017, 0.0077] | [−0.003, 0.0075] | |
| Central Herd—Winter | [−0.0059, −0.0015] | [−0.0014, 0.0031] | [−0.0036, 0.0043] | |
| Northern Herd—Summer | [0.0017, 0.0077] | [−0.0014, 0.0031] | [−0.0005, 0.0022] | |
| Northern Herd—Winter | [−0.003, 0.0075] | [−0.0036, 0.0043] | [−0.0005, 0.0022] | |
Principal component analysis (PCA) plot of the summer breeding herds of Yellowstone bison using 725 SNPs (a) or 24 microsatellite markers (b). Breeding herds are identified by color with ellipses of 95% confidence intervals.
Clustering plots of the Yellowstone bison breeding herds for K = 1 to K = 6 using fastStructure with 725 SNPs (a) or Structure with 24 microsatellite markers (b). Each vertical line represents an individual sample and horizontal bars represent a run for each K. Plots were generated using pophelper v2.3.1 (Francis 2017).
Estimates of FST were used to measure genetic differentiation among each collection group using both the SNP and microsatellite data sets (Table 1). All FST estimates were less than or equal to 0.0050, with the confidence intervals for most comparisons overlapping and including zero. The central herd summer and winter samples had the lowest measures of differentiation in both the SNP and microsatellite comparisons, FST = −0.007 and −0.0043 respectively. In the SNP comparison, the highest level of differentiation was between the central winter and northern summer samples (FST = 0.0050) though the 95% confidence interval overlapped with the central and northern summer comparison (FST = 0.0049). However, in the microsatellite comparison, the highest differentiation was between the central and northern summer samples (FST = 0.0044) though the 95% confidence interval overlapped with all comparisons except for the central herd summer and winter comparison.
To characterize the genetic diversity of each collection group (central herd—summer, central herd—winter, northern herd—summer, and northern herd—winter) estimates of mean observed heterozygosity (HO), mean gene diversity within populations (HS), mean FIS, and mean allelic richness (AR) were calculated across each group for 725 autosomal SNPs and 24 microsatellites (Supplementary Table 2). Estimates of genetic diversity had a larger range when calculated with the microsatellite data set compared to the SNP data set. In the SNP data set, the central herd had on average higher observed heterozygosity than the northern herd. Both central groups were above the overall mean (HO = 0.4239), whereas both northern groups were below the overall mean. Though in the microsatellite data set, the central summer was above the overall mean (HO = 0.5906) and the other groups were below the overall mean with the central winter group estimate the lowest. All 95% confidence level range estimates of mean FIS calculations overlapped with each other and the overall mean 95% confidence interval in both data sets.
All Yellowstone animals had bison-derived mitochondrial haplotypes based on ten interspecific mtDNA SNPs (Ward et al. 1999; Forgacs et al. 2016; Forgacs 2019). No sample had evidence of domestic cattle mitochondrial DNA. Additionally, three SNPs distinguished the major intraspecific mitochondrial clades found in bison (Forgacs et al. 2016). All 272 samples in the final SNP data set were assigned a mitochondrial haplotype based on the criteria outlined in Stroupe and Derr (2024). Of the two bison-specific mitochondrial haplotype clades, Clade II had the highest frequency of 65% in overall assignment (Fig. 4). Within each group, Clade II was also the predominant mitochondrial haplotype. However, the central herd winter samples had the highest frequency of Clade II followed by central summer, northern winter, and northern summer with frequencies of 78%, 70%, 65%, and 57%, respectively.
Percentages of samples that were assigned to each mitochondrial haplotype, Clade I or II, for each grouping of Yellowstone bison, as defined by collection location and time of year. Bison mitochondrial clades are described by Forgacs et al. (2016).
Discussion
No evidence supports the hypothesis that the bison in Yellowstone National Park are currently composed of genetically distinct and independently breeding subpopulations. The presented genetic analyses using both microsatellite and SNP markers did not reveal substantial differentiation between bison sampled in the northern and central ranges during the summer breeding season. Moreover, analyses showed clear support for considering the bison in Yellowstone as one interbreeding population, composed of two subunits. However, there is undeniable evidence of multiple genetic lineages contributing to the current genetic diversity as can be seen from historic documentation and the presence of multiple maternal lineages.
The sets of microsatellites and SNPs used in this study have both been used in previous population genetic studies and have revealed genetic distinctions between related populations (Halbert and Derr 2008; Stroupe and Derr 2024). Therefore, it is unlikely the lack of population subdivision is due to a lack of sensitivity in the selected genetic markers. Both previous studies revealed an observational difference in multiple measures of genetic diversity between Yellowstone National Park and other federal bison populations. Based on estimates using the same set of SNP markers, the genetic differentiation between the Yellowstone National Park breeding herds is lower than between closely related federal bison populations such as Theodore Roosevelt National Park’s North and South units (FST = 0.0049 vs 0.0884). Furthermore, the level of genetic differentiation in the Yellowstone breeding herds is more comparable to two Turner Enterprise Inc. bison populations, Snowcrest Ranch and Vermejo Park Ranch, that are derived from the same source and have been separated less than 8 yr (FST = 0.0049 vs 0.0030) (Stroupe and Derr 2024).
Although the SNP and microsatellite analyses both revealed the lack of subdivision among Yellowstone bison, there were some observed differences. PCA and admixture analysis using SNPs seemed to reveal a higher sensitivity of subtle differences. PCA of SNPs exposed the central herd as a subset of the whole, whereas the microsatellite PCA showed no observable differences in the herds. Admixture analysis of the SNPs showed variable levels of admixture among samples, whereas the microsatellite data was congruent across samples. Differences could be due to the SNP data set being more sensitive to individual signals of admixture derived from multiple genetic lineages compared to microsatellites.
Due to anthropogenic movements of bison since the population bottleneck in the late 19th century, all modern bison populations are derived from multiple historic lineages, and all have evidence of domestic cattle introgression (Coder 1975; Stroupe et al. 2022). Therefore, when engaging in tasks such as gauging population divergence and migration rates, a comprehensive understanding of a population’s historical trajectory is imperative for accurate interpretation of findings.
During the population bottleneck in the late 19th century, bison in the Greater Yellowstone Area consisted of a dwindling population of indigenous bison that had persisted since prehistoric times. However, when the population reached an estimated low of 23 individuals in 1902, genetically distinct animals were brought in from Montana (Pablo-Allard herd) and Texas (Goodnight herd) to form an introduced population in the northern range (Meagher 1973). Therefore, the bison in Yellowstone would be better described as initially two subpopulations of genetically distinct lineages that have become a single interbreeding population through gene flow, population growth, range expansion, response to environmental pressures, and migration instead of a historic divergence then convergence of a single source population (Meagher 1989, 2002; White and Wallen 2012).
Although previous studies provided evidence of genetic subdivision in Yellowstone bison, there is no evidence this persists in the modern population (Halbert et al. 2012b). No conclusions at the time could be drawn as to whether the identified genetic population substructure extended to the breeding groups due to limited location data during the summer breeding season. Additionally, two decades (1997–2003 vs 2019–2021) separated when the population was sampled and showed subdivision compared to the presented study, which did not. It is not clear if the observed differences represent a change in population dynamics over time, management actions, increase in apex predators, sampling strategy and study design, or a combination of these factors.
Our estimates of genetic differentiation between the central and northern breeding herds were much lower (Microsatellite FST = 0.0044, SNP FST = 0.0049) compared to previous estimates (FST = 0.0321) (Halbert et al. 2012b). This difference could represent a change in population structure over time or could be due to grouping samples by genetic cluster rather than geographic distribution or a combination thereof. Moreover, the previously estimated migration rate of two bison per generation between the two identified clusters is likely underestimated due to the separation of samples based on genetic clustering rather than location and lack of samples collected during the breeding season (Halbert et al. 2012b).
Throughout the years, there have been many changes in conditions and management that could alter the behavior and movement patterns of bison (White et al. 2015). Yellowstone National Park is one of the only places where bison population size and behavior freely respond to environmental changes such as predators, resource limitations, and climate. Furthermore, the Yellowstone bison population also fluctuates due to culling bison during the winter at the parks boundaries to reduce numbers and the risk of brucellosis transmission to domestic cattle (White et al. 2015; Geremia 2022). In some years, this management practice has resulted in the removal of over 1,000 animals from this population, with disproportionate removals from the central and northern herds (White et al. 2011; Geremia 2022).
In addition, we were able to evaluate maternal lineages and mitochondrial haplotype distribution within the Yellowstone bison population. Similar to the findings of Forgacs et al. (2016), mitochondrial haplotype clades I and II were distributed across the central and northern ranges in both the summer (breeding) and winter. These two mitochondrial haplotype clades likely reflect the mixed lineage history of indigenous and introduced animals that comprise the current Yellowstone National Park bison herd. Overall mitochondrial haplotypes from Clade II were found in higher frequency than Clade I haplotypes (65% for Clade II and 35% for Clade I) and Clade II was also found in higher frequency in both the Central and Northern populations in summer and winter collections (Fig. 4).
This reevaluation of the Yellowstone bison population to delineate the current population structure was able to improve upon previous studies by evaluating the two major summer breeding populations with previously used methods and new SNP-based technologies. The SNP panel has proved effective in population differentiation among closely related bison populations and estimating measures of genetic diversity (Stroupe and Derr 2024). The SNP-based evaluation of this important population provides a point-in-time observation of the Yellowstone bison. With the increased efficiency of SNP-based platforms due to large number of loci, consistency, easier automation (Anderson and Garza 2006), lower mutation rates (Amorim and Pereira 2005), reduced influence of inbreeding (Fernández et al. 2013), ease of data sharing (Forcina and Leonard 2020), and an increase in precision and power for population genetics analyses over microsatellites (Zimmerman et al. 2020), there is more opportunity for utilization in long-term management. Moreover, samples collected during the breeding season are imperative for establishing the current structure since that is when genetic exchange between populations occurs and migration patterns can overlap in the winter ranges. Though the addition of samples collected in the winter in our analyses showed similar results.
Here we present evidence, developed from multiple analyses, that indicate bison within Yellowstone National Park represent a single interbreeding population. Even though there are multiple breeding herds and clear evidence of historical bison lineages, it appears substantial gene flow is occurring throughout the population. Thus, there is no way to confidently assign individuals outside of the summer breeding season to their respective breeding herd without tracking individual movement histories. The bison at Yellowstone National Park are an important biological resource that is essential to the long-term conservation of this species. Continued genetic monitoring is imperative to track genetic diversity indices of both nuclear and mitochondrial DNA in order to maintain stewardship of this iconic bison resource.
Supplementary Material
Supplementary material can be found at http://www.jhered.oxfordjournals.org/.
Acknowledgments
Portions of this research were conducted with the advanced computing resources provided by Texas A&M High Performance Research Computing.
Funding
This work was supported by the Department of Interior, National Park Service [grant M1802772]. S.S. support was provided by a Throlson American Bison Foundation Scholarship, The Houston Safari Club, Jim Womack Endowed Fund in Animal Genomics, and the School of Veterinary Medicine and Biomedical Sciences, Texas A&M University.
Data Avaliability
The data underlying this article are available in the article and in its online supplementary material. Additional data will be shared on reasonable request to the corresponding author.
