Remotely Sensed Spatiotemporal Variation in Crude Protein of Shortgrass Steppe Forage

In the Great Plains of central North America, sustainable livestock production is dependent on matching the timing of forage availability and quality with animal intake demands. Advances in remote sensing technology provide accurate information for forage quantity. However, similar efforts for forage quality are lacking. Crude protein (CP) content is one of the most relevant forage quality determinants of individual animal intake, especially below an 8% threshold for growing animals. In a set of shortgrass steppe paddocks with contrasting botanical composition, we (1) modeled the spatiotemporal variation in field estimates of CP content against seven spectral MODIS bands, and (2) used the model to assess the risk of reaching the 8% CP content threshold during the grazing season for paddocks with light, moderate, or heavy grazing intensities for the last 22 years (2000–2021). Our calibrated model explained up to 69% of the spatiotemporal variation in CP content. Different from previous investigations, our model was partially independent of NDVI, as it included the green and red portions of the spectrum as direct predictors of CP content. From 2000 to 2021, the model predicted that CP content was a limiting factor for growth of yearling cattle in 80% of the years for about 60% of the mid-May to October grazing season. The risk of forage quality being below the CP content threshold increases as the grazing season progresses, suggesting that ranchers across this rangeland region could benefit from remotely sensed CP content to proactively remove yearling cattle earlier than the traditional October date or to strategically provide supplemental protein sources to


Introduction
In arid and semi-arid rangeland ecosystems, sustainable livestock production is dependent on matching the timing of forage availability and associated quality with animal intake demands [1][2][3][4]. Ranchers often develop grazing plans that use established guides for initial stocking rates [5,6] and diverse management strategies to accomplish operational goals and desired outcomes [7]. Both intra-and inter-annual variations in forage availability, however, provide ranchers substantial challenges in matching with animal demand [8][9][10]. These challenges have major impacts on ranching operation economics, gains can increase economic returns by 10-25% for ranchers by removing yearling grazing cattle from grazing rangelands earlier than is typical [13]. This earlier removal of yearlings would be most effective if implemented when forage CP content declines below the 8% threshold. It is unclear, however, how often the 8% CP content threshold is reached during the grazing season [25][26][27][28].
Our main objective was to quantify the level of risk experienced by farmers within the shortgrass steppe due to inadequate CP content of forage during the growing season. To reach this main objective, we proposed two specific objectives. First, we evaluated the associations between field estimates of CP content from paddocks with contrasting forage availability and botanical composition with seven different spectral bands acquired through MODIS satellite imagery for a semi-arid shortgrass steppe in North America. Second, using a calibrated model developed from these associations, we then assessed risk levels of reaching the 8% CP content threshold during the grazing season for paddocks with longterm differences in their grazing intensity (light, moderate, or heavy grazing intensities).

Site Description
Research was conducted on the USDA-Agricultural Research Service Central Plains Experimental Range in north-central Colorado, USA (40 • 49 N, 107 • 46 W), a Long-Term Agroecosystem Research (LTAR) network site. Mean annual precipitation is 340 mm, with 71% of this occurring during the growing season (April-August). Dominant ecological sites are Loamy Plains (R067BY002CO) and Sandy Plains (R06BW024CO) [47,48] with paddocks either being dominated by one or mixed. Ecological sites are distinct types of land with specific soils that produce distinct types and amount of vegetation. The C4, shortgrass blue grama (Bouteloua gracilis) is the dominant plant species for the Loams Plains ecological site. The C3, mid-height grasses needle and thread (Hesperostipa comata) and western wheatgass (Pascopyrum smithii) dominate the more productive plant community (>75% more forage [6,49]) for the Sandy Plains.

Model Calibration
CP content was estimated from fecal samples of steers in paddocks with different botanical composition. Fecal samples were collected weekly from five individual yearling steers that were kept in three different herds of 22-26 individuals each. Each herd was allocated to a paddock, of 130 ha, either dominated by Loamy or Sandy Plains, and the third was represented by a mix of both ecological sites ( Figure 1). These differences in the ecological type provided a range of variation in the abundance of C3 (Sandy Plains) and C4 (Loamy Plains) grass species. During the 2015 and 2016 grazing seasons (mid-May through October) samples from the individual steers were composited by pasture, frozen, and shipped to the Grazingland Animal Nutrition laboratory (https://cnrit.tamu.edu/ index.php/frequently-asked-questions/ accessed on 10 January 2022) for analyses of CP content using Near Infrared Reflectance Spectroscopy (NIRS). This resulted in 115 different CP content field estimations: three paddocks × two grazing seasons, and 19-20 weeks per grazing season.
MODIS reflectance data were obtained for the same paddocks where the fecal samples were collected. Reflectance from seven spectral bands was obtained from pixels that occupied 80% or more of the analyzed paddocks ( Figure 1). Information from the MOD09A1 synthetic product was downloaded through the Google Earth Engine platform [50]. Specifically, 8-day reflectance data were obtained for seven spectral bands, along with associated metrics of the quality of the spectral measurements. The quality of satellite information was assessed according to the MODLAND QA bits (all the bands presented good quality). To match field estimates of CP content, the 8-day MOD09A1 data was converted to a daily step by linear interpolation between two successive dates.
A two-step approach was used to evaluate the association between field estimates of CP content with the seven spectral bands. First, we tested linear models between field estimates of CP content with all possible combinations of the seven spectral bands, through multiple regression analyses and simple regressions of the spectral bands combined as the normalized difference between them (band 1 − band 2/band 1 + band 2). Second, the predictive accuracy for the best model was assessed by a bootstrapping procedure (70% of the data used to calibrate and 30% to evaluate) using 50 iterations.
An important aspect of this calibration approach is the use of fecal samples as ground truth for CP content. Within remote sensing there are, at least, two main mechanisms for calibrations. The first approach considers sampling within the specific area of a pixel [51,52]. The second considers sampling within different areas of a paddock and, then, through data fusion, combines different satellites with contrasting spatial resolutions [20]. Here, we used the second approach: we combined fecal samples from each entire paddock and related them with remote sensing data that fully covered that paddock ( Figure 1). Through this approach there is a great advantage associated with the ability of cattle to sense the entire pasture. At the same time, this sensing ability may be blurred due to stocking density differences. In addition, topography differences may interfere in the signal captured by satellite not associated with CP content variations. The first source of uncertainty is tackled using contrasting stocking rates. The second source, the topography effects on satellite signal, will be part of the noise in our model. good quality). To match field estimates of CP content, the 8-day MOD09A1 data was converted to a daily step by linear interpolation between two successive dates. A two-step approach was used to evaluate the association between field estimates of CP content with the seven spectral bands. First, we tested linear models between field estimates of CP content with all possible combinations of the seven spectral bands, through multiple regression analyses and simple regressions of the spectral bands combined as the normalized difference between them (band 1 − band 2/band 1 + band 2). Second, the predictive accuracy for the best model was assessed by a bootstrapping procedure (70% of the data used to calibrate and 30% to evaluate) using 50 iterations.
An important aspect of this calibration approach is the use of fecal samples as ground truth for CP content. Within remote sensing there are, at least, two main mechanisms for calibrations. The first approach considers sampling within the specific area of a pixel [51,52]. The second considers sampling within different areas of a paddock and, then, through data fusion, combines different satellites with contrasting spatial resolutions [20]. Here, we used the second approach: we combined fecal samples from each entire paddock and related them with remote sensing data that fully covered that paddock ( Figure 1). Through this approach there is a great advantage associated with the ability of cattle to sense the entire pasture. At the same time, this sensing ability may be blurred due to stocking density differences. In addition, topography differences may interfere in the signal captured by satellite not associated with CP content variations. The first source of uncertainty is tackled using contrasting stocking rates. The second source, the topography effects on satellite signal, will be part of the noise in our model.

CP Content Risk Assessment
At CPER, three paddocks were assigned to the long-term (since 1939) grazing intensity treatments ( Figure 1). All paddocks were dominated by the Loamy Plains ecological sites, and targeted 20%, 40%, and 60% utilization of peak growing-season biomass for light, moderate, and heavy, respectively [54]. Across the three paddocks, C4 grass biomass increases and C3 graminoid biomass decreases as grazing intensity increases [55]. Plant community composition is similar between the heavy and moderate grazing intensity treatments but differs under the light grazing intensity treatment [56]. For these three

CP Content Risk Assessment
At CPER, three paddocks were assigned to the long-term (since 1939) grazing intensity treatments ( Figure 1). All paddocks were dominated by the Loamy Plains ecological sites, and targeted 20%, 40%, and 60% utilization of peak growing-season biomass for light, moderate, and heavy, respectively [54]. Across the three paddocks, C4 grass biomass increases and C3 graminoid biomass decreases as grazing intensity increases [55]. Plant community composition is similar between the heavy and moderate grazing intensity treatments but differs under the light grazing intensity treatment [56]. For these three paddocks we used the calibrated model to assess the risk levels of CP content across the early growing (April) and grazing season (mid-May to October).
To understand the risks to which ranchers can be exposed during the grazing season, we first described seasonal dynamics of the CP content for each grazing treatment across the 22 years (2000-2021) of observations from MODIS data. Second, risk assessment by grazing intensity level and observed date by MODIS was estimated as the ratio between Remote Sens. 2022, 14, 854 5 of 12 the years where the CP was ≤8% over the total 22 observed years. Third, we quantified the number of consecutive days within each of the 22 grazing seasons where CP content was >8%. Through this quantification we were able to extract three features from each grazing season: the initial and final date not limited by the 8% CP content threshold, and from the difference between these two dates we estimated the total number of days not limited by the CP content. Finally, we used trend analyses to determine if the initial and/or final date, and the total number of days not limited by CP content, changed in the last 22 years. To do so, we performed a linear regression analysis where the dependent variable was either the initial or final date, expressed as day of the year (DOY), and the independent variable was time expressed in years. paddocks we used the calibrated model to assess the risk levels of CP content across the early growing (April) and grazing season (mid-May to October). To understand the risks to which ranchers can be exposed during the grazing season, we first described seasonal dynamics of the CP content for each grazing treatment across the 22 years (2000-2021) of observations from MODIS data. Second, risk assessment by grazing intensity level and observed date by MODIS was estimated as the ratio between the years where the CP was ≤8% over the total 22 observed years. Third, we quantified the number of consecutive days within each of the 22 grazing seasons where CP content was >8%. Through this quantification we were able to extract three features from each grazing season: the initial and final date not limited by the 8% CP content threshold, and from the difference between these two dates we estimated the total number of days not limited by the CP content. Finally, we used trend analyses to determine if the initial and/or final date, and the total number of days not limited by CP content, changed in the last 22 years. To do so, we performed a linear regression analysis where the dependent variable was either the initial or final date, expressed as day of the year (DOY), and the independent variable was time expressed in years.

Crude Protein Content Calibration with Remotely Sensed Data
Over the two evaluated approaches, multiple regression or normalized difference indices, the first prevailed ( Figure 3 and Table S1 Supplementary Material). The most parsimonious model included bands one and four of the MOD09A1 synthetic product, representing the red and green portions of the spectrum, respectively: CP% = 10.46 − 0.025b1 + 0.026b4. These two bands explained 69% of spatial and temporal variation within and across grazing seasons with RMSD of 1.54 ± 0.15 (n = 115) (Figure 3). No differences were observed in the most parsimonious model for paddocks dominated by different proportions of C 3 and C 4 plant functional types (Supplementary Material Table S1). dices, the first prevailed ( Figure 3 and Table S1 Supplementary Material). The most parsimonious model included bands one and four of the MOD09A1 synthetic product, representing the red and green portions of the spectrum, respectively: CP% = 10.46 − 0.025b1 + 0.026b4. These two bands explained 69% of spatial and temporal variation within and across grazing seasons with RMSD of 1.54 ± 0.15 (n = 115) (Figure 3). No differences were observed in the most parsimonious model for paddocks dominated by different proportions of C3 and C4 plant functional types (Supplementary Material Table S1).

Figure 3.
Relation between CP content (%) derived from remotely sensed and the observed CP content. Remotely sensed data were obtained from the MODIS-Terra satellite through the MOD09A1 synthetic product using the green (b4) and red (b1) portions of the spectrum. The continuous line represents the fitted model between the observed and predicted values. The dashed line at 8% CP content represents the threshold for individual intake.

Risk Assessment of the Lack of Crude Protein Content
Clear patterns within the grazing season were observed across the 22 years for CP content, with maximum values occurring at the beginning of June and then declining as the grazing season progressed (Figure 4, top panel). The mean value for CP content across the grazing season was 9.6%, with a relative inter-annual variation of 21%. We did observe that the 8% CP threshold could occur at any time during the grazing season.

Risk Assessment of the Lack of Crude Protein Content
Clear patterns within the grazing season were observed across the 22 years for CP content, with maximum values occurring at the beginning of June and then declining as the grazing season progressed (Figure 4, top panel). The mean value for CP content across the grazing season was 9.6%, with a relative inter-annual variation of 21%. We did observe that the 8% CP threshold could occur at any time during the grazing season.
The risk of CP content being ≤8% was inverse to the seasonal pattern for CP content. The risk of being below the threshold was approximately one in five years (20%) at the beginning of June, but increased to three in four years (75%) by early September and 100% of the years for the three grazing intensity treatments by October (Figure 4, lower panel). Few differences were observed for risk among the grazing intensity treatments, suggesting that the observed patterns were primarily driven by plant phenology. Finally, during the two major droughts experienced in 2002 and 2012, the 8% threshold was not surpassed during the entire grazing season. The risk of CP content being ≤8% was inverse to the seasonal pattern for CP content. The risk of being below the threshold was approximately one in five years (20%) at the beginning of June, but increased to three in four years (75%) by early September and 100% of the years for the three grazing intensity treatments by October (Figure 4, lower panel). Few differences were observed for risk among the grazing intensity treatments, suggesting that the observed patterns were primarily driven by plant phenology. Finally, during the two major droughts experienced in 2002 and 2012, the 8% threshold was not surpassed during the entire grazing season.
The 8% CP content threshold at the end of the season was reached increasingly earlier across the 2000-2021 period. If ranchers had to end the grazing season based on that date,

Discussion
Our study provides the first complete description of the association between field estimates of crude protein content (CP content) and its remotely sensed estimation using spectral bands partially independent of NDVI for one of the main semi-arid rangelands in North America, the shortgrass steppe. We determined that CP content was positively associated with the green portion of spectrum and negatively with the red portion from MOD09A1. A calibrated model based on this association allowed us to assess for the first time the long-term risk that CP content of forage will fall below the 8% threshold during the traditional grazing season, and thus expose ranchers to significant declines in weight gains of yearling steers. Surprisingly, we found that CP content is expected to be a limiting factor for yearling growth in four of five years (80%) during about 60% of the mid-May to October grazing season. We also observed a weak negative temporal trend for the last day in the grazing season when content remains above the 8% threshold, suggesting that the number of grazing season days not limited by CP content decreased over the past 22 years (2000-2021).
To date, we lacked a full picture of the association between CP content and portions of the spectrum independent of NDVI for long-term remote sensing data. Moreover,

Discussion
Our study provides the first complete description of the association between field estimates of crude protein content (CP content) and its remotely sensed estimation using spectral bands partially independent of NDVI for one of the main semi-arid rangelands in North America, the shortgrass steppe. We determined that CP content was positively associated with the green portion of spectrum and negatively with the red portion from MOD09A1. A calibrated model based on this association allowed us to assess for the first time the long-term risk that CP content of forage will fall below the 8% threshold during the traditional grazing season, and thus expose ranchers to significant declines in weight gains of yearling steers. Surprisingly, we found that CP content is expected to be a limiting factor for yearling growth in four of five years (80%) during about 60% of the mid-May to October grazing season. We also observed a weak negative temporal trend for the last day in the grazing season when content remains above the 8% threshold, suggesting that the number of grazing season days not limited by CP content decreased over the past 22 years (2000-2021).
To date, we lacked a full picture of the association between CP content and portions of the spectrum independent of NDVI for long-term remote sensing data. Moreover, NDVI has shown limited sensitivity to moderate to high chlorophyll content or high quantities of leaf layers under erectophile canopies because, under such a canopy arrangement, more radiation is scattered to the bottom leaves, thus decreasing NIR reflectance [57]. Under low vegetation cover, soil reflectance also contributes to NIR reflectance. Specifically, soil moisture causes decreases in NIR reflectance, affecting NDVI signal in sites where green vegetation cover is far from reaching 95% interception. Finally, NIR reflectance is species dependent because it is closely associated with leaf and canopy structure [58]. Our calibrated model provides new insight as previous works focused on the association between NDVI and CP content [38,39,59]. However, the phenology derived from NDVI also explained a high variation in CP content [57]. Three key aspects from this work should advance additional efforts. First, we confirmed that CP content is closely negatively related to the red portion of spectrum, indicating a direct association with photosynthetic activity [59][60][61]. Second, we found a positive association between CP content and the green portion of the spectrum. This represents a seasonal and independent estimation of CP content from either growth rate of biomass [20,57] or the total amount of biomass [62]. These latter efforts use strong associations with the red and different wavebands within the infrared portion of spectrum for their remote sensing approaches. Finally, our model seems to be independent of the botanical composition and the growth rate, providing a robust model to extrapolate across different grazing intensity levels, as recently shown in savannas of northern Australia [40].
Our findings provide novel insights to understand intra-and inter-seasonal variations in forage quality in semi-arid environments, and linkages to assessing risks of grazing strategies. In the shortgrass steppe of the western part of the North American Great Plains, forage production is closely associated with current [63,64] and previous year [65] spring precipitation. Climate changes in this system predict increases in spring precipitation variability, increasing temperatures, and rising atmospheric CO 2 [66,67]. Increasing precipitation variability during the spring suggests enterprise changes to ranching, with more yearlings added and reduced numbers of cow-calf pairs [12]. This anticipated increase in the number of yearlings grazing in this rangeland ecosystem, combined with expected increases in forage production but with reduced CP content arising from increasing atmospheric CO 2 [19], showcases the need for robust, remotely sensed estimations of CP content for ranchers. Obtaining CP content estimates at temporal and spatial scales relevant for paddock and ranch-level decision making will be needed to determine when grazing animals may need supplementation of protein or when to remove them from extensive rangelands to confined finishing operations. Being proactive for both scenarios should be economically advantageous for ranchers and contribute to long-term resilience of this rangeland ecosystem.

Conclusions
Ranchers and land managers have access to remotely sensed biomass estimates that can be used for paddock and ranch-level grazing management (e.g., https://grasscast.unl.edu/ accessed on 10 January 2022). However, lacking are corresponding forage quality estimates that would assist in predicting near-real-time trade-offs between forage quantity and quality that impact the growth rates of grazing animals. For example, can remote sensing provide robust and accurate estimates of CP content independently of growth rate estimations associated directly with NDVI [20,60], and thereby enhance opportunities for ranchers to adaptively manage grazing lands to optimize the weight gain of yearling animals? The model predictions developed here can be used with emerging technological advances in the virtual herding of animals to match high-quality forage locations more effectively on the landscape with growing animals, without needing the botanical composition to estimate the growth rate [20]. Combining these technologies may reduce the intrinsic risks of yearling grazing forages that do not provide a diet with sufficient CP for growth, and increase optimization of weight gains across the grazing season to enhance the efficiency of the yearling grazing enterprise for ranchers in this highly variable semi-arid rangeland environment.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/rs14040854/s1, Table S1: R 2 values of the linear relationship between field estimations of CP content with all possible normalized difference indices. Highlighted are those that represented the top association and similar fitted values.