Climate mitigation potential and soil microbial response of cyanobacteria‐fertilized bioenergy crops in a cool semi‐arid cropland

Bioenergy carbon capture and storage (BECCS) systems can serve as decarbonization pathways for climate mitigation. Perennial grasses are a promising second‐generation lignocellulosic bioenergy feedstock for BECCS expansion, but optimizing their sustainability, productivity, and climate mitigation potential requires an evaluation of how nitrogen (N) fertilizer strategies interact with greenhouse gas (GHG) and soil organic carbon (SOC) dynamics. Furthermore, crop and fertilizer choice can affect the soil microbiome which is critical to soil organic matter turnover, nutrient cycling, and sustaining crop productivity but these feedbacks are poorly understood due to the paucity of data from certain agroecosystems. Here, we examine the climate mitigation potential and soil microbiome response to establishing two functionally different perennial grasses, switchgrass (Panicum virgatum, C4) and tall wheatgrass (Thinopyrum ponticum, C3), in a cool semi‐arid agroecosystem under two fertilizer applications, a novel cyanobacterial biofertilizer (CBF) and urea. We find that in contrast to the C4 grass, the C3 grass achieved 98% greater productivity and had a higher N use efficiency when fertilized. For both crops, the CBF produced the same biomass enhancement as urea. Non‐CO2 GHG fluxes across all treatments were low and we observed a 3‐year net loss of SOC under the C4 crop and a net gain under the C3 crop at a 0–30 cm soil depth regardless of fertilization. Finally, we detected crop‐specific changes in the soil microbiome, including an increased relative abundance of arbuscular mycorrhizal fungi under the C3, and potentially pathogenic fungi in the C4 grass. Taken together, these findings highlight the potential of CBF‐fertilized C3 crops as a second‐generation bioenergy feedstock in semi‐arid regions as a part of a climate mitigation strategy.


| INTRODUCTION
Developing sustainable bioenergy and carbon capture systems (BECCS) is a critical component of global climate change mitigation (IPCC, 2018), but the environmental implications have yet to be thoroughly examined at field scales. The climate mitigation potential of sustainable bioenergy systems largely hinges on understanding how these systems affect carbon (C) and nitrogen (N) cycles (Robertson et al., , 2017. Given geographical variations in soil and climate, it is necessary to identify bioenergy crop species and fertilizer regimes that maximize crop yield while minimizing negative environmental impacts, including greenhouse gas (GHG) emissions (e.g., CH 4 , CO 2 , N 2 O) and loss of soil organic carbon (SOC). One promising avenue for BECCS pairs second-generation lignocellulosic perennial grass with cyanobacteria-based Nbiofertilizer in semi-arid climate . However, the interactions between C and N cycling, and the potential trade-offs between biomass production and environmental impacts have yet to be evaluated in cyanobacteria-fertilized bioenergy systems.
Semi-arid ecosystems cover 45% of the earth's land surface, are critical reservoirs for carbon sequestration (Lal, 2009;Poulter et al., 2014), and often contain vast tracts of marginal land unsuitable for traditional cultivation due to environmental constraints over crop production from factors such as low precipitation and soil organic matter. Thus, areas in semi-arid climates are likely candidate regions for BECCS expansion with low land-use conflicts (Ahmed et al., 2021;Stoy et al., 2018). However, it remains unclear how BECCS expansion would influence ecosystem C, especially in response to different N fertilization regimes (Dolan et al., 2020). The Upper Missouri River Basin (UMRB) is a promising candidate region that serves as a model system for evaluating BECCS in semi-arid regions (Stoy et al., 2018). The UMRB spans ~800,000 km 2 across Montana, Wyoming, North and South Dakota and has over ~150,000 km 2 of land classified as marginal (Dolan et al., 2020). Additionally, soils across the UMRB have a high C soil sink capacity based on relatively slow decomposition rates and currently low SOC saturation levels (Georgiou et al., 2022;Lal, 2009;Sanderman et al., 2017). Despite these favorable factors, only one study in the western UMRB has reported empirical data on soil and gas C and N dynamics in a dedicated bioenergy study, and only for switchgrass (Panicum virgatum; Schmer et al., 2012). Schmer et al. (2012) found that N 2 O emissions were higher under modest synthetic N fertilization (67 kg N ha −1 ) compared to a non-fertilized control and that overall CO 2 -equivalent GHG flux was dominated by CO 2. A recent regional modeling study across the entire UMRB found that switchgrass expansion and associated management practices led to net losses of ecosystem C (Dolan et al., 2020), but was not parameterized with other potential bioenergy species that are likely to have different outcomes.
A critical aspect of developing a BECCS system is identifying crop species that are well adapted to the regional climate and soil, achieve high levels of aboveground biomass (AGB) and belowground biomass (BGB) production, and have traits (e.g., biomass chemical composition) that support the efficient transfer of biomass-C to SOC (Carvalho et al., 2017;Córdova et al., 2018). Perennial grasses have been widely identified as an environmentally advantageous feedstock in bioenergy systems due to their low nutrient and irrigation requirements and ability to limit net GHG losses and leaching losses of fertilizer-N (Abraha et al., 2019;Ledo et al., 2020;Robertson et al., 2000). Switchgrass has been identified as a model bioenergy crop in the United States based on its favorable biological, agronomic, and physiological characteristics (Wright & Turhollow, 2010). However, despite the agronomic success of switchgrass in warmer and wetter regions of the US, the establishment and productivity of switchgrass west of 98°W latitude remains unclear based on potential climate constraints (Berdahl et al., 2005). Although C4 grasses like switchgrass are generally thought to be more water-use efficient, their productivity is strongly limited by extended periods of cold temperatures and, therefore, may have reduced growth in cold semi-arid regions (Epstein et al., 1997;Monson et al., 1983;Paruelo & Lauenroth, 1996). Comparatively, C3 species may be better matched to the current climate of the western UMRB given they have greater cold tolerance, the ability to maximize photosynthetic capacity during the cool and moist spring, and can competitively acquire N in semi-arid climates (Karp & Shield, 2008;Luo et al., 2018).
Across much of the UMRB, synthetically produced urea is the principal form in which N is applied to crops, and its mismanagement has resulted in sizeable reactive N losses (e.g., ammonia volatilization) from the plant-soil system (Engel et al., 2011;John et al., 2017). Cyanobacterial biofertilizer (CBF) is an organic, microbial-based N-fertilizer that has the potential to replace petroleum-based fertilizer in perennial bioenergy systems . Across field and greenhouse studies, CBF has shown the ability to promote SOC accumulation, increase nutrient availability, and increase soil N retention while maintaining crop productivity comparable to a urea fertilizer (Alvarez, Weyers, Goemann, et al., 2021;Alvarez, Weyers, Johnson, et al., 2021;Castro et al., 2017;Goemann et al., 2021). Moreover, the production of algal biofertilizers has a significantly lower carbon cost in comparison to urea (Rupawalla et al., 2021;Shi et al., 2020). In addition to C costs from production processes, field application trials are critical to assess the contribution of soil GHG emissions to the entire cradle-to-field life cycle assessment of C and the CO 2 -equivalence. Yet, no previous field study has evaluated how the application of CBF over multiple growing seasons affects GHG emissions, soil carbon stocks, and soil microbial communities.
The cultivation of a beneficial soil microbiome is an emerging management strategy to enhance the sustainable production of bioenergy biomass, especially in challenging environments such as marginal or semi-arid lands (Zhalnina et al., 2021). Beneficial plant-microbe interactions in agricultural settings include increases in yield via nutrient acquisition through arbuscular mycorrhizal fungi colonization (AMF; Gao et al., 2020;Smith and Read, 2008), the alleviation of abiotic stresses (e.g., drought), pathogenic protection, and enhanced C sequestration through persistent SOM formation (Mendes et al., 2011;Sher et al., 2020;Six et al., 2006). However, mechanisms for these feedbacks remain poorly constrained, and more field data are required to better assess how soil microbial communities respond to crop and fertilizer interactions in sustainable agricultural operations (Lekberg & Helgason, 2018;Rillig et al., 2016). Thus, to optimize beneficial plant-microbe interactions, it is necessary to first quantify the magnitude and direction of microbial community shifts to best evaluate their influence on crop production and C and N processes important to climate mitigation goals.
To address these gaps in field-based knowledge of perennial bioenergy cropping systems and biofertilizers in the UMRB, we evaluate the production and potential biogeochemical tradeoffs of establishing two perennial types of grass, switchgrass (C4) and tall wheatgrass (Thinopyrum ponticum; C3) and compare the efficacy of a novel cyanobacteria biofertilizer to conventionally used urea fertilizer. Our objectives were to (1) quantify AGB and BGB production, (2) quantify changes in SOC stocks at two soil depths, (3) assess GHG flux of CO 2 , CH 4 , and N 2 O during the growing seasons, and (4) characterize shifts in soil microbial composition for crop and fertilizer treatments.

| Site description and experimental design
The study site is located at Montana State University's (MSU) Arthur H. Post Farm in Bozeman, Montana (45.66°N, −111.15°W; Figure 1) at 1450 m above sea level. Mean annual temperature is 6.5°C and mean annual precipitation is 408 mm, most of which falls during the spring (Miller & Holmes, 2012). Field soils at the Post Farm are classified as Amsterdam silt loam, a fine-silty, mixed, superactive, frigid Typic Haplustoll. Experimental plots were established in April 2018 and seeded at 16 kg ha −1 with a no-till disk seeder. The previous crop was proso millet (Panicum milaceum) that was fertilized with 67 kg N ha in spring 2017. Preceding 2017 and since 1999, this site had been cropped annually with a wide range of alternative crops using low-disturbance no-till seeding. The plot locations have low background soil fertility (SOC 1% ± 0.03, total nitrogen [TN] 0.11% ± 0.03; Oldfield et al., 2019) and no irrigation. The experiment was established as a F I G U R E 1 Aboveground biomass yields from reported switchgrass and tall wheatgrass studies from within and around the Upper Missouri River Basin (UMRB) between 2008 and 2020. All studies are non-irrigated but vary in nitrogen fertilization rates (Table S1). Grey areas depicted are major bodies of water. The star in the inset marks the location of the field site for this study (Bozeman, MT). complete randomized split-plot (each split-plot is 3.34 m 2 ) design with three fertilizer treatments (urea, CBF, and unfertilized) for each grass species (switchgrass; Dakotah variety) and (tall wheatgrass; Alkar variety) that was replicated five times ( Figure S1). Fertilizer treatment was considered the main plot while grass species was subplot. Switchgrass is a well-studied C4 species favored for second generation bioenergy production because of its chemical composition for post-processing, high ethanol yield, and significant biomass potential (Wright & Turhollow, 2010). Conversely, tall wheatgrass is a less characterized C 3 species (Martyniak et al., 2017), but has been shown to have an optimal chemical composition for efficient combustion (Xue et al., 2011), is a relatively salt-tolerant (Bleby et al., 1997). The experimental design, agronomic management, and initial soil sampling are described in detail in Goemann et al. (2021).

| Growing degree days, rain use efficiency, and soil pH
Growing degree days were used to evaluate the effect of temperature on plant-specific growth rates and were calculated for each crop by subtracting the base temperature (tall wheatgrass = 0°C, switchgrass = 10°C; Clifton-Brown et al., 2011;Miller et al., 2001) from the daily mean temperature measured at the Post Farm in Bozeman, MT where the base temperature is defined as the lower developmental threshold of each species. Rain use efficiency (RUE) is measured as crop production per unit rainfall (kg dry matter ha −1 mm −1 ) across the agricultural year (Le Houérou, 1984). Here, the agricultural year for calculating RUE was October 2019 to August 2020, as this encompasses the groundwater recharge period from the previous fall and precipitation during the primary growing season months of 2020. Soil pH was measured in the laboratory on air-dried soil from 0 to 15 cm and 15 to 30 cm soil depths and analyzed using a 2:1 milli-Q deionized water to soil solution and glass electrode.

| Fertilization
The CBF utilized was a mixture of live Nostoc sp. isolated from soil directly adjacent to the field plots and commercially available dried Spirulina. Due to limitations on culturing capacity, we used 20:80 Nostoc:Spirulina ratio to meet the N fertilization goal. The Nostoc sp. was isolated by culturing a slurry of moist field soil in BG-11 0 (nitrogenfree) media and plating the slurry on sterile petri plates of BG-11 0 agar. Repeated streaking was conducted until visual isolation of a single morphology was evident via light microscopy at which point the culture was scaled up into increasingly larger volumes of sterile liquid media. CBF cultivation was conducted as described in Goemann et al. (2021). Urea was hand-spread across the plots while the CBF was applied as a slurry from a backpack sprayer and was applied either immediately before or after precipitation events to increase infiltration into the soil. In 2018, 100 kg ha −1 of fertilizer was applied to maximize the harvestable yield of barley (Hordeum vulgare, Hayse forage variety) as an intercrop. During the subsequent two growing seasons (2019 and 2020), 50 kg N ha −1 of urea and CBF were applied equally across the spring months due to production constraints of CBF. In total, there were five replicate plots of urea, CBF, and unfertilized treatments. In 2020, plots were fertilized over three application periods May 13, May 27, and June 5. After the first year, the fertilization rate was reduced from 100 to 50 kg N ha −1 to reduce potentially negative environmental effects (Lai et al., 2018).

| Soil and sampling and analysis
Pre-treatment soil samples were collected before seeding during April 2018, and soils were resampled at the end of the growing season during September 2020 to assess the 3-year change in SOC. During September 2020, two fixed depth intact soil samples were retrieved with an auger (1.9 cm inner diameter) from each split-plot and composited into 0-5 cm, 5-15 cm, and 15-30 cm depth increments for chemical and microbial analyses. In each split-plot, two soil samples were taken (one near the base of a plant and one from an inter-row area to capture plot heterogeneity) and pooled together. Since baseline soil samples were only collected at 0-15 cm and 15-30 cm depths in 2018, a weighted average was calculated for the 2020 samples at the 0-15 cm depth (Methods S1). Before chemical analyses, soil samples were oven dried at 105°C for 48 h before being sieved at 2 mm and then milled in a slowrolling jar tumbler for 12 h. Prior to combustion, for SOC and TN quantification, samples were acidified using 1 N and then 2 N HCl to remove the inorganic C fraction following the MSU Environmental Analytical Laboratory's in-capsule aqueous acidification protocol (Methods S3). After 48 h, % SOC and % TN were determined via dry combustion using an ESC 4010 elemental analyzer (Costech) at MSU's Environmental Analytical Laboratory. Duplicates were run for ~10% of the samples to determine analytical precision. Bulk density was quantified using a slide hammer to reduce soil compaction by taking soil cores (5 cm diameter and 15 cm height) at 0-15 cm and 15-30 cm adjacent to each replicate row in 2019 (n = 10; Figure S1). Cores were then composited by depth before drying at 105°C to a constant mass and sieved at 2 mm. Bulk densities of 1.28 g cm −3 at the 0-15 cm soil depth and 1.23 g cm −3 at the 15-30 cm depth were calculated and used to calculate SOC pools from 2018 and 2020.

| Aboveground and belowground productivity and biomass chemistry
Harvestable AGB was estimated by hand-harvesting a 1 × 1 m section 10 cm above the ground from the center of each plot on September 9, 2020. We acknowledge that this method underestimates total annual net primary productivity because of the exclusion of the plant stubble. Other estimates from tall wheatgrass plots at the Post Farm suggest this may underestimate total AGB by approximately 15% such that our measurements are a conservative estimate of AGB. AGB samples were air-dried at 40°C until they reached a constant mass and milled to a fine powder. To assess cumulative belowground root biomass accumulation, we collected two soil cores from each plot by first clearing litter material and then using a root auger (7.5 cm diameter) to collect a 15 cm deep soil core. To calculate the mean annual production of BGB, we divided the total biomass of each plot by the 3 years since crop establishment. Since the majority of perennial grass root biomass is usually concentrated in the surface-soil layer we used 0-15 cm as a root coring depth (Sainju et al., 2017). In total, three out of the five reps from both bioenergy crops (n = 18) were used to collect soil root cores. After collection, rootsoil cores were gently rinsed underwater in a fine-mesh sieve until all soil was removed. Retained root material was then oven-dried at 40°C until it reached a constant dry mass before being hand-milled to a fine powder. Total C and N concentrations for all AGB and BGB samples were determined by dry combustion on a Truspec CN analyzer (LECO) at MSU. C:N ratios are reported by mass.

| GHG flux collection
GHG fluxes (CO 2 , CH 4 , and N 2 O) were measured over the 2020 growing season (March-October) using the closed static chamber technique. Gas measurements were taken via a portable cavity ring down spectroscopy (CRDS) gas analyzer (G2508, Picarro Inc.) with gas concentration data output recorded every 2 s in ppm and ppb. Dry mole fractions of CO 2 , CH 4 , and N 2 O are reported to account for the cross-influence of H 2 O during sampling. Capturing soil gas flux spatial and temporal heterogeneity are two challenging aspects of assessing ecosystem C and N dynamics. As such, sampling times and dates were scheduled to capture a range of field conditions and specific events such as fertilizer application, heavy rainfall, and freeze-thaw cycles. Measurements were made over 24 separate sampling dates from March to October 2020 for a total of 204 individual chamber measurements and were made to pair crops as closely together as possible, usually a day apart. Moreover, to best capture the high temporal variability of N 2 O fluxes sampling frequency was increased earlier in the growing season during periods when larger fluxes were expected due to higher soil moisture, and we decreased sampling frequency later in the growing season when fluxes were expected to be low (Barton et al., 2015). Daily mean fluctuations of soil GHG flux are largely driven by temperature and to account for diurnal variability measurements we made measurements in midmorning to best correspond to daily average temperatures. As such, we started the 2 h measurement cycle as close to 10:00 local standard time but never exceeded 15:30 (Jian et al., 2018), and GHG measurements follows the Greenhouse gas Reduction through Agricultural Carbon Enhancement network (GRACEnet) protocol established by the US Department of Agriculture-Agricultural Research Services (USDA-ARS).
Collars were inserted to a 5 cm depth in early March 2020 (one per plot) and left in place for the growing season to minimize effects of soil disturbance. Collars were installed to include a single plant from each plot, and each plant was trimmed to ensure a tight seal when the chamber was attached. During each sampling event, fluxes were measured in three out of the five replicates, and the sequence by which the plots were measured was randomized to avoid bias. Flux measurements were made by connecting the gas analyzer via 2 m Teflon tubing (1/8″ inner diameter) to a vented and white polyvinyl chloride (PVC) chamber attached to the PVC collar inserted in the soil. Together, the chamber and collar have a diameter of 0.032 m 2 and volume of 4.86 L. Chamber caps were left on the collar for 5 min, gas concentrations were measured automatically every 2 s, and data output was recorded in ppmv and ppbv. The chamber lid was flushed via sweeping through the air three times or until a stable atmospheric baseline was met between sampling each plot. Soil temperatures were recorded with a digital thermometer for each chamber measurement and collected from a shaded location directly outside the PVC collar at a 5 cm depth. During each GHG flux measurement, small soil cores (2.5 cm diameter) to a 5 cm soil depth were collected directly outside the collar, placed in an airtight bag in a cooler, and transported directly back to the laboratory to determine gravimetric soil water content. Water filled pore space (WFPS) was calculated via its gravimetric soil water content after drying soil samples at 105°C for at least 24 h and adjusted for bulk density and soil porosity (Methods S1).

| Statistical analysis
Soil organic C and N concentrations from each soil depth were converted to a mass per unit area (Mg ha −1 ) by multiplying the concentrations by bulk density and the soil depth increment. Changes in SOC between sampling times were calculated by subtracting initial values from values after 3 years within the same plot (Liebig et al., 2008). SOC change was then evaluated within plot by the fixed effects of soil depth, crop, and fertilizer treatment using a linear mixed model with plot nested in replicate as random effects to account for unmeasured site variability and the split-plot design using the lmer package (Bates et al., 2015) in R. Linear mixed models were also used to assess the effect of crop type and fertilizer treatments on biomass production and tissue chemistry also with plot nested by replicate as the random effect. To evaluate if the fixed effects (and their interactions) accounted for variation in response variables, we performed type II ANOVA using the "Anova" function in the lmerTest package (Kuznetsova et al., 2017). When fixed effects accounted for variation, we compared treatments using post-hoc pairwise Tukey HSD and report the model means and ± SE. We used Levene's test to evaluate the plot-level variance of SOC concentrations from 2018 to 2020 using the leveneTest in the car package (Fox & Weisberg, 2019). The soil chemistry and biomass results are reported as model mean values (±SE). Soil pH is reported as the mean of hydrogen ion concentrations followed by a log 10 back-transformation with uncertainty reported as a 95% confidence interval (CI). Nitrogen use efficiency (NUE crop ) was calculated following Equation (1) and represents the fraction of fertilizer N that was utilized and allocated to biomass N (Congreves et al., 2021): To calculate GHG fluxes, we quantified whether each chamber gas measurement exhibited a linear increase or decrease in their concentrations over the 5 min enclosure time ( Figure S4, Methods S2). We fit linear mixed models to assess the fixed effects of crop, fertilizer treatment, and time (month), soil temperature, and WFPS on daily mean GHG fluxes and included the random effects of replicate crossed with sample day to account for repeated measures. To examine the temporal differences of crop-specific fluxes based on different plant phenology and in response to environmental factors, we examined possible interactions between month and crop. WFPS and soil temperature were assessed for inclusion in all statistical flux models. To assess whether the fixed effects (and their interactions) accounted for variation in daily mean and cumulative fluxes, we performed type II ANOVA using the "Anova" function. When fixed effects accounted for variation, we compared treatments using post-hoc using pairwise Tukey HSD and report the model means and ±SE (emmeans package; Lenth, 2021). Flux uncertainty is reported as ±SE of the model means across the three replicates measured during each sampling event. All model residuals were assessed for normality, and logarithmic transformations were applied to the AGB, biomass %N (AGB and BGB), and CO 2 flux data before analysis to achieve a normal distribution. All analyses were performed in R (R Core Team, 2022). Treatment significance was set a priori to α = 0.05 for all analyses.

| Soil microbial community analysis
rRNA gene amplicon sequencing was performed to determine shifts in the composition and diversity of soil microbial communities under crop and fertilizer treatments. The soil used for sequencing was flash-frozen in an ethanol-dry ice bath in the field, and the samples were collected as subsets of the cores extracted for soil chemistry analysis to examine microbial relationships directly.
Genomic DNA was extracted from 0.5 g frozen soil samples using the Fast DNA SPIN Kit for soil (MP Biomedicals) according to the manufacturer's instructions. We targeted the V4 region of the 16S rRNA gene using f515a and r806b (GTGYCAGCMGCCGCGGTAA and GGACTACVSGGGTATCTAAT, respectively) from the Earth Microbiome Project (Thompson et al., 2017) and the D2 hypervariable region of the fungal LSU using LR22r and LR3 (Mueller et al., 2016). Details of PCR reactions and sequencing are provided in Goemann et al. (2021).
The two libraries were pooled at 4 nM and loaded onto an in-house Illumina MiSeq (Illumina). Reads were merged, trimmed, and dereplicated with USEARCH, and zero-radius operational taxonomic units (ZOTUs) were identified with UNOISE3 (v.11.0.667;Edgar, 2010). Representative 16S and LSU ZOTUs were classified against their respective databases using the online Bayesian classifier of the Ribosomal Database Project (Set 18 for bacteria, Set 11 for fungi; Wang et al., 2015). All reads classified to the chloroplast (16S), or animalia, protozoa, or viridiplantae (LSU) were removed before analysis, and datasets were randomly subsampled according to the sample with the lowest read number for statistical analysis using the R package vegan (Oksanen et al., 2022). The data are accessible at MG-RAST under accession numbers: mgs875156 (LSU) and mgs875153 (16S). Differences in community composition were assessed with PERMANOVA using the vegan package and relative abundances calculated with the microbiome package (Lahti & Sudarshan, 2017). The rstatix package (Kassambara, 2021) was used to calculate (1) NUE crop = Yield N f Fertilizer N pairwise t-tests with multiple test corrections for differences in relative abundances between crops and compared to the 2018 samples.

| Soil WFPS, climate, and pH results
WFPS (0-5 cm) was highest in the spring and lower in the summer with monthly means of 47% in March and 12% in August (p = 0.02; Figure S10b, Table S5). WFPS was highest under the control plots (p < 0.0001) regardless of crop or fertilizer type. Soil temperature increased throughout the growing season from 4.9 ± 2.1°C in March to 22.7 ± 2.1°C in August across treatments ( Figure S10a). Weather conditions in 2020 favored tall wheatgrass (C3) growth as there were 2.8 times more growing degree days (2845 vs. 1019; Figure S4) and greater RUE compared to switchgrass (24.7 vs. 8.7 kg dry matter mm −1 rainfall; Figure S4). There were no differences in soil pH across crop type (p = 0.67) or treatment (p = 0.67) but note a higher pH at the 15-30 cm depth (8.4, CI: 7.9-8.8; p = 0.0001) compared to the 0-15 cm depth (8.0, CI: 7.5-8.4), due to the greater carbonate presence at this depth (Engel et al., 2017).

| SOC change from 2018 to 2020
The direction and magnitude of SOC change over the 3year establishment period differed across bioenergy crop type (p = 0.009; Figure 3) and soil depth (p < 0.0001), while fertilizer type did not affect SOC dynamics over the duration of this study. We observed a net gain of 0.75 ± 0.44 Mg ha −1 of SOC (0.25 Mg ha −1 year −1 accumulation rate; Table S1) in the tall wheatgrass plots and a net loss of −0.23 ± 0.44 Mg SOC ha −1 (−0.08 Mg SOC ha −1 year −1 ) in the switchgrass plots across a 0-30 cm soil depth. Tall wheatgrass plots had a ~ 3-fold greater SOC increase at the 15-30 cm depth over 3 years (1.58 ± 0.48 Mg SOC ha −1 vs 0.61 ± 0.5 SOC Mg ha −1 ). At the 0-15 cm depth, there was no change in SOC over 3 years in the tall wheatgrass (p = 0.32) and a 3-year net loss of 1.09 ± 0.5 Mg SOC ha −1 in the switchgrass plots (p < 0.0001). We find decreasing soil sample variation of SOC concentrations from 2018 to 2020 at the 0-15 cm soil depth, whereas variation increased at the 15-30 cm soil depth in both crops ( Figure S9 and Table S10). Over F I G U R E 2 Least squared mean values of above (AGB) and belowground (BGB) root biomass (0-15 cm) production in 2020 for switchgrass (SG) and tall wheatgrass (TWG). *Denotes significance from Tukey's HSD post-hoc (α = 0.05) between N fertilizer treatments and the unfertilized control while † denotes differences between crop type biomass pools, irrespective of fertilizer treatment. the 3 years, no differences in bulk soil N were observed across depth, crop, or fertilizer treatments (p = 0.53, p = 0.2, p = 0.84). Thus, C:N ratios increased from a background of 7.9 ± 0.35 under tall wheatgrass by 0.51 ± 0.35 and decreased under switchgrass by 0.54 ± 0.35 across the 0-30 cm soil depth (p = 0.008).

| Soil GHG flux
Soil efflux for CO 2 displayed strong seasonal variability by crop type (p < 0.0001) while CH 4 and N 2 O did not (p = 0.97 and p = 0.39, respectively; Figure 4). Monthly CO 2 fluxes were greater in the tall wheatgrass plots in March and April (p = 0.03 and p = 0.002; Table 1) with some evidence for greater fluxes in May (p = 0.07; Table 1). This shifted in June when CO 2 fluxes converged and remained similar between crops through October. Daily CO 2 fluxes were affected by both crop type (p < 0.0001, Table 2) and fertilizer treatment (p = 0.0007, Table 3). Tall wheatgrass plots had higher mean CO 2 fluxes compared to switchgrass (1.7 ± 0.19 vs. 1.2 ± 0.13 μmol CO 2 m −2 s −1 , p = 0.0118) irrespective of fertilizer treatment. CO 2 fluxes were greatest under the CBF fertilization (Table 3; p = 0.0007). There was little evidence for differences in CH 4 uptake between crop (p = 0.1) or fertilizer type (p = 0.18). All methane fluxes except one measurement in June were negative across the sampling period. Mean N 2 O emission rates were not different between fertilizer treatments or crops (p = 0.66, p = 0.18). We note that N 2 O emissions were generally small; the largest N 2 O flux was captured in March when WFPS was highest (2.1 μmol N 2 O m −2 hour −1 ; Table S5). CO 2 fluxes followed a general positive nonlinear relationship with soil temperature, and CO 2 fluxes under switchgrass were more sensitive to temperature than tall wheatgrass (p = 0.05). There was no observed effect of WFPS on CO 2 flux, nor was there an observed effect of temperature on WFPS or CH 4 fluxes. N 2 O flux increased linearly with WFPS (p = 0.01) but was not affected by soil temperature.

| AGB and root biomass chemistry
The %C in switchgrass AGB and root biomass was 2.9% and 4.7% greater compared to tall wheatgrass (p < 0.0001, p < 0.0001; Table 4) and the %N in switchgrass root biomass was 26% greater than tall wheatgrass (p = 0.05). However, we find little evidence for differences in mean root C:N ratios (53.9 vs 63.5; p = 0.14) or AGB C:N ratios (72.1 vs. 77.4 p = 0.44) between switchgrass and tall wheatgrass. %C was greater in fertilizer treatments relative to controls irrespective of crop type (p = 0.02). A greater amount of biomass-N was lost from harvest (31.7 vs. 13.6 kg N ha −1 ; p < 0.0001) in the tall wheatgrass plots, suggestive of a higher nitrogen use efficiency (Equation 1; NUE crop ). If we assume all plant N was acquired from a fertilizer source, 63% of fertilizer N was recovered in harvestable tall wheatgrass yield, while only 27% was recovered from switchgrass.

| Microbial biodiversity shifts
Switchgrass and tall wheatgrass plot bacterial/archaeal and fungal community compositions shifted between 2018 and 2020 (Figure 5a,b, Tables S6 and S7). Soil depth contributed to the dissimilarity in the fungal community, with greater shifts observed in the 15-30 cm depth (p = 0.001, F I G U R E 3 Box plots (median, 5th and 95th percentiles) of the 3-year SOC change, with fertilizer treatment controlled for, at soil depth increments of 0-15 (a) and 15-30 (b) and 0-30 cm (c) across tall wheatgrass (TWG) and switchgrass (SG) plots. Black dots are least squared means and * denotes significance from Tukey's HSD post-hoc test (α = 0.05) between the two crop treatments. Large black circles depict mean values. Table S5). Fertilizer-driven treatment effects on community composition were detected in the fungi (p = 0.01) but not in the bacteria/archaea (p = 0.28). Fungal but not bacterial/archaeal alpha diversity also increased over time at the 15-30 cm soil depth, particularly in the tall wheatgrass plots ( Figure S7).
To examine the differences in community composition between crop types and their shift from the 2018 pre-establishment conditions, we compared the relative abundances of major bacterial/archaeal phyla and fungal classes. In the bacterial/archaeal community, Planctomycetes relative abundance increased in both crops while Actinobacteria increased in only the tall wheatgrass and Thaumarchaeota in the switchgrass plots (Figure 5c). Gemmatimondetes and Proteobacteria decreased from 2018 to 2020 in both crops. Due to the observed differences in CH 4 between crops, we further investigated clades of known methanogens and methanotrophs. Euryarchaeota, a clade primarily composed of methanogens, was enriched in both switchgrass and tall wheatgrass compared to the background soil and had the greatest abundance in switchgrass ( Figure S6).

F I G U R E 4
Time series of soil fluxes of CO 2 , CH 4 , and N 2 O across tall wheatgrass and switchgrass fertilizer treatments. Dotted black lines identify fertilizer application dates, and the solid black line represents the September harvest of aboveground biomass. Data points represent the mean (±SE) of the three reps measured during each sampling date. *Significance from ANOVA (α = 0.05) of the fixed effects.
T A B L E 1 Differences in CO 2 fluxes (μmol m −2 s −1 ) between months. Values reported are modeled means (±SE) with fertilizer type and soil temperature controlled for as a fixed effect. In the fungal community, while Agaricomycetes increased in relative abundance in both crops, only tall wheatgrass plots exhibited significantly increased Glomeromycetes while the Sordariomycetes were enriched in the switchgrass plots (Figure 5d). A decrease in relative abundance in Tremellomycetes was observed in both crops. The fertilizer effects detected by ordination were primarily due to a higher relative abundance of Exobasidiomycetes, a class of broadly pathogenic fungi, in urea-fertilizer plots, particularly under the switchgrass. Furthermore, we detected a greater relative abundance of Ustilaginomycetes, another class of broadly pathogenic fungi, in switchgrass compared to tall wheatgrass (Table S7), particularly under urea fertilization. blast analysis of the Exobasidiomycete and Ustilaginomycete ZOTUs revealed that the species detected are all documented fungal pathogens of grass species (Table S8).

| Tall wheatgrass (C3) was more productive than switchgrass (C4)
Fertilized tall wheatgrass achieved high levels of AGB production in 2020 (6.7 Mg ha −1 ), rivaling the production of switchgrass and tall wheatgrass grown in warmer and wetter regions in the eastern UMRB ( Figure 1, Table S1) despite low soil fertility, modest N-fertilization application (50 kg N ha −1 ), and no irrigation. Switchgrass showed minimal potential as a bioenergy crop in the study area, regardless of fertilization. We contend that the large differences in aboveground production between switchgrass and tall wheatgrass likely reflect an important regional climate constraint on C4 species in the cool, early-season precipitation dominated, and semi-arid climate of the western UMRB and highlight the potential of C3 species for bioenergy production despite opportunities for C4 bioenergy production in future climates (Dolan et al., 2020).
In the semi-arid western UMRB, annual productivity in dryland cropping systems is dependent on precipitation inputs and adaption to a shorter growing season. Previous studies have shown that efficient use of rainwater is critical for bioenergy crop sustainability in semi-arid climates (Hendrickson et al., 2013;Peterson et al., 1996). Our finding of nearly three times more growing degree days and almost three times greater RUE for the C3 crop supports its potential as a bioenergy crop in this area. In addition to T A B L E 2 Mean growing season CO 2 , CH 4 , and N 2 O flux by crop type with fertilizer treatment controlled for as a fixed effect. Values are reported as modeled means (±SE) and a negative flux indicates soil or plant uptake. climate constraints, we acknowledge that other biophysical factors could have compounded the poor establishment of switchgrass. For example, under the switchgrass plots, we found a greater increase in the relative abundance of Ustilaginomycetes and Exobasidiomyctes, both pathogenic species, although disease symptoms were not assessed.

| CBF matches urea crop production and did not increase non-CO 2 GHG emissions
We determined that a fertilization with CBF achieved the same biomass production as urea N-fertilizer at equivalent application rates and importantly, the field use of CBF did not induce appreciable non-CO 2 GHG emissions compared to urea fertilization. Notably, CBF's lack of effect on N 2 O emissions is important as the global warming potential (GWP) of bioenergy systems is often a function of the N 2 O emissions since the 100-year GWP of N 2 O (298) is much higher than CH 4 (36) and CO 2 (1) (Lawrence et al., 2021;Robertson et al., 2000). In addition, the Nostoc species is capable of N-fixation, and it is possible that post-application N-fixation occurred, although unquantified in this study. Assessing biofertilizer N-fixation post-application is an important next step forward for this field of research. CH 4 is also a potent GHG, and quantifying how different management practices affect soil CH 4 fluxes is critical in assessing regional agricultural mitigation potential and calculating C budgets (Levine et al., 2011). In this study, CH 4 fluxes remained negative throughout the entire growing season across all crop and fertilizer treatments and were greatest under tall wheatgrass. This contrasts with Schmer et al. (2012) that reported positive CH 4 fluxes under switchgrass late in the growing season after harvest in September. Interestingly, we find a greater relative abundance of methanogens (CH 4 producing microbes) under switchgrass plots ( Figure S6), and this coupled with findings from Schmer et al. (2012) may warrant future exploration of upland CH 4 dynamics under switchgrass, especially for irrigated systems in the UMRB.
CO 2 fluxes drove GHG differences between crop and fertilizer treatments. The temporal trends in CO 2 flux between crops highlight plant phenology as a critical driver of C fluxes and an important consideration for evaluating seasonal GHG budgets in this system. Since the CO 2 measured in our chambers is a composite of autotrophic respiration and heterotrophic respiration, it is difficult to infer the relative contributions from each in these data. However, another study in a mixed grassland reported root respiration covarying with the season, ranging from 31% to 51% of total soil respiration (Wang et al., 2005). The greater CO 2 fluxes under tall wheatgrass from March to May likely reflect a larger fraction of autotrophic respiration from early root growth from the cold-adapted metabolism of C3 species (Hanson et al., 2000). The similar CO 2 fluxes between crops during June-October suggest a seasonal increase in the ratio of heterotrophic to autotrophic respiration under warming temperatures. The similar signatures of CO 2 flux from June to October suggest a larger fraction of heterotrophic CO 2 production under switchgrass given the substantially smaller plant inputs compared to tall wheatgrass. Taken together, the divergent pattern in CO 2 flux between crops early in the season indicates a strong autotrophic C3 signature that appears to shift to a dominant imprint of heterotrophic respiration as soil temperatures rise. This pattern is important for understanding the pattern of SOC loss under switchgrass although future experimental validation would be needed to confirm this.
Despite greater CO 2 emissions under the CBF fertilizer, rates were within published values for similar agronomic systems (Oertel et al., 2016), and we identified no effect of fertilizer on SOC over the 3 years. Possible explanations for greater CO 2 emissions under the CBF fertilized include: direct increases in respiration from the cyanobacterial inoculants, a decrease in microbial carbon use efficiency, or increased soil microbial respiration from the addition of a labile C substrate or a "priming effect" (Kuzyakov et al., 2000;Soong et al., 2020). Moving forward, partitioning soil respiration sources under both crop types and CBF will be critical for evaluating net ecosystem productivity and C budgets, especially heterotrophic respiration as this is a critical soil C loss vector. Furthermore, upstream production C costs are needed to develop a full carbon lifecycle assessment of CBF and assess its climate mitigation potential compared to conventional fertilizer practices.

C3 and loss under C4 crops partitioned by soil depths
The establishment of tall wheatgrass led to a net gain of SOC over 3 years, while the establishment of switchgrass led to a net loss of SOC irrespective of fertilization. While this finding is likely in part due to the lower C inputs from switchgrass it may be compounded by the fact that plant tissue from C4 plants decomposes faster than C3 plants (Ye & Hall, 2020). Patterns in SOC gain and loss were partitioned by soil depth in the different crops. There was a ~3 times greater increase in SOC under tall wheatgrass compared to switchgrass at the deeper 15-30 cm depth, while at the shallower 0-15 cm soil depth SOC losses were observed under switchgrass with no change under tall wheatgrass (Figure 3a,b). The strong directionality and magnitude of SOC loss under switchgrass suggest the two crops had differing effects on SOC in the topsoil. This pattern of change in SOC aligns with studies that found significant losses of SOC in topsoil, but not subsoil due to increases in labile root exudation that stimulated decomposition of older SOC through the priming effect (de Graaff et al., 2014;Salomé et al., 2010). Together, these findings align with the idea that different mechanisms govern SOC dynamics in the topsoil and subsoil (Fierer et al., 2003) and therefore, the magnitude and rate of soil carbon sequestration will vary with soil depth.
The rate of SOC accumulation across 0-30 cm depth under tall wheatgrass (0.25 Mg SOC ha −1 year −1 ) in this study is similar to SOC changes under perennial grasses in North Dakota (0.28 Mg SOC ha −1 year −1 ; Liebig et al., 2008). Our observed SOC accumulation rate is also similar to a previous study on a mix of alfalfa and three native grasses adjacent to this study (0.24 Mg SOC year −1 ; Engel et al., 2017). Since our study focused on bulk SOC, we cannot exclude that proportional fractions of particulate or mineral associated SOC were affected by crop type or N addition, which have important implications for the long-term stability of SOC (Cotrufo et al., 2019). Despite many observations that the variability of SOC is of greater importance than bulk density for detecting changes in SOC pool changes, our SOC pool results need to be interpreted carefully as we did not measure the 3-year change in bulk density ( Figure S9; Schrumpf et al., 2011). Important questions remain, including whether the SOC increases under tall wheatgrass were due to an increase in the mineral or particulate-associated carbon pools and how these crops influence the depth partitioning of SOC.

| Soil microbial shifts are crop specific and relate to plant productivity
Our results show that in both the fungal and bacterial/ archaeal communities, the establishment of the two perennial crops altered the soil microbiome composition and that crop type was a stronger driver of these shifts than the fertilization regime. In the fungal community, Glomeromycetes (AMF) increased in relative abundance at the 15-30 cm depth under the tall wheatgrass aligning with the observed increases in SOC at the same depth. AMF can contribute significantly to soil C stocks through the production of extraradical hyphal networks while also increasing N and phosphorus availability to roots in nutrient-limited (i.e., marginal) soils (Emery et al., 2018). Thus, our observation of increased root biomass under tall wheatgrass and greater AMF colonization could represent a positive feedback loop that led to increased accumulation of SOC possibly through mineral protection from enhanced aggregation (Wilson et al., 2009). However, this hypothesis is under the assumption that our measurements of root biomass at the 0-15 cm depth extend to the 15-30 cm and warrants future investigation and targeted experimentation to evaluate the role of AMF in tall wheatgrass productivity (Rillig et al., 2019). Switchgrass is highly dependent on mycorrhizal nutrient acquisition, so the lack of change in AMF relative abundance in the switchgrass plots further signals the imprint of low productivity or lack of their establishment (Dirks & Jackson, 2020;Hetrick et al., 2011). We also found relative increases in Agaricomycetes under the tall wheatgrass; these are a broad group of fungi that include taxa with lignin-degrading capabilities and have been found to be critical regulators of C flow into soil (Lodato et al., 2021). This increase is likely due to greater litter biomass inputs and the lower substrate quality compared to switchgrass (i.e., higher C:N ratio). Total fungal richness increased in both bioenergy crops irrespective of N fertilization type, which has previously been documented during perennial crop establishment due to decreases in soil disturbance from tillage (Bowles et al., 2017;Helgason et al., 2010;Schmidt et al., 2019).
In contrast to the fungal community, the bacterial/ archaeal community richness exhibited undetectable changes over the 3 years of this study. However, cropspecific changes in community composition were detected, including an increased relative abundance of Actinobacteria and Planctomycetes in the tall wheatgrass plots which play an important role in plant biomass and soil carbon turnover in low fertility soils (Bao et al., 2021;Hartman et al., 2017;Wieczorek et al., 2019). We also report an increase in the relative abundance of Thaumarchaeota in the switchgrass plots, all of which belonged to the genus Nitrososphaera. Nitrososphaera are ammonia-oxidizing archaea (AOA) that perform the rate-limiting step in nitrification, a critical loss pathway since it transforms the reduced form of N (NH 4 + ) to an oxidized and more mobile form (NO 3 − ). Interestingly, we observed the inverse response in Proteobacteria under switchgrass, which could be due to competition amongst ammonia-oxidizers as some members of the Beta-and Gammaproteobacteria are ammonia-oxidizers (AOB). Although N 2 O fluxes were extremely low in this study determining the factors governing AOA vs. AOB abundance in soils will be important for future studies as AOB have been observed to produce more N 2 O than AOA in some soil types (Banning et al., 2015;Prosser et al., 2020). Lastly, we observed a decrease in Gemmatimonadetes in the switchgrass and tall wheatgrass plots. Gemmatimonadetes are overall poorly characterized in soil systems, but their mirrored decline may signal they have a negative N fertilization response as observed by Cederlund et al. (2014).

| Conclusion
Our findings provide the first multi-year evaluation of second-generation BECCS perennial grass feedstock candidates in the western UMRB and establish a baseline for C3 and C4 yields for bioenergy candidate species for these cold semi-arid regions. However, our study is limited in that it was conducted at only one location within the UMRB. An important next step is to extend field trials to different locations in the URMB with different environmental conditions, previous management strategies, and original land uses. Despite this limitation, we conclude that the combination of tall wheatgrass (a cool-season C3) and a cyanobacteria biofertilizer showed strong potential as a productive bioenergy crop-fertilizer pairing as we observed greater overall productivity and fertilizer-N recovery in the C3 tall wheatgrass compared to the C4 switchgrass crop. CBF matched urea AGB production and did not stimulate increases in non-CO 2 GHG emissions compared to urea. The soil microbiome changed over time as we observed increases in fungal diversity and crop-specific shifts in bacterial/archaeal and fungal communities linked to the regulation of soil C turnover and nitrification. These microbial shifts may have contributed to the success of the C3 crop due to the emergence of a beneficial microbiome not observed under the C4 crop (e.g., AMF and pathogenic defense). Despite ongoing debate over the potential of agricultural soils to sequester C (Schlesinger, 2022), our data suggest that C3 perennial grasses may contribute to a climate mitigation strategy in cool semi-arid regions that is in part due to their ability to promote soil C accumulation. Finally, our findings are important for integrating and developing coupled microbialbiogeochemical models parameterized for agroecosystems (e.g., DayCent; Berardi et al., 2020) that may help in management decisions aimed at climate mitigation goals.