Research was conducted in continuous cotton systems at the Texas A&M AgriLife Research Center in Lubbock, TX (33.687°, -101.827°) where conservation systems were implemented in Fall 2014 and sampling was conducted from 2018-2020. Annual precipitation is 481 mm and the 30-year mean annual temperature for Lubbock is 16.1°C (30-year average, 1991-2020) (27). The soil series at the research site is an Acuff loam described as a fine-loamy, mixed, superactive, thermic Aridic Paleustoll (28).

The study was designed as a split plot with the main plot consisting of different conservation practices (tillage systems), and the split plot consisting of N fertilizer treatments differing in their time of application. Conservation practices included no-tillage with a winter wheat cover crop (NTW) and no-tillage winter fallow (NT), compared to a conventional tillage winter fallow (CT) system. Both the NTW and NT systems were introduced in November 2015. The field was in conventional tillage for at least the previous 60 years. Split plot N fertilizer application timings include the following: no-added N (control); 100% of N applied pre-plant (PP); and 40% of N applied PP and 60% applied mid-season (SPLIT). The study consisted of 3 blocks (replicates) for a total of 27 plots, each measuring 15 m in length and four rows wide with 1 m row spacing. Full agronomic information for the three-year study and a timeline for management and sampling in 2020 is included in the supplemental material ( Supplementary Material 1 : Methods).

Agronomic information for the final year is briefly presented here. Furrow irrigation of the control plots and all treatment plots was conducted on 03 June, 30 June, 10 August, and 26 August with each event applying about 152 mm of water. Pre-plant N fertilizer was applied 6 May 2020. Mid-season N fertilizer applications were made 21 July 2020 near the initiation of reproductive growth. Application of N fertilizer was made at rate of 168 kg ha^-1 as urea-ammonium-nitrate (UAN) and was conducted using a knife and coulter fertilizer applicator mounted to a tractor. Nitrogen fertilizer was applied 10 – 15 cm from the cotton row at a depth of 15 cm. Nitrate added in the irrigation water amounted to 47 kg N ha^-1 in 2020, calculated from the concentration of nitrate-N in the irrigation water, the rate of irrigation applied, and the number of irrigation events (29). Cotton (DP 1820) was planted on 1 June 2020 and harvested on 20 November 2020. Cover crops were planted in November 2019 and terminated on 15 April 2020.

Soil sampling was conducted at three key cotton growth periods: early-season (June, vegetative growth), mid-season (August, peak plant production), and late-season (October, reproductive growth) in 2018, 2019, and 2020. However, only 2020 is presented in the main text due to the inclusion of soil gas and microbial genetic potential measurements being made in conjunction with soil sampling in the final year. In 2020, samples were specifically collected on 24 June, 05 August, and 10 October. Soil samples were collected using 2.5 x 40 cm hand probes to a depth of 20 cm separated into 0-10 and 10-20 cm depths. Six soil cores per plot were collected, homogenized, and stored on ice until a ~25 g subsample could be collected and frozen at -80°C for DNA analysis. The remaining soil was weighed and then dried at 60°C for 7 days to determine gravimetric water content. Soil samples were extracted for soil nitrate-N ( NO 3 − -N) and ammonium-N ( NH 4 + -N) using a 1 M KCl (1:10 soil to extractant ratio, 4 g of soil) and analyzed for NO 3 − -N using flow injection spectrometry (FIAlab 2600, FIAlab Instruments Inc., Bellevue, WA) following cadmium reduction to NO 2 − and for NH 4 + -N by the Berthelot reaction involving salicylate (30). In addition, soil samples were analyzed for mineralizable C (Cmin) by a 3-day rewetting incubation with 40 g of air-dried soil sieved at 4 mm (31). Soil pH was determined by a 1:2 dilution of 10 g soil with DI water (32).

Microbial analyses for genetic potential were conducted by quantitative polymerase chain reaction (qPCR) of major bacterial nitrification and denitrification genes and the genetic marker for bacterial abundance ( Table 1 ) for samples collected in the final year of the study (33–39). DNA was extracted using a Quick-DNA Fecal/Soil Microbe 96 Magbead Kit (Zymo Research, Irvine, CA, USA) with a 30-min, room-temperature drying period prior to final DNA elution. Quantitative polymerase chain reaction was conducted on a Thermo Fisher Quant Studio 5 (Thermo Fisher Scientific Inc., Waltham Mass, US) with gBlock gene fragments specific for each gene assayed ( Supplementary Table 1 , Integrated DNA Technologies, Inc., Coralville, IA, USA). In addition, quality control qPCR was conducted on the Quant Studio 5 by spiking a known concentration qPCR reaction with extracted DNA from this study and evaluating if any inhibition occurred in a process similar to Hartman et al. (40), where extracted DNA was added to a qPCR assay to quantify the abundance of Vibrio alginolyticus with gyrB as the gene target (41). A spiked sample with a threshold cycle (C_t) value within 3 standard deviations of the quality control C_t mean was determined to not contain inhibitors. A small number (<25) of samples were determined to inhibit PCR reactions and were diluted 1:10 with sterile water to reduce inhibition and re-analyzed with qPCR as described above. Initial extractions and dilutions were used to conduct qPCR evaluation of the soil for the 16S rRNA gene, nosZ clade I and clade II (N₂O reductases). Due to elution volume limitations, additional extractions were conducted, and diluted 1:10 with sterile water, as described above for nirK (copper dependent nitrite reductase), nirS (cytochrome cd nitrite reductase), and amoA (ammonia monooxygenase) qPCR evaluations from the same soil subsamples as the previous extractions.

Nitrous oxide emissions and in-ground concentrations of N₂O and NO were collected during the final year of the study using a Gasmet DX4040 portable FTIR gas analyzer (Gasmet Technologies, Helsinki, Finland) integrated with a LiCor 8100A opaque chamber for soil surface emissions. The Gasmet DX4040 was also connected to in-ground pore-space probes (PSPs) to measure in-ground concentrations (42). Soil N₂O emissions and pore-space samplings were conducted within 48 hours of soil sample collection. Soil NO emissions were measured but were not determined to be significant (flux rate below measurable limit) and thus are not included. Pore-space concentration measurements were conducted according to Brummell et al. (43) and McDonald et al. (42). Briefly, PSPs were installed at the beginning of the growing season and then connected to reservoir bottles 24 hours before measurement. On the day(s) of sampling, PSPs with reservoir bottles were integrated into a closed loop system with the Gasmet analyzer for soil N₂O and NO concentration determination. The reservoir bottle increased the system volume to greater than that of the sample cell within the Gasmet, thus eliminating any vacuum being placed on soil pores.

Analyses of variance (ANOVA) for soil chemical parameters was conducted across the three growing seasons within each depth using the lmer function within the lme4 package in R (44) where the initial analysis of variance model included sampling period, conservation system, and N treatment as fixed effects with evaluation of their main effects as well as their combined interactions. Random effects specified in the initial model were rep and the interaction of rep and sampling period to account for repeated measurement of the same research plots over time. The two sampling depths were parsed prior to ANOVA due to the well-known effects of depth on soil parameters, and the desire to evaluate the data in the context of the area from which soil pore-space gases will be measured. In addition, the application of N fertilizer at a depth of 15 cm created a gradient of inorganic-N between the two depths which may have affected the production of greenhouse gases at depth. Consistent over the evaluation of soil parameters was the significant effect of sampling period as a main effect, while only the repeated measure of NO 3 − -N was affected with an interaction of sampling period and N timing. Due to the strong observed effect of sampling period as a main factor, analyses for the final year of the study were split by depth and sampling period to understand within seasonal effects without the variability of the distinct climatic and soil conditions at each sampling period. Analyses of variance of soil chemical parameters for the three-year dataset and final year are presented in the supplementary material ( Supplementary Material 2 : Soil Chemical Results). Microbial genetic potential for N cycling was also evaluated by ANOVA within sampling period and depth for the final year of the study ( Supplementary Material 3 : Microbial abundance in the final year).

Principal components analysis of the soil parameters was conducted using the prcomp function in the stats package in R (45), where an initial evaluation of the relationship of soil parameters was conducted across all three years of the study, and then separately for the final year to support the use of the final year as a reasonable representative of the variability in soil parameters over the length of the study. Soil parameters were normalized to a range of 0-1 prior to analyses to evaluate the variability between samples without the vast differences in scale using the preProcess function of the caret package in R (46). As mentioned above, soil gas measurements were only collected in the final year of the study due to the use of the first two years as test periods for the soil pore-space probes and subsequent mechanical issues with the FTIR analyzer that prevented full-season evaluation of pore-space gases and surface emissions, further necessitating the final-year analyses. Microbial abundance data was subsequently limited to analysis during the final year of the study as well due to the research team’s decision to specifically focus on the final year for the holistic evaluation of the soil biochemical effects of the implemented treatments on gas production. Principal components analysis of microbial gene abundances was conducted using the same process as the analysis for soil chemical parameters detailed above. However, soil chemical parameter vectors were fit to the ordination of microbial gene data using the envfit function in the vegan package in R (47). Principal components ordinations and contributions to axes were plotted using ggplot2 (48)

LASSO regression analysis was conducted using the glmnet and cv.glmnet functions within the glmnet package in R (49). The data was split into training and testing datasets at random with a 75%, 25% split resulting in 114 samples for training and 38 samples for testing (46). The model for each gas was initially developed using leave-one-out cross validation with the cv.glmnet function, after which the best model was refined using the glmnet function and the lambda generated by the initial model. The model was then used to predict gas concentration for the testing dataset using the predict function, and the R² of the fit was calculated to determine predictive capacity of the model for the new dataset. Lastly, the measured vs. the predicted concentrations for the testing dataset were plotted to visualize the fit of the data to the linear model using ggplot2 (48).

Soil chemical analyses were conducted for the combined effect of the implemented treatments and time over the course of the growing season across all three years of the study and are presented in the supplementary material ( Supplementary Material 2 : Soil chemical results). Due to significant effects of sampling period as a main effect and minimal interaction with soil chemical parameters, soil chemical characteristics were evaluated within sampling period and depth for the final year of the study where soil chemical analyses were combined with microbial gene abundance and gas sampling. The results of ANOVA for the final year are also presented in the supplementary material (see Supplementary Material 2 .2 Final Year Analyses).

The relationships between soil chemical parameters as they are affected by the implementation of conservation system and altered N timings was addressed using PCA across the entire dataset ( Supplementary Figure 2 ) and within the final year of the study ( Figure 1 ) across both depths. It is clear that the data in the final year are representative of the pattern of variance present in the three-year dataset due to the similar pattern in variance for the full dataset and the final year alone and the similar total variance explained by PCA of soil chemical parameters between the final year (66.8%, Figure 1 ) and the full three-year dataset (59.2%, Supplemental Figure 2 ). For the final year, PC1 explained 38.7% of the variance in the data where significant contributions to the variability explained included soil NO 3 − -N, pH, and NH 4 + -N which each contributed greater than 20% of the variability ( Figure 1 ). PC2 explained 28.1% of the variance in soil chemical parameters where significant contributions included GWC and Cmin content ( Figure 1 ), Lastly, PC3 also explained a moderate proportion of variability, accounting for 14.6% of the variability with both Cmin content and soil pH significantly contributing to the variance explained ( Figure 1 ).

Correlations between the significant contributors to variance explained for each PC are displayed in Figure 1 . Both NO 3 − -N and NH 4 + -N were negatively correlated with PC1 and soil pH was positively correlated with PC1 ( Figure 1 ); Cmin content was positive correlated with PC2 and GWC was negatively correlated with PC2; and both Cmin content and pH were negatively correlated with PC3. As mentioned above, PCA analyses were conducted across depths and depth was not determined to cluster within the ordination of either soil chemical dataset.

Analyses of variance for genetic abundances are reported in the supplemental material (see Supplemental Material 5 : Microbial abundance in the final year). In brief, bacterial amoA (AOB) abundance was affected at all three sampling periods, while nirS and nirK abundances were unaffected at any of the three sampling periods. In general, bacterial N cycle gene abundances increased across the growing season, with a distinct increase in AOB abundance from the early- to mid-season sampling. As with soil chemical characteristics, microbial genetic abundances were evaluated by PCA to determine overall relationships of implemented treatments and sampling periods on the microbial genetic potential for N cycle activities ( Figure 2 ). In addition, soil chemical characteristics were fit to the PCA ordination of the microbial gene abundances, and their vectors are presented on Figure 2 . It should be noted that vector coordinates were multiplied by 4 for visualization on the figure, but the relationships to the ordination of the microbial gene abundances remain unchanged. In addition, some clustering was observed by depth for microbial gene abundances and ellipses are fit to these clusters. Principal components analysis of microbial gene abundances revealed no distinct clustering by sampling period, conservation system, or N timing. The data was moderately clustered by depth, where the 0-10 cm depth contained much more variability than the 10-20 cm depth ( Figure 2 ). The total variability explained by PCA for gene abundances was 67.2 for PC1 and PC2 combined, where PC1 explained 50.4% of the variance and PC2 explained 16.8% of the variance. In addition, PC3 moderately contributed to variance explained (13.8%). Significant contributors to the variability explained by PC1 included the 16S gene, nosZ clade I, nirK, and nirS where all genes were negatively correlated with PC1 ( Figure 2 ). Only nosZ clade II significantly contributed to the variance explained by PC2, where clade II was positive correlated with PC2 ( Figure 2 ). Lastly, only AOB significantly contributed to the variance explained by PC3, where AOB was positively correlated with the axis ( Figure 2 ).

For the analysis fitting soil chemical parameters to the ordination of genes in Figure 2 , only Cmin content and GWC were determined to be significant, with nearly antithetical relationships to the ordination data ( Figure 2 ).

Nitrous oxide emissions were evaluated at each soil sampling event to assess the impacts of conservation system and N timing on production of N₂O from the soil-atmosphere interface. At the early-season sampling, N timing significantly affected N₂O emissions (p=0.002, Figure 3 ) where emissions from the PP treatment (8071 µg N₂O-N m^-2 d^-1) were greater than from the SPLIT treatment (1697 µg N₂O-N m^-2 d^-1) and the control (-218 µg N₂O-N m^-2 d^-1).

At the mid-season sampling, N timing affected N₂O emissions (p=0.002, Figure 3 ) where emissions from the SPLIT treatment (2307 µg N₂O-N m^-2 d^-1) were greater than from the PP treatment (-106 µg N₂O-N m^-2 d^-1) and the control (-362 µg N₂O-N m^-2 d^-1).

No net production or consumption of N₂O was determined at the late-season sampling ( Figure 3 and Supplementary Table 12,13 ).

Pore-space concentrations of N₂O and NO were analyzed for each sample date at a depth of 7.5 cm and 15 cm. At the early-season sampling at the 7.5 cm depth, pore-space concentration of N₂O was affected by N timing (p=0.041, Figure 4 and Supplementary Table 13 ) where the PP treatment had greater N₂O concentration (1.2 µL N₂O L^-1) than the control (0.88 µL N₂O L^-1), but not the SPLIT treatment (1.0 µL N₂O L^-1). No treatment effects were determined for NO at 7.5 cm or for either N₂O or NO at 15 cm at the early season sampling ( Figure 4 and Supplementary Table 12 ). Pore-space concentrations of NO and N₂O were unaffected at both 7.5 cm or 15 cm depth for the mid- and late-season samplings ( Figure 4 and Supplementary Tables 12, 13 ). Across the growing season, pore-space N₂O concentrations consistently decreased at both the measured depths, while NO concentrations increased from the early- to mid-season samplings and then decreased at the late season sampling at the 7.5 cm depth and were unchanged from early- to mid-season, but decreased at the late-season sampling at the 15 cm depth.

Least absolute shrinkage and selection operator (LASSO) regression modeling was conducted for pore-space concentrations of N₂O and NO to determine the best representation of the drivers of each gas’ production in these semi-arid soils. The models developed for each gas represent the most parsimonious combination of soil chemical parameters and microbial gene abundances that can predict the concentration of N₂O and NO in the soil pore-space. Model coefficients for each gas within each sampling period are reported in Table 2 , and predictive capacity of the model is presented in Figure 5 .

The model generated for the prediction of pore-space N₂O concentration performed moderately, with an R² of 0.46 for the measured vs. predicted concentration ( Figure 5 ). The model coefficients for pore-space N₂O included GWC, pH, and NO 3 − -N ( Table 2 ). The model generated for prediction of pore-space NO performed poorly compared to N₂O with an R² of 0.20 ( Figure 5 ). The model coefficients for pore-space NO included GWC and AOB abundance ( Table 2 ).

Soil chemical characteristics were analyzed across the 3-year study and separately within the final year, where variability in the final year was determined to be similar to that of the entire three-year dataset. Principal components analysis confirmed the strong effect of sampling period on soil chemical conditions, where the early- and mid-season samplings clustered together, with the late-season sampling period moderately separated. Sampling period separation was most apparent along PC2, where the extremely dry soil conditions separated the late-season sampling from the two samplings earlier in the growing season. There was no clear clustering by conservation system or N treatment, likely due to the strong effect of climate and treatment implementation which resulted in temporal clustering by sampling period. We hypothesized that soil management practices would impact soil chemical parameters and water content of the soil, which was confirmed through ANOVA, but did not result in clear patterns of separation with PCA, suggesting a larger overall effect of either climate, or another unmeasured factor.

Principal components analysis of genetic potential generally did not reveal any distinct clustering of implemented treatments or patterns in gene abundance through the growing season for the N cycling community as a whole. Interestingly, the contributions to each PC of the microbial gene abundance PCA were distinct. PC1 was associated with the overall bacterial population and the population of typical denitrifiers (nirS/K and nosZ clade II), PC2 was highly driven by the abundance of nosZ clade II which is known to be highly abundant in the soils in this study (23), and PC3 was affected by AOB abundance. The application of environmental variables to the ordination revealed a significant relationship between the variability in genetic potential and Cmin and GWC. Although no treatments clustered with these parameters, the significant relationships suggest some distal control of C and water on the genetic potential for N cycling in this semi-arid soil. However, there was clearly an additional unmeasured factor driving the variability in genetic potential observed in these semi-arid soils, although it must be noted that the relationship between soil chemical changes and microbial gene abundance may not be linear due to genetic potential including active, dormant, or even relic DNA.

A strong relationship between pore-space N₂O and NO and emissions from the soil surface was expected, as was previously reported in forest soils (26); however, no significant relationship was found in our study due to the sporadic nature of significant positive emissions of N₂O. A reduction in N₂O concentrations in the soil pore-space observed with the progression of the growing season at both depths was likely related to diminishing N resources. The same reduction was also oberved for N₂O emissions but emissions and pore-space concentrations were not correlated (analysis not reported). Compared to previous measurements in agricultural fields (24) and forests (25), the concentration of NO in the soil pore-space was much greater for this semi-arid agricultural soil. We suspect that these high NO concentrations are a consequence of greater N fertilizer input than previous studies. Increases in gas production aligned with reported trends in N gas production in fertilized systems (10–12), and with seasonal variation in gas production (25). Increased NO concentrations at the mid-season sampling likely also reflect general increases in microbial activity due to increased temperatures where NO concentration wasn’t affected by the application of N fertilizer.

It was hypothesized that denitrification was driving the production of N gas production in these semi-arid soils due to the large abundance of nosZ clade I and II genes previously observed (23). LASSO regression analysis produced a model supporting this hypothesis for N₂O pore-space concentration where the model featured positive coefficients for GWC and NO 3 − -N, and a negative coefficient for pH. Both GWC and NO 3 − are known drivers of N₂O production (9–12), and their inclusion in the model as positive coefficients supports the idea that the pore-space N₂O measured is resulting from denitrification activities. When considering the timing of N₂O emissions and the method of application (liquid UAN knifed into the soil), it is not surprising that denitrification is the main driver of excess gas production that translates to enrichment of N₂O in the pore-space and non-directly to emission of N₂O from the soil surface. Furthermore, the inclusion of pH in the model is key, as the effects of pH on N₂O reduction capacity are well known (10). The regression model produced for N₂O pore-space concentrations was moderately predictive, resulting in an R² of 0.46 for the testing dataset. Basing the prediction of pore-space gas on such a small set of field measurements (and only five total soil chemical parameters) likely resulted in reduced predictive power. However, the moderate predictivity of the model does not dimmish the lack of genetic potential contribution to the model for N₂O in the pore-space.

In contrast to the denitrification driven concentration of N₂O in the soil pore-space, the LASSO regression model for NO concentration featured indicators of nitrification. Increased ammonia oxidation by AOB could have resulted in pooling of NO in the soil pore-space if it overwhelmed the oxidizing and reducing capacity of the N cycling population (13). The addition of GWC to this model indicates the controlling nature of water in this semi-arid soil. However, it must be noted that the model was poorly predictive, producing an R² of 0.2 for the testing dataset comparison of measured and predicted pore-space NO concentrations. Due to the poor performance of the model, we cannot clearly define the production of NO as a nitrification activity, especially without the inclusion of NH 4 + in the model. Despite the lack of connection to genetic potential these models represent early indications of distinct N-cycle pathways producing each of these important GHGs and further represent the clear differences between N gas production in semi-arid lands and mesic systems where carbon plays a stronger role in N gas production (6, 7, 15).

To further consider the interactions of soil chemistry, genetic potential, and N gas production we will examine the PP timing of N fertilizer application. Emissions and pore-space concentrations of N₂O peak during the early-season sampling for the PP timing due to the recent application of N fertilizer (10–12) supporting our hypothesis that N fertilizer application would relieve the poor N content of the soil, jumpstart N cycling, and overwhelm the N₂O reducing capacity of the soil. By the mid-season sampling, net consumption of N₂O was recorded for the PP timing and was likely related to the loss of soil NO 3 − -N by plant and microbial consumption. Low inorganic-N concentrations are known to be associated with N₂O consumption (50–52), and by mid-season, NO 3 − -N was greatly reduced compared to the early-season.

Despite the reduction in emissions, and pore-space concentrations, genetic potential for these pathways increased from the early- to mid-season samplings in the PP timing. We suspect there are two potential reasons for this seemingly antithetical outcome. First, we suspect the genetic abundance increase measured at the mid-season point may reflect the past increase in N cycle gene abundance and subsequent N cycle activity which resulted in the loss of NO 3 − from the soil, and reduced emissions at the mid-season sampling period. This is supported by the abundance of nosZ clade II, which increased from early- to mid-season reflecting a need for N₂O-reductase in the preceding months where high N cycle activity produced more N₂O and thus increased the need for this gene. The high level of N cycle activity would also thus be reflected in greater abundances of other N cycle organisms from the early- to mid-season sampling periods which we also observed.

A second potential reason for this outcome follows similar logic to the first where a significant loss of soil NO 3 − -N occurred between the early-and mid-season sampling, related to plant and microbial consumption. However, the loss of soil N would have selected for non-traditional N₂O reducers as has been seen in several environments (26, 51, 52). In this scenario, the reduced NO 3 − is indicative of an active population of N cycling organisms and an active crop that is consuming soil inorganic-N. With continually reducing N resources, the plant may be selecting for organisms that can mineralize organic matter and supplement the N supply. For AOB, this would likely be a highly plausible reason for the increase in abundance as there would be fresh NH 4 + continually being released from the soil matrix. However, for the other genetic abundance increases, this is less plausible than the reflection of past activity and environmental selection resulting in the measured community abundance as explained above. Furthermore, due to environmental selection and plant selection, one would expect to see a decrease in genetic potential for N cycle activities between the mid- and late-season samplings, which was not reflected in this study.

The lack of reduction in genetic potential over the growing season further reflects the difficulty in using gene abundance to predict microbial activity. When evaluating microbial processes that aren’t consistently active in the study system, the stable nature of DNA in the environment and the potential for dormant or inactive microbial species, clouds the connection between genetic potential and actual environmental outcomes.