Seasonal patterns of bison diet across climate gradients in North America

North American plains bison (Bison bison) have been reintroduced across their former range, yet we know too little about their current diet to understand what drove their past migrations as well as observed continental-scale variation in weight gain and reproduction. In order to better understand the seasonal diets of bison at the continental scale, bison fecal material was collected monthly from April to September in 2019 across 45 sites throughout the conterminous United States. Fecal material was analyzed for dietary quality using near infrared spectroscopy and dietary composition with DNA metabarcoding. As observed in previous research, dietary quality peaked in June and was on average greatest for sites with cold, wet climates. Yet, in April, dietary quality was highest in warmer regions, likely reflecting earlier phenology of plants in southern than northern regions. Independent of climate and season, bison that consumed more warm-season grasses had lower dietary protein concentrations. Interpreting the relative abundance of sequences from different plant species as the relative intake of protein from those species, only 38% of bison protein intake came from grasses. An equal amount of dietary protein came from legumes (38%) and 22% from non-leguminous forbs. Seasonal shifts in bison diet were also clear, in part, following the phenology of functional groups. For example, cool-season grass protein intake was highest in May, while legume protein intake was highest in August. Comparing data taken in June and September 2018 in a previous study with corresponding data in 2019, on average, June [CP] was 20% higher in 2019 than 2018, while September [CP] did not differ between years. Dietary functional group composition was generally similar in amounts and relationships with climate between years, yet in September 2019, legumes contributed 20% more protein and warm-season grasses 14% less than in September 2018. In all, this research demonstrates that bison consistently rely on eudicots for protein with the functional group composition of their diet in some ways consistent across space and time, but also spatially and temporally variable. The early-season inversion of plant quality gradients would have been a strong driver of migratory behavior for large numbers of bison optimizing protein intake. As most bison currently experience protein deficiency, optimizing protein intake under current non-migratory conditions will require increasing the relative abundance of high-protein species such as N2-fixing species.

www.nature.com/scientificreports/ but not in easily predictable ways 11 . 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 summer 4,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 grasses [13][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 eudicots 18 . In Yellowstone National Park and Alaska, bison regularly browse woody species 19,20 , something also observed for the closely related European bison 21,22 . DNA metabarcoding of plains bison fecal material in Kansas, found that only 39% of the protein intake of bison were from grasses 4 . Even at the northernmost edge of range (Manitoba), less than half of protein bison consumed was from graminoids 23 . In order to begin to assemble continental patterns of bison diet, 50 bison herds were sampled in June and September 2018 across the United States 10 . 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 northward 24 . 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 limitation 10,25 . Fecal material was also analyzed for dietary composition with metabarcoding from Next Generation Sequencing 26 . 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).
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 quality 27 and applied to bison here as previously 4 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.667 30 . The forward primer (5′-CGA AAT CGG TAG ACG CTA CG-3′) and reverse primer (5′-CCA TTG AGT CTC TGC ACC TATC-3′) were removed using Cutadapt v1.18 31 . 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 5 33 . 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 data 34 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 C 3 photosynthesis as well as Carex and Equisetum species), warm-season grasses (all C 4 grasses), legumes (Fabaceae species), herbaceous eudicots, i.e. forbs, and woody eudicots. 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 warmseason 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  Table S1, Fig. 2).
Examining monthly patterns across sites, C 3 grass protein intake was highest in May (40.8 ± 5.0%) and lowest in September (15.8 ± 2.8%) while C 4 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). 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 warmseason grasses were in the diet, independent of climate and season (Supplementary Table S4 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).  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 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.
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).  (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    www.nature.com/scientificreports/ 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). 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  www.nature.com/scientificreports/ 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.

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 cattle 8,9 as well as those inferred from bison weight gain 7 , 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. www.nature.com/scientificreports/ 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-present 36 . 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 N 2 -fixing plant intake varied by 15% with a similar magnitude in variation in peak midsummer warm-season grass intake 4 . 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 studies 4,37 , but still differs from assumptions of bison diet [13][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 leaves 38 . 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 trials 39 , comparisons of diet reconstructions from microhistology and fecal DNA 40 , and isotopic estimates of the relative contribution of C 3 and C 4 species 4,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 species 42 .
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 Hypothesis 43,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 caribou 45 and saiga 46 , 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 general 47,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 N 2 -fixing species, could be promoted.

Data availability
The datasets generated during and analysed for the current study are available in the Dryad repository, https:// doi. org/ 10. 5061/ dryad. k3j9k d56m. www.nature.com/scientificreports/