Inter-Annual Climate Differences Supersede Grazing Effects on the Anatomy and Physiology of a Dominant Grass Species.

Grassland ecosystems are historically shaped by climate, re, and grazing as essential ecological drivers. These grassland drivers inuence morphology and productivity via physiological processes, resulting in unique water and carbon use strategies among species and populations. Leaf-level physiological responses in plants are framed by the underlying microanatomy, previously shown to reect patterns of carbon assimilation and water-use in leaf tissues. However, the magnitude to which microanatomy and physiology are impacted by grassland drivers, remains unstudied. To address this knowledge gap, we sampled from three locations along a latitudinal gradient in the mesic grassland region of the central Great Plains, USA during the 2018 and 2019 growing seasons. We measured annual biomass and forage quality at the plot level, while collecting physiological and microanatomical traits at the leaf-level in cattle grazed and ungrazed locations at each site. Leaf-level measurements were focused on the dominant grass species Andropogon gerardii (big bluestem) because of its high abundance, continental-scale distribution, and forage value. The two sampling seasons received markedly different levels of precipitation: drought conditions in 2018 and excessive early season precipitation in 2019. Ambient drought conditions negatively impacted A. gerardii physiology and drastically reduced productivity regardless of grazing. Leaf-level microanatomical traits, particularly those associated with water-use, varied within and across locations and between years. Our results highlight how trait plasticity can serve as an important tool for predicting future grassland responses to climate change and variable disturbances. Specically, climate played a stronger role than grazing in shaping above-ground processes in microanatomy and physiology.


Introduction
The Great Plains is the largest expanse of grasslands in North America, reaching from Saskatchewan through Texas (Robinson et al. 2019;Jones et al. 2020). The community composition and productivity of native grasses in the Great Plains varies as a result of the precipitation (longitudinally) and temperature (latitudinally) gradients (Teeri and Stowe 1976;Sala et al. 1988; Lura et al. 2019). The impacts of these gradients are re ected in the grassland ecotones of the Great Plains (arid to mesic) that separate regions of shortgrass, mixed-grass, and tallgrass prairies (DeLuca and Zabinski 2011; Dixon et al. 2014). Each of these prairie systems are dominated by a few grass species that account for a majority of annual production. For example, Bouteloua gracilis, commonly found throughout the shortgrass steppe, accounts for upwards of 90% of total biomass in that system (Sasaki and Lauenroth 2011;. Similarly, the northern mixed-grass prairie is dominated by B. gracilis, Heterostipa comata, and Pascopyrum smithii which accounts for the majority of the annual production (Lura et al. 2019). Andropogon gerardii can account for over 70% of annual biomass in the tallgrass prairie when ample rainfall is received (Weaver 1968;Smith and Knapp 2003). These dominant grasses thrive in their native habitats because each has evolved specialized functional traits as mechanisms of persistence within each region's disturbance regimes (Anderson 2006; ; Jardine et al. 2021). These adaptations include but are not limited to: 1) large shallow rooting systems comprised of ne roots that quickly absorb water (Nippert and Knapp 2007;Nippert et al. 2012); 2) belowground meristematic tissues ("bud banks") which provide new growth after senescence, re, and grazing (Dalgleish and Hartnett 2006; Ott and Hartnett 2015; Ott et al. 2019); and 3) specialized leaf morphology and anatomy to maximize light capture and minimize water loss to combat drought (Hameed et al. 2012; Nunes et al. 2020). While these functional traits aid our understanding as to the continued success of grass species in their respective region, less is understood about how these traits vary, and morphological and physiological adaptations offer the same advantages across sites with varying local climates, yet all within the same broad geographic region. For instance, which anatomical and physiological traits confer persistence locally (tallgrass prairie) and do these traits express the same relationship across different locations and climates (Great Plains).
Research investigating local adaptations of dominant C 4 grass species like A. gerardii have been primarily focused on assessing productivity, whole-leaf economics, or genomics/phenotypes ( importance of this research should not be overlooked due to its signi cance in feeding a growing global population, these data are collected from controlled environments with tightly controlled environments and abundant resources. Under real-world conditions, resources for native species are not always in steady supply -typically they are variable and often limiting. In addition, data for annual agricultural species does not always translate to regions like the Great Plains, which are comprised of native perennial grasses that must maintain structural integrity and investment beyond a single annual reproductive cycle (Benson et al. 2004;Benson and Hartnett 2006).
Genetic differentiation and physiological plasticity in a dominant and ubiquitous species have been used to explain responses to changes in water availability (Avolio and Smith 2013;McAllister et al. 2015). The evolved responses within populations are due to environmental differences that occur at a single location, not across regional scales (Valladares et al. 2014;). Therefore, understanding traits within populations exposed to different intra-annual climate within a single site will provide a better understanding for how a single species respond to climate variability at various locations. To our knowledge, a multi-location and multi-year investigation across a climate gradient in conjunction with grazing effects on leaf-level anatomy and physiology has not been done for a native grass species.
Here, we propose an alternate approach by investigating naturally occurring populations in their home environments but under a range of environmental conditions. This allows for an assessment of responses to climate variability within a site and then comparisons of variability across sites. Because A. gerardii has evolved within the Great Plains, it's important to consider the responses to climate variability in conjunction with other ecosystem drivers (i.e., re and grazing). For example, the effects of grazing on native grass species can be exacerbated by periodic droughts, resulting in reduced ecosystem services and productivity (Koerner and Collins 2014;Souther et al. 2020). While intraspeci c responses to grazing and climate may re ect genetic differences, the mechanisms that underlie physiological responses are found at the microanatomical level (Christin et al. 2013;Bellasio and Lundgren 2016;Guha et al. 2018). This study aims to provide a mechanistic understanding of how varying climate and grazing impacts a dominant species' physiological and microanatomical traits across a latitudinal gradient in the Great Plains. We hypothesized that: 1) due to site-level differences in climate histories across the latitudinal gradient, and contrasting growing season conditions in 2018 and 2019, there would be signi cant differences in mean and variability (measured here as the coe cient of variation) of leaf-level nutrient content, microanatomical traits, and instantaneous physiological responses across sites; 2) because microanatomical traits constrain/frame physiological responses to water availability, the existing trait relationships will show signi cant differences between years sampled due to the disparity in precipitation received; and 3) due to the stress of compensatory growth and reallocation of resources, grazing will accentuate leaf-level microanatomy and nutrient content differences between treatments and across locations.

Leaf physiology and anatomy
Gas exchange rates were measured using a Li-COR 6400XT (Li-COR Biosciences, Lincoln, NE, USA) equipped with an LED light source (maintained at 2000 µmol m -2 s -1 ), CO 2 concentration at 400 ppm, and maintained relative humidity in the chamber between 40-60%. This instrument was used between 10:00 and 14:00 CDT to collect photosynthetic rates (A n ), stomatal conductance (g s ), and transpiration rates (E) during two periods (June and August) in the 2018 and 2019 growing seasons. At each sampling period, leaves from three individual A. gerardii grasses were measured in each plot. Measurements were recorded when gas exchange levels remained stable for ~ 2 minutes. These same individual leaves were also used to determine nutrient content and microanatomical traits.
Following physiological gas exchange measurements, the previously measured leaf tissues were then clipped (~ 30 mm) and immediately placed into FAA (10 % formalin/5 % glacial acetic acid/50 % ethanol (95 % EtOH)/35 % DI water) for vacuum in ltration to analyze microanatomical traits. Leaf tissues were then cross sectioned to a 4 µm thickness with a Leica RM2135 microtome (Leica Biosystems, Newcastle, UK), stained with Safranin-O and Fast Green (Ruzin 2000), and imaged at 100X and 200X on an Olympus BH-2 compound microscope (Olympus America Inc, Melville, NY) (Fig. 1). We then quanti ed microanatomical traits by using IMAGEJ software (Rasband 1997) and the established procedure detailed by . The selected microanatomical traits included: the total crosssectional area measured (TMA), bundle sheath cell area (BS A, ), mesophyll area (MS A ), bundle sheath: mesophyll area (BS:MS), bulliform area (B A ), xylem area (X A ), and xylem reinforcement (t/b): the ratio of xylem wall thickness (t) with xylem diameter (b). The following traits were measured on an area basis (as a percentage of TMA): BS A , MS A , B A , and V A . In addition, due to the small size of minor veins in the sampled leaf tissue, xylem characteristics were restricted to the major vascular bundles.
Leaf stoichiometry and biomass Carbon (C) and nitrogen (N) content were measured on the same leaves used for gas exchange. These leaves were dried and ground for elemental composition of carbon and nitrogen per plot (protocol outlined in Connell et al. 2020). Aboveground biomass was determined by clipping herbaceous tissues in one 0.1 x 0.1 m frame per plot at the conclusion of each growing season. This biomass was sorted to exclude dead biomass (when necessary) and then dried at 60°C for 48 hours and weighed to determine dry mass.

Statistical analyses
All analyses were completed in the statistical program R V3.5.3 (R Core Team 2020). For all analyses, we evaluated homogeneity of variances by examining residuals vs tted and also examined normality using qq-plots and when necessary, a Shapiro-wilk's test. TMA and t were the only traits that required nonparametric analyses via Kruskal-Wallace test accompanied with a post hoc pairwise Wilcox test. To assess the effects of grazing and climate differences between locations, we utilized repeated measures mixed-effects model ANOVAs with separate models for each physiological, microanatomical, and nutrient trait as the response variables, and location, treatment, and year sampled as predictor variables, and plot as the random effects. Tests were performed using the "lmer" function within the "lmerTest" package (Kuznetsova et al. 2017) using R blah blah blah citation. .

Leaf-level physiological traits
Leaf-level physiological traits in A. gerardii varied by location (P < 0.005). However, grazing had no effect on gas exchange rates (P > 0.40) (Fig. 2; Table 1). E was statistically similar across the years sampled (P > 0.05; Fig. 2C; Table 1). Grasses at PRP had the highest gas exchange rates in 2018, while FHPP displayed the highest rates in 2019 (Fig. 2). A n and g s increased between 2018 and 2019 (P < 0.001), most notably at FHPP (74% and 156% respectively) and KPBS (119% and 150% respectively) ( Fig. 2A, B). In addition, there was an interaction between location and year sampled for both A n and g s (P < 0.001; Table 1). No statistically signi cant latitudinal trend was discernible in 2018 (Fig. 2B, C). Table 1 ANOVA results, reported as F-values for leaf-level physiological, microanatomical, stoichiometric traits, and biomass. Subscript text in parentheses refers to data transformation necessary to meet assumptions of normality. ˆ P < 0.10, *P < 0.05, **P < 0.01, ***P < 0.001.  Table 1). In addition, C:N ratios, MS A and biomass were the only traits that were affected by the grazing treatment, but only within KPBS in 2019 (P < 0.05; Table 1). Overall, MS A did not change between years nor among locations (P > 0.05), maintaining 40% of TMA. The ratio of bundle sheath area and mesophyll area (BS:MS) displayed signi cant effects from location, year, and their interaction (P < 0.03, P < 0.0001, P < 0.0001; Table 1). V A varied signi cantly between years (P < 0.05) but was not affected by treatment or location sampled (P > 0.05, P = 0.056; Table 1). V A at FHPP and PRP increased from 2018 to 2019; in contrast, V A at KPBS decreased (Table   S2). Tissues within V A were consistently between 12-18% of TMA (Table S2). B A did not vary across locations (P = 0.96) or treatment (P = 0.59; excluding KPBS in 2019), but signi cantly decreased from 2018 to 2019 in all locations except KPBS (P = 0.25; Fig. 3B; Table 1). In addition, TMA consisted of ~ 20-30% B A across each location, year, and treatment (Table S2). X A also increased across years sampled ( Fig. 3; P < 0.005) but remained statistically similar across locations and treatment (P = 0.36, P = 0.69 respectively; Table 1). FHPP was the only location that exhibited a large difference in X A between control and grazing treatment in both 2018 and 2019 (P < 0.02, P < 0.005 respectively), while grazing only impacted X A at KPBS in 2018 (P < 0.05; Fig. 3A; Table 1). Lastly, xylem reinforcement (t/b) followed a similar pattern to X A , resulted in signi cant decreases across years sampled (Table 1; P < 0.005), but were similar across locations and unaffected by grazing (P > 0.05; Table 1).

Stoichiometry and Productivity
Carbon and nitrogen content in A. gerardii leaves varied according to year and location, but C:N was the only stoichiometric measurement affected by the grazing treatment (P > 0.05; Table 1). Lea nitrogen content was consistently higher in 2019 than in 2018 (P < 0.0001; Table 1); Grass leaves at PRP had the highest nitrogen content, regardless of year (Table S1). In addition, C:N ratios were higher in 2018 compared to 2019 and also varied by location sampled and treatment (P < 0.05; Table 1). The C:N ratio was higher at both FHPP and KPBS relative to PRP in both years sampled, regardless of treatment (Table  S1). Aboveground biomass varied by location, year, and not surprisingly-treatment (P < 0.05; Table 1). PRP was the most productive location in both 2018 and 2019, in both grazed and control plots (Table  S1).

Trait relationships and variation
While traits did show relationships with average climate parameters (MAP and MAT) for the three sites, grazing had little effect on most traits and relationships (Fig S1, 2). However, higher temperatures were associated with lower N content and higher C:N ratios; while gas exchange data patterns were similar to leaf N patterns ( Figure S1). Trait data collected across locations, years, and treatment displayed considerable, yet statistically signi cant variation (Fig S1, 2). Speci cally, increased photosynthetic rates (A n ) were correlated with increased levels of N (P < 0.001), while increasing C:N ratios was shown to decrease A n (Fig S1; P < 0.05). In addition, individuals with increased amounts of BS A were also observed to have a positive relationship with A n , (Fig S1; P < 0.001). The mean coe cient of variation (CV) in physiological traits (A n , g s , and E) was signi cantly higher than the mean CV in microanatomical traits (Fig. 5A, B). However, water usage/storage traits (X A , t/b, and B A ), were responsible for the majority of microanatomical variation (Fig. 5C). In addition, slight changes in microanatomical CV were observed between years and locations, while physiology displayed signi cantly higher CV in 2018 than 2019 (Fig. 5). Weak relationships between functional trait CV were visible in several traits, however many of these relationships were dependent on climate, irrespective of grazing (Fig. 6).

Discussion
For the grasslands in the Great Plains, the east-west precipitation gradient results in three distinct grassland types: tallgrass prairie from Illinois to Kansas, associated with a region where rainfall amounts exceed evaporation losses, shortgrass steppe in the west limited by rainfall and growing season temperatures, and the mixed-grass prairie in the central potion as a transition between the wetter and drier across multiple years and locations with distinct climate histories (contrasting precipitation and temperature) (Fig. 1). Our data illustrates that patterns of response in this widespread species vary across locations, but perhaps most importantly, that this pattern of variation in response to wet/dry years is not uniform within a single latitude. Signi cant differences in leaf level physiology, microanatomy, stoichiometry, and biomass were observed across sites and between years in this study. The long-term climate histories of each location were responsible for shaping functional traits of local populations (Fig. 1), allowing for site-speci c leaf-level anatomy and physiology (Fig. 2, 3; Table 1; Table S1, 2) (Hoffman and Smith 2020; Bachle and Nippert 2021). In addition, decreased soil moisture availability reduces carbon assimilation, decreases nutrient uptake, and leads to reduced productivity (Lemoine et al. 2018;Jardine et al. 2021). Our data illustrate similar patterns, at the FHPP and KPBS sites, which received signi cantly less rainfall in the 2018 growing season than the subsequent year (Fig. 1). The drought conditions at both locations resulted in signi cantly reduced photosynthetic rates, stomatal conductance, and leaf nitrogen content ( Fig. 2; Table 1; Table S1). Increasing water stress decreases stomatal aperture, allowing for reduced water loss.
However, long durations of water stress can lead to carbon starvation (Lawson and Matthews 2020; Nunes et al. 2020). Similarly, reductions in X A and increased B A were also observed in 2018 (Fig. 3), re ecting changes in water-use strategies. Previous research showed that increased X A allows for greater water transport, but it also increases the likelihood of cavitation during droughts or when the water column is under high tension (Olson et al. 2020).
Intraspeci c trait variability (CV) was statistically different between years, but relatively similar across locations (Fig. 5). The greatest variation was reported for gas exchange measurements (A n , g s , E) in 2018, which were ~ 2 times higher than the following year (at both FHPP and KPBS) (Fig. 5A). While high variability may be inherent to the instantaneous nature of gas exchange measurements, the CV of physiological responses in 2019 was similar to all microanatomical traits regardless of function (Fig. 5B,  C). This decrease in physiological CV may indicate a baseline physiology and associated physiological plasticity of A. gerardii, when water is less limiting. While mean microanatomical traits varied signi cantly between 2018 and 2019, there was little change in variability (CV) across years (Fig. 5B, C).
In fact, most microanatomical variation resulted from water-speci c traits (X A , t/b, B A ) (Fig. 5C) These results emphasize how disparate climates across years can result in dissimilar relationships among traits and between traits and climate variables ( Fig. 4; Fig. 6; Figure S1, 2). A. gerardii photosynthetic rates correlated positively with increasing leaf nitrogen content (Fig. 4A) when analyzed between years. However, this seemingly tight relationship broke down when analyzing each year and treatment separately (Fig. 4B). Several mean trait relationships in physiological and microanatomical traits displayed opposing trends between 2018 and 2019 ( Figure S1, 2), including BS A against gas exchange traits (A n , g s , and E). In addition, the timing of precipitation has also been known to impact grassland productivity (Nippert et al. 2006;Craine et al. 2012), which is a result of altered microanatomy and physiology (Fay et al. 2002;Wang et al. 2016;Lemoine et al. 2018). For example, early season rainfall (coinciding with tissue development) allows for the production of larger vessel areas for greater transport potentials, while early season droughts constrain development, which results in smaller vessel areas (Mauseth 1988 Elson and Hartnett 2017). More recently however, the majority of grazing is accomplished by non-native grazers like cattle. Similar to climate variability and re, responses to grazing are typically examined at the community or ecosystem levels, while less is understood about the physiological and microanatomical mechanisms responsible for those responses (O'Keefe and Nippert 2017). However, grazing and other forms of herbivory can increase gas exchange rates in order to compensate for the loss of tissue (Pinkard et al. 2011;O'Connor et al. 2020). While this allows for greater carbon assimilation, it requires increased stomatal conductance which inherently leads to greater water loss (Bertolino et al. 2019). During drought conditions, this compensatory response of recently grazed tissues would negatively impact grass physiology, thereby decreasing carbon assimilation and future productivity (Feller 2016;Souther et al. 2020). However, gas exchange rates within grazed locations in this study were nearly identical to the control (Table 1; Fig. 2), even during the dry 2018 growing season. In addition, only three functional traits were impacted by the grazing treatment: MS A , C:N ratios, and biomass production ( Table 1). The grazing treatment at KPBS was responsible for most MS A variation, in both 2018 and 2019 (Table S2). In 2018, grazing increased C:N ratios in leaf tissues from FHPP and PRP (Table 1; Table S1). While grazing did impact functional trait variability, it was only observed during the 2018 growing season and only in physiological and water-use microanatomical trait CV (Fig. 5). The lack of treatment response may be due to several factors including: 1) stocking rates at each location may not be conducive to re ect substantive grazing pressure; 2) the experimental design may not have adequately covered/represented each site and subsequent treatment; 3) due to the evolutionary history of A. gerardii in the Great Plains, a heightened grazing intensity may be necessary to induce alternative physiological responses.
These results highlight how trait plasticity can serve as an important tool for understanding the anatomical and physiological mechanisms that facilitate wide distributions of a dominant grass species.
This research was completed during the 2018 and 2019 growing seasons which had signi cantly different water availability among years. Drought conditions in 2018 resulted in decreased gas exchange rates and subsequent biomass production, irrespective of grazing. However, increased water availability in 2019 facilitated high gas exchange rates and the doubling of aboveground biomass. In addition, there was signi cant variation in microanatomical traits across locations and between sampling years.
Together, these results indicate that there are speci c leaf construction strategies based on intra-annual climate conditions across the Great Plains. Such leaf construction strategies frame instantaneous physiological responses to climate variability, and also other grassland drivers (i.e., grazing and re). Results from this study underlie the importance of collecting multiple years of data from native species in natural environments. Our data also emphasizes the need for increased microanatomical research, as we clearly demonstrate site and climate-speci c leaf construction strategies are important for understanding and contextualizing physiological responses in a dominant grass species.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. BachleNippertSuppFinal.docx