Estimating integrated measures of forage quality for herbivores by fusing optical and structural remote sensing data

Northern herbivore ranges are expanding in response to a warming climate. Forage quality also influences herbivore distributions, but less is known about the effects of climate change on plant biochemical properties. Remote sensing could enable landscape-scale estimations of forage quality, which is of interest to wildlife managers. Despite the importance of integrated forage quality metrics like digestible protein (DP) and digestible dry matter (DDM), few studies investigate remote sensing approaches to estimate these characteristics. We evaluated how well DP and DDM could be estimated using hyperspectral remote sensing and assessed whether incorporating shrub structural metrics affected by browsing would improve our ability to predict DP and DDM. We collected canopy-level spectra, destructive-vegetation samples, and flew unoccupied aerial vehicles (UAVs) in willow (Salix spp.) dominated areas in north central Alaska in July 2019. We derived vegetation canopy structural metrics from 3D point cloud data obtained from UAV imagery using structure-from-motion photogrammetry. The best performing model for DP included a spectral vegetation index (SVI) that used a red-edge and shortwave infrared band, and shrub height variability (hvar; Nagelkerke R 2 = 0.81, root mean square error RMSE = 1.42%, cross validation ρ = 0.88). DDM’s best model included a SVI with a blue and a red band, the normalized difference red-edge index, and hvar (adjusted R 2 = 0.73, RMSE = 4.16%, cross validation ρ = 0.80). Results from our study demonstrate that integrated forage quality metrics may be successfully quantified using hyperspectral remote sensing data, and that models based on those data may be improved by incorporating additional shrub structural metrics such as height variability. Modern airborne sensor platforms such as Goddard’s LiDAR, Hyperspectral & Thermal Imager provide opportunities to fuse data streams from both structural and optical data, which may enhance our ability to estimate and scale important foliar properties.


Introduction
The ranges of some northern herbivores are expanding in response to increased forage biomass resulting from a warming climate (Tape et al 2016, Zhou et al 2017. Herbivore range expansions may impact nutrient cycling (Doiron et al 2014, Zamin et al 2017, Schmitz et al 2018. For example, high levels of herbivory may accelerate the successional transition of palatable forage species such as willow (Salix spp.) to unpalatable species such as alder (Alnus spp.) or conifers (Pastor et al 1988, Kielland and Bryant 1998, Christie et al 2015. Such transitions in species composition can alter ecosystem carbon and nitrogen dynamics (Schmitz et al 2018). For instance, nutrient turnover rates may decrease cellulose, and defensive chemicals are higher in less palatable species (Pastor et al 1993). In addition to forage biomass, herbivore distributions have also been linked to variation in forage quality (Ball et al 2000, Wu et al 2019. Forage quality is influenced by the concentration of chemical constituents and has important bottom-up effects on herbivore life-history traits such as maternal body condition, pregnancy rates, and survival (Parker et al 2009). Characterizing forage quality for herbivores is complex. Fiber, crude protein, and defensive chemical concentrations all play a role in defining forage quality (Mirik et al 2005, McArt et al 2009, DeGabriel et al 2014. Despite the importance of these individual chemical constituents, previous studies have called for increased attention to integrated measures of forage quality such as digestible protein (DP) and digestible dry matter (DDM;Moore 2005, McArt et al 2009), because simpler metrics of forage quality such as crude protein do not capture the full range of characteristics that can influence palatability and fitness.
Integrated forage quality metrics such as DP are strongly influenced by the presence of tannins (Robbins et al 1987a, 1987b, Hanley et al 1992, an important class of plant secondary compounds that significantly reduce protein digestion for herbivores by binding to plant proteins (Spalinger et al 2010). DP is particularly important in northern landscapes because plant available nitrogen is often limited (Sponseller et al 2016), which in turn limits nitrogen available to herbivores (White 1993, McArt et al 2009. Because DP estimates reflect a chemical relationship between the protein precipitating capacity of tannins and the total concentration of protein it can be used as a metric of forage quality for numerous herbivore species. However, DP does not account for the overall digestibility of forage, which varies according to herbivore species.
Body size exerts a strong influence on the digestive capabilities of herbivores (van Soest 1996, Barboza andBowyer 2000). Smaller body size facilitates selective feeding choices, whereas larger body size increases the overall fraction of digestible forage retained by herbivores (Jarman 1974, Senft et al 1987, van Soest 1996. For instance, moose (Alces alces) are largebodied ruminants whose large size enables even very poor-quality forage to be digested. DDM is an important integrated forage quality metric that quantifies the portion of plant matter that is digestible by an herbivore. DDM estimates vary depending on the herbivore in question and can be estimated in two ways. First, in vitro DDM can be estimated directly with fistulated animals, which provide an opening to the animal's stomach. DDM can also be estimated mathematically for ruminants using developed equations that account for the concentrations of DP and fiber present in plant samples (Robbins et al 1987a, Hanley et al 1992. Estimates of forage quality have classically relied on laboratory analyses, or on direct assessments from fistulated animals. However, in the past two decades remote sensing has emerged as a viable method for assessing forage quality across broad geographic extents (Mirik et  Remote sensing enables landscape-scale monitoring of forage quality, which is of great interest to wildlife managers (Walton et al 2013, Vance et al 2016. Optical remote sensing approaches use reflected light from the visible (400-700 nm) to the shortwave infrared (SWIR; 1400-2500 nm) region and have been used to detect variation in foliar chemistry. To date, much of the optical remote sensing research has focused on detecting individual components of forage quality such as crude protein, fiber, or defense chemicals like condensed tannins (Mirik et al 2005, Ferwerda et al 2006, Skidmore et al 2010, Thulin et al 2012, Jennewein et al 2020. However, Youngentob et al's (2012) pioneering work showed that DP and DDM could be successfully estimated across the landscape in Eucalyptus trees using hyperspectral remote sensing, which samples reflected light in very narrow (3-10 nm), contiguous spectral bands (Goetz 2009). More recent work demonstrated that DDM can also be estimated in grassland and pasture ecosystems using multispectral imagery obtained from unoccupied aerial vehicles (UAV) (Insua et al 2019, Michez et al 2020). Similarly, Wu et al (2019) provided the first example of successful DP estimation using the multispectral WorldView-3 satellite, which contains much broader spectral bands (30-180 nm) compared to hyperspectral data. Despite the importance of integrated forage quality metrics like DP and DDM, few studies apply remote sensing approaches to map these characteristics across the landscape and to our knowledge no studies have assessed them in Arcticboreal regions that are undergoing rapid changes due to warming (Serreze et al 2000, Verbyla 2008, Wolken et al 2011. Therefore, our first hypothesis (H1) was that integrated forage quality metrics-DP and DDMin palatable willow shrubs can be predicted using hyperspectral remote sensing. We used hyperspectral remote sensing (as opposed to multispectral) because high spectral resolution data can enable direct linkage to absorption and scattering features of foliar properties known to influence forage quality (Curran 1989, Elvidge 1990, Kokaly et al 2009. Additionally, we used hyperspectral remote sensing as previous work indicates that common multispectral vegetation indices-such as the normalized difference vegetation index (NDVI)-can have mixed results when predicting forage quality metrics (Johnson et al 2018). We focused on willow species because they are the preferred forage resource for many vertebrate herbivores in Alaska. We predicted that red-edge (680-730 nm) spectral bands would be important for both DP and DDM predictions as red-edge indices are sensitive to plant chlorophyll and nitrogen content (Eitel et al 2007, Ramoelo et al 2012, Wang et al 2012. Additionally, we predicted that the SWIR region would be important for DDM predictions as its longer spectral wavelengths are sensitive to plant fibers (  Additionally, browsing intensity can drastically alter plant canopy architecture (Christie et al 2014) and the concentrations of foliar chemical properties that influence forage quality (Bryant 1981, Bryant andChapin 1986). For instance, moderate browsing stimulates compensatory growth, which in turn creates 'bushier' shrubs that are frequently rebrowsed (Stouter 2008). In contrast, heavy browsing stunts growth, decreases shrub height, and increases canopy openness Bryant 1998, Christie et al 2014). Although plant canopy architecture can be strongly influenced by browsing intensity, ground-based assessments traditionally use only three categories to classify browsing historyunbrowsed, browsed, and broomed (figure 1). Yet variation in broomed architecture can be pronounced (figures 1(A) and (D)) and may indicate distinct functional differences such as added nitrogen from herbivore excreta Bryant 1998, Butler andKielland 2008).
Recently, remote sensing technologies such as lidar have shown utility in assessing vegetation structure for wildlife applications (e.g. Vierling et al 2008, Lone et al 2014, Melin et al 2016. For example, studies in Europe's boreal forests have demonstrated that metrics derived from aerial lidar can successfully detect insect defoliation (Solberg et al 2006, Vastaranta et al 2013 and heavy moose browsing (Melin et al 2016) on young Scots pine (Pinus sylvestris) stands. There is also great potential for fusing lidar with optical remote sensing to improve the characterization of ecosystems (Asner et al 2012, Torabzadeh et al 2014, Luo et al 2017. However, the collection of aerial lidar data may be cost-prohibitive, and 'structure from motion' (SfM) data acquired from UAVs are considered a viable alternative to aerial lidar (Wallace et al 2016).
Our second hypothesis (H2) was that statistical models for estimating DP and DDM using hyperspectral data would be improved when shrub structural metrics obtained from UAV SfM point clouds were incorporated because herbivores can drastically alter canopy structure (figure 1; Christie et al 2014). Additionally, because sunlit and shaded leaves in Alaska often differ in their concentrations of important foliar chemicals such as fiber (Molvar et al 1993) and tannins (Bryant and Chapin 1986, Klein 1990, Thompson and Barboza 2014 that influence the palatability of forage species, we also hypothesized that including the cumulative irradiance (W m −2 ) incident on a shrub in a growing season would improve models (H3). We predicted that as cumulative irradiance declined, DP and DDM would increase because nitrogen concentrations increase and fiber decreases in shaded plants (Molvar et al 1993, Lenart et al 2002. Similarly, topographic attributes such as aspect and slope influence the amount of solar radiation received by plants. Additionally, elevational gradients influence plant phenology, where higher-elevation plants often have a delayed onset of budburst, which in turn influences migrant herbivore behaviors as they move to 'surf the green wave' (Bischof et al 2012, Mysterud et al 2017). Thus, our fourth hypothesis (H4) was that including topographic attributes would improve our ability to remotely monitor forage quality because previous work demonstrates that combining hyperspectral data with topographic features such as aspect, slope, and elevation can improve model predictions of forage quality (Knox et al 2012, Pullanagari et al 2018.
Despite the importance of integrated forage quality metrics like DP and DDM, few studies to date investigated remote sensing approaches to estimate these characteristics, and none have occurred in Arctic-boreal landscapes. Thus, the overarching objective of this study was to assess the suitability of optical and structural remote sensing data to estimate DP and DDM in north central Alaska. Assessing remotely sensed measures of DP and DDM is an important first step toward mapping nutritional landscapes in high northern latitude regions subject to ongoing environmental change.

Study area
The upper Koyukuk River drainage in north central Alaska (figure 2) contains a wide range of terrain and vegetation types, including boreal forest dominated by black spruce (Picea mariana), alpine tundra and shrubs such as alders, willows, and dwarf birch (Betula glandulosa), and muskegs and other riparian vegetation such as white spruce (Picea glauca) and poplar (Populus spp.). Located in the southern end of the Brooks Range, topography is rugged and varies from 80 to 2250 m. The region experiences continental climate patterns. In winter, the average temperature ranges from −40 • C to −22 • C, with snow depths exceeding 60 cm most winters. In summer, the average temperature ranges from 3 • C to 16 • C, but can also reach >30 • C.

Field-based forage-quality assessment study
We collected vegetation spectra and destructivevegetation samples from willow shrubs along a latitudinal transect in the Koyukuk River drainage (n = 45 in July 2019; figure 2). We stratified sites across soil types, elevation, and burn history. We used targeted sampling instead of a random sample because we were limited in our ability to select shrub locations a priori as no land cover product in this region accurately differentiates between preferred forage shrubs (e.g. willows) and other shrub species (e.g. alder and dwarf birch). Additionally, we purposefully selected willows that represented the three classes of browsing intensity (i.e. unbrowsed, browsed, and broomed; figure 1) to ensure we captured a full range of possible digestibility as herbivores influence plant biochemical properties and structure (Kielland and Bryant 1998, Stouter 2008, Christie et al 2014. We collected spectral information using a Field-Spec Pro Full Range Spectroradiometer (Analytical Spectral Devices, Incorporated), which ranged from 350 to 2500 nm. This instrument has a full-width half-max spectral resolution of 3 nm in the visible and near infrared (NIR) range (i.e. 350-1050 nm), and 10-12 nm in both the NIR and short-wave infrared (SWIR) regions of the electromagnetic spectrum (i.e. 1050-2500 nm). We collected canopy-level spectra under low-cloud (<20%) conditions to minimize atmospheric interference and between 11:00 Canopy-level spectra were collected on sunexposed leaves in each of the four cardinal directions and averaged into a single representative spectral collection to eliminate canopy-level variation in nutrient distribution. Prior to sampling in each direction, white reference measures were obtained using a white reference panel (Spectralon Labsphere, Inc., North Sutton, New Hampshire, United States). Dark current-systematic noise from the instrument-was also recorded prior to collection. We then calculated spectral vegetation indices (SVIs) using all possible band combinations in the simple ratio-type vegetation index (Band A/Band B) and normalized difference-type vegetation index (Band A − Band B)/(Band A + Band B) formats and related them to calculated DP and DDM (H1).
We collected destructive vegetation samples, which we dried at 30 • C -40 • C for 12 h. Each sample was then analyzed for percent: (a) crude protein, (b) neutral detergent fiber, (c) acid detergent fiber, (d) acid detergent lignin, (e) acid insoluble ash, and (f) tannins using the CBB-BSA (2000) methodology. Integrated measures of forage quality were calculated using the Robbins et al (1987aRobbins et al ( , 1987b) equations for DP and DDM. Estimates of DDM and DP were quantified on a percent dry matter basis.

Landscape and shrub structural variables
We collected UAV data prior to destructive vegetation sampling using a DJI Phantom 4 Pro (Los Angeles, California, USA). Flight elevation ranged from 20 to 25 m above ground level with a frontal and side overlap of 80% resulting in a spatial resolution of 1 cm. To minimize atmospheric interference and minimize the effects of viewing geometry, flights occurred on sunny days between 11:00 and 15:00. Three-dimensional (3D) structural information was obtained from UAV We also interpolated DSMs (10 cm resolution) for the plots from our UAV SfM point clouds, which were used to model the light environment of the outer portion of sampled willow canopies (H3) using the 'insol' package (Corripio 2015). This package estimates the instantaneous irradiance (W m −2 ) for a given location using a DSM and atmospheric (i.e. relative humidity, ozone, visibility, air temperature) and surface reflectance (i.e. albedo) properties, which we acquired through the Env-DATA annotation service (Dodge et al 2013) for each sampled shrub location. Based on solar geometry, local topography, surface reflectance, and atmospheric properties we estimated diffuse and direct canopy irradiance for every minute of the summer, which we summed into 'total irradiance' for each shrub location. We modeled the total irradiance experienced by a shrub one week, two weeks, and one month prior to harvest. These estimated insolation values should be interpreted as the maximum possible estimation of insolation rather than a direct measure, as true insolation is likely lower than our estimates due to cloud cover.
Finally, topographic attributes including elevation, transformed aspect (TRASP) (Roberts and Cooper 1989), and slope were sourced from the Arc-ticDEM (Porter et al 2018) and included as additional landscape metrics to predict DP and DDM (H4). We also used the ArcticDEM to calculate a topographic wetness index (TWI), which uses slope and the upstream contributing area to determine topographic effects on hydrological processes. TWI has been shown to influence nitrate concentrations (Ogawa et al 2006), which directly impacts the nitrogen available for plant uptake and hence plant protein levels.

Data analyses
All data analyses were conducted in R statistical software (R Core Team 2020). The best performing SVIs were identified using Akaike information criterion for small sample sizes (AICc; Cavanaugh 1997). We also assessed the role of known spectral indices for detecting plant structure and function such as water content (the normalized difference water index (NDWI; Gao 1996)), plant pigments (chlorophyll carotenoid index (CCI) (Springer et al 2017), red-edge indices (Eitel et al 2007, Ramoelo et al 2012) and biomass (NDVI (Shippert et al 1995, Jia et al 2003 appendix A). Analysis of semivariograms for both DDM and DP indicated the potential for spatial autocorrelation in one or both of those response variables ( figure B1). Accordingly, we used generalized least squares (GLSs) regression in the 'nlme' package in R (Pinheiro et al 2017) to determine (a) whether accounting for spatial autocorrelation improved model fit, and (b) if so, which spatial correlation structure was optimal for our data (table 2). We evaluated competing models using (a) AICc, (b) root mean square error (RMSE), and (c) either adjusted R 2 (for linear models) or two pseudo R 2 values, the McFadden and Nagelkerke (for GLSs models). McFadden R 2 is often used to compare nested models (McFadden and Zarembka 1974), but values are less comparable to adjusted R 2 from linear Next, we assessed the benefit of adding shrub structure, topographic attributes, and irradiance to models by sequentially adding one metric at a time.
We also calculated AIC weights, which sum to one with values ranging from 0 to 1 and may be interpreted as the probability that a given model is the best model (Symonds and Moussalli 2011) in the set of candidate models for both DP and DDM. Additionally, we conducted leave-one-out cross validation (LOOCV) by sequentially leaving out one willow at a time and using the remaining observations to predict the excluded one. We compared LOOCV slopes (observed values plotted against predicted values) to assess bias and Spearman rank coefficients (ρ) to quantify predictive ability of these models. The residuals of the models met regression assumptions (i.e. homogeneity of variance, normality of residuals). Finally, we included all predictors that improved model fit (decreases in AICc > 2) into a single model for both DP and DDM. Predictors were only combined in the same model if collinearity between them was <0.70 (Dormann et al 2013).

Results
DP ranged from 1.57% to 13.37% of dry matter, while DDM ranged from 22.73% to 61.76% of dry matter. Most models under predicted DP and DDM (LOOCV slopes < 1; tables 3 and 4). However, we found that hyperspectral SVIs successfully predicted DP (adjusted R 2 = 0.75, RMSE = 1.42%, AICc = 166.27; appendix C) and DDM (adjusted R 2 = 0.65, RMSE = 5.11%, AICc = 330.97; appendix C), which supported our first hypothesis. The best performing SVI for DP included a rededge and a SWIR band in the normalized difference format ((R 703nm − R 1719nm )/(R 703nm + R 1719nm )). The best performing SVI for DDM included a blue and a red band in the simple ratio format (R 483nm /R 657nm ). Many viable simple ratio SVIs were found for both DP and DDM (figure 4). Of the irradiance options Table 3. Results from DP GLS models included a Gaussian correlation structure. The SVI was regressed against percent DP. Shrub structural and landscape metrics included: median canopy height (hmed), mean canopy height (hmean), 90th (h90p) and 99th (h99p) height percentiles, coefficient of variation of height (hcv), variance of heights (hvar), kurtosis (hkur) and skew (hskew) of height distributions, median canopy height * canopy fractional cover (can_med_frac), fractional cover 1-3 m (can_fcov_1_3) and >3 (can_fcov_>3) height range, TWI, TRASP, slope, elevation, and total irradiance the week before sample harvest and were sequentially added to models and evaluated based on two pseudo R 2 s: McFadden and Nagelkerke. RMSE quantified error. We used AICc to assess model fit and AICc weights to assess the probability of a given model being the best model. We used slopes and Spearman Rank correlation coefficients (ρ) from LOOCV to assess bias and predictive ability of models. The best model for DP was bolded for easy visualization.  modeled-total irradiance experienced by a shrub one week, two weeks, and one month prior to harvest-the one-week cumulative irradiance before sample harvest produced the best model fit (not shown). Thus, we only present results from this model.

Models
Model performance for predicting DP did not improve when well-known existing spectral indices for detecting plant characteristics (e.g. NDVI) were included (appendix A). The exponential correlation structure had the lowest AICc for DP (AICc = 158.87; table 2) but LOOCV models did not converge using this structure. Thus, we used the Gaussian correlation structure (AICc = 159.24). According to the DP semivariogram (appendix B), spatial correlation for DP content continued to show spatial autocorrelation throughout the range of distances present in our sample (i.e. ∼2 km). The best performing SVI showed very strong predictive power (LOOCV ρ = 0.88, slope = 0.76) and a low error estimate (RMSE = 1.42%; table 3) without the inclusion of any structural metrics (H1). Model fit improved when some structural metrics (H2)hmed (AICc = 154.30), hmean (AICc = 154.91), h90p (AICc = 155.72), hvar (AICc = 153.92)-were added to the base model (table 3). However, these metrics were all highly correlated (r > 0.7) and thus could not be considered in a single model. Thus, we identified hvar and SVI model as the best performing model (AIC weight = 0.32; figure 5). However, this model also demonstrated a slight decrease in predictive power when compared to the base model using only the SVI (∆LOOCV ρ = − 3%).
DMD model performance improved when NDVI and NDRE were included in addition to the SVI (appendix A). However, NDVI and NDRE were highly correlated (r > 0.7) and therefore only NDRE was included as it had the lowest AICc value. No correlation structures improved model fit for DDM (table 2). Thus, we used the simple linear model (AICc = 259.10). According to DDM semivariogram (appendix B), no spatial autocorrelation was present in our samples. Using only the SVI and NDRE, the model showed good predictive power (LOOCV ρ = 0.77, slope = 0.67) with a low error estimate (RMSE  (table 4). The best performing combined model for DDM included SVI, NDRE, and hvar (table 4; figure 4; adjusted R 2 = 0.73, AIC weight = 0.25, LOOCV ρ = 0.80). Increased error estimates in DDM compared to DP may be related to the relative range of values associated with these estimates, 1.57%-13.37% for DP and 22.73%-61.76% for DDM. We observed limited evidence that including the light environment (H3) nor topographic attributes (H4) produced model improvements.

Discussion
We assessed how well integrated measures of forage quality-DP and DDM-could be predicted using a fusion of hyperspectral SVIs, shrub structural metrics, topographic attributes, and the light environment. Results from our study demonstrated that DP and DDM could be successfully estimated using hyperspectral remote sensing (supporting our H1). Remotely-sensed estimates of DP showed a strong correlation with observed estimates and low error (LOOCV ρ = 0.88, RMSE = 1.42%). DDM estimates also were highly correlated with measured values but with higher error (LOOCV ρ = 0.77, RMSE = 4.61%). These findings were consistent with previous work in Australian Eucalyptus forests where DP (R 2 = 0.64) and DDM (R 2 = 0.78) were estimated with high accuracy using hyperspectral remote sensing approaches (Youngentob et al 2012), though our results suggested DP in northern willows is more easily predicted than DDM.
Although wavelengths used in the SVIs in this study did not correspond exactly to existing absorption features previously identified for tannins, protein, or fibers (Curran 1989, Elvidge 1990, Ferwerda et al 2006, they were within 25 nm. The SVI for DP contained one wavelength (703 nm) in the rededge portion of the spectrum, which is known to be sensitive to chlorophyll and has been used as a proxy for nitrogen content (Eitel et al 2007, Ramoelo et al 2012, Wang et al 2012. The second wavelength (1719 nm) used in the SVI for DP was directly adjacent to a known absorption feature of hemicellulose at 1720 nm (Elvide 1990), though numerous tannin absorption features can be found in the SWIR region (Ferwerda et al 2006). In contrast to DP, the wavelengths in DDM's SVI were both in the visible portion of the spectrum (483 and 657 nm) and near known spectral absorption features of chlorophyll pigments (Curran 1989, Ben-Dor et al 1997, which are often related to plant nitrogen concentrations. This was somewhat surprising because DDM estimates incorporated fiber concentrations that usually have absorption in the SWIR region (Curran 1989, Elvidge 1990. For instance, the SWIR region was shown to be sensitive to both DP and DDM in Australian Eucalyptus trees (Youngentob et al 2012). In our case, we found several viable vegetation indices for DP and DDM (figure 4), many of which also included SWIR wavelengths.
We also observed some improvement in model fit by incorporating shrub structural metrics such as hvar (H2). Model fit and AIC weights indicated that the best model for DP included hvar in addition to the SVI (∆AICc = 5.19; AIC weight = 0.32; table 3), although this addition slightly reduced predictive power (∆LOOCV = −3%). The best model for DDM included hvar and NDRE, which somewhat enhanced predictive power (∆LOOCV = +3%, AIC weight = 0.25; table 4). Herbivores strongly influence plant canopy architecture of palatable species such as willow. Previous work showed browsing influenced shrub height, canopy openness, and branching structure (Kielland and Bryant 1998, Christie et al 2014, 2015, which in turn can affect the palatability of forage species. To our knowledge, this study was the first to incorporate remotely-sensed shrub structural metrics as a proxy for browsing history in models to predict forage quality. Our results indicated that incorporating shrub structure may be an important, and often unconsidered, aspect of remotely sensed forage quality metrics. Based on these findings we suggest that future work should consider shrub structure when using remote sensing to study forage quality metrics. Variation in canopy structure increases the complexity of the three-dimensional space where photons interact (Asner 1998, Knyazikhin et al 2013; hence, SVIs calculated from reflectance spectroscopy are influenced by plant structural characteristics Table 4. Results from DDM linear models. The SVI and the NDRE was regressed against percent DDM. Shrub structural and landscape metrics included: median canopy height (hmed), mean canopy height (hmean), 90th (h90p) and 99th (h99p) height percentiles, coefficient of variation of height (hcv), variance of heights (hvar), kurtosis (hkur) and skew (hskew) of height distributions, median canopy height * canopy fractional cover (can_med_frac), fractional cover 1-3 m (can_fcov_1_3) and >3 (can_fcov_>3) height range, TWI, TRASP, slope, elevation, and total irradiance the week before sample harvest and were sequentially added to models and evaluated based adjusted R 2 , RMSE, AICc, and AICc weights to assess the probability of a given model being the best model. We used slopes and Spearman Rank correlation coefficients (ρ) from LOOCV to assess bias and predictive ability of models. The best models for DDM were bolded for easy visualization.  (Chen and Cihlar 1996, Turner et al 1999, Eitel et al 2008. This, coupled with the impacts of browsing, suggests a growing need to incorporate structure into remotely sensed models of forage quality to better estimate and map these metrics across the landscape. Modern airborne sensor platforms such as Goddard's LiDAR, Hyperspectral & Thermal Imager (G-LiHT; 1 m pixels, with 6 lidar pulses per m 2 ), the Airborne Visible-Infrared Imaging Spectrometer -Next Generation (AVIRIS-NG; 0.3-4.0 m pixels) and the land, vegetation, and ice sensor (LVIS; 5-25 m spots) provide opportunities to fuse data streams from both structural and optical data, which may enhance our ability to estimate and scale important foliar properties such as DP and DDM. Strategies for estimating spatiotemporal variation in forage quality metrics are needed because northern ecosystems are rapidly changing. The range of northern herbivores is expanding as the quantity of forage resources increases (Tape et al 2016, Zhou et al 2017. However, the impact of climate warming on forage quality is less clear and will likely vary depending on region and species (Lenart et al 2002, Hansen et al 2006, Turunen et al 2009, Elmendorf et al 2012. Since forage quality strongly influences herbivore lifehistory traits like maternal body condition, pregnancy rates, and survival (Parker et al 2009), monitoring 'nutritional landscapes' (sensu Merems et al 2020) that include integrated metrics of forage quality-like DP and DDM-is urgently needed. In addition to the importance of forage quality on fitness, secondary effects related to changes in herbivore populations can have cascading effects on ecosystem structure and function.

Models
We did not see any improvement in model fit or predictive power when we included the light environment (H3), which contrasted with previous work indicating that light conditions influenced fiber and nitrogen concentrations as well as DDM (Molvar et al 1993, Lenart et al 2002. This may in part be because the light modeling employed in this study did not account for cloud cover or variation in illumination throughout the canopy (i.e. we only modeled the surface foliage of the shrub canopy). Indeed, one study compared various techniques for quantifying the light environment of Salix pulchra and found that only lidar-based techniques captured photosynthetic partitioning of nitrogen and chlorophyll in canopies (Magney et al 2016).
Moreover, our findings may be related to the relatively coarse spatial scale of the atmospheric variables (with a spatial resolution of 32 km) included within the models of solar irradiance employed in this study. Future work may benefit from applying the approaches used in Magney et al (2016) or ground-based sensors that estimate the instantaneous irradiance at the location of shrubs to determine how solar energy influences DP and DDM. Similarly, including topographic attributes produced no model improvements (H4). This is surprising considering that elevation, slope and aspect have been shown to influence forage quality (Knox et al 2012, Pullanagari et al 2018.
One limitation of our study was that we did not quantify shrub biomass, because we did not have the resources locally to do so. Yet, the total energy acquired through browsing is a function of both forage quantity and quality, and thus both of these metrics are needed to truly map nutritional landscapes. Additionally, forage quality is normally the highest after budburst and then declines throughout the remainder of the growing season as fiber content increases and digestibility and nitrogen decrease (Klein 1990, McArt et al 2009, Shively et al 2019. Our study demonstrated that remote sensing of DP and DDM is possible during peak biomass. However, future work should consider the seasonal dynamics of these integrated metrics. Finally, future work should investigate the possibility of scaling these plot-level assessments to airborne and satellite platforms. Although we advocate for future work to include both optical and structural data streams, we also observed strong relationships between hyperspectral SVIs and DP and DDM without additional shrub structural inputs (tables 3 and 4; appendix C; figure C1). Therefore, we anticipate that as new hyperspectral satellite platforms-such as the Environmental Mapping and Analysis Program (EnMAP; Guanter et al 2015) and PRecursore IperSpettrale della Missione Applicativa (PRISMA; Loizzo et al 2018)-become more widely available, our ability to monitor integrated forage quality metrics seasonally and interannually at the landscape scale will be enhanced. We also suggest that future work that employs satellite data should couple with finer scale structural remotely sensed data that would help characterize the uncertainty in canopy structure attributes as the spatial resolution will be much coarser than our shrub-level assessments. Structural data may be sourced from aerial lidar transects, which can provide estimates of plant height variability (hvar) but would not provide wallto-wall coverage. One alternative to aerial lidar in high latitude regions would be canopy height models estimated from the ArcticDEM (Meddens et al 2018) to estimate canopy hvar, which we found to be an important predictor of both DP and DDM. Assessments would be needed for aerial lidar (∼1 m grid cells), satellite imagery, and canopy height models (∼5 m grid cells) to determine if these data sources are sufficiently fine-scale to provide useful information for canopy hvar, as our UAV SfM point clouds were 1 cm spatial resolution.

Conclusion
Results from our study demonstrate that integrated forage quality metrics like DP and DDM can be successfully quantified using hyperspectral remote sensing data, and that models based on those data can be improved by incorporating shrub structural metrics. Mapping DP and DDM to create a spatially explicit representation of the nutritional landscape available to herbivores may assist in management decisions in the face of ongoing environmental change. Mapping nutritional landscapes will become a more important tool for understanding wildlife ecology in the rapidly changing Arctic.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors. Table A1. Results from plant structure and function spectral index analysis. Each spectral index was added to the best performing index identified in this study.

Model
Adjusted R