Introduction

The Bison genus (bison) is native to Asia and central Europe, but migrated to the American continent via the steppe bison (Bison priscus) and the giant bison (Bison latifrons). Current populations of American bison (AB; Bison bison) are the product of adaptation, evolution and natural selection; there are two allopatric subspecies, the plains bison (Bison bison bison) and the mountain bison (Bison bison athabascae). Historical and archaeological data suggest that the AB developed on the North American prairies, with estimated populations as high as 60 million individuals. During the 19^th Century, bison hunting for food and hides decimated the population, bringing it near extinction^1,2,3,4.

In Mexico, there are historical accounts of AB in the states of Chihuahua, Coahuila, Durango and Sonora; the Janos-Hidalgo herd was a transboundary herd that moved between Chihuahua and New Mexico^5,6. There is currently a conservation herd at El Uno Ranch, in the Janos Biosphere Reserve (Chihuahua, Mexico), which was created with 23 individuals from Wind Cave National Park in the United States⁷. As a genetic resource, the AB exhibits the time and space components, as well as use and option values of biodiversity. The time and space components are determined by evolution and changes in species richness, relative abundance and dominance. The use value consists of the benefits provided by the resource, and the option value is defined by a genetic resource’s role in or contribution to ecosystem stability^8,9.

Population biodiversity is the product of adaptation to and integration into ecosystems driven by evolutionary forces and population genetics (e.g., natural selection, genetic drift and migration). Genetic diversity is a component of biodiversity and comprises differences in heritable genetic material. Genetic variability is a measure of genotype differentiation as a function of population size and the criteria used to define inheritance of genetic material. Determined by its evolutionary history, a population’s genetic structure expresses the genetic diversity it harbors and this is distributed within the population. Loss of genetic diversity is the main challenge in at-risk populations, and is therefore a vital concept in the design of conservation schemes^10,11. The present study objective was to analyze the genetic structure and variability of the AB herd at El Uno Ranch using simple nucleotide polymorphism (SNP) genetic markers.

Material and methods

The AB herd at El Uno Ranch exists in a wild environment with almost no human contact, and is isolated and protected from populations of other bovids or other species that could alter its normal development. All ranch personnel are specialized and facilities are exclusively for managing the herd. A herd census and identification is done annually. For the present study, 174 animals (80 % of the total herd) were sampled: 102 females and 72 males born in 2012. Blood deposited on specialized cards in the GeneSeek Laboratory of Neogen^® Corporation was used for DNA extraction. Analyses were done of 42,366 SNP genotypes distributed in 29 chromosomes and defined in the GGP Bovine 50K chip. During editing, loci were discarded if they had a usable percentage (UP) <90 %, were monomorphic, had a minor allele frequency (MAF) <0.01 and/or were in Hardy-Weinberg disequilibrium (HW; P<0.05).

After editing, the SNPs panel was used to estimate six genetic variability indicators: expected heterozygosity (He); observed heterozygosity (Ho); polymorphic information content (PIC); the fixation index (F_ST); the Shannon index (SI) and linkage disequilibrium (LD)^12,13. The LD was evaluated based on the correlation (r²) between haplotype frequencies through loci¹⁴. Correlation (r²) values range from zero to one, with values near zero indicating an absence of LD and independent segregation and those near 1 indicating non-random association between loci. The kinship relationship (Rij; %) was estimated using all the sampled individuals and the effective population size (Ne) based on adjusted average r² via the Waples method¹⁴. Estimates of r² were done using the FSTAT program¹⁵; the GenAlex program¹⁶ was used to estimate He, Ho, PIC and the F_ST; the LDNE program¹⁷ was used to analyze Ne; and ML-Relate¹⁸ was applied to estimate Rij.

Genetic structure was elucidated with the Structure genetic analysis program¹⁹. This uses Bayesian grouping to infer the number of lines or genomes (k) within a population by using genetic markers for genotype analysis. The procedure assumes that individuals are of pure ancestry (k= 1) vs. ancestry of two or more lines (k ≠ 1), and proportionally assigns a genome to each line. Use of Bayesian clustering to infer k is derived from the a posteriori probability distribution generated by the Markov Chain-based Monte Carlo method. Five possible lines were evaluated in the present study, and individuals were assigned to them probabilistically. The number of lines (k) that provides the best fit is derived from the logarithmic likelihood of each sampling step, and the maximum or optimal value was obtained with the approach of Evanno et al.²⁰ and the Structure Harvester program²¹.

Results and discussion

Editing produced a panel with 2,135 identified and selected SNPs (5.04 % yield versus total number of evaluated SNPs), with an average of 74 SNPs per chromosome. A total of 40,231 SNPs were discarded: 84.5 % were monomorphic; 8.5 % by UP<90 %, 6.3 % by MAF<0.01, and 0.70 % by HW in disequilibrium. The panel of selected SNPs had a SI of 0.30, a Ho of 0.187, a He of 0.182, a F_ST of -0.029, and a PIC of 0.152 (Table 1). The He, Ho and PIC values determine genetic marker viability in genetic variability studies. All SI estimates were nearer zero than one, which is associated with homogeneity in the population and reduces uncertainty when predicting the probability of assignment of an individual to a population. In all the chromosomes F_ST had values ranging from -0.002 to -0.062. This indicator measures levels of heterozygosis and homozygosis, and produces values between -1 and 1. Positive values indicate a heterozygote deficit and negative ones an excess. Values near zero are a sign of stability in the homozygous/heterozygous relationship. The average Rij value was 7.61 % based on 15,051 estimates from 174 individuals. Within the 0 to 100 % range of this indicator, the present results could be classified into five strata: 45.1 % of the estimates were equal to zero; 14.5 % were from 0.01 to 4.9; 27.3 % from 5.0 to 19.9; 11.8 % from 20.0 to 49.9; and 1.3 % were equal to or greater than 50.0. An individual’s consanguinity (F) is half the Rij of its parents. Estimates of Rij can therefore be used to select subpopulations for reproduction and conservation schemes, with a view to maintaining F levels.

In the present results Ne was 12.5, and overall average r² was 0.099, with a range of 0.019 to 0.504 (Table 1). Genetic structure analysis using five lines identified three lines in the study population, with proportions of 0.730, 0.157 and 0.113 (Figure 1). Based on Ne, the possible change or increase in inbreeding levels per generation is 4.0 % (ΔF = 1/(2*Ne)). In small populations managed for conservation, increases in F levels indicate loss of genetic variability. This can drive consanguineous depression which can affect population viability, survival, reproduction, disease resistance, and environmental stress, among other factors^22,23. A Ne ≥50 is recommended for populations under conservation management¹⁰, with the aim of keeping any increase in inbreeding at or below 1 % per generation. For example, a study of a European bison population reported estimated Ne values of 7.0 to 28.4 through five generations²⁴. However, population increases did not result in higher Ne values, highlighting the fact that Ne may be influenced by founding population size and that low Ne levels may be associated with genetic drift and greater loss of diversity. A population’s evolutionary potential depends on its genetic variability and Ne values; if Ne is low then genetic drift is strong and may negatively impact its evolutionary potential²⁵.

Each vertical line represents an individual and segment color represents the proportion of each group.

Figure 1 Structure and composition of El Uno Ranch American bison population based on five lines (k= 1, 2, 3, 4, 5). 

In similar studies, the Bovine SNP50K chip was used in three bison populations (one European and two American), producing SNP percentages of 1.8, 2.6 and 2.9, and He estimates of 0.135, 0.197, and 0.199^26,27. A SNP percentage in the same range (2.8 %) was reported for B. bonasus²⁸, although higher values (9.35 %) have also been reported for European bison, with an accompanying Ho of 0.306 and He of 0.250²⁹. Another study of European bison identified 1,536 SNPs, distributed at 8 to 136 SNPs per chromosome with an average of 51.2³⁰.

The current AB population in the United States of America is derived from a genetic bottleneck process with significant variability and genetic structure^31,32. Three genetic lines were defined for the present study population. Given the origins of the El Uno Ranch herd, the line corresponding to the highest proportion probably corresponds to plains bison. A complementary line may be a contribution of the mountain bison and a third was likely generated by separation and development of the studied population. Any differentiation in the study population may have been caused by the genotype-environment interaction, although its adaptation and contribution to the ecosystem may also have had an effect. Genetic isolation between subpopulations affects some demographic and evolutionary processes; the consequent reduced gene flow can lead to accumulation of genetic differences between subpopulations³³. Overall, differences within populations can derive from the genetic diversity of the founding ancestors and their relative contributions, as well as Ne and its evolution over time³². Genetic substructure does not always coincide with obvious morphological or geographic differences between subpopulations. Data from genetic markers and complementary analyses are required to draw contrasts between populations, identify possible sources of genetic material, and/or, where appropriate, define any possible differentiation. For example, a study of eleven bison populations identified genetic differentiation grouped into eight clusters³², while another study identified two genetically distinct subpopulations within the Yellowstone National Park herd³³. Finally, an analysis of genetic structure in twelve bison herds identified three lines or genetic substructures (average constitution = 0.412, 0.303 and 0.285)³⁴.

Conclusions and implications

Three genetic lines were identified within the El Uno Ranch American bison herd. Two may be associated with the source populations while a third is probably linked to the separation process and the effects of natural selection or genetic drift. The present results highlight the need to consider genetic variability and parentage levels when designing reproduction and conservation plans.