Introduction

Nutrition is one of the most important keys to understanding the ecology of North American plains bison (Bison bison)1,2,3,4. Many factors such as mineral concentrations and secondary compounds affect grazer nutrition, but protein and energy concentrations of the plants bison eat determine their performance and reproduction, regulate their fine-scale movements as well as long migrations, and structure their interaction with vegetation3,5,6.

Some of the distal determinants and consequences of variation of bison nutrition have begun to be quantified. For example, at the continental scale, bison weight gain is greater in cool, wet climates than in hot, dry climates7. 6.5-year-old male bison in Ordway Prairie, South Dakota were 260 kg heavier than those in the hotter and drier Wichita Mountains, Oklahoma7. This implied that plant protein concentrations aggregated over the year are highest in cool, wet climates. This assertion was subsequently supported by quantifying cattle and bison dietary quality along these climatic gradients4,8,9,10.

Complementary to geographic patterns of bison nutrition, interannual studies of bison performance demonstrated that weight gain and reproduction varies among years11, with reproduction dependent on the mother’s mass12. For example, in Kansas, a female bison that weighed 320 kg in the fall would have an approximately 20% probability or reproduction the next spring, while one that weighed 520 kg would have > 80% probability. In part, the interannual variation in bison performance was tied to variation in timing and amounts of precipitation, but not in easily predictable ways11. For example, in Kansas, years with greater precipitation did not necessarily lead to greater weight gain. Greater precipitation early in the summer led to less weight gain, while precipitation later in the summer lead to greater weight gain. Bison dietary [CP] varies over the year and tends to peak in early summer4,11, which is synchronous with period of greatest nutritional demands for lactating bison.

One of the last links in linking geographic or temporal patterns in bison performance with nutrition is to understand the plant species composition of bison diet. Variation in dietary quality is driven by the nutritional quality of individual plant species, but also the relative intake of different plant species. Bison were once thought to consume predominantly grasses13,14,15,16,17, but recent studies that have quantified bison diet have shown that it to be much broader. For example, in Oklahoma, about half the seeds found in bison fecal material were from eudicots18. In Yellowstone National Park and Alaska, bison regularly browse woody species19,20, something also observed for the closely related European bison21,22. DNA metabarcoding of plains bison fecal material in Kansas, found that only 39% of the protein intake of bison were from grasses4. Even at the northernmost edge of range (Manitoba), less than half of protein bison consumed was from graminoids23. In order to begin to assemble continental patterns of bison diet, 50 bison herds were sampled in June and September 2018 across the United States10. In that study, forbs and legumes contributed over half the protein to bison diet across their range.

Despite these findings and all that has been learned about geographic and temporal patterns of bison performance, dietary quality, and dietary composition, there is still uncertainty about a number of links among these factors. For example, the geographic patterns of bison weight gain suggest that in aggregate bison nutritional quality over the year is greater in cool, wet climates. Despite the 2018 geographic survey of bison dietary quality, it is still unknown whether the greater dietary quality in cold, wet climates observed in June also held earlier in the season. As a general rule, spring onset east of the Rockies Mountains should occur approximately 4 days later per degree of latitude northward24. This would equate to a difference of roughly 52 days between the North Dakota-South Dakota border and central Texas. As such, it is possible that April and May dietary quality could be higher in warmer than cooler climates. It is also unknown what types of plant species bison are relying on early in the season. Similarly, the 2018 project provided no data on dietary quality and composition patterns midsummer, when temperatures are the hottest and nutritional stress is still high from lactation and mating. Regarding interannual patterns, despite the multiple years of data on bison diet for bison in Kansas, we have little understanding of interannual variation of dietary quality and composition for bison across their range. It is unknown whether regional variation in weather could promote some species over others in bison diet, which would potentially explain variation in dietary quality and performance.

Given this background, in order to better understand the seasonal and interannual patterns of bison diet, we collected bison fecal material monthly (April–September) from 45 herds across the US. Fecal material was analyzed with near infrared spectroscopy (NIRS) to quantify dietary crude protein concentrations ([CP]) and digestible organic matter concentrations ([DOM])25. The ratio of [DOM] to [CP] is considered an index of energy to protein limitation with ratios above 7 representing strong protein limitation and below 4 strong energy limitation10,25. Fecal material was also analyzed for dietary composition with metabarcoding from Next Generation Sequencing26. With these data, we sought to (1) understand seasonal patterns of dietary quality and dietary composition, (2) examine how these patterns vary across climate gradients across the six months, and (3) begin to understand interannual variation in dietary quality and composition across the bison’ range with comparisons between the current set of samples and those taken in June and September 2018 across the same climate gradients.

Methods

Sample collection

The goal of selecting bison herds for the study was to include herds where bison grazed on grasslands representative of a given region without any supplementation of protein or energy. To recruit managers of bison herds for the study, participants of the 2018 study were contacted and emails were sent out to additional members of the National Bison Association requesting their participation. Each registered participant was sent a fecal collection kit, which included a six sampling cups, instructions on sampling, a cooler, gloves, disposable spoons, and an icepack. Participants were instructed to combine small amounts of fresh bison fecal samples from each of 10 fecal pats the first week of each month from April to September. Monthly samples were then frozen and stored until all samples were collected and then sent in to the Grazingland Animal Nutrition Lab (GANLab) at Texas A&M. In all, a total of 260 samples were received as some sites did not collect samples every month. Data for 2 sites from CA were excluded as they fell outside of the continental climate envelope we were investigating. Another site was excluded having had supplemented their animals with protein. A total of 248 samples remained in the data set after exclusion (Fig. 1).

Figure 1
figure 1

Map of sampling locations for herds sampled in 2018 (red) and 2019 (blue). Colors are semi-transparent such that herds sampled in both years appear purple. Map generated with maps package in R 3.5.2.

Full size image

Dietary quality

Samples received at GANLab were subsequently dried at 60 °C, ground in a Udy mill to pass a 1-mm screen, and analyzed using Near Infrared Reflectance Spectroscopy (NIRS) to assess dietary quality parameters that included crude protein concentrations [CP] and digestible organic matter concentrations [DOM]. Spectra (400–2500 nm) were collected on a Foss NIRS 6500 scanning monochrometer with spinning cup attachment. Calibration curves were derived from NIRS spectra of cattle fecal material and directly measured forage quality27 and applied to bison here as previously4,10,28.

Dietary composition

After dietary quality was assessed, dried samples were sent to the Jonah Ventures laboratory in Boulder, Colorado for DNA metabarcoding using the c-h primers of the trnL intron in plant chloroplast4,29. Genomic DNA from samples was extracted using the MoBio PowerSoil htp-96 well Isolation Kit (Cat#12,955–4) according to the manufacturer’s protocol. Genomic DNA was eluted into 100 µl and frozen at − 20 °C. Each 25 µL PCR reaction was mixed according to the Promega PCR Master Mix specifications (Promega catalog # M5133, Madison, WI) which included 0.4 µM of each primer and 1 µl of gDNA. Both forward and reverse primers also contained a 5′ adaptor sequence to allow for subsequent indexing and Illumina sequencing. DNA was PCR amplified using the following conditions: initial denaturation at 94 °C for 3 min, followed by 40 cycles of 30 s at 94 °C, 30 s at 55 °C, and 1 min at 72 °C, and a final elongation at 72° C for 10 min. Amplicons were then cleaned by incubating amplicons with Exo1/SAP for 30 min at 37 °C following by inactivation at 95 °C for 5 min and stored at − 20 °C. A second round of PCR was performed to give each sample a unique 12-nucleotide index sequence on the forward read. The indexing PCR included Promega Master mix, 0.5 µM of each primer and 2 µl of template DNA (cleaned amplicon from the first PCR reaction) and consisted of an initial denaturation of 95 °C for 3 min followed by 8 cycles of 95 °C for 30 s, 55 °C for 30 s and 72 °C for 30 s. 5 µl of indexing PCR product of each sample were visualized on a 2% agarose gel to ensure the success of the barcoding PCR. Final indexed amplicons from each sample were cleaned and normalized using SequalPrep Normalization Plates (Life Technologies, Carlsbad, CA). 25 µl of PCR amplicon is purified and normalized using the Life Technologies SequalPrep Normalization kit (cat#A10510-01) according to the manufacturer’s protocol. Samples are then pooled together by adding 5 µl of each normalized sample to the pool. Sequencing occurred on an Illumina MiSeq (San Diego, CA) running the 2 × 150 bp chemistry with a v2 300-cycle kit.

For bioinformatic processing, sequencing success and read quality was verified using FastQC v0.11.8, and reads were demultiplexed by using Illumina-utils v2.6 (iu-demultiplex) using default settings. Sequences of each sample were then merged using the -fastq_mergepairs option in Usearch v11.0.66730. The forward primer (5′-CGAAATCGGTAGACGCTACG-3′) and reverse primer (5′-CCATTGAGTCTCTGCACCTATC-3′) were removed using Cutadapt v1.1831. Cutadapt is also used to discard sequences with length below 108 bp. Expected error filtering as implemented in Usearch is then used to discard low quality reads (max_ee = 0.5)32. Instead of OTU clustering, reads affected by sequencing and PCR errors are then removed using the unoise3 algorithm with an alpha value of 533. This denoising is applied to each individual sample, and Exact Sequence Variants (ESV) compiled in an ESV table including sequences and read counts for each sample. Taxonomy is assigned to each ESV by mapping them against a GenBank reference data34 as well as Jonah Ventures voucher sequences records, using usearch_global with –maxaccepts 0 and –maxrejects 0 to ensure mapping accuracy. Consensus taxonomy is generated from the hit tables, by first considering 100% matches, and then going down in 1% steps until hits are present for each ESV. For final statistical analyses, the top 200 ESVs across all sites were used for analyses. The top 200 ESVs represent 89.89% of all reads and no single ESV had less than 10% of the reads for any one sample. Relative abundances for a sample were standardized so the sum of all reads of the top 200 ESVs was equaled 100%. Each ESV was assigned a representative genus if species from multiple genera matched the ESV.

A functional group assignment for each ESV was then made based on the taxonomy assignment of the ESV. Functional group classifications included: cool-season graminoids (grasses with C3 photosynthesis as well as Carex and Equisetum species), warm-season grasses (all C4 grasses), legumes (Fabaceae species), herbaceous eudicots, i.e. forbs, and woody eudicots.

Climate data

Site location was used to obtain 30-year climate normal (1981–2010) from the Oregon State University PRISM Climate Group data explorer database (PRISM Climate Group).

Statistical analyses

All statistical analyses and figure preparations were conducted in R v. 3.5.235. Beyond standard univariate analyses, a standard least squares model was used to understand the relationship between diet quality and climate overall and for each month. In the models, [CP], [DOM], [DOM]: [CP] were predicted with fixed effects of MAT, MAP, the identity of the month, and all pairwise interactions. To understand whether functional group composition affected dietary quality independent of climate and season, another set of identical models were run but each also included the relative abundance of one of the functional groups. The only functional group that was significant was warm-season grasses and the results of the model including warm-season grasses are reported. In addition to dietary quality, the climate and seasonal influences on functional group composition were tested with models that had the relative abundance of individual functional groups as the response and predictors of MAT, MAP, the identity of the month, and all pairwise interactions. As in no case was the interaction between one of the climate variables and month identity significant, these interactions were removed from the final model. To compare whether dietary quality and composition values in 2018 and 2019 differed, linear models were run for each parameter for each month (June and September) that included just the identity of the year. To test whether relationships between parameters and site climate differed between the years additional linear models were run that included MAP, MAT, identity of the year, the identity of the month, and interactions between 1) MAP and MAT, 2) year and MAT, 3) year and MAP, and 4) year and month identity. Maps of CP and DOM for each month in 2019 and functional composition of diet for June, 2019 were generated utilizing the raster package for reclassifying and masking pixels. Forested areas were masked from the map based on data from the Commission for Environmental Cooperation, http://www.cec.org.

Results

General patterns of dietary quality and composition

Averaged across the months, site-level [CP] ranged from 62.6 to 147.3 mg g−1 and averaged 96.7 ± 3.3 mg g−1. [DOM] ranged from 571.3 to 643.6 mg g−1 and averaged 605.8 ± 2.4 mg g−1 while [DOM]: [CP] ranged from 4.4 to 10.3 with an average of 6.9 ± 0.2.

Examining seasonal patterns averaged across all sites, [CP] averaged 82.5 ± 5.4 mg g−1 in April, peaked at 113.2 ± 3.3 mg g−1 in June, and declined to 80.6 ± 4.5 mg g−1 in September (Supplementary Table S1, Fig. 2). Similarly, [DOM] averaged 583.0 ± 6.1 mg g−1 in April, peaked at 630.2 ± 3.4 mg g−1 in June, and declined to 585.1 ± 4.2 mg g−1 by September (Supplementary Table S1, Fig. 2). [DOM]: [CP] averaged 8.1 ± 0.5 in April, declined to 5.8 ± 0.2 in June, and increased to 8.1 ± 0.4 by September (Supplementary Table S1, Fig. 2).

Figure 2
figure 2

2019 monthly dietary quality metrics averaged (± S.E.) across all sites. Shown are (a) crude protein concentrations ([CP]), (b) digestible organic matter concentrations ([DOM]), and (c) the ratio between [DOM] and [CP].

Full size image

Examining dietary functional group composition across sites, on average 38.2 ± 2.6% of dietary protein intake came from grasses, ranging from 12.7 to 82.2% across sites. Examining patterns for the two grass functional groups, 26.5 ± 2.9% (1.3–82.0%) of protein intake was derived from cool-season (C3) graminoids and 11.7 ± 1.6% (0–36.7%) came from warm-season (C4) grasses. Of the Eudicots, protein intake from legumes averaged 37.8 ± 2.8% (1.3–70.2%), from non-leguminous forbs averaged 21.5 ± 2.0% (0.2–68.0%), and 2.5 ± 0.8% (0–31.6%) from woody species.

Examining monthly patterns across sites, C3 grass protein intake was highest in May (40.8 ± 5.0%) and lowest in September (15.8 ± 2.8%) while C4 grass protein intake peaked in September (16.2 ± 2.9%) and was lowest in July (10.2 ± 2.4%) (Supplementary Table S2, Fig. 3). Legume protein intake was highest in August at 55.8 ± 5.3% and was lowest in May (20.0 ± 4.0%) (Supplementary Table S2, Fig. 3). Forb protein intake peaked in June at 27.6 ± 3.5% and was lowest in August (14.2 ± 2.8%) (Supplementary Table S2, Fig. 3).

Figure 3
figure 3

Average contributions of different functional groups to 2019 dietary protein intake for bison across sites for each month from April to September. Data based on the 200 most abundant Exact Sequence Variants (ESVs), which were assigned a consensus taxonomy and then a function classification. Cool-season graminoids includes C3 grasses, sedges, and Equisetum. Warm-season grass was derived exclusively from grasses. Legumes only included species from Fabaceae.

Full size image

Climate and dietary quality

Examining patterns of [CP] across sites, [CP] was highest in cool, wet climates (Fig. 4, Supplementary Table S3). For example, bison in June at a site with 1200 mm of rain and 6 °C MAT would have [CP] of 186.3 mg g−1. Bison in June at a site with 400 mm of rain and 18 °C MAT would have [CP] of just 67.8 mg g−1. There was also a statistical interaction between the identity of the month and MAT on [CP] (Fig. 4, Supplementary Table S3). In April, [CP] tended to increase with increasing MAT (2.51 ± 1.46 mg g−1 °C−1), as warm-climate sites tended to have higher [CP] than cold-climate sites (Fig. 4, Supplementary Table S3). By May, colder sites tended to have higher [CP] (− 2.01 ± 1.81 mg g−1 °C−1), which became stronger by June (− 5.34 ± 1.83 mg g−1 °C−1) and lasted through September (Fig. 4, Supplementary Table S3). Examining patterns of functional group composition on top of climate relationships, [CP] was lower when greater amounts of warm-season grasses were in the diet, independent of climate and season (Supplementary Table S4).

Figure 4
figure 4

Maps of dietary crude protein and digestible organic matter concentrations for non-forested areas. Maps are modeled based on Supplementary Table S3. Units are mg g−1. Maps generated utilizing the raster package for reclassifying and masking pixels in R 3.5.2.

Full size image

Like [CP], [DOM] was also highest in cool, wet climates (Fig. 4, Supplementary Table S3). For example, bison at a site with 1200 mm of rain and 6 °C MAT in June would have [CP] of 670.4 mg g−1 as opposed to bison at a site with 400 mm of rain and 18 °C MAT in June, which would have [CP] of just 580.4 mg g−1. In April, bison in warm-climate sites tended to consume a diet that was higher in [DOM] than bison in colder-climate sites (2.04 ± 1.73 mg g−1 °C−1) (Fig. 4, Supplementary Table S3). By May, [DOM] was invariant across MAT gradients (− 0.08 ± 1.98 mg g−1 °C−1) and by June, cold-climate sites had higher [DOM] (5.51 ± 2.00 mg g−1 °C−1) (Fig. 4, Supplementary Table S3).

Examining patterns of DOM:CP reveals the seasonal and spatial patterns of protein limitation. In general, bison in cool, wet climates had the lowest [DOM]: [CP] (Supplementary Table S3). Yet, the statistical interaction between MAT and month shows that in April, bison in cold climates were more protein limited than in warm climates (− 0.16 ± 0.11 °C−1) (Supplementary Table S3). By May, [DOM]: [CP] tended to become higher in warm climates, with the gradient becoming strongest by August (0.61 ± 0.15 °C−1) (Supplementary Table S3).

Climate and dietary composition

Examining cross-site relationships between the relative amounts of different functional groups in diet and the predictors of climate and season, the proportion of cool-season grass in the diet was highest in cool sites (decreasing at a rate of 3.8 ± 0.5% °C−1), with no significant influence of MAP (P = 0.16) (Fig. 5, Supplementary Table S5). The proportion of cool-season grass in the diet varied among months, but there was no difference among months in the relationship with MAT (P = 0.54). In contrast, the proportion of warm-season grass in diet was greatest in hot, dry regions (Fig. 5, Supplementary Table S5). For example, bison at a site with 400 mm of rain and 18 °C MAT would have 42.9% of their dietary protein from warm-season grasses as opposed to a site with 1200 mm of rain and 6 °C MAT, where warm-season grasses would provide 0% of their dietary protein.

Figure 5
figure 5

Maps of functional group composition of bison diet (excluding woody plants, which had low abundance and little pattern with climate) for non-forested areas. As there were no interactions between either MAT or MAP and the identity of the month, only maps for June are shown here. See Fig. 2 for how functional group composition of the diet changes over the six six months. Maps generated utilizing the raster package for reclassifying and masking pixels in R 3.5.2.

Full size image

Among Eudicots, there was no seasonal variation in the percentage of forbs in the diet, but the percentage of dietary protein from forbs was highest in hot, dry regions (Fig. 5, Supplementary Table S5). For example, bison at a site with 400 mm of rain and 18 °C MAT in June would have 54.2% of their dietary protein from forbs as opposed to a site with 1200 mm of rain and 6 °C MAT in June, where it would be just 9.6%. Legumes had strong seasonal variation in its dietary contribution, but little variation with climate, increasing at a rate of 2.0 ± 0.8% per 100 mm MAP (P = 0.02) (Fig. 5, Supplementary Table S5).

Comparison with 2018 data

Comparing [CP] between years, June [CP] was higher on average in 2019 than 2018 (113.2 ± 3.3 vs. 93.6 ± 3.2 mg g−1; P < 0.001). In contrast, September [CP] was not different between the two years (80.6 ± 4.5 mg g−1 in 2019 vs. 80.6 ± 4.5 mg g−1 in 2018; P = 0.98). Like [CP], June [DOM] was higher in 2019 than 2018 (630.2 ± 3.4 mg g−1 vs. 617.6 ± 3.7 mg g−1; P = 0.01) while September [DOM] did not differ between years (585.1 ± 4.2 mg g−1 in 2019 vs. 591.7 + 4.4 mg g−1 in 2018; P = 0.28). June [DOM]: [CP] was lower in 2019 vs. 2018 (5.8 ± 0.2 vs. 6.9 ± 0.2; P < 0.001) while September [DOM]: [CP] did not differ between years (8.1 ± 0.4 vs. 8.7 ± 0.6; P = 0.36).

Comparing monthly functional group composition of bison diet between 2018 and 2019, in almost all cases, patterns in June 2019 were somewhat similar to June 2018 (Table 1). For example, June protein intake from cool-season grasses was approximately 31% in both years (32.5% ± 4.4% in 2019 vs. 30.1% ± 3.9% in 2018; P = 0.68; Table 1). In 2019, cool-season graminoids, warm-season grasses, and legumes differed significantly from 2018. In 2019, 45.6% ± 5.1% of September dietary protein intake came from legumes (Table 1). In 2018, it was 24.6% ± 3.5% (P < 0.001). September dietary protein intake from cool-season graminoids was 10% lower in 2019 than 2018 (15.8 ± 2.7% vs. 26.0 ± 3.5%; P = 0.027). For warm-season grasses, it was 14% lower in 2019 than 2018 (30.2% ± 4.1% vs. 16.2% ± 2.9%; P = 0.009) (Table 1).

Table 1 Comparison between years of mean dietary quality and functional group abundance in diet.
Full size table

Comparing relationships between [CP] with climate between the two years, although [CP] was higher on average in 2019 after standardizing for climate (P < 0.001), there was no difference in the relationships between [CP] and either MAT or MAP between years (P > 0.05 for both) (Table 2). For [DOM], 2019 [DOM] was 12.3 mg g−1 higher after standardizing for climate (P = 0.02) (Table 2). The relationships between [DOM] and MAP were similar between years (P = 0.22), with [DOM] decreasing at a faster rate with increasing MAT in 2019 than 2018 (− 2.6 mg g−1 °C−1 higher in 2019 than 2018; P = 0.05) (Table 2). After controlling for site climate, [DOM]: [CP] was 1.0 ± 0.5 lower in 2019 than 2018 (P = 0.02), but the relationships with climate did not differ between years (P > 0.1 for both) (Table 2).

Table 2 Model results for comparing dietary quality for June and September between 2018 and 2019.
Full size table

Controlling for site climate, which would take into account potential differences in the distribution of sites in climate space between the years, there was no difference in the protein contribution of cool-season grasses to bison diet between years (P = 0.77) and no difference in relationships between climate and the protein contribution of cool-season grasses (P > 0.7 for both MAT and MAP) (Table 3). The same was true for warm-season grasses (Table 3). Despite significant differences in legume contribution to diet between years for the different months (P = 0.002), controlling for climate, legume contribution to diet was similar between years (P = 0.34) with no differences in the relationship between legume contribution and climate (P > 0.3 for both MAT and MAP) (Table 3). There were also no significant relationships with climate for forb or woody contribution to diet between years (Table 3). Despite many of the similarities between years, comparing dietary functional group composition in September between years, legumes contributed 20% more protein and warm-season grasses 14% less in 2019 than 2018.

Table 3 Model results for comparing functional group abundance in bison diet for June and September between 2018 and 2019.
Full size table

Discussion

Comparing the relationships between climate and dietary quality between 2018 and 2019 reveals that cooler, wetter sites generally have higher forage quality for bison than warmer, drier sites. Although not invariant, this pattern further strengthens patterns observed for cattle8,9 as well as those inferred from bison weight gain7, where bison from cool, wet climates have the highest weight gain. In contrast to previous work, these patterns did not hold for samples collected in April, where warmer sites had higher forage quality. More than likely, this is due to the earlier phenology of plants in warmer sites. Although warmer sites are more likely have a higher relative abundance of warm-season grasses than cooler sites, which can narrow phenological differences along broad temperature gradients, in April, many northern sites would have just emerged from being covered in snow while southern sites are closer to peak forage quality. This inversion is short-lived, but still would have been a driver of animal migrations assuming that bison could migrate in concert with phenological waves.

Along with the consistency in midsummer relationships between forage quality and climate between the years, functional group composition of bison diet was similar in June 2018 and June 2019. Yet, in September 2019, bison on average consumed almost twice as much protein from legumes in 2018 and correspondingly less warm- and cool-season graminoids. To better understand whether broad-scale differences in weather were behind this pattern, for each of the 45 sites measured in 2019, daily meteorological data were extracted from the Daymet surface weather database, which are provided on a 1-km grid for North America 1980-present36. From these data, mean August (DOY 214–244) temperature and total precipitation were calculated for 2018 and 2019. Examining August temperatures in 2018 and 2019, which presumably would influence differences in early September phenology, August temperatures were 0.1 °C higher in 2019 than 2018 (23.1 °C vs. 23.2 °C) with just 20 more mm of precipitation (109 vs. 89 mm) precipitation. Yet the similarity in average temperature masked latitudinal patterns where southern sites were approximately 1 °C warmer in 2019 and northern sites 1 °C cooler, compared to 2018 (y = 4.6 − 0.117*Latitude, P < 0.001). Single-site repeated measures of diet reveals similar magnitude of interannual variation in dietary composition. For example, across 3 years in a Kansas grassland, peak spring N2-fixing plant intake varied by 15% with a similar magnitude in variation in peak midsummer warm-season grass intake4. Given that, whether the differences among years are driven by differences in production between years, quality of individual species, or selection by bison is a topic for further study.

Although we have largely focused on the broad patterns here, there are still lessons to be learned from individual sites. For example, the diet of bison at Kankakee Sands in Indiana was unique. Kankakee Sands is located on a sand plain and largely consists of restored prairies dominated by warm-season grasses such as Sorghastrum and Schizachyrium. Despite the predominance of the warm-season grasses, the protein intake of bison at Kankakee Sands was dominated by Equisetum and Salix. For example, in June, > 35% of the protein intake was from Equisetum. Salix was > 85% of the protein intake in May. Overall, dietary protein concentrations were still high for the herd, reaching as much as 122 mg g−1 in June. Despite these high values, it appears that in a sea of low-quality grasses, bison concentrate browsing on a few species. What the dietary quality would be without Equisetum and Salix, or whether this diet is sustainable is unknown. This may be a case where management needs to begin promoting cool-season grasses and legumes to maintain high quality diets in the future in case populations of Equisetum and Salix are reduced either through over-consumption or other factors.

Overall, the predominance of eudicots for protein intake for this network of sites across 2018 and 2019 is congruent with past studies4,37, but still differs from assumptions of bison diet13,14,15,16,17. Part of this discrepancy could be explained by initial biases in how dietary composition was assessed, the locations of the studies, and/or the difference in focus on mass intake in earlier studies vs. protein intake in the metabarcoding studies. Still, even if eudicots had triple the protein concentration of grasses, they still would constitute approximately 40% of the mass intake. Before accepting this importance of eudicots in bison diet, there is still some further testing required. For example, the linkage of relative read abundance and relative protein intake requires scaling between chloroplast DNA concentrations of chloroplasts, chloroplast density, and protein concentrations. If some species have more copies of chloroplast DNA per chloroplast they would be “overrepresented” in the sequence counts relative to their protein contribution. Little is known about the patterns of chloroplast DNA concentrations across species and how they may change over time, but they are not constant over the maturity of leaves38. Setting aside the detailed molecular work on chloroplast abundance and cpDNA copies, a simple experiment that measured relationships between protein concentrations and cpDNA copies across leaves of different species would be sufficient to determine whether the generalization of scaling between the relative abundance of cpDNA among species and protein intake was a valid approximation. To date, feeding trials39, comparisons of diet reconstructions from microhistology and fecal DNA40, and isotopic estimates of the relative contribution of C3 and C4 species4,41 all support the validity of fecal DNA for diet reconstructions. One feeding trial did not, but that was conducted with dried hay that likely would have caused cpDNA degradation among species42.

In all, the research presented here illuminates one of the reasons that bison might have migrated long distances in the Great Plains, similar to the Green Wave Hypothesis43,44. Bison that began the spring in southern ranges would have experienced higher protein concentrations than those in northern ranges. Assuming they could have migrated fast enough to follow phenological development, this would have provided them with higher total protein intake than those that did not migrate. A migration rate of ~ 20 km d−1, which is similar to the migration rate of caribou45 and saiga46, would be sufficient to cover the distance between central Texas and southern Nebraska over a 2-month period. Given the interannual variation in dietary quality observed between 2018 and 2019 in June and September, it is likely that this benefit to migration would have varied among years, although more years of monitoring with data covering the entire growing season is required to more fully evaluate this question. That being said, grasslands are different now than they were even a century ago. Forage quality has declined in general47,48 and many northern grasslands are dominated by high-quality, introduced plant species that likely provide more protein than those they replaced. Regardless, it is a simple lesson to learn from this work that improving bison nutrition can be attained by providing access to high-protein plant species. Bison that experience protein stress during the growing season could be supplemented with protein, rangelands could be fertilized, or plants with high protein, such as N2-fixing species, could be promoted.