Environmental factors function as constraints on soil nitrous oxide fluxes in bioenergy feedstock cropping systems

Nitrous oxide (N2O) is a potent greenhouse gas and major component of the net global warming potential of bioenergy feedstock cropping systems. Numerous environmental factors influence soil N2O production, making direct correlation difficult to any one factor of N2O fluxes under field conditions. We instead employed quantile regression to evaluate whether soil temperature, water‐filled pore space (WFPS), and concentrations of soil nitrate ( NO3− ) and ammonium ( NH4+ ) determined upper bounds for soil N2O flux magnitudes. We collected data over 6 years from a range of bioenergy feedstock cropping systems including no‐till grain crops, perennial warm‐season grasses, hybrid poplar, and polycultures of tallgrass prairie species each with and without nitrogen (N) addition grown at two sites. The upper bounds for soil N2O fluxes had a significant and positive correlation with all four environmental factors, although relatively large fluxes were still possible at minimal values for nearly all factors. The correlation with NH4+ was generally weaker, suggesting it is less important than NO3− in driving large fluxes. Quantile regression slopes were generally lower for unfertilized perennials than for other systems, but this may have resulted from a perpetual state of nitrogen limitation, which prevented other factors from being clear constraints. This framework suggests efforts to reduce concentrations of NO3− in the soil may be effective at reducing high‐intensity periods—”hot moments”—of N2O production.

fossil fuel displacement for many bioenergy feedstock cropping systems (Crutzen, Mosier, Smith, & Winiwarter, 2008;Robertson, Paul, & Harwood, 2000), making their management and mitigation of N 2 O production a major part of assessing net long-term environmental impact (Gelfand & Robertson, 2015). This is especially important for dedicated biomass production cropping systems such as miscanthus plantations or harvested native grass fields (Gelfand et al., ), whose appeal rests heavily on their anticipated positive environmental benefits (Landis et al., 2017;Werling et al., 2014). Many studies report markedly different N 2 O emissions levels among cropping systems, with lower emissions frequently observed in perennial, species-rich, or minimally fertilized systems (Gelfand, Shcherbak, Millar, Kravchenko, & Robertson, 2016;Niklaus, Wardle, & Tate, 2006;Oates et al., 2016;Stehfest & Bouwman, 2006). It is less clear whether these differences primarily derive from environmental conditions (e.g., higher soil N) or how agroecosystems respond to those conditions, with major implications for predicting their behavior under novel conditions.
Studies linking environmental factors directly to N 2 O production have found these factors vary substantially among regions (Dechow & Freibauer, 2011) and cropping systems . Similarly, studies such as a Bayesian recalibration of the nitrous oxide emission (NOE) module of the agroecosystem model CERES by Lehuger et al. (2009) and a meta-analysis of 14 published models by Surendran Nair et al. (2012) show responses to environmental factors differing among sites or cropping systems. Thus, while the biology of N 2 O production is well understood, predicting and modeling this process in a particular agroecological context remains challenging.
Much of the difficulty in monitoring, modeling, and managing soil N 2 O arises from "hot moments," high-intensity short-duration bursts of activity occurring when multiple environmental factors align to create near-ideal conditions for denitrification (Groffman et al., 2009). For example, fertilizer application followed by heavy rain and high temperatures creates a substrate-rich, oxygen-limited environment with high microbial activity. These events produce flux orders of magnitude greater than is typical of the system and can contribute 25%-50% of cumulative annual N 2 O emissions (Flesch et al., 2018;Saha et al., 2017;Zenone et al., 2016). This dynamic likely contributes to the high interannual variability of cumulative annual N 2 O emissions  and makes it challenging to accurately characterize the probable range of annual emissions without multiple years of measurement.
The prevalence of hot moments suggests N 2 O-generating processes may be subject to the ecological law of the minimum (Hiddink & Kaiser, 2005), wherein the rate of a process separately constrained by multiple factors is determined by the single most limiting factor. Conventional regression, which evaluates the central tendency of the process rate for a given value of a predictor (the conditional mean), would perform poorly under such conditions, as the response of each observation would be independent of all but the single most limiting factor. In contrast, quantile regression can be used to evaluate the maximum observed value for the process rate for a given value of the predictor (the conditional upper limit; Cade & Noon, 2003). If the upper limit for N 2 O flux magnitudes changes with the value of an environmental factor, it would suggest that factor serves as a constraint on N 2 O production.
We correlated individual N 2 O flux measurements with paired environmental measurements collected from seven bioenergy feedstock cropping systems at two sites over a 6 year period . These cropping systems covered multiple dimensions of agroecological intensity including perenniality, plant species richness, and nitrogen addition, allowing us to explore relationships with environmental drivers across a wide range of relevant field conditions. Using quantile regression, we evaluated two hypotheses: (a) that soil temperature, water-filled pore space, and concentrations of NO − 3 and NH + 4 all constrained the upper limits of N 2 O flux measurements, and (b) that the nature of these constraints differed among the cropping systems in our study. This work provides a framework for an alternative approach to relating N 2 O fluxes to environmental measurements in field studies.

| Experimental design and study sites
We conducted this study on the DOE Great Lakes Bioenergy Research Center's Bioenergy Cropping Systems Experiment  Table S1). The remaining treatments were in place throughout the measurement period and consisted of continuous no-till corn without a cover crop, monocultures of switchgrass (Panicum virgatum L.), miscanthus (Miscanthus × giganteus), and hybrid poplar (Populus nigra × P. maximowiczii "NM6"), a five species native grass mix, an early successional field recruited from the pre-existing seedbank, and an 18 species restored tallgrass prairie. Species and variety information is presented in Supporting Information Table S1 of Oates et al. (2016).
All treatments were planted in 27 × 43 m plots in a fivereplicate randomized complete block design (n = 5 blocks) and managed with field-scale equipment. Annual grain systems were managed following recommendations from local university extension programs. The poplar system was fertilized in 2010 (210 kg N/ha as 34-0-0 granular ammonium nitrate) and harvested by coppicing during the 2013-2014 winter. Microplots (10 × 43 m) were established in all other systems to test effects of N addition. The restored prairie microplot and main plots of all other systems received annual spring N addition (56 kg N/ha as 34-0-0 granular ammonium nitrate). The main restored prairie plot and microplots of all other systems received no added N. Aboveground biomass in these systems was harvested to 10 cm residual stubble height following the first frost event in the fall. N-addition dates for all systems are given in Supporting Information Table S2, and full details on agronomic management were presented in Sanford et al. (2016).
Soils at KBS are primarily Kalamazoo loam (USDA soil classification: Fine-Loamy, Mixed, Semiactive, Mesic Typic Hapludalfs). Mean annual temperature from 1981 to 2010 was 9.9°C, and mean annual precipitation was 1,027 mm (MSCO, 2013). Prior to BCSE establishment, the field was planted to alfalfa (Medicago sativa L.) and a corn-soybean rotation. The switchgrass, native grass mix, and restored prairie treatments at KBS suffered seed loss following flooding in 2008 and were reseeded in 2009. Soils at ARL are predominantly Plano silt loam (USDA classification: Fine-Silty, Mixed, Superactive, Mesic Typic Arguidolls). Mean annual temperature from 1981 to 2010 was 6.8°C, and mean annual precipitation was 869 mm (NWS, 2013). Pre-BCSE land use differed among blocks: the corn phase of a corn-soybean rotation (blocks A1-A3) or an alfalfa-orchardgrass (Dactylis glomerata L.) hay mixture (blocks A4-A5). We replanted miscanthus at ARL in 2010 following stand loss during the 2008/2009 winter.

| Data generation
Nitrous oxide was measured biweekly during the growing season, with additional sampling following fertilization and major precipitation events. Sampling frequency was reduced during the winter, particularly at KBS and in the early years of the study. Static chambers were used to estimate trace gas emissions, with one chamber per plot/microplot. Chambers were cylindrical (28.5 cm diameter, ~17 cm effective height, ~10 L volume) and inserted to a soil depth of ~5 cm. We ensured adequate headspace mixing by selecting an appropriately low volume:surface area ratio and by keeping vegetation trimmed within chambers (Livingston, Hutchinson, & Spartalian, 2006). Chamber lids were fitted with a septum for gas extraction and a 2 mm diameter vent tube for pressure equilibration. For each chamber, four headspace gas samples of 30 ml were collected: immediately upon chamber closure then subsequently at three ~20 min intervals. Samples were placed in glass 5.9 ml Exetainer vials (Labco Limited, Buckinghamshire, UK), using 20 ml of sample to flush the vial before overpressurizing with another 10 ml. Following gas chromatography, CO 2 concentration was detected using an infrared gas analyzer (IRGA, LiCor 820, Lincoln, NE, USA) and N 2 O concentration was detected using an electron capture detector (micro-ECD, Agilent 7890A GC System, Santa Clara, CA, USA). We avoided ECD cross-sensitivity issues by using an argon-methane carrier gas and setting the detector temperature to 350°C (Wang, Wang, & Ling, 2010).
Prior to estimation of N 2 O fluxes, CO 2 accumulation curves were visually inspected for outlier samples indicating compromised vial integrity or other mechanical errors, removing these outliers. In time series with four valid measurements, nonlinearity of fluxes was evaluated using the HMR package (v0.3.1, Pedersen, 2015) in the R statistical environment (v3.5.0, R Core Team, 2018). Following this classification, time series received a secondary visual inspection focused on identifying outlier samples in N 2 O concentrations, particularly those that might drive a nonlinear fit. Nonlinear flux estimates from the HMR() function were used for samples that passed this secondary inspection without any data removal and whose nonlinear estimate was outside the 95% confidence interval for the linear flux estimate. All other samples used linear flux estimates. All flux observations were used for analysis.
From 2010 onward, soil cores (3.7 cm diameter, 15 cm depth at ARL, 25 cm depth at KBS) were collected concurrently with trace gas sampling. Inorganic soil N was extracted from a 10 g field-moist subsample using 2 M KCl following Robertson, Sollins, Ellis, and Lajtha (1999). Ammonium (NH + 4 ) and nitrate (NO − 3 ) concentrations were determined using a Flow Solution 3100 segmented flow injection analyzer (OI Analytical, College Station, TX, USA), using USEPA methods 27200110 and 27190110, respectively. The instrument has a detection limit of 0.05 μg N/g soil, which we used as a floor for concentrations.
Soil temperature was measured at the time of trace gas sampling using a 15 cm temperature probe (Checktemp 1C, Hanna Instruments, Smithfield, RI, USA). Soil moisture was measured at KBS by determining gravimetric water content (GWC) for the soil N samples. At ARL, moisture was measured as volumetric water content (VWC) within 1 m of the static chamber using a time domain reflectometer with 20 cm rods (FieldScout 300, Spectrum Technologies, Plainfield, IL, USA). Bulk density was measured for all plots in 2008 and 2013. We estimated annual changes in bulk density by linear interpolation, calculating mean values for groups of cropping systems (three groups: annual grain crops, poplar, and all other systems) and sets of blocks (four sets: A1 and A3, A2, A4 and A5, and all KBS blocks) with similar distributions of measurements. Water-filled pore space (WFPS) was estimated from bulk density (Bd) and soil particle density (Pd, assumed to be a constant 2.65 g/cm 3 ): For our analyses, we constrained WFPS values to be ≤100%, as deviations between local bulk density and the

| Data analysis
Nitrous oxide fluxes, NH + 4 , and NO − 3 values all varied over orders of magnitude and exhibited a strong right skew; to mitigate this, observations were inverse hyperbolic sine (IHS) transformed prior to analysis. This transformation commonly used in the social sciences to handle overdispersed variables (Burbidge, Magee, & Robb, 1988) but can also be used in the natural sciences when values near or below 0 are relevant (Sekhon et al., ). This allowed us to include negative N 2 O fluxes, which are periodically observed (Molodovskaya et al., 2012) and to avoid amplifying measurement errors for values close to the detection limit of our instruments. Following this transformation, median values were centered between the 25th and 75th percentiles (Figure 1 and Supporting Information Figure S1).
All analyses were conducted in the R statistical environment. Graphics were generated using the ggplot2 package (Wickham, 2009). Boxplot quantiles (Figure 1, Supporting Information Figures S1 and S2) and LOESS fits (Supporting Information Figure S3) were generated using the default settings of geom_boxplot() and stat_smooth(). Quantile regression was carried out using the rq() function in the quantreg package (Koenker, 2018) using τ = 0.95 and estimating standard errors using the kernel method. Briefly, quantile regression operates similarly to regular regression but with asymmetrically weighted errors with the parameter tau indicating the quantile used to calculate weights. In our case, τ = 0.95 meant 5% of observations would be above the regression line and their errors would be weighted to match the 95% of observations below the regression. Analogous to regular regression, a slope of 0 signifies the τth quantile of the response is the same for all values of the predictor, while a significant slope indicates the maximum values observed in the response will depend on the predictor. The statistical significance of differences between nested model structures was evaluated with anova.rq(), using a Wald test. Statistical significances of differences among cropping systems slopes were evaluated by setting the annual fertilized systems (described in Results) as the baseline and assessing the significance of interaction terms.

| Evaluating the range and depth of N 2 O flux and environmental factor measurements
Our dataset consisted of 9,542 individual N 2 O flux measurements with at least one paired measurement of NH + 4 , NO − 3 , WFPS, or soil temperature, with 4,273 observations having all four factors (Table 1). Soil inorganic N data were most limited, with no data in 2009 and reduced data collection frequency from 2013 onward. We analyzed the annual cropping systems as a group (rotations and rotational phases described in Supporting Information Table S1), as they had nearly identical distributions within a site for all five measured variables (Supporting Information Figure S1).
Within cropping systems, N 2 O flux measurements varied over multiple orders of magnitude (Figure 1). Fertilized perennial systems at both sites overlapped considerably with the annual systems, although annual systems had a much higher prevalence of high fluxes. Similarly, unfertilized perennial systems tended to have lower and less variable distributions of fluxes than their fertilized counterparts. Similar patterns were visible in NO − 3 concentrations, although with a smaller range. The effects of fertilization at ARL were more pronounced, as were differences between annuals and perennials at both sites. In contrast, NH + 4 concentrations were less dynamic across cropping systems and fertilization levels, although values tended to be higher and less variable at KBS than at ARL. Cropping system differences in soil temperature and WFPS were minimal relative to within-system variability, although WFPS was much higher at ARL than at KBS (Supporting Information Figure S2). At low concentrations, NO − 3 and NH + 4 were minimally correlated, although at high levels they were somewhat correlated in annual systems (Supporting Information Figure S3).

| Environmental factors correlate to upper bounds of N 2 O fluxes
We used quantile regression to independently evaluate the correlation between the four environmental factors and the upper limit of N 2 O fluxes (Figure 2). Modeling separate slope and intercept terms for each site significantly improved fits for NH + 4 (F 2,4497 = 10.9, p < 0.05), WFPS (Wald F 2,9178 = 321.6, p < 0.05), and soil temperature (Wald F 2,9417 = 22.1, p < 0.05) but not NO − 3 (Wald F 2,4497 = 1.2, p = 0.30). In all cases, regression slopes were positive and significant at p < 0.05, signifying that the maximum observed N 2 O flux magnitude increased at higher levels of all four environmental factors. Minimal fluxes were observed at high levels of all factors, suggesting individual factors were not sufficient to drive high fluxes. There may have been insufficient observations at the upper range of NH + 4 and NO − 3 concentrations to appropriately support quantile regression. Even after IHS transformation, NH + 4 measurements had a long rightward tail. Observations with these extremely high NH + 4 values had fluxes below the upper limit predicted by quantile regression, suggesting a factor other than NH + 4 limited N 2 O production. NO − 3 was more evenly distributed across its range, providing greater support for its role as a limiting factor. The clearest difference between sites was observed with WFPS, with the coarse texture of the soil at KBS restricting the range of values that could be observed. The high concentration of soil temperature values near 0°C resulted from higher sampling intensity outside of the growing season at that site (Table 1).

| Environmental factor constraints differ broadly among cropping systems
Fitting separate quantile regression slopes and intercepts for each cropping system-fertilization combination significantly improved model performance (all Wald test p < 0.05). Based on the slope and intercept parameters by cropping system and site (Supporting Information Figures S4 and S5), we aggregated cropping systems into three groups: annual crops, fertilized perennials, and unfertilized perennials. The unfertilized systems were the only ones that did not have significantly positive slopes for all environmental factors, with a neutral slope for soil temperature at both sites and a negative one for NH + 4 at ARL (Figure 3). At ARL, perennial systems had lower slopes than annual systems in their response to NH + 4 and soil temperature. There was more inconsistency at F I G U R E 2 Quantile regression between soil N 2 O fluxes and environmental factors. Lines indicate the quantile regression at τ = 0.95, representing the 95th percentile of fluxes conditional on each environmental factor. Regressions were calculated independently for each environmental parameter. Flux data and soil inorganic N concentration data are presented on an inverse hyperbolic sine scale. All slope and intercept terms were significant at p < 0.05 KBS, with unfertilized perennials responding more strongly to NH + 4 , fertilized perennials responding more strongly to WFPS, and both responding less strongly to soil temperature. Intercepts were all positive except for the response of unfertilized perennials at KBS to NH + 4 (Supporting Information Figure S6).

| Environmental factors correlate to upper bound of soil N 2 O fluxes
Using quantile regression (Cade & Noon, 2003), we found that NO − 3 , NH + 4 , WFPS, and soil temperature were each correlated to the upper quantiles of soil N 2 O flux measurements collected from a broad range of potential bioenergy cropping systems. This supports our hypothesis that these environmental factors function as constraints on soil N 2 O flux magnitudes. None of these factors drove high fluxes by itself. The largest measured fluxes occurring when WFPS, soil temperature, and NO − 3 concentrations were all high, but minimal fluxes were observed at the highest levels of all factors. The importance of these factors fits with the consensus understanding of their importance to rates of N 2 O-generating processes (Robertson & Groffman, 2015).
Quantile regression requires a large amount of data across the dynamic range (Cade & Noon, 2003). We had the greatest number of observations and relatively even distributions of observations for soil temperature and WFPS. Most models evaluating the impacts of WFPS use more complex relationships than linear regression (Heinen, 2006); although as shown in Figure 1, this seems to be a reasonable approximation at ARL. The low WFPS values observed at KBS are due to the low water holding capacity of the site's sandy soils. Castellano et al. (2010) found that matric potential was a better predictor of N 2 O production in soils, explaining why a sandy soil could generate high N 2 O fluxes at low WFPS. The limited effect of temperature at ARL likely results from the large number of observations taken from near-frozen soils, many of which had substantial fluxes. Near-frozen soils can generate substantive N 2 O fluxes, particularly in association with freeze-thaw events (Teepe, Brumme, & Beese, 2000;Wagner-Riddle et al., 2017). The different distribution at KBS almost certainly resulted from reduced sampling outside of the growing season, as previous work at that site found sizeable wintertime N 2 O fluxes (Ruan & Robertson, 2017). Interestingly, these findings suggest the impact of soil temperature may differ during the winter, as restricting the analysis to soil temperatures >0°C increases the slope of the soil temperature constraint by 16% at KBS and 71% at ARL.
We observed a much stronger effect for NO − 3 than NH + 4 at both sites, although large N 2 O fluxes could occur at low levels of both compounds. Separate processes use the two compounds, and while denitrification of NO − 3 drives the largest N 2 O fluxes, nitrification of NH 3 /NH + 4 is still relevant under certain conditions (Gelfand & Yakir, 2008;Mathieu et al., 2006). We observed a correlation between NO − 3 and NH + 4 only when both were at high concentrations, shortly after fertilizer application, while at lower concentrations, the two were uncorrelated. Thus, in a system in which denitrification-derived N 2 O production was limited by NO − 3 , nitrification-derived N 2 O might not be similarly constrained (and vice versa), weakening the apparent importance of both substances at lower concentrations. While statistically significant, the effect of NH + 4 is certainly the weakest we observed, which is consistent with its role as a N 2 O-generating process of secondary importance.

| Limited evidence that system properties, rather than environmental conditions, drive differences in N 2 O production
We hypothesized that cropping systems would differ in their response to environmental conditions, driving the observed differences in N 2 O production. The literature provides ample reasons to expect this to be the case. Oxygen sensitivities for key denitrification enzymes differ among locations (Cavigelli & Robertson, 2001), which could shift the response to soil moisture conditions. The response to soil inorganic N concentrations also differs among systems (Lehuger et al., 2009), with some work suggesting nonagricultural systems may respond less strongly to N additions (Lu et al., 2011). Soil carbon availability and microbial biomass influence denitrification potential (Heinen, 2006), which in turn would determine how much N 2 O flux would increase with the lifting of an environmental constraint. In our study, ARL had substantially higher microbial biomass and soil carbon than KBS (Liang et al., 2012), suggesting it should have higher denitrification potential and be more responsive to reductions in constraints. Despite these potential mechanisms, our dataset provides limited support for this hypothesis.
Quantile regression depends on individual measurements limited primarily by the factor being tested across its dynamic range, requiring far more data than more conventional regression approaches. Quantile regression coefficients were highly variable across individual cropping systems, with extremely large standard errors in some cases, and some systems had highly implausible negative coefficients. Aggregating systems into annual, perennial fertilized, and perennial unfertilized groups reined in much of that variability, suggesting there were insufficient data at the level of individual cropping systems for an accurate analysis.
More problematically, unfertilized systems almost always had low concentrations of NO − 3 and NH + 4 . In this light, the low regression slopes calculated for WFPS and soil temperature in the group of unfertilized systems seems to reflect that these factors were rarely limiting, rather than a different response. This illustrates to a key challenge of this approach in a field setting; it is highly dependent on the coincidence of permissive levels of multiple environmental factors, which may not happen very frequently. The specificity of cropping system responses to environmental factors remains a highly relevant question that is extremely important for designing management approaches aimed at minimizing N 2 O fluxes.

| Constraint framework suggests steps to avoid hot moments
Conceptualizing N 2 O production in terms of hot moments and environmental constraints provides a useful framework for identifying areas for reducing both the magnitude and variability of cumulative annual N 2 O production. Managing NO − 3 concentrations in the soil may be particularly effective at minimizing the impact of hot moments when other factors are highly conducive to denitrification and is more easily accomplished than managing WFPS or soil temperature. Land managers are unlikely to eschew nitrogen fertilization altogether, but there may be alternative means to reduce NO − 3 concentrations. The use of enhanced efficiency fertilizers (EEFs) such as those containing nitrification inhibitors can promote a more gradual release of NO − 3 over time (e.g., Akiyama, Yan, & Yagi, 2010). Also, there is considerable evidence that increasing the diversity and perenniality of cropping systems reduces the amount of available NO − 3 and NH + 4 in the soil (Duran, Duncan, Oates, Kucharik, & Jackson, 2016;Lu et al., 2011;Stehfest & Bouwman, 2006), likely constraining the magnitude of fluxes that could occur even with optimal temperature and moisture. Diversity and perenniality extend the range of spatial, temporal, and functional niche that can be exploited, reducing high resource concentrations (Oelmann et al., 2007;Palmborg et al., 2005). Further reductions may be possible by addressing processes that increase mineralization of organic nitrogen, such as the increase in freeze-thaw cycles caused by snow removal (Ruan & Robertson, 2017). While these recommendations are largely sound on their own merit, thinking of how they interact with other constraining environmental factors could lead to a more nuanced set of management approaches.