Industrialization and agricultural intensification have increased the supplies of reactive forms of nitrogen (N) and phosphorus (P) across the globe, creating a myriad of ecological consequences. Nitrogen fixation from industrial and agricultural human activities now exceeds all reactive N creation from natural biological N fixation and lightning (Schlesinger 2009;Fowler and others 2013), and human activities have increased P deposition by 200% in the past 200 years, particularly in areas with high agricultural fertilizer application (Filippelli 2019). Increasing nutrient supply can have important effects on ecosystem carbon (C) stocks and cycling. For example, in grasslands, which cover 40% of the earth's ice-free surface and store 30% of global soil C (Foley and others 2005), increased supplies of biologically limiting nutrients (for example, N and P) might increase productivity and ecosystem C storage (Bai and Cotrufo 2022;Georgiou and others 2022;Pellegrini and others 2023). However, the strength of these effects may vary widely and depend on the types and quantities of nutrients added and on the interaction of these nutrients with important biotic (for example, plant and microbial community composition and productivity) and abiotic (for example, soil chemistry and climate) factors.
Nutrient addition may increase ecosystem C pools by increasing plant productivity and organic inputs to soils. Globally, inorganic N addition tends to stimulate aboveground plant biomass in many grasslands (LeBauer and Treseder 2008;Stevens and others 2015), and this increase can be enhanced by the addition of other nutrients, like P and potassium (K; Elser and others 2007;Fay and others 2015), although continued inputs of nutrients over time can lead to diminishing returns (Isbell and others 2013a). In soils, N alone often has weak effects on C pools, which depend on the supply of other nutrients and site characteristics, such as soil texture and climate (Crowther and others 2019;Keller and others 2021). However, plant and soil C pools represent the balance of multiple processes, and the addition of nutrients alone and in combination may affect these C pools via changes in C fluxes. Because of their differing physiological roles, N, P, and K could each, independently or jointly, limit the rate of photosynthesis and plant respiration (Marschner 1995), ultimately affecting plant net C uptake and storage (Terry and Ulrich 1973;Reich andothers 2001, 2009;Saito and others 2008;Bracken and others 2015). By measuring the degree of co-limi-tation of a range of interacting plant C fluxes and pools we gain a more mechanistic understanding of the ways in which nutrients influence pools and fluxes of C.
Nutrient-induced changes in soil C fluxes and pools are influenced by changes in soil chemistry. Nitrogen addition acidifies soil (for example, Liu and others 2014;Liu and Zhang 2019), which can decrease soil bacterial alpha diversity and shift bacterial community composition (Li and others 2016). The acidifying effect of N addition can also decrease arbuscular mycorrhizal abundance (Pan and others 2020) which may induce P-limitation of plant growth in areas with low P soils. Nitrogen addition often, but not always, leads to reduced soil microbial respiration (Riggs and others 2015;Riggs and Hobbie 2016;Du and others 2018;Gill and others 2021), which can enhance soil C storage; however, low rates of N addition can stimulate both soil (Wilcots and others 2022) and ecosystem respiration (Peng and others 2017;Wilcots and others 2022). Cations, such as magnesium (Mg 2+ ) and calcium (Ca 2+ ), are important for stabilizing soil organic matter (SOM) against microbial decomposition via cation bridging to mineral surfaces (Hobbie and others 2007), but N-induced acidification may reduce concentrations of these and other cations. The effects of N alone versus N in combination with other nutrients on rates of soil C fluxes and C storage have not been well-documented. For example, the well-known acidifying effect of N can be counterbalanced by addition of cations (for example, Seabloom and others 2021a), which increases soil C, likely through increases in productivity and SOM stabilization (Whittinghill and Hobbie 2012). Additionally, effects of elevated N supply on soil pH may change through time as soils accumulate hydrogen ions and lose base cations from multiple years of N addition and nitrate leaching (Fang and others 2012;Du and others 2018), so the effects of N may change from C storage to C loss as the negative consequences of acidification on soil C outweigh the fertilization effects of N on C gain.
Nitrogen-induced changes in plant community composition also may affect ecosystem C fluxes and pools. Increased nutrient supply can lead to species losses over time (Suding and others 2005;Clark and Tilman 2008;Bobbink and others 2010;Midolo and others 2019), shifts in plant resource allocation between above-and belowground biomass production (Chapin 1980;Chapin and others 1987;Cleland and others 2019;Wilcots and others 2022), and changes in total ecosystem C (Wedin and Tilman 1996;Crowther and others 2019). Theoretical and experimental evidence has been interpreted to indicate that increased N availability in grasslands favors fast-growing poor-N competitors, like C 3 grasses, to the detriment of N-conservative C 4 grasses and long-lived perennials (Wedin and Tilman 1996). Nutrient addition also typically results in plant diversity losses, reduced belowground-and increased aboveground-biomass allocation, and increased soil C storage (Tilman 1985;Suding and others 2005;Borer and others 2014;Eskelinen and others 2022). Plant compositional shifts that alter above-and belowground allocation can also influence C fluxes (Fornara and Tilman 2008). For example, a shift from a community composed of species with higher root/shoot ratios, like native C4 prairie grasses such as Andropogon gerardii or Schizachyrium scoparium, to one with lower root/shoot ratios, like non-native grasses like Poa pratensis or Elymus repens can reduce belowground autotrophic respiration and C inputs, which may decrease microbial respiration; communities with higher aboveground allocation may also have higher total rates of photosynthesis, which could lead to higher net carbon uptake (for example, Reich and others 2001;Craine and others 2002;Fargione and Tilman 2005;Tjoelker and others 2005;Dijkstra and others 2006). Furthermore, nutrient-stimulated changes in plant functional group dominance can increase ecosystem respiration, and thus C loss, via changes in litter quality and photosynthetic capacity (Du and others 2018). Quantifying changes in plant composition in response to single-and multiplenutrient addition could thus help predict changes in ecosystem C pools.
Research on nutrient fertilization in terrestrial ecosystems has historically emphasized the role of N as the most limiting nutrient (for example, Vitousek and Howarth 1991;LeBauer and Treseder 2008). As such, many experiments designed to alleviate plant N limitation of plant communities have either added large quantities of N such that limitation would definitely be overcome (Suding and others 2005), or added N at a range of addition rates in conjunction with addition of P, K, and other nutrients at rates to overcome potential limitation by other elements (for example, Tilman 1982Tilman , 1987)). However, work over the past two decades has emphasized the importance of co-limitation by N, P, and other nutrients across all ecosystems (Elser and others 2007;Harpole andothers 2011, 2016;Borer and others 2014;Bracken and others 2015;Fay and others 2015). Given the independent changes in cycling of reactive N and other elements due to human activities, isolating the effect of N on C dynamics from that of other elements is critical to improving our understanding of anthropogenic changes on ecosystem C balance.
In this study, we measured the response of grassland C pools and fluxes to N added alone versus in combination with other nutrients. Our two main objectives were to:
1. Determine the effects of long-term multiplenutrient addition on plant aboveground biomass, and 2. Assess the relationship between C pools (plant biomass and soil C) and C fluxes (net ecosystem exchange, gross primary productivity, and ecosystem respiration) in four focal treatments to isolate the effect of N addition on C dynamics.
To accomplish our objectives, we quantified whether limitation by different nutrients changed through time in a 13-year factorial nutrient experiment by examining trends in plant aboveground biomass and key soil abiotic factors, such as pH and base cations. We then measured responses of aboveground biomass and plant functional groups (C 3 vs C 4 grasses), C fluxes, and soil organic matter at the end of the 13-year period to disentangle the effects of N versus other potentially limiting nutrients on C flux rates and pools. Throughout the duration of the experiment, we predicted that N added alone would increase plant aboveground biomass and C uptake compared to Control, but aboveground biomass and plant C uptake in plots with N, P, and K plus other micronutrients (µ) added concurrently would increase more than in the N alone treatment (for example, Fay and others 2015), a relationship referred to as 'serial co-limitation' (as opposed to 'simultaneous co-limitation' where an effect would only be seen when adding all nutrients) (Harpole and others 2011). We tested the expectation of serial co-limitation by comparing Control versus N, PKµ versus NPKµ, and N versus NPKµ using both the time series and final year of data, whereby a change in biomass in the NPKµ and either the N or PKµ treatments would indicate serial co-limitation but a change in biomass in only the NPKµ treatment would indicate simultaneous co-limitation.
Using the final year of the experiment, we tested whether nutrient effects on soil organic matter were mediated by changes to plant biomass, ecosystem C fluxes, plant functional group composition, or soil abiotic conditions. We used structural equation modeling (SEM) to determine the relationships among soil organic matter and these mediating factors, both through time and in the final year of sampling. Given our prediction that changes in aboveground biomass would be driven by changes in plant functional group composition, we tested whether changes in C fluxes across the growing season in the final year of sampling.
All data were collected at Cedar Creek Ecosystem Science Reserve (CCESR) in East Bethel, MN, USA, a Long-Term Ecological Research site (45.4020°N, 93.1994°W, mean annual temperature=6 degrees C, mean annual precipitation=750 mm). CCESR is located on the Anoka Sand Plain, characterized by sandy soils (>88% sand) with low nutrient content (Grigal and others 1974). The experiment was established in 2007 in an old field, fifty-seven years after the field was abandoned from agriculture, where soils are Typic Udipsamments (Grigal and others 1974). Background N deposition at this site ranged from 0.7 to 0.9 g N m -2 year -1 during the experimental years (equivalent to 7-9 kg N ha -1 year -1 ; National Atmospheric Deposition Program), and P deposition was 0.05 g P m -2 year -1 (Twaroski and others 2007). The experiment was a completely randomized block design composed of five replicate blocks of eight 25-m 2 plots each with a factorial design of N, P, and Kµ addition treatments, at a rate of 10 g of N, P, or K m -2 year -1 . The treatments are as follows: Control (no nutrients added), N, P, Kµ (see below for details on micronutrient additions), NKµ, NP, PKµ, and NPKµ. Nitrogen was added to the plots as time-released urea (CH₄N₂O), which has been found to be an effective mimic of N deposition (Jia 2020), despite the majority of N deposition being composed of NH x and NO y compounds (Bebber 2021;National Atmospheric Deposition Program). Phosphorus was added as triple super phosphate (Ca (H 2 PO 4 ) 2 •H 2 O, and K was added as potassium sulfate (K 2 SO 4 ). Micronutrients (µ) were added in the plots with added K as 100 g m -2 yr -1 of a 6% Ca, 3% Mg, 12% S, 0.1% B, 1% Cu, 17% Fe, 2.5% Mn, 0.05% Mo, and 1% Zn mix in Year 1 only to prevent toxicity (see Borer and others 2013). Fertilizer was added at the end of each growing season, and herbivores (mainly whitetailed deer, Odocoileus virginiana, and pocket gophers, Geomys bursarius) were not excluded from any of the plots.
To meet Objective 1, aboveground biomass was sampled in mid-August 2007-2020 as described in Borer and others (2013). Briefly, we measured standing aboveground biomass (live and dead) within two 0.1-m 2 strips (10 cm9100 cm). All live plant biomass was clipped to ground level and sorted, dried at 60℃ to a constant mass and weighed to the nearest 0.01 g. The locations of the clipped strips were moved every year to adjacent, but previously unsampled, locations, and the data collected in 2007 were pre-treatment data.
Linking Composition, Biomass, and C Fluxes, 2020
In 2020, we measured plant community composition, growing season ecosystem carbon fluxes (NEP, GPP, and ER), and soil chemistry parameters in the Control, N, NPKµ, and PKµ treatments only to meet Objective 2, and specifically assess serial colimitation and mechanisms underlying responses to N (for example, Harpole and others 2011).
We measured plant community composition in the Control, N, NPKµ, and PKµ treatments by visually estimating the percent cover to the nearest 1% of each plant species within a permanent 1 m91 m subplot each month from April through September 2020 by two observers. Percent cover has been estimated in the same subplot of each experimental plot in August of each year since 2007; an analysis of the first ten years of plant community composition change in this experiment can be found in Wilcots and others (2021). Because cover for each species was assessed independently, total percent cover may sum to over 100%. To investigate compositional change, we summed plant cover by plant functional groups (C 3 grasses, C 4 grasses, legumes, and non-leguminous forbs (forbs, hereafter)). Although photosynthetic pathway can be a relatively coarse categorization (Donnelly and others 2023), at our site, belowground traits are correlated with photosynthetic pathway (that is, warm-season, native C4 grasses are more deeprooted, while cool-season, non-native C3 grasses are more shallow-rooted) (Craine and others 2002). Importantly, shifts in C4 and C3 grass abundance due to fertilization can strongly impact soil C and aboveground biomass (Wedin and Tilman 1996). Here, we present only C 3 and C 4 grass cover data, as these two functional groups made up over 80% of biomass in each plot; additional analyses of forb cover in these plots can be found in Wilcots and others (2021).
In 2020, we made 10 measurements of net ecosystem productivity (NEP), gross primary pro-ductivity (GPP), and ecosystem respiration (ER) between April and October in the Control, N, NPKµ, and PKµ treatments only. We measured NEP by measuring the change in CO 2 concentration in the headspace of a chamber (a 1 m 3 PVC frame covered with 6 mil clear plastic sheeting with 30 cm flaps that lay on ground) placed over a plot, using a Licor 850-3 (Licor, Lincoln, NE) connected to a Lenovo tablet running the "Flux Puppy" software (Carbone and others 2019). The chamber was sealed to the ground using two heavy chains placed on the ground flaps, and 3 small fans were used inside the chamber to ensure air mixing. Light conditions inside the chamber were measured using an Apogee MQ-200X PAR sensor (Apogee Instruments). Chamber CO 2 concentrations were measured for at least two minutes or until they stabilized (up to five minutes). To standardize NEP values across light levels, we used garden shade cloth and tarps to create a gradient of decreasing light environments within the chamber. Twominute measurements were conducted at 50% ambient light, 25% ambient light, and 0% ambient light (which yielded measures of ecosystem respiration, ER). NEP was then calculated using equation 1:
where ρ is the air density (mol air m -3 ) calculated as P/RT where P is the average atmospheric pressure (Pa), R is the ideal gas constant (8.314 m 3 Pa mol air -1 K -1 ), T is the average temperature in units of Kelvin (K), and dC/dt is the change in CO 2 -C concentration over time (mol C mol air -1 s -1 ). The volume of the cube (V) and ground area (A) sampled were 1 m 3 and 1 m 2 , respectively. From the four light measurements, we fit linear and hyperbolic curves to predict NEP at a standard level of photosynthetically active radiation (PAR= 800 µmol m -2 s -1 ). We used the R 2 values to assess the fit of each curve, and the hyperbolic estimation was always the better fit. Gross primary productivity (GPP) was estimated using equation 2 (Kohli and others 2021; Wilcots and others 2022):
Biomass-specific fluxes were calculated by dividing GPP, ER, or NEE by the total dry mass per m 2 of aboveground biomass in each plot.
To assess the effect of N and PKµ addition on soil chemistry and C stocks, we collected 20 cm deep soil cores in each plot in late July 2020. Soils at CCESR have negligible carbonates (Grigal and others 1974) so all C in these soils is organic. Soil was sieved to 2 mm, dried, ground, and weighed before being analyzed for pH (1:1 v/w in water), soil organic matter (loss on ignition for two hours at 400°C), and common exchangeable nutrients (base cations; K + , Mg 2+ , and Ca 2+ ) using the Mehlich-3 extraction method (Mehlich 1984; Salehi and others 2011) (Waypoint Analytical, Tennessee). Soil organic matter content, soil pH, and base cations were also measured at 10 cm depths using the same methods in 2008, 2012, 2016, and 2019 (see Seabloom and others 2021a). In this long-term experiment, soils are collected every 3-4 years to avoid damaging the plots and potentially confounding the effect of nutrient addition with the effect of disturbance from taking soil cores.
All statistical analyses were performed in R (v.4.1.1) and RStudio. For comparing changes in biomass (logged), SOM, and pH through time in all treatments, we used linear mixed-effect models (R package 'nlme'; Pinheiro and others 2018) with treatment and year (logged) as fixed effects, and year as a random effect (pre-treatment data were excluded from this model). In this way, the model accounts for any linear change in mean over time while also accounting for nonlinear among-year variation arising from other sources (for example, among-year variation in precipitation). We used least-square means (R package 'lsmeans') to compare pairwise treatment effects. We used planned orthogonal contrasts to assess the difference in response between Control and N, PKµ and NPKµ, and N and NPKµ in the final year of sampling. For repeated measures sampling (that is, seasonal C flux data and yearly biomass data), we used linear mixed-effect models (R package 'nlme'; Pinheiro and others 2018) with treatment (Control, PKµ, N, and NPKµ) added as a fixed effect and block and plot as nested random effects. We constructed the structural equation models using the 'piece-wiseSEM' package in R (Rosseel 2012) and assessed fit using Fisher's C to assess the impact of direct and indirect variables on aboveground biomass, soil organic matter, and C fluxes. A detailed description of the SEMs can be found in Appendix S2.
Through time, all treatments with multiple added nutrients (that is, NKµ, NP, PKµ, and NPKµ) had higher biomass than the Control (Figure 1a). Of the treatments with only one nutrient added (that is, N, P, and Kµ), only N led to more biomass than Control (Appendix S1: Table S1). The data provide evidence for serial co-limitation, with less biomass in the N treatment compared to both NP (t-ratio -5.008, p<0.0001) and NPKµ treatment (t-ratio=-4.674, p<0.0001; Figure 1b) and more biomass in the PKµ treatment compared to the P and Kµ alone treatments (t-ratio -3.128 and -5.033 and p-values 0.0391 and<0.0001, respectively; Figure 1b). To further investigate these trends, we focused our analyses of C pools and fluxes between four focal treatments: Control, N, PKµ, and NPKµ.
Using structural equation modeling (SEM), we found no effect of N alone on aboveground biomass (AGB) using all thirteen years of data from the experiment (Fisher's C=1.13, p=0.569, d.f.=2; Figure 2). In contrast, N alone consistently reduced soil pH, which in turn was negatively associated with SOM. Both PKµ and NPKµ treatments increased base cation concentration and AGB, and NPKµ reduced soil pH. Base cations were positively associated with SOM. Regressed against time, pH decreased across all treatments, but most strongly in the N alone treatment; SOM and base cations did not significantly change through time (Appendix S1: Figure S1 and Table S2).
In the final year of the study, we quantified concurrent changes in plant biomass, SOM, and C fluxes in four focal treatments: Control, N, NPKµ, and PKµ. Plant AGB in the NPKµ treatment was 84% higher than in the N alone treatment (F 1,16 = 10.66, p=0.0049, Figure 3a); in contrast to earlier in the experiment, AGB in the N alone treatment was no longer significantly different from in the Control, contrary to our expectations given past work in this area (for example, Isbell and others 2013a) (F 1,16 =0.090, p=0.77; Figure 3a). We did not find evidence of any differences between AGB the PKµ and NPKµ treatments (F 1,16 =0.118, p= 0.74; Figure 3a).
All plot-level fluxes, measured in the final year of the study (NEP, GPP, and ER), exhibited seasonal trends, with values peaking in mid-August (Figure 4a-c). However, contrary to our expectations, the treatments' flux rates did not differ at any point during the growing season (Appendix S1: Table S1). At peak biomass, biomass-specific GPP and biomass-specific ER (GPP or ER per gram of biomass) were higher in the N alone treatment compared to NPKµ (Figure 4d-f, F1,16=6.121, p= 0.025 and F 1,16 =4.579, p=0.048, respectively). Biomass-specific GPP and ER rates did not differ between the PKµ and NPKµ treatments (F 1,16 = 0.28, p=0.60 and F 1,16 =0.027, p=0.87, respec- tively), or between the Control and N alone treatments (F 1,16 =0.00, p=0.99 and F 1,16 =2.233, p= 0.15, respectively). We found no effect of treatment for biomass-specific NEP (F 3,16 =0.32, p=0.32).
The relative abundance of plant functional groups was insensitive to N alone but responded to addition of NPKµ, providing evidence of a multinutrient influence on C 3 and C 4 grasses, similar in shape to simultaneous (rather than serial) co-limitation (Figure 3b). C 3 grass cover more than doubled in the NPKµ treatment compared to the N alone treatment (F 1,16 =6.182, p=0.11). Compared to the N alone treatment, mean C 4 grass cover was reduced by nearly 87% with NPKµ addition (F 1,16 = 2.891, p=0.11), and by 90% compared to the PKµ treatment (F 1,16 =17.281, p=0.0007). Monthly NEP was best predicted by C 4 cover (r=0.34, p=0.0001; Appendix S1: Figure S2). When evaluated using an SEM, however, we did not find evidence that C 3 or C 4 grass cover influenced SOM or biomass-specific fluxes (Appendix S1: Figures S3 and 4). In fact, we found evidence that biomass-specific GPP might decline with increasing base cation abundance. Using only the final year of data, we found strong direct effects of both N and NPKµ addition on soil acidification, which may have reduced SOM, an interpretation supported by the SEM using all 13 years of data.
The long-term experimental work presented here demonstrates that biomass production in a northtemperate grassland was not limited solely by N throughout 13 years of fertilization. Rather, multiple lines of evidence indicated the importance of different forms of co-limitation of grassland biomass by N and PKµ. Using the framework laid out by Harpole and colleagues (Harpole and others 2011), we found evidence for simultaneous colimitation of N, P, and Kµ (whereby only the simultaneous addition of all three nutrient treatments had an effect) for aboveground biomass and functional group composition in the final year of treatment, but evidence of serial co-limitation of N and PKµ (where by each treatment had an effect when added alone and a positive synergistic effect when added together) when assessing biomass trends through time. Plant AGB in the N alone and PKµ treatment was greater than in the Control, and lower than in both the NP and NPKµ treatments. However, by the final treatment year, AGB in the N alone treatment was similar to that in the Control.
Similarly, C 3 cover increased, and C 4 cover decreased, in response to NPKµ, but not to N alone, similar to findings from the first ten years in this experiment (Wilcots and others 2021). These findings indicate that effect of N with P, K, and µ is important for shaping plant communities and determining production at this site, and that the effect of N alone at this site changes with continued N addition and is less pronounced than has been emphasized in previous investigations and their interpretations (Wedin and Tilman 1996;Clark and Tilman 2008).
There was significantly more plant biomass in the NPKµ treatment compared to the N alone treatment, despite similar levels of soil acidification. Past studies have found that adding N and P together can increase soil acidity which can reduce mycor- ). There were no differences between treatments for mass-specific NEP, but mass-specific GPP (e) was higher in the N alone treatment compared to the NPKµ treatment (f). Transparent points represent plot-level fluxes, and solid points represent treatment average flux rates. Error bars are ±1 SE.
rhizal abundance (Lekberg and others 2021), cation exchange capacity (for example, Pan and others 2021), and bacteria/fungi ratios in soil (Li and others 2016). If the degree of soil acidification was enough to negatively impact arbuscular mycorrhizal fungi infection or abundance, we would expect the community to demonstrate P-limitation, given that arbuscular mycorrhizal fungi primarily help plants acquire P (Smith and others 2011). Thus, the consistently greater biomass in the NPKµ treatment compared N alone treatment may have had reduced capacity for P or K uptake, or higher N:P or N:K ratios, leading to lower rates of biomass gain compared to the NPKµ that had high rates of P and K supplied annually (for example, Pen ˜uelas and others 2019). This mechanism would be consistent with recent work from this same field, which found increased plant AGB in treatments with low, not high, rates of N addition, and no change in soil pH between Control and low N addition treatments (Wilcots and others 2022), and would also explain why initial biomass gains in the N alone treatment disappeared through time in this study.
Furthermore, we found that nutrient-induced changes in soil pH and base cations impacted both SOM stocks and biomass-specific flux rates. Specifically, N alone led to acidification, and pH was negatively correlated with SOM. Both the PKµ and NPKµ treatments had direct positive effects on base cations, and the positive indirect effect of NPKµ on SOM via base cation abundance outweighed the negative effect of NPKµ on SOM via decreases in pH. Thus, previously reported increases in soil C with multiple-nutrient addition (for example, Seabloom and others 2021b) may be due to the positive influence of K and cations on SOM (Crowther and others 2019;Osburn and others 2023). Interestingly, the effects of base cations on SOM did not appear to be driven by increases in organic matter inputs to soils, at least from aboveground biomass, as the aboveground biomass-SOM linkage was weak. Thus, base cations may have promoted stabilization of SOM into aggregates or onto mineral surfaces (for example, Hobbie and others 2007). Increasing soil acidity, however, often negatively impacts soil microbial communities (Li and others 2016;Shi and others 2016) and reduces soil microbial respiration (Du and others 2018;Wilcots and others 2022), which would lead to less C loss. Soil organic matter pools were similar across all treatments despite large differences in soil pH, suggesting that other mechanisms, such as increased root respiration (Peng and others 2017), changing patterns in below-ground allocation (for example, Cleland and others 2019), or increased root turnover (Wilcots and others 2022), may lead to no net SOM gain in treatments with enhanced productivity. Despite large differences in aboveground biomass and functional group abundances induced by the treatments, we found no treatment-induced differences in ecosystem fluxes. The N alone treatment had higher rates of biomass-specific GPP compared to the NPKµ treatment. Given that the N alone treatment had more C 4 grass cover and less C 3 cover compared to NPKµ and that C 4 grasses are typically more resource-efficient (particularly Nefficient) and compared to C 3 grasses (Chapin and others 1987; Augustine and others 2011), the difference in biomass-specific GPP could be due to changes in plant composition, especially given recent evidence that N addition can reduce rates of photosynthesis in Elymus repens, a common C 3 grass at our site (Pastore and others 2020). Although we did not find significant pathways between C 4 or C 3 cover and biomass-specific GPP in our SEM, other traits of common C 3 or C 4 grasses at our study site could explain the difference in biomass-specific GPP rates between the C 3 -dominated NPKµ treatment and the C 4 -dominated N treatment. For example, a common C 3 grass at our site, Elymus repens, has been found to increase litter accumulation in plots where it is dominant, and thus litter buildup and light extinction could explain the lower rates of mass-specific GPP in the NPKµ treatment (Isbell and others 2013b).
Despite strong theoretical and experimental evidence of co-limitation in plant communities since the 1980s (for example, Chapin 1980), many studies still focus on single-nutrient limitation. Our findings reinforce the idea that the singlenutrient (often N) focus of terrestrial ecological research, first proposed by Liebig and reinforced by reviews such as LeBauer and Treseder (2008), does not always capture the complicated dynamics between plant and soil communities and added nutrients. Our study joins a growing body of studies from global grasslands that show the importance of co-limitation of plant biomass production by N, P, and other nutrients (Harpole and others 2011;Bracken and others 2015;Fay and others 2015). Our results build from this to demonstrate that nutrient effects on GPP and ER are decoupled from nutrient effects on biomasswhereas N alone did not affect C fluxes compared to NPKµ, we did find an increase in plant AGB in NPKµ compared to N alone. This mismatch is important globally, because anthropogenic nutrient supplies are decoupled from one other and are occurring against varied background soil fertility levels, so effects of nutrient inputs on ecosystem C pools and fluxes in areas with increasing N deposition will depend on P or K supply or availability. Where N inputs are rising but P, K and micronutrients remain low, systems may be buffered from overall C losses. We found evidence, at a rate of N addition that exceeds N deposition in most locations in the world, that N addition alone still may not trigger changes in functional group composition of the community. However, areas receiving multiple-nutrient pollutants, such as those with high agricultural use, may experience large shifts in plant functional group composition at elevated N input levels (Borer and Stevens 2022). Furthermore, areas receiving low amounts of N deposition may exhibit characteristics of N single-limitation if soil base cation concentration or pH remains high (Wilcots and others 2022).
Our findings provide insight into high-profile results from a neighboring field at Cedar Creek, where multiple levels of N addition, in conjunction with P, K, and other micronutrients, caused large increases in aboveground biomass (Tilman 1987;Isbell and others 2013a), shifts from C 4 to C 3 dominance (Wedin and Tilman 1996), and increases in soil C (Seabloom and others 2021b). These studies have been widely interpreted as evidence that N is the critical element limiting biomass and controlling the composition of plant communities (but see Seabloom and others 2021b). Here, by disentangling the effect of N alone versus N with other nutrients, we found strong evidence for serial limitation of N and PKµ throughout the duration of the experiment. These results indicate that while N is important for biomass accumulation and productivity, this effect may be dependent on temporal changes in soil base cation abundance and soil pH. Our findings also clarify past results regarding changes in plant community composition. As in Wedin and Tilman (1996), we found large increases in C 3 plant cover and decreases in C 4 grass cover, but only when N was added with PKµ. This demonstrates that the loss of C 4 species is not solely due to elevated N, but is driven by interactive effects between N and P, K, or µ that simultaneously co-limit these dominant functional groups.
The response of C storage and plant growth to ongoing global change depends on interconnected changes in plant and microbial community composition and soil chemistry. Here, we present evidence that co-limitation by N and PKµ of total biomass and relative abundance of grassland functional groups. Soil organic matter was not strongly influenced by changes in plant aboveground biomass but was sensitive to changes in soil pH and base cations. Past work including treatments from this experiment found no change in soil C stocks after more than a decade of N or P enrichment (Keller and others 2021;Wilcots and others 2022), but studies that included treatments with N with PKµ found increases in total soil C under many site conditions (Crowther and others 2019;Seabloom and others 2021b). Taken in context with previous studies, our results show that both biotic (shifts in plant functional groups) and abiotic (changes in soil base cation concentration and pH) effects of nutrient enrichment can be important for C accumulation (Pu ¨spo ¨k and others 2023), even in extremely low-N grassland soils, and that measuring pools and fluxes of C concurrently can give a clearer understanding of how grassland C accumulation may change with changing elemental supply rates due to anthropogenic activities.