Increasing atmospheric CO₂ is changing the dynamics of tropical savanna vegetation. C₃ trees and grasses are known to experience CO₂ fertilization, whereas responses to CO₂ by C₄ grasses are more ambiguous. Here, we sample stable carbon isotope trends in herbarium collections of South African C₄ and C₃ grasses to reconstruct ¹³C discrimination. We found that C₃ grasses showed no trends in ¹³C discrimination over the past century but that C₄ grasses increased their ¹³C discrimination through time, especially since 1950. These changes were most strongly linked to changes in atmospheric CO₂ rather than to trends in rainfall climatology or temperature. Combined with previously published evidence that grass biomass has increased in C₄-dominated savannas, these trends suggest that increasing water use efficiency due to CO₂ fertilization may be changing C₄ plant-water relations. CO₂ fertilization of C₄ grasses may thus be a neglected pathway for anthropogenic global change in tropical savanna ecosystems.

Specimen collection. Samples for analysis were collected from herbaria at Skukuza Biological Reference Collection in Kruger National Park and the South African National Biodiversity Institute (SANBI). In total, we sampled 344 grass specimens of four species: two C₄ (Digitaria eriantha [N = 121], Hyparrhenia hirta [N = 56]) and two C₃ species (Pentameris pallida [N = 84], Koeleria capensis [N = 73]). To minimize variation based on underlying physiology, both C₄ species were chosen to have the same C₄ photosynthetic subtype (using a nicotimamide-adenine dinucleotide phosphate-depending malic enzyme, commonly referred to as the NADP-ME subtype), with species selected depending on herbarium specimen availability. Sample collection dates spanned 110 years (1890s-2000s), with a mode in the early to mid 1900s (Table S2).

Each herbarium specimen is associated with a quarter-degree square, commonly used to geolocate specimens. Most of our samples, especially pre-1950, only have locations accurate to the quarter-degree square. We prioritized quarter-degree squares with more than five records spanning at least 40 years. In total, our sampling covers 35 quarter-degree squares focused in the Cape Floristic Region, Drakensberg and Highveld grasslands, and Lowveld savannas of South Africa, and covering considerable variation in rainfall regimes (Figure 1).

Isotope Analyses. Specimens were sampled for approximately 10 mg of leaf tissue, avoiding tissue with obvious deterioration or surface contamination. We did not sample specimens that had been treated with chemical fumigants as part of their mounting process due to contamination risk. As an additional precaution against surface contamination, we swabbed sample material with ethanol immediately after sampling. In preparation for analysis, we ground samples to integrate carbon across leaf parts. After samples were ground, they were re-dried at 60^oC for 3+ days before weighing for analysis.

Chemical analyses of C and N concentrations, 𝛿¹³C, and 𝛿¹⁵N were done at the Yale Analytical Stable Isotope Center using standard methods on a Thermo DeltaPlus Advantage mass spectrometer with a Costech ECS 4010 Elemental Analyzer with a Conflo III interface. Precision for our isotope values was estimated by the standard deviations of measurements of known working standards. Samples were run in randomized order, and all batch quality control standards were measured to within one standard deviation or ± 0.2‰ against standard sample 𝛿¹³C and 𝛿¹⁵N.

Atmospheric 𝛿¹³C and [CO₂]. Carbon isotope discrimination (∆¹³C) can be calculated with known plant tissue and atmospheric 𝛿¹³C via Equation 2 using estimated annual atmospheric 𝛿¹³C for 𝛿_a. Estimates were derived by merging the Law Dome Ice Core record (Rubino et al., 2013), which extends from the pre-Industrial to the early 2000s, and monthly carbon concentrations measured by CSIRO (Krummel et al., 2018) at the global atmospheric watch station in the South Pole which extends from 1991-2019. The 15-year overlap between records was used to standardize across datasets, and annual estimates were interpolated with a cubic-spline fitted from 1700 to 2018 on the merged datasets. Annual estimates of atmospheric CO₂ concentrations were derived in a similar fashion using a splined dataset from the Law Dome ice-core (pre-industrial-1996) (Rubino et al., 2013) merged with data from NOAA’s South Pole Observatory (1975-2019) (NOAA ESRL GML CCGG Group 2019). 

Rainfall and Temperature. Rainfall variables were calculated from station data, via the Global Historical Climatology Network (Menne et al., 2012) and from Kruger National Park weather stations. We calculated annual wet-season rainfall (mm year^-1), the length of the wet season (days), and the magnitude of the storm depth (amount of rainfall in mm per rainfall event during the wet season) via methods described by (Liebmann et al., 2012).

Station data were available within each quarter degree square sampled, so mean rainfall metrics across years always represent local estimates. However, annual rainfall metrics were not always available in the year and cell that a sample was collected; in these cases, we took the closest available estimate within 50 km from the center of the cell when one was available and otherwise excluded the sample. We included statistical models for both data subsets: all 334 samples with mean-annual rainfall metrics and the 251 samples for which annual rainfall metrics were available. Sample subsets and sizes used for each analysis are specified alongside results in Tables S2-S5.

Temperature estimates were extracted from NOAA’s Merged Land-Ocean Surface Temperature (Martínez-Yrízar, 1995), a global raster dataset of monthly temperature anomalies estimated at a resolution of 5°. We were interested in growing-season temperature and so calculated wet-season temperature anomalies, selecting relevant rasters via annualized or averaged wet-season duration estimates as described above. Anomalies were calculated as the difference from the long-term mean growing-season temperature.