Effects of Fire and Large Herbivores on Canopy Nitrogen in a Tallgrass Prairie

This study analyzed the spatial heterogeneity of grassland canopy nitrogen in a tallgrass prairie with different treatments of fire and ungulate grazing (long-term bison grazing vs. recent cattle grazing). Variogram analysis was applied to continuous remotely sensed canopy nitrogen images to examine the spatial variability in grassland canopies. Heterogeneity metrics (e.g., the interspersion/juxtaposition index) were calculated from the categorical canopy nitrogen maps and compared among fire and grazing treatments. Results showed that watersheds burned within one year had higher canopy nitrogen content and lower interspersions of high-nitrogen content patches than watersheds with longer fire intervals, suggesting an immediate and transient fire effect on grassland vegetation. In watersheds burned within one year, high-intensity grazing reduced vegetation density, but promoted grassland heterogeneity, as indicated by lower canopy nitrogen concentrations and greater interspersions of high-nitrogen content patches at the grazed sites than at the ungrazed sites. Variogram analyses across watersheds with different grazing histories showed that long-term bison grazing created greater spatial variability of canopy nitrogen than recent grazing by cattle. This comparison between bison and cattle is novel, as few field experiments have evaluated the role of grazing history in driving grassland heterogeneity. Our analyses extend previous research of effects from pyric herbivory on grassland heterogeneity by highlighting the role of grazing history in modulating the spatial and temporal distribution of aboveground nitrogen content in tallgrass prairie vegetation using a remote sensing approach. The comparison of canopy nitrogen properties and the variogram analysis of canopy nitrogen distribution provided by our study are useful for further mapping grassland canopy features and modeling grassland dynamics involving interplays among fire, large grazers, and vegetation communities.


Introduction
Fire and large herbivores have historically played indispensable roles in determining and maintaining the canopy structure and species composition of vegetation communities in grassland and savanna ecosystems [1]. Effects of periodically recurring fires on grassland vegetation cover are complex, which involve multiple above-and below-ground processes [2][3][4]. In general, fire stimulates plant growth through removing the standing dead litter aboveground and increasing the light interception in the canopy [5][6][7], while simultaneously eliminating substantial sources of foliar nitrogen discrete categorical levels, in which each level is considered a patch type. In the resulting categorical imagery, spatial heterogeneity of vegetation characteristics can be quantified by a variety of metrics, such as contagion, evenness, and patchiness [51]. The spatial heterogeneity metrics can be calculated either for a single specified patch type or for all the patch types as a whole, describing the composition and spatial configuration of patches from multiple aspects. In terms of grassland heterogeneity, the spatial distribution of patches with high levels of vegetation quality or quantity is of special interest [46]. The quality and quantity of vegetation resources are important factors that influence the foraging pattern by large herbivores at either a fine feeding-station scale or a coarse landscape scale [54,55].
American Plains Bison (Bison bison) and domestic cattle (Bos taurus) are two major ungulate grazer species in the Great Plains of North America [31,56]. Comparison between bison and cattle in the grazing effects on grassland evolution can provide conservation implications for ecosystem managers. The objectives of our study were to (i) evaluate the effects of pyric herbivory on the spatial heterogeneity of canopy nitrogen in tallgrass prairie and (ii) compare bison and cattle from the perspective of grazing history (long-term bison grazing since 1992 vs. recent cattle grazing since 2010). In our study, variogram analysis was applied to continuous remotely sensed canopy nitrogen images to examine the spatial structure and variability in grassland canopy nitrogen distribution. Spatial heterogeneity metrics were calculated from the categorical canopy nitrogen maps and compared across the different fire and ungulate grazing treatments. Investigating the spatial distribution of canopy nitrogen associated with fire and grazing activities throws light on critical questions concerning interactions between the forage nutrient pattern and resulting grazer behaviors [31,56]. In addition, comparison in canopy nitrogen distributions among experimental sites grazed by bison, grazed by cattle, and not grazed in several decades provides new information on the role of grazing history in modulating grassland resource heterogeneity.

Study Area
This study was conducted at Konza Prairie Biological Station (KPBS, Figure 1), a tallgrass prairie of 34.87 km 2 in the Flint Hills near Manhattan, Kansas, USA (39 • 05 N, 96 • 35 W). The site is dominated by a continental climate with warm, wet springs, hot summers, and dry, cold winters. Average monthly temperature ranges from −2.7 • C in January to 26.6 • C in July. Average annual total precipitation is 835 mm with 75% falling during the seasons from April to September. The vegetation at the site is dominated by C4 perennial grasses, intermingled with a minor proportion of C3 forbs and a few C3 grasses [57]. The soil conditions vary with topographic positions. In general, the lowland soils are thick, silty clay loams. The upland and hillside soils are much shallower and rockier [58]. Topography influences the plant community structure at Konza Prairie. Lowland prairie has deeper, moister soil with higher productivity [59] compared to the shallow upland soil [60][61][62].
The KPBS is divided into more than 50 watersheds, each subjected to a treatment of either fire or ungulate grazing, or both [63]. Fire treatments vary with fire frequencies (1, 2, 4, 10, and 20 years). Grazing treatments include bison, cattle, and no grazing. At KPBS, a 970 ha enclosure (9 watersheds) is permanently stocked (year-round) with bison at a density of one individual per 4.5 ha. Bison have had full access to this enclosure since 1992 [41]. After a more than 30-year absence of grazing KPBS, a cattle grazing study was initiated in 2010. Cow-calf pairs at a density of one pair per 3.2 ha graze four watersheds for a 5-month grazing season from early May to early October each year [64,65]. Twenty-one watersheds (with an area of 17.13 km 2 ) were selected for data analysis in this study given the coverage of the aerial images. All these watersheds had not been grazed since the late 1960s and had been subjected to their prescribed burning intervals since the mid-1970s [41]; therefore, we evaluated the spatial heterogeneity of 2011 canopy nitrogen maps for (i) ungrazed watersheds that have not been grazed for more than 30 years (n = 8 watersheds), (ii) bison-grazed watersheds that have been grazed year-round since 1992, after a~30-year (1960s-1992) absence of grazing (n = 9 watersheds), and (iii) cattle-grazed watersheds that have been grazed during the growing season since 2010, after a more than 30-year (1960s-2010) absence of grazing (n = 4 watersheds).
Remote Sens. 2019, 11, x 4 of 21 growing season since 2010, after a more than 30-year (1960s-2010) absence of grazing (n = 4 watersheds). Twenty-one out of in total more than 50 KPBS watersheds (labeled by the watershed names) were selected for analysis in this study, including 8 watersheds without ungulate grazing, 9 watersheds with long-term bison grazing, and 4 watersheds with short-term cattle grazing. In the watershed names, "N" indicates bison grazing treatment, "C" denotes cattle grazing treatment, and all other watersheds are not grazed by large herbivores.

Data Collection and Preprocessing
The fire history records and the digital elevation model (DEM) of the study area were provided by the National Science Foundation Long Term Ecological Research Program at KPBS [58]. Hyperspectral imagery covering the 21 watersheds was collected on four dates (26 May, 29 June, 2 August, and 26 September) spanning the 2011 growing season, using an AISA Eagle camera (Spectral Imaging Ltd., Teknologiantie 6D, Fin-90750 Oulu, Finland) mounted on a Piper Warrior aircraft operated by the Center for Advanced Land Management Information Technology (CALMIT) of the University of Nebraska-Lincoln. The spatial resolution of the aerial imagery is 2 m. The spectral resolution varied between ~3 nm for the May-August flights, and ~10 nm for the September flight, with a spectral range from 435 to 950 nm.
The fire interval for each watershed was calculated by the days between the dates of the previous burn and the image capture. Based on the DEM and the derived slope data, each watershed was divided into four topographic types, including drainage bottoms, lowlands, uplands, and hillslopes. The hillslopes were defined as the transition zones between the lowlands and the uplands with the slope value >10° determined by a natural break [66]. The gallery forest pixels on imagery were removed before analysis ( Figure 2). The continuous canopy nitrogen maps (Figures 3a-d) were retrieved from hyperspectral imaging statistically. The empirical models used for Twenty-one out of in total more than 50 KPBS watersheds (labeled by the watershed names) were selected for analysis in this study, including 8 watersheds without ungulate grazing, 9 watersheds with long-term bison grazing, and 4 watersheds with short-term cattle grazing. In the watershed names, "N" indicates bison grazing treatment, "C" denotes cattle grazing treatment, and all other watersheds are not grazed by large herbivores.

Data Collection and Preprocessing
The fire history records and the digital elevation model (DEM) of the study area were provided by the National Science Foundation Long Term Ecological Research Program at KPBS [58]. Hyperspectral imagery covering the 21 watersheds was collected on four dates (26 May, 29 June, 2 August, and 26 September) spanning the 2011 growing season, using an AISA Eagle camera (Spectral Imaging Ltd., Teknologiantie 6D, Fin-90750 Oulu, Finland) mounted on a Piper Warrior aircraft operated by the Center for Advanced Land Management Information Technology (CALMIT) of the University of Nebraska-Lincoln. The spatial resolution of the aerial imagery is 2 m. The spectral resolution varied between~3 nm for the May-August flights, and~10 nm for the September flight, with a spectral range from 435 to 950 nm.
The fire interval for each watershed was calculated by the days between the dates of the previous burn and the image capture. Based on the DEM and the derived slope data, each watershed was divided into four topographic types, including drainage bottoms, lowlands, uplands, and hillslopes. The hillslopes were defined as the transition zones between the lowlands and the uplands with the slope value >10 • determined by a natural break [66]. The gallery forest pixels on imagery were removed before analysis ( Figure 2). The continuous canopy nitrogen maps (Figure 3a-d) were retrieved from hyperspectral imaging statistically. The empirical models used for retrieving canopy nitrogen content were developed (n = 108, R 2 = 0.86) and validated (n = 54, R 2 = 0.92, RMSE = 0.24 g/m 2 ) by 162 field sample measurements collected from May to September from three bison-grazed watersheds (N1B, N4D, and N20B) with different fire frequencies (1 year, 4 years, and 20 years). These field samples characterized a variety of canopy structures and plant species compositions affected by fire and grazing treatments over seasons. The estimated canopy nitrogen content reflected the combined information of the green vegetation amount (the leaf area index) and the leaf-level nitrogen concentration, while the nitrogen in the dead materials was not included in analysis [57]. The categorical canopy nitrogen maps (Figure 3e-h) were derived from the continuous maps by classifying the numerical canopy nitrogen values into five groups from low to high levels using natural breaks in ArcGIS (Environmental Systems Research Institute, Inc., Redlands, CA, USA).
Remote Sens. 2019, 11, x 5 of 21 retrieving canopy nitrogen content were developed (n = 108, R 2 = 0.86) and validated (n = 54, R 2 = 0.92, RMSE = 0.24 g/m 2 ) by 162 field sample measurements collected from May to September from three bison-grazed watersheds (N1B, N4D, and N20B) with different fire frequencies (1 year, 4 years, and 20 years). These field samples characterized a variety of canopy structures and plant species compositions affected by fire and grazing treatments over seasons. The estimated canopy nitrogen content reflected the combined information of the green vegetation amount (the leaf area index) and the leaf-level nitrogen concentration, while the nitrogen in the dead materials was not included in analysis [57]. The categorical canopy nitrogen maps (Figures 3e-h) were derived from the continuous maps by classifying the numerical canopy nitrogen values into five groups from low to high levels using natural breaks in ArcGIS (Environmental Systems Research Institute, Inc., Redlands, CA, USA).

Watershed-Level Descriptive Statistics and Spatial Heterogeneity Metrics
Based on the continuous canopy nitrogen imagery, the mean and standard deviation of canopy nitrogen were calculated for the zones defined by the topographic types in each watershed. On the categorical canopy nitrogen maps, the spatial heterogeneity for each topographic type in a watershed was analyzed using the FRAGSTATS program [67], in which the spatial heterogeneity metrics including the contagion index, Simpson's diversity index, and the

Watershed-Level Descriptive Statistics and Spatial Heterogeneity Metrics
Based on the continuous canopy nitrogen imagery, the mean and standard deviation of canopy nitrogen were calculated for the zones defined by the topographic types in each watershed. On the categorical canopy nitrogen maps, the spatial heterogeneity for each topographic type in a watershed was analyzed using the FRAGSTATS program [67], in which the spatial heterogeneity metrics including the contagion index, Simpson's diversity index, and the interspersion/juxtaposition Remote Sens. 2019, 11, 1364 7 of 21 index (IJI) were calculated (Table 1). Contagion and Simpson's diversity indices were calculated for all the patch types with five canopy nitrogen classes from low to high levels. The IJI was calculated for the high-nitrogen patches. The calculated canopy nitrogen properties were then compared among the fire and grazing treatments across topographic positions and seasons.  To study the effects of fire separately, the data for the ungrazed watersheds were gathered. The fire treatments were classified into two levels, including the fire intervals within one year (burned in spring) and those greater than one year (unburned in spring of the current year). The differences in the canopy nitrogen properties between the two fire interval levels were analyzed using the two-sample t-test. For the combined treatments of fire and grazing, the data for different topographic divisions were analyzed separately. Factorial analysis of variance (ANOVA) was used to examine the presence of interactions of fire and ungulate grazing effects on canopy nitrogen properties. If an interaction occurred, the two-sample t-test was used to compare the differences in canopy nitrogen properties between sites with grazing and non-grazing treatments for each of the two fire interval levels in the post-hoc tests. In addition, within the watersheds burned in spring, the canopy nitrogen properties between sites with bison and cattle were compared using two-sample t-tests at a significance level of 0.05 to reveal differences in grassland spatial heterogeneity in watersheds with different grazing histories.

Variograms for Different Fire and Grazing Treatments
Variogram analysis was used to determine the spatial structure of the continuous canopy nitrogen data. Theoretically, a variogram levels off after a certain distance ( Figure 4). This distance is called the range, within which the samples are spatially autocorrelated. The value of the variogram at the range is called the sill. The nugget is the height of the variogram at the zero separation distance, which represents the measurement errors or the variation at a distance smaller than the sampling interval. The combination of the sill and range indicates the separation distance in which the spatial variability of the characteristic is noted. The range provides insight into the spatial scale of the sampled characteristic; the sill and variance proportion (variance proportion = (sill − nugget)/sill) indicate the contrast of the characteristic [52,68]. In our study, the range reflected the vegetation patch size within which the canopy nitrogen was homogeneous and varied in a small range. The sill reflected the variance of the canopy nitrogen among vegetation patches that were not spatially autocorrelated.
Variogram analysis is sensitive to the direction, whereas the topographic divisions of this study manifested highly anisotropic patterns. To reduce the effects from the topographic complexity, zones with different topographic types were analyzed separately. For each topographic division in a watershed, three patches, each enveloped by a convex hull, were randomly selected. The patch size allowed a variogram analysis at a scale of 2-50 m. All the pixels in the selected patch were included in analysis. Results of the patches with the same treatment of fire and ungulate grazing were pooled division in a watershed, three patches, each enveloped by a convex hull, were randomly selected. The patch size allowed a variogram analysis at a scale of 2-50 m. All the pixels in the selected patch were included in analysis. Results of the patches with the same treatment of fire and ungulate grazing were pooled and averaged. Ranges and sills of the variograms were examined and compared among different fire and grazing treatments [69].

Figure 4.
Variogram analysis for spatial variability. In our research, the range reflected the vegetation patch size within which the canopy nitrogen was homogeneous and varied in a small range. The sill reflected the variance of the canopy nitrogen among vegetation patches that were not spatially autocorrelated.

Fire Effects on Canopy Nitrogen in Ungrazed Watersheds
Canopy nitrogen concentrations were higher in the watersheds with shorter fire intervals ( Figure 5). This indicated that the spring fires promoted plant productivity and forage quality; the fire stimulus to plant growth was most evident in the first year following the burn. Differences associated with the topographic positions were not evident in May, June, or September ( Figures  5a,b,d), while in August, drainage bottoms supported more canopy nitrogen content than upland locations (Figure 5c). In September, there was no apparent tendency observed in the variation of canopy nitrogen related to the fire intervals or topography.

Fire Effects on Canopy Nitrogen in Ungrazed Watersheds
Canopy nitrogen concentrations were higher in the watersheds with shorter fire intervals ( Figure 5). This indicated that the spring fires promoted plant productivity and forage quality; the fire stimulus to plant growth was most evident in the first year following the burn. Differences associated with the topographic positions were not evident in May, June, or September (Figure 5a,b,d), while in August, drainage bottoms supported more canopy nitrogen content than upland locations (Figure 5c). In September, there was no apparent tendency observed in the variation of canopy nitrogen related to the fire intervals or topography. Factorial ANOVA showed that the interactive effects of fire intervals and topography were watersheds K2A, SuB, K4A, and K20A had fire intervals greater than one year. In September, watersheds K2A and K20A were not included in analysis, given that the two watersheds were burnt shortly before the flyover in a small wildfire, resulting in a fire interval of only one week.
Factorial ANOVA showed that the interactive effects of fire intervals and topography were seldom present on the canopy nitrogen properties ( Table 2). This result suggested the independence of the two factors, namely fire intervals and topography, in affecting the spatial heterogeneity of canopy nitrogen in most cases. Effects of fire intervals on the canopy nitrogen properties were then examined using the two-sample t-tests for differences between sites with the two different fire interval levels (burned vs. not burned in spring of the current year). Results showed that during May-August the watersheds burned in spring (with fire intervals ≤ 1 yr) had higher canopy nitrogen concentrations than the unburned watersheds (with fire intervals > 1 yr, see Figure 6a); the high-nitrogen patches had higher proportions ( Figure 6e) and lower levels of spatial interspersion (Figure 6f) in the burned watersheds than in the unburned watersheds. In September, differences between the burned and unburned watersheds were insignificant in the mean of canopy nitrogen concentrations (Figure 6a) and the distribution of the high-nitrogen patches (Figure 6f).
The standard deviation of canopy nitrogen content (Figure 6b) was significantly lower in the burned watersheds than the unburned watersheds in May, which indicated a relative uniformity of canopies shortly after the fire treatments. The difference in canopy nitrogen standard deviation between burned and unburned sites became insignificant as the canopy matured in June. Canopy nitrogen standard deviation in burned watersheds became significantly higher than in unburned watersheds in August-September. Differences in contagion ( Figure 6c) and Simpson's diversity (Figure 6d) between the burned and unburned sites were statistically significant in May-June. The burned sites had higher contagion and lower diversity among the patch types than the unburned sites.
One-way ANOVA across various topographic positions showed that differences in canopy nitrogen were not statistically significant during May-August, whereas differences in the mean of canopy nitrogen concentrations and the IJI of the high-nitrogen patches were present in September (Figure 7). Tukey's honest significant difference (HSD) post-hoc test for the data in September showed that the mean of canopy nitrogen concentrations in the drainage bottoms was higher than that in the uplands (difference = 0.177, p-adj = 0.023); the IJI of the high-nitrogen patches in the drainage bottoms was significantly lower than that in the uplands and hillslopes (difference between drainage bottoms and uplands = −30.789, p-adj = 0.042; difference between drainage bottoms and hillslopes = −31.166, p-adj = 0.029).
Remote Sens. 2019, 11, x 10 of 21 unburned watersheds were insignificant in the mean of canopy nitrogen concentrations (Figure 6a) and the distribution of the high-nitrogen patches (Figure 6f). The standard deviation of canopy nitrogen content (Figure 6b) was significantly lower in the burned watersheds than the unburned watersheds in May, which indicated a relative uniformity of canopies shortly after the fire treatments. The difference in canopy nitrogen standard deviation between burned and unburned sites became insignificant as the canopy matured in June. Canopy nitrogen standard deviation in burned watersheds became significantly higher than in unburned watersheds in August-September. Differences in contagion (Figure 6c) and Simpson's diversity (Figure 6d) between the burned and unburned sites were statistically significant in May-June. The burned sites had higher contagion and lower diversity among the patch types than the unburned sites.
One-way ANOVA across various topographic positions showed that differences in canopy nitrogen were not statistically significant during May-August, whereas differences in the mean of canopy nitrogen concentrations and the IJI of the high-nitrogen patches were present in September (Figure 7). Tukey's honest significant difference (HSD) post-hoc test for the data in September showed that the mean of canopy nitrogen concentrations in the drainage bottoms was higher than

Interactive Effects of Fire and Large Herbivores
Results from factorial ANOVA showed that the interactive effects of fire and ungulate grazing on the mean of canopy nitrogen concentrations and the IJI of the high-nitrogen patches were present in all of the three types of topographic positions (i.e., drainage bottoms, lowlands, and hillslopes) during May-August. In September, the interactive effects were seldom present (Table 3).

Interactive Effects of Fire and Large Herbivores
Results from factorial ANOVA showed that the interactive effects of fire and ungulate grazing on the mean of canopy nitrogen concentrations and the IJI of the high-nitrogen patches were present in all of the three types of topographic positions (i.e., drainage bottoms, lowlands, and hillslopes) during May-August. In September, the interactive effects were seldom present (Table 3). Table 3. Analysis of variance F statistics for the canopy nitrogen properties responding to the interaction of the fire interval and the ungulate grazing type across different topographic positions. The uplands were not included in analysis given the insufficient samples of non-grazing treatments with fire intervals greater than one year.  Post-hoc analysis showed that during May-August the grazed watersheds had generally lower canopy nitrogen concentrations (Figure 8a-c) and a higher IJI of the high-nitrogen patches (Figure 8e-g) than the ungrazed watersheds when the fire intervals were within one year; in September, however, lower nitrogen concentrations ( Figure 8d) and a higher IJI of the high-nitrogen patches (Figure 8h) were observed in the ungrazed watersheds. When the fire intervals were greater than one year (Figure 8i-p), the differences between grazed and ungrazed watersheds were generally insignificant during June-September. However, at the beginning of the growing season in May, the grazed watersheds had generally higher nitrogen concentrations ( Figure 8i) and a lower IJI of the high-nitrogen patches (Figure 8m) than the ungrazed watersheds. The lack of contrast in canopy nitrogen between grazed and ungrazed sites with fire intervals greater than one year suggested a less influential effect of the grazing disturbance.
Within the burned watersheds, where the grassland canopies were evidently influenced by ungulate grazing disturbances, the canopy nitrogen properties with bison and cattle grazing were compared (Figure 9). In May, bison-grazed watersheds had significantly lower canopy nitrogen concentrations (Figure 9a) and more interspersed high-nitrogen patches (Figure 9e) than the cattle-grazed watersheds, whereas the cattle-grazed watersheds were not significantly different from the ungrazed watersheds. During June-August (Figure 9b,c,f,g), bison and cattle grazing did not differentially affect canopy nitrogen content or patchiness with the exception of the uplands, where cattle grazing had a less evident influence on the grassland canopy than year-long grazing by bison. In September, there was no significant difference detected across different topographic positions.
insignificant during June-September. However, at the beginning of the growing season in May, the grazed watersheds had generally higher nitrogen concentrations ( Figure 8i) and a lower IJI of the high-nitrogen patches (Figure 8m) than the ungrazed watersheds. The lack of contrast in canopy nitrogen between grazed and ungrazed sites with fire intervals greater than one year suggested a less influential effect of the grazing disturbance.  bottom, L = lowland, S = hillslope. *The difference was statistically significant at the 95% confidence level.
Within the burned watersheds, where the grassland canopies were evidently influenced by ungulate grazing disturbances, the canopy nitrogen properties with bison and cattle grazing were compared (Figure 9). In May, bison-grazed watersheds had significantly lower canopy nitrogen concentrations ( Figure 9a) and more interspersed high-nitrogen patches (Figure 9e) than the cattle-grazed watersheds, whereas the cattle-grazed watersheds were not significantly different from the ungrazed watersheds. During June-August (Figures 9b,c,f,g), bison and cattle grazing did not differentially affect canopy nitrogen content or patchiness with the exception of the uplands, where cattle grazing had a less evident influence on the grassland canopy than year-long grazing by bison. In September, there was no significant difference detected across different topographic positions.

Spatial Structure of Canopy Nitrogen Distrubtion
Results from the variogram analysis ( Figure 10) revealed distinctive spatial structures of the canopy nitrogen across seasons. Generally, the sills were at higher levels in May (Figure 10a-d) and August (Figure 10i-l), and dropped to lower levels in June (Figure 10e-h) and September (Figure 10m-p). Differences in the variograms across topographic positions were evident. The sills were typically high in the hillslopes (Figure 10c,g,k,o) and low in the uplands (Figure 10d,h,l,p). The ranges were generally greater in the drainage bottoms (Figure 10a,e,i,m) and lowlands (Figure 10b,f,j,n) than in the hillslopes and uplands. There was a steeper, faster rising shape of the variograms for the hillslopes, whereas the variogram curves for the lowlands were gently rising at a longer distance before leveling off. This result suggested a greater contrast in the canopy nitrogen and a smaller patch size at the hillslope sites than at the lowland sites. The disparity in terms of different fire and ungulate grazing treatments was most evident at the lowland sites in May (Figure 10b), where the sill and range were relatively high in the unburned watersheds (with fire intervals >1 yr) compared to that in the burned watersheds (with fire intervals ≤ 1 yr).
Remote Sens. 2019, 11, x 14 of 21 range were determined ( Table 4). The sills and variance proportions were noticeably high for the bison-grazing treatments across seasons and topographic positions. This indicated that more heterogeneous canopy occurred in bison-grazed watersheds, leading to greater canopy nitrogen variance than the cattle-grazed and ungrazed watersheds. The sills and variance proportions for the cattle-grazed sites were greater than that for the ungrazed sites in May (Figures 10b,d and Table 4). During June-August, the difference associated with topographic positions was observed. In the lowlands (Figures 10f,j and Table 4), the cattle-grazing sites had lower sills and variance proportions than the ungrazed sites, whereas in the uplands (Figures 10h,l and Table 4), the cattle-grazed sites had higher sills and variance proportions than the ungrazed sites. In September, the difference between cattle-grazing and non-grazing treatments was more evident in the lowland than in the uplands (Figures 10n,p and Table 4). It was found that the ranges for the lowlands (Figures 10b,f,j and Table 4) were greater than those for the uplands (Figures 10d,h,l and Table 4) across different grazing treatments during May-August. The ranges for the grazed sites were generally lower than that for the ungrazed sites across For the lowlands (Figure 10b,f,j,n) and uplands (Figure 10d,h,l,p) burned within one year, the variograms were approximated using spherical models [68], through which the nugget, sill, and range were determined ( Table 4). The sills and variance proportions were noticeably high for the bison-grazing treatments across seasons and topographic positions. This indicated that more heterogeneous canopy occurred in bison-grazed watersheds, leading to greater canopy nitrogen variance than the cattle-grazed and ungrazed watersheds. Table 4. Variogram analysis of the canopy nitrogen in lowlands and uplands burnt within one year, where the effects of non-grazing, cattle, and bison on the spatial structures of canopy nitrogen were compared. The sills and variance proportions for the cattle-grazed sites were greater than that for the ungrazed sites in May (Figure 10b,d and Table 4). During June-August, the difference associated with topographic positions was observed. In the lowlands (Figure 10f,j and Table 4), the cattle-grazing sites had lower sills and variance proportions than the ungrazed sites, whereas in the uplands (Figure 10h,l and Table 4), the cattle-grazed sites had higher sills and variance proportions than the ungrazed sites. In September, the difference between cattle-grazing and non-grazing treatments was more evident in the lowland than in the uplands (Figure 10n,p and Table 4).

Topography
It was found that the ranges for the lowlands (Figure 10b,f,j and Table 4) were greater than those for the uplands (Figure 10d,h,l and Table 4) across different grazing treatments during May-August. The ranges for the grazed sites were generally lower than that for the ungrazed sites across topographic positions during June-August (Table 4).

Fire, Ungulate Grazing, and Vegetation Canopy
In this study, we found that during the growing season from May to August, the watersheds burned in spring apparently had high canopy nitrogen concentrations of herbaceous vegetation and large proportions of high-nitrogen patches. At the start of the growing season, significantly higher contagion and lower diversity indices at the sites with shorter fire intervals suggested that the canopy nitrogen patches were more aggregated with uneven proportional distributions in the more recently burned watersheds. It can be implied that (i) the spring fires had an immediate effect on the plant species composition in the canopy, which decreased the plant species diversity dramatically and led to large uniform patches; (ii) the fire stimulus on the plant growth was most evident within the first year following the burn and in decay shortly thereafter as the plant species diversity increased and the system approached equilibrium [15,53]. As for the role of fire as a component in the nitrogen cycling of a grassland system, a large amount of nitrogen within grassland is in fact removed by fire and released in the atmosphere [70,71]. Although burning can volatilize nitrogen [72], no significant declines in total soil nitrogen stocks have been reported from annual burning at Konza Prairie [10]. The greater canopy nitrogen concentrations of herbaceous vegetation in recently burned watersheds than in unburned watersheds observed in our study are a reflection of the fire stimulus to plant growth, which is transient as the fire stimulus weakens.
In addition, our results from the ungrazed watersheds suggested that the influences of topography on the vegetation canopies were most evident in the senescent season and less so in the early growing season. The timing of senescence in the drainage bottoms was potentially later than that in the uplands due to the different conditions of soil moisture, depth, and temperature found in small-scale plot-based studies [73]. This landscape-level result extends our understanding of topographic heterogeneity in providing vegetation resources at key times of year when high-quality forage resources are limited [74,75].
Fire and grazing are historically important determinants of the structure and functioning of tallgrass prairie, and are widely used practices to manage grasslands [35]. For the first time, we highlight the role of pyric herbivory and seasonality in driving the landscape ecology of aboveground nitrogen distribution in tallgrass prairie vegetation. The interactive effects of fire and ungulate grazing on the canopy nitrogen were observed in May-August, suggesting that grazing effects on the vegetation canopy were dependent on the fire interval in the growing season. The ungulate grazing preference for the new, burned patches in the growing season, referred to as pyric herbivory, has been observed repeatedly in previous studies [29,36]. As the vegetation canopy was senescent overall in September, the interactive effects of fire and ungulate grazing were weakened. It can be implied that the ungulate forage was more evenly distributed and may shift to the unburned areas in the dormant season [75,76].
In the burned watersheds, the canopy nitrogen at the ungrazed sites was significantly higher in May-August but lower in September, compared with that at the grazed sites. The contrast between the growing season and the dormant season suggested that the effects of ungulate grazing on the vegetation canopy varied across seasons. In the watersheds with fire intervals greater than one year, the grazed sites had significantly higher canopy nitrogen than the ungrazed sites in May, contrary to the situation with shorter fire intervals. This result may be associated with grazing intensity. In the early growing season, the higher grazing intensity concentrated at the sites with shorter fire intervals reduced the canopy density, whereas the low level of grazing disturbances at the sites with longer fire intervals improved the vegetation production [22,25].
The role of fire and grazing at different intensities on aboveground canopy nitrogen distribution has not been illustrated through examination of experimental watersheds with long-term disturbance treatments [77]. Here, comparison between watersheds with recent grazing by cattle and watersheds with persistent, long-term grazing by bison showed that the sites with higher grazing pressure (year-round bison grazing from 1992 to 2011) had more interspersed high-nitrogen patches than the sites with more recent (since 2010) and short-term (May-October) cattle grazing in May across different topographic positions. This result was consistent with the variogram analysis, which indicated that the bison-grazed sites had greater grassland resource heterogeneity than the cattle-grazed sites across topographic positions and across the growing season. The differences in the vegetation canopies between sites with bison and cattle may be associated with differences in the time each ungulate species has been present at our study site. Year-round bison grazing for two decades has likely resulted in greater patchiness of high-quality vegetation on the landscape (particularly in upland portions) than the more recent growing-season cattle grazing, which had similar levels of high nitrogen patchiness to ungrazed watersheds (Figure 9e,f).
Variogram analysis and comparison in canopy nitrogen properties across the grazing treatments also showed evidence for the role of ungulate grazing history in relation to topography and grassland phenology. The variogram analysis of the canopy nitrogen distribution at the sites with recent cattle grazing indicated that the cattle preferentially selected the uplands in June-August and potentially shifted to the lowlands in September. Comparison in the canopy nitrogen content and interspersion of the high-nitrogen patches between sites with long-term bison grazing and recent cattle grazing indicated that greater grazing pressure by bison had an even more evident influence on the canopy nitrogen distribution. These findings extended the previous studies of interactions between large herbivores and grassland resource heterogeneity across topographic positions and seasons [45,56,78] from the perspective of the grazing histories. Theoretical work has identified grazing history as a driver of spatial heterogeneity in grassland vegetation [79], but few field experiments have been attempted [80,81]. Here, using a remote sensing approach, for the first time, we quantify canopy nitrogen heterogeneity at the landscape-scale. We feel this approach has much potential in elucidating ecological processes that underpin aboveground grassland structure and function including the interactive effects of fire and ungulate grazing, which have been shown to be driven by canopy nitrogen content [29]. Furthermore, like other grazed landscapes, where compositional shifts in plant communities from one state to another can take decades [82,83], our findings suggest that high landscape-level patchiness of high foliar quality may also require long-term sustained grazing to occur.

Use of Remote Sensing Imagery
Because collecting appropriate field data is so difficult, research on the development of grazing and vegetation patterns over time has previously relied on simulation models [84]. In this study, we demonstrate that remote sensing provides quantitative and continuous data of the landscape over an extensive area, which aid in our understanding of grassland pattern and process relationships. Results of the variogram analyses from numeric remote sensing imagery at a spatial scale of 2-50 m in our study were comparable with field measurements at a coarse scale of 5-30 m taken by Augustine and Frank [45], who concluded that the influences of ungulate grazers on plant communities were associated with a topographic gradient. The analysis of spatial heterogeneity in our study was limited by the 2 m spatial resolution of the remote sensor.
Analysis of spatial heterogeneity metrics from the categorical maps in our study suggested that the interspersion/juxtaposition index of high-nitrogen content patches can be used to reveal differences in the vegetation resource distribution affected by fire and ungulate grazing factors. Our findings extended the previous study from field measurements by Wallace, Turner, Romme, Oneill, and Wu [46], which compared the distributions of high biomass points at a scale of 2-30 m between burned and unburned sites.

Conclusions
Our study analyzed the spatial heterogeneity of grassland canopy nitrogen concentrations over the growing season in relation to the effects from fire, ungulate grazing, topography, and phenology at a coarse watershed scale from remote sensing imagery. Results showed an apparent and transient fire stimulus to grassland plant growth, which resulted in large uniform herbaceous vegetation patches across the growing season. The interactive effects of fire and ungulate herbivores were observed in the watersheds burned in spring, where fire-promoted herbivory reduced the canopy density and increased interspersions of high-nitrogen vegetation patches. Moreover, ungulate grazing effects on the grassland heterogeneity varied with topographic positions and site grazer history. The differences in the spatial distribution of canopy nitrogen concentrations observed across topographic positions and grazing treatments suggested preference for the uplands by the ungulate grazers in the growing season. The differences observed between sites with long-term bison grazing and recent cattle grazing indicated that high intensity bison grazing for decades apparently promoted the patchiness of grassland canopies, which revealed the role of grazer history in establishing grassland resource heterogeneity.
The results of remote sensing imagery in our study are comparable to some analyses from field measurements [45,46]. However, our approach extends previous research on the effects of fire [1,15] and large herbivores [45,56] on grassland heterogeneity across topographic positions [85] by highlighting the role of grazing history in modulating the landscape-level distribution of herbaceous plant quality. The comparison of canopy nitrogen properties and the variogram analysis of canopy nitrogen distribution provided by our study are useful for further mapping grassland canopy features and modeling grassland dynamics involving interplays among fire, large grazers, and vegetation communities.