Strong growing season soil methane uptake observed in a restored prairie-savanna

Humid tallgrass prairies in the Upper Midwest have been largely converted to agriculture, creating large regional sources of greenhouse gases through soil oxidation and nutrient enrichment. Although prairie restoration is expanding, greenhouse gas fluxes from restored ecosystems remain poorly quantified. We measured growing season soil carbon dioxide, methane, and nitrous oxide fluxes across restored oak woodland, open tallgrass prairie, and wet prairie at Nachusa Grasslands, Illinois, and tested how fluxes were controlled by soil temperature, moisture, cover type, and elevation, and associations with surface soil microbial communities. Across all cover types, soils were consistent carbon dioxide sources (median 2.29 μmol m-2 s-1), methane sinks (median −0.36 nmol m-2 s-1), and neutral for nitrous oxide. Woodland soils were exceptionally strong methane sinks, while open prairie was often near neutral. Notably, net methane emission was absent even in wet prairie. Decision tree models explained 50%–70% of variability in carbon dioxide and methane fluxes, with soil temperature as the dominant control on carbon dioxide and cover type plus soil moisture most important for methane uptake. Inferred microbial community composition was dominated by putative aerobic taxa, including decomposers, nitrogen cyclers, and methane oxidizers. These results indicate that restored Midwest prairie-savanna ecosystems likely function as net greenhouse gas sinks, with particularly high rates of methane uptake in woodlands. Future work should incorporate non-growing season fluxes and vegetation fluxes to develop full ecosystem greenhouse gas budgets and better quantify the climate mitigation potential of prairie restoration.

Soils regulate the Earth’s climate by functioning as atmospheric sources or sinks of the three most important well-mixed greenhouse gases (GHG): carbon dioxide (CO 2), methane (CH 4), and nitrous oxide (N 2 O) (Singh et al., 2010). Soil GHG emissions in natural, well-drained upland soils are dominated by microbial respiration of CO 2 during decomposition of organic matter, and microbial CH 4 uptake, both via aerobic metabolisms, while typically having little to zero net N 2 O flux (Yang et al., 2021). As soils become wetter, CO 2 is still produced and emitted primarily from decomposition, but increasingly via anaerobic metabolisms, which also increases N 2 O emissions and, if soils are wet for long enough, fermentation and methanogenesis dominates decomposition leading to higher CH 4 emissions (McNicol et al., 2014). For each gas, upland and wetland soils combined accounts for a large share of their respective global GHG budgets (Friedlingstein et al., 2026; Saunois et al., 2025; Tian et al., 2020).

Though in situ GHG flux data coverage across U.S. ecosystems has expanded rapidly in the first 2 decades of the 21st century, exemplified by the growth of the AmeriFlux network (Chu et al., 2023), there are still relatively few data for all three GHGs in native prairie ecosystems, and even fewer for soil fluxes in humid Upper Midwest. The COSORE community database for continuous soil-atmosphere greenhouse gas flux data reports only 8% of 8.1M fluxes as from croplands or grasslands (Bond-Lamberty et al., 2020). Moreover, within restored or native U.S. tallgrass prairie, most GHG fluxes are concentrated at a few sites in the central and southern Great Plains, such as at the LTER and NEON sites at Konza Prairie Biological Station, and only recently expanded to cover the western Great Plains (Koyama et al., 2024).

More data are thus needed to fill observation gaps in the humid tallgrass prairie region of the United States, encompassing Illinois, Iowa, Missouri, and eastern Kansas. In northern tallgrass prairie regions, deep and carbon-rich Mollisols (Liu et al., 2012) underwent vast agricultural transformation in the 19th and 20th centuries, which led to the almost total loss (∼99%) of native ecosystem cover and contributes substantially to global atmospheric CO 2 accumulation via soil oxidation (Nuzzo, 1986; Samson and Knopf, 1994; Sanderman et al., 2017). Restoration of these ecosystems could therefore begin to rebuild native Mollisols and their large belowground carbon stocks (Matamala et al., 2008; Mielnick and Dugas, 2000). In addition to the mitigation of agricultural CH 4 sources, woodland restoration within the humid prairie-savanna ecotone could also strengthen CH 4 sinks (Dataur and Verchot, 2007). This may however be offset by net CH 4 emission in wet prairie in landscape depressions if saturated, low-redox conditions persist long enough to support methanogen growth (Jerman et al., 2009; Cui et al., 2024). A key uncertainty in this context is role of sedge radial oxygen loss in wet prairie, which can reduce or negate CH 4 emissions in other wetlands (Le Mer and Roger, 2021; Vroom et al., 2022). Finally, reductions or complete stoppage of fertilizer application and irrigation could mitigate or reverse large Corn Belt N 2 O emissions (Lawrence et al., 2021).

Characterizing associations between soil GHG fluxes and underlying key soil biogeochemical controls including temperature, moisture, and substrate availability is important to build accurate soil GHG flux models that extend to tallgrass prairie regions (Yang et al., 2021). Soil microbial community associations are increasingly valued to link redox-sensitive metabolic processes that produce and consume soil GHGs with microbial ecological compositional patterns and dynamics (Orr et al., 2025). Though establishing the link between microbial community composition and function is challenging, guild-based frameworks that aggregate taxa based on preferred or obligated metabolic pathway (e.g., methanogens vs. methanotrophs) may be able to simplify the challenge of high taxonomic diversity and partially obviate the problem of high functional redundancy (Hartman et al., 2024).

Here we begin to fill the observation gap for GHG fluxes of CO 2, CH 4, and N 2 O from the humid tallgrass prairie ecoregion of the U.S. in restored native oak woodland, open prairie, and wet prairie cover types. We test hypotheses that seasonality in temperature regulates the temporal pattern in GHG fluxes, while ecosystem cover type and moisture, mediated by elevation, regulate spatial variation in fluxes.

Site description

The study was conducted in the Tellabs West management unit at Nachusa Grasslands, a ∼4,000-acre Nature Conservancy preserve in Lee and Ogle Counties, northern Illinois (approx. 41.8944°N, 89.3793°W; Figure 1). The site lies within the Prairie-Forest Border Ecoregion at the northern boundary of Grand Prairie and mixed woodland/forests region to the north and features a mosaic of remnant and restored tallgrass prairie, oak savanna, and woodland on a complex of Mollisols (dominant in open prairie areas, including Nachusa silt loam series) and transitional soils on woodland knolls. The regional climate is humid continental, with mean annual temperature of approximately 8.5 °C–10 °C and mean annual precipitation of 940–1,016 mm based on nearby Dixon, IL NOAA station normals (1991–2020). Open prairie communities are dominated by warm-season grasses including big bluestem (Andropogon gerardii) and Indiangrass (Sorghastrum nutans), and diverse forbs including smooth blue aster (Symphyotrichum laeve), showy goldenrod (Solidago speciosa), and white wild indigo (Baptisia alba), while woodland knolls support open-grown bur oak (Quercus macrocarpa) with a prairie/savanna understory, maintained through prescribed fire and restoration management.

Study design

In 2022, we measured growing season GHG fluxes approximately monthly in July, August, September, and October in three native tallgrass prairie cover types that span a ∼30 m elevation gradient from oak woodland knolls (highest elevation), through open tallgrass prairie (mid-elevations), to depressional wet prairie (lowest elevation). Sites for each cover type were selected within adjacent woodlands and prairies. In open prairie, one wet praire (Wt-1), and woodlands, static soil flux chambers were installed evenly at the same elevation across a ∼20 m transect, while avoiding tree stumps and rocky outcrops, running perpendicular to the slope. In the lower elevation wet prairie cover site (Wt-2), given only slight elevation changes, the same overall design was used (3 × 4 chamber transects), but here micro-topography rather than a sloping gradient was captured. Flux chamber bases were installed in June, and soil was allowed to settle for 1 month before sampling began.

Soil flux chamber sampling

On each sampling date, rates of dark soil GHG fluxes were measured. Briefly, opaque chamber lids fitted with a small internal fan were placed over the chamber base. Vegetation was trimmed with scissors as needed to enable a good fit, and the connection was made air-tight using duct tape. Over 30 min, five 15 mL air samples were collected at 0, 5, 10, 20, and 30 min intervals after lid closure with a 20 mL syringe and over-pressurized in pre-evacuated 12 mL vials (Exetainers, Labco Limited, United Kingdom). Pressure changes were prevented with a narrow-diameter stainless steel pigtail tube fitting and air temperature within the chamber was measured at the end of the 30-min measurement period. This enabled use of the Ideal Gas Law to convert measured concentrations in parts-per-million volume (ppm(v)) into molar quantities for flux estimation. Alongside each sampling point, percent vegetation cover was visually estimated and recorded (5% bins from 0%–100%), a photograph of the soil and vegetation enclosed by the chamber was taken vertically from above the chamber to enable later diagnosis of unusual flux behavior, and 0–10 cm soil temperature and moisture were recorded with handheld probes.

Laboratory analyses

Gas samples were analyzed for CO 2, CH 4, and N 2 O concentration on a Shimadzu 2030 Ga Chromatograph (GC) fitted with thermal conductivity and electron-capture detectors housed at the University of Illinois Chicago. Replicate one-point reference standard gas (CO 2 991 ppm(v); CH 4 10.1 ppm(v); and N 2 O (1 ppm(v)) run alongside unknowns enabled as precision of at least 1%–3% for CO 2 and CH 4, and 5% for CH 4 and N 2 O. Non-linear electron-capture detector responses for N 2 O were not observed at the ambient atmospheric concentration range observed in this study, enabling the use of a one-point reference standard. Chromatograph peak heights were exported and converted to ppm(v) before being combined with measured chamber volumes and air temperatures to compute soil fluxes in standard units (CO 2, μmol m -2 s -1; CH 4 |N 2 O, nmol m -2 s -1).

Soil sampling and DNA extraction

At each soil sampling location, a small sample of soil was collected from the interior of an intact soil core collected at 0–10 cm, beneath the root and litter layer, using a stainless-steel spoon sterilized via flaming with 70% isopropyl alcohol. Sub-samples were gently homogenized by hand in sterile Whirl-Pak bags, transported on ice to the University of Illinois Chicago, and stored at −16 °C until processing. DNA was extracted from 0.25 g of soil using the DNeasy PowerLyzer PowerSoil Kit (Qiagen, Germantown, MD, USA) following the manufacturer’s protocol. Briefly, soil was added to a PowerBead Tube and subjected to mechanical lysis by bead beating. Inhibitors were removed, after which DNA was bound to a silica spin column. The column was washed, and DNA was eluted in a low-salt solution.

16S rRNA gene amplicon sequencing and bioinformatics

Amplicon libraries were generated from environmental DNA using the 515F/806R primer pair (with linkers) targeting the V4 region of the 16S rRNA gene (∼291 bp amplicon). Sequencing was performed at the University of Illinois Chicago Research Informatics Core on an Illumina MiniSeq platform using paired-end 2 × 150 bp chemistry. Library preparation and sequencing followed standard protocols for soil-derived amplicons, with basic data processing (including quality assessment) provided by the core. Raw 16S rRNA gene amplicon sequences have been deposited in the NCBI Sequence Read Archive (SRA) under BioProject accession number PRJNA1490676 (soil Targeted loci environmental).

Raw Illumina paired-end 16S rRNA gene amplicon sequences were processed using the DADA2 pipeline (v1.40.0; Callahan et al., 2016) in R. Forward (R1) and reverse (R2) reads were first inspected using plotQualityProfile(). Reads were then filtered and trimmed with filterAndTrim() using the parameters trimLeft = c(19, 20), truncLen = c(150, 120), maxN = 0, maxEE = c(2, 2), and truncQ = 2. Read depth and filtering metrics are provided in Supplementary Table S3. Due to consistently low merging success (<1% in most samples), only forward reads were used for downstream analysis, which is a common and accepted approach when reverse-read quality is insufficient for reliable merging (Callahan et al., 2016; DADA2 documentation). However, as there is a reduction in amplicon length without merging leaving fewer variable positions available for classification, we limit our main taxonomic analysis and inference to the genus and family level since these broader patterns are less sensitive to loss of resolution at the species level.

Error rates were learned separately for forward and reverse reads using learnErrors() with nbases = 1e8, multithread = TRUE, and randomize = TRUE. Sample inference was performed with dada() (pool = TRUE). Paired-end reads were merged using mergePairs(), but merging success was extremely low (<1% in most samples). Given the poor overlap and quality of reverse reads, downstream analyses proceeded using forward reads only, which maintained high quality and sufficient depth.

Chimeric sequences were removed using removeBimeraDenovo() with the consensus method. Taxonomy was assigned to amplicon sequence variants (ASVs) using IdTaxa() from the DECIPHER package (v2.30.0; Murali et al., 2018) with the SILVA SSU r138.2 training set (Quast et al. 2013). A phyloseq object (McMurdie and Holmes, 2013) was constructed containing the ASV table, taxonomy table, and sample metadata. All analyses were performed in R (v4.4+). The code for processing and analyzing the sequencing data are in the Github repo: https://github.com/McNicol-Lab/2022-nachusa-dna.

Functional guild analysis

To assess greenhouse gas (GHG)-relevant microbial guilds, ASVs were aggregated at the genus level and relative abundances were calculated. Taxa were assigned to six functional guilds based on established literature associations. CH 4 production (methanogens): Key methanogenic Archaea (e.g., Methanobacterium, Methanosaeta, Methanosarcina).

CH 4 consumption (methanotrophs): Aerobic methane-oxidizing bacteria (e.g., Methylococcus, Methylosinus, Methylocystis).

CO 2 aerobic: Aerobic heterotrophic and respiratory taxa (e.g., Mycobacterium, Pseudomonas, Sphingomonas).

CO 2 anaerobic: Anaerobic fermentative and syntrophic taxa (e.g., Clostridium, Bacteroides, Desulfovibrio).

N 2 O nitrification + denitrification: Nitrifiers (e.g., Nitrospira, Nitrosomonas, Nitrobacter, Denitratisoma).

N 2 O fixation: Nitrogen-fixing and ammonifying taxa (e.g., Bradyrhizobium, Rhizobium).

Guild relative abundances were visualized using stacked bar plots (‘ggplot2’ package) and summarized by sample.

Data analysis and statistics

All data processing, statistical analyses, and visualizations were conducted in the R statistical computing environment (R Core Team, 2025). Data wrangling, exploratory analysis, summary statistics, and figure generation were performed using the tidyverse collection of packages (Wickham et al., 2019). Linear models were fitted with base R functions; the lme4 package (Bates et al., 2015) was loaded for potential extensions to linear mixed-effects models. Decision tree (CART) analyses were implemented with the rpart package (Therneau and Atkinson, 2026) and visualized with rpart.plot (Milborrow, 2026).

Ecosystem codes in the raw 2022 Nachusa flux dataset were recoded into three functional cover types: open prairie (Pr; combining big prairie BP and small prairie SP sites), oak woodland (Wd; east woodland EW and west woodland WW sites), and wet prairie (Wt; north wetland NW and south wetland SW sites), with two replicate sites per cover type. Retained variables included sampling date, transect, replicate, gas type (CO 2, CH 4, N 2 O), flux value, p-value from the linear regression fit to chamber concentration time series, soil moisture (%), and soil temperature (°C). A final tidy dataset incorporating site geographic coordinates was exported for reproducibility (to be deposited in Dryad) and associated codes are stored on Github (https://github.com/McNicol-Lab/nachusa-2022-ghg).

Flux quality control was applied prior to analysis. If CO 2 concentration changes from a given chamber measurement resulted in a poor linear fit (p-value >0.05) then all three GHGs fluxes were set to missing (‘NA’), because a non-significant p-value indicated that the flux estimate was unreliable due to sample collection, handling, or GC analysis procedures. In addition, GHG fluxes from the same chamber measurements as anomalously high CO 2 fluxes (>9 μmol m -2 s -1) recorded at woodland sites were excluded (set to ‘NA’) because these outliers were attributed to observed ant-burrow openings rather than being representative baseline woodland soil respiration. When CO 2 concentration changes from a given chamber measurement resulted in a good linear fit (p-value <0.05) but CH 4 or N 2 O linear fits were poor (p-value >0.05) then CH 4 and N 2 O fluxes were set to zero, given no net flux is plausible, and the lowest analytically detectable fluxes were ∼0.6 nmol m -2 s -1 for CH 4 and ∼0.05 nmol m -2 s -1 for N 2 O. These quality control steps removed nine measured fluxes, and remaining complete-case data (n = 91–95) were used for descriptive statistics, modeling, and visualization.

Potential differences between the two replicate sites within each cover type were tested with ordinary least-squares linear models fitted separately for each gas (flux ∼ cover + factor(site)). Model summaries provided t-statistics and p-values for the site term, allowing evaluation of within-cover spatial variability.

To identify key predictors of flux without assuming linearity or normality, regression decision trees were grown independently for each gas using the rpart package (ANOVA method; Therneau and Atkinson, 2026). Predictors were soil temperature, soil moisture, and cover type. Trees were initialized with cp = 0.01, minsplit = 10, and minbucket = 5, then pruned to the complexity parameter that minimized cross-validated error (xerror). Performance was assessed with training-set R 2, which is the squared correlation of observed vs. predicted values, and 10-fold cross-validation R 2 (1 – minimum xerror), which iteratively trains only on 90% of the data, leaving out a distinct 10% of the data each time, then aggregates the observed vs. predicted value errors. Uncertainty was evaluated and is reported for each terminal node prediction using the within-node standard error. Relative variable importance scores were extracted from the pruned trees. For CO 2, temperature sensitivity (Q 10) was additionally quantified from flux ratios across temperature bins (10 °C–15 °C, 15 °C–20 °C, 20 °C–25 °C) and via linear regression of log-transformed binned fluxes against midpoint temperature.

GHG flux by cover type

Greenhouse gas fluxes for the prairie-savanna growing season pooled across all cover types were, for CO 2, all positive values and unimodal with a right skew, for CH 4, all negative values and bimodal with a left skew but with a main peak near zero, and, for N 2 O, spanned positive to negative values, and were unimodal around zero (Supplementary Figure S1). When flux data were parsed by prairie-savanna cover type (Figure 2), the distribution of CO 2 fluxes was bimodal in prairie and woodland cover types, and unimodal with a long right tail for woodlands. The distribution of CH 4 fluxes was bimodal with a main peak at zero for all cover types, and a longer left tail, and more negative secondary peak, for woodland cover. The distribution for N 2 O fluxes was unimodal and centered on zero for all cover types. Thus, as only N 2 O flux data followed a normal distribution when pooled at the prairie-savanna site-level or parsed by cover type, data is summarized as min./max. ranges, medians, and interquartile ranges (IQR) for all GHGs to ensure comparability (Table 1; Supplementary Table S1).

Across the prairie-savanna growing season, soils were consistent sources of CO 2 across all cover types (open tallgrass prairie (Pr), oak woodland (Wd), wet prairie (Wt)) with a median CO 2 flux of 2.29 (IQR 1.52–3.55) μmol m -2 s -1. When parsed by cover type, CO 2 emissions ranged from a min. of 0.43 μmol m -2 s -1 observed in Wd, to a max. of 6.87 μmol m -2 s -1 observed in Wt. Overall, CO 2 emissions were ∼30% lower on average in Wd (p < 0.001) than in Wt or Pr, which did not differ significantly (Table 1).

In contrast, soils were always either sinks (negative) or neutral (zero flux) with respect to CH 4 fluxes across all cover types during the growing season, with a median CH 4 flux of −0.36 (IQR -1.26–0) nmol m -2 s -1. When parsed by cover type, the highest CH 4 uptake rates of −2.59 nmol m -2 s -1 observed in Wd, while CH 4 neutrality (no uptake) was observed at times in all cover types. Overall prairie-savanna seasonal median CH 4 sink-strength was greatest in Wd at −1.46 (IQR -1.67 to −1.21) nmol m -2 s -1, followed by Wt at −0.43 (IQR -0.85–0) nmol m -2 s -1, while median Pr CH 4 flux was zero (p < 0.001).

All cover types, and therefore also the prairie-savanna ecosystem overall during the growing season, displayed neutral N 2 O fluxes with a median and IQR of 0 nmol m -2 s -1. Maximal N 2 O sink-strengths observed ranged from −0.39 nmol m -2 s -1 in Wd, to −0.19 nmol m -2 s -1 in Wt, to −0.07 nmol m -2 s -1 in Pr. Maximal N 2 O source-strengths showed the reverse pattern, with an observed range from 0.14 nmol m -2 s -1 in Pr, to 0.12 nmol m -2 s -1 in Wt, to 0.09 nmol m -2 s -1 in Wd. Thus, the overall pattern of N 2 O fluxes had a slight trend moving from a slight source in Pr (p < 0.05), to neutral in Wt and Wd, which did not differ significantly (Table 1).

Soil temperature and moisture variation

Soil temperature steadily declined (Figure 3, left panel) from an observed peak in August (site mean ± st. dev., 22.8 °C ± 2.7 °C) through September (site mean 16.6 °C ± 1.2 °C), to a minimum in October (site mean 12.0 °C ± 2.4 °C). Open tallgrass prairie (Pr) and Wd soils were slightly warmer (+1.3 to +1.4 °C) than Wt soils in August (Supplementary Table S2), while the same Pr and Wd soils were slightly cooler than Wt soils in September (−1.4 to −1.1 °C). In October, Wd soils remained warmer than Pr (+2.0 to +2.2 °C), likely due to the south-facing slope aspect and lower sun angle.

Soil moisture (Figure 3, right panel) reached a maximum across all cover types in September (site mean ± st. dev., 13.7% ± 9.9%) and was lower in August (11.5% ± 9.2%) and October (11.3% ± 10.2%). Wet prairie (Wt) soil moisture was higher throughout the sampling period (Supplementary Table S2) but also showed greater spatial (within-month) variability, peaking in September at 19.1% ± 13.8%. Open prairie (Pr) soil moisture was intermediate ranging from to 11.9% ± 3.8% in October to 15.7% ± 3.6% in September. Woodland (Wd) soil moisture was consistently lowest and ranged from 3.9% ± 0.8% in October to 6.4% ± 2.6% in September.

GHG flux modeling

Decision tree models in which growing season soil GHG flux was related to soil temperature and moisture, and cover type, explained 71% of the variation (measured as five-fold cross-validation R 2) in CO 2 fluxes and 54% of the variation in CH 4 fluxes but did not explain any variation in N 2 O fluxes (Table 2). Therefore, modeling results are only presented for CO 2 and CH 4 (Table 2; Figures 4, 5).

For CO 2, growing season emissions increased from low-to-higher soil temperature, except for Wd cover at the highest soil temperatures (Figure 4). Data were spread relatively evenly across five terminal (leaf) nodes in the fitted decision tree, which had four branching points. Node mean (±1 standard error) CO 2 fluxes ranged from 1.1 ± 0.04 μmol m -2 s -1 for the coolest soil temperature conditions (<15 °C; 30% of data), to 4.1 ± 0.12 μmol m -2 s -1 for the warmest soil temperature conditions outside Wd (>20 °C; 22% of data). Most observations were associated with two leaf nodes for intermediate soil temperatures (15 °C–20 °C; 39% of data) for fluxes ranging from 1.8–2.5 μmol m -2 s -1. Notably, in Wd soils, the highest soil temperatures (>20 °C) showed an unexpected decline in CO 2 emissions, producing an intermediate leaf node at 2.8 ± 0.19 μmol m -2 s -1 (8% of data). Soil temperature was by far the most important variable determining CO 2 flux (Table 2), though soil moisture and cover played a role, albeit 10-fold smaller than temperature.

For CH 4, growing season sink-strength increased from Pr and Wt cover types to Wd cover (Figure 5) and increased within Pr cover when moisture was low (<7%). Data were grouped less evenly across three terminal nodes resulting from only two branching points, from −1.3 ± 0.06 nmol m -2 s -1 for Wd soils in general (33% data), to −0.75 ± 0.10 nmol m -2 s -1 for drier conditions in Pr cover types (11% of data), to −0.17 ± 0.03 nmol m -2 s -1 in moister soil conditions for the same Pr cover types (56% of data). Cover type was the most important variable followed by moisture (67% importance compared to cover), and soil temperature was less important (10-fold smaller importance) (Table 2).

Soil microbial community survey

Total relative abundance of the 10 most abundance genera and families for each sample are shown in Supplementary Tables S4, S5 but are not reported in detail given the functional guild analysis approach. Microbial functional guild composition at the genus level varied slightly across the six open prairie cover type samples (Figure 6) and constituted on average 3%–5% of the total relative abundance (95%–97% were “Other”). Neither putative methanogen-associated or methanotroph-associated taxa were detected at the genus level. Taxa commonly associated with aerobic heterotrophy were the most abundant among assigned guilds, with Candidatus Solibacter and Mycobacterium consistently dominant among taxa identified. Anaerobic taxa were generally rare and dominated by Anaeromyxobacter. Nitrogen-cycling guilds showed the strongest site-specific patterns. Notably, denitrifying genera were not identified at any site. In samples from one woodland and one wet prairie site (Wd-1, Wt-1), Nitrospira, a genus commonly associated with nitrification, was relatively abundant. Samples from the same wet prairie (Wt-1) and one open prairie site (Pr-1) exhibited genera commonly associated with N-fixation (e.g., Bradyrhizobium). In contrast, samples from the second woodland and open prairie site (Wd-2 and Pr-2) displayed comparatively few assigned taxa from N-cycling guilds.

At the family level, microbial guild composition varied more substantially than by genera across the six prairie-savanna surface soil samples (Figure 7) and constituted 10%–22% of the total relative abundance (78%–90% were “Other”). Aerobic taxa dominated community composition in all samples, accounting for 30%–48% of the non-other total relative abundance, with the highest proportions observed in Wd-2 and Pr-2 samples. In contrast to the genus level analysis, families commonly assigned to methane-oxidizing guilds were detected at low relative abundance in drier woodland and prairie (Wd-1/2, and Pr-2) samples. Families often linked to detnitrification were identified, but again scarcely present, whereas families often linked to nitrification were relatively more abundant in Pr-1/2, Wd-1, and Wt-2. Families commonly associated with anaerobic respiration or fermentation (e.g., Geobacteraceae and Syntrophaceae), remained low across all samples, except Wt-2. No families commonly assigned to methanogenesis were detected.

Climate benefits of prairie restoration

Our results provide initial empirical evidence in support of the claim that tallgrass prairie restoration in the Upper Midwest could provide measurable climate benefits compared to agriculturally converted land uses, which are well-documented to be large regional-scale sources of the GHGs CO 2, CH 4, and N 2 O (e.g., Lawrence et al., 2021; Sanderman et al., 2017).

Median growing season soil respiration varied from 1.91 to 2.78 μmol m -2 s -1, which scales annually to 273–398 gC m -2 y -1 or 10.0–14.6 kgCO 2 ha -1 yr -1. Although we did not measure gross primary production, literature values for tallgrass prairie in the Midwest and Great Plains region range from 916 g C m -2 yr -1 (reconstructed prairie in central Iowa; Dold et al., 2017) to 1,372–1,789 g C m -2 yr -1 (native tallgrass prairie in the Southern Plains; Suyker and Verma, 2001; Bajgain et al., 2018). Even under conservative assumptions, these GPP rates substantially exceed measured soil respiration, suggesting that restored prairie systems in our study region are most likely to function as net CO 2 sinks, potentially strong sinks, consistent with reported positive soil carbon accrual rates of ∼43 gC m -2 yr -1 at the nearby Fermi National Accelerator Laboratory restoration (Matamala et al., 2008). Although we did not measure non-growing season soil CO 2 emissions, these are likely to be much lower than growing season emissions due to the strong temperature sensitivity of CO 2 production in soils (Davidson et al., 2012), and the lower autotrophic respiration rates from senesced vegetation.

The median rate of growing season CH 4 uptake in oak woodland soils was −1.46 (IQR -1.67 to −1.21) nmol m -2 s -1, which scales to an annual uptake of −0.56 (IQR -0.63 to −0.46) gC-CH 4 m -2 yr -1 or -7.5 (IQR -8.4 to −6.1) kgCH 4 ha -1 yr -1, assuming seasonal variability is low. The woodland soils therefore displayed higher rates of CH 4 uptake than any sites monitored in a recent Great Plains grasslands CH 4 flux survey across seven sites, in which maximum annual median rates of −5.4 (IQR -7.7 to −4.5) kgCH 4 ha -1 yr -1 were observed at a tallgrass prairie at Nine Mile, Nebraska (Koyama et al., 2023; Koyama et al., 2024). The global range in mean soil CH 4 uptake rates have been reported to span 1.1 and 1.5 kgCH 4 ha -1 yr -1 in desert and tundra ecosystems, respectively, to 3.3 and 5.7 kgCH 4 ha -1 yr -1 in tropical and temperate forests, respectively (Dutaur and Verchot, 2007). Thus, although upland forest soils are well understood to be the strongest sinks for CH 4 globally, measured rates for oak woodlands in this study are very high in comparison to published rates. We acknowledge, however, the limitation of our estimates being drawn only from growing season flux data, and that cooler, wetter, and snow-covered soil conditions during the non-growing season and winter in the Upper Midwest would likely reduce annual average uptake rates due to kinetic and thermodynamic constraints (Davidson et al., 2012).

Despite the large statistical uncertainties arising from extrapolations of our limited soil GHG flux dataset to larger spatio-temporal scales, the relatively high soil CH 4 uptake rates imply a potential climate benefit of oak woodland ecosystem restoration, which might apply to other sandy, well-drained oak woodland soils in the Midwest and Mississippi Valley. For instance, to assess a plausible upper bound on this potential, we can use the growing season median Wd-2 CH 4 uptake rate of −7.5 kgCH 4 ha -1 yr -1 and assume a 100-year Sustained Global Cooling Potential for CH 4 of 203 (Neubauer et al., 2015), yielding a CO 2 -equivalent GHG flux of −1.5 × 10 3 kgCO 2 e ha -1 yr -1. This potential GHG sink associated with oak woodland CH 4 uptake is thus an order of magnitude greater than upper estimates of soil respiration, underscoring a high likelihood of a net GHG sink for the prairie-savanna overall. Moreover, given soil and cover-type patterns observed in this study, the high CH 4 uptake rates are likely linked to the coarse textured and well-drained sandy soils on the oak woodland knolls formed from St. Peter Sandstone outcrops (Ketner et al., 1978), which covers 10 4 –10 5 ha of soils today across NW Illinois, SE Minnesota, and NE Iowa, and could support a substantial regional CH 4 sink.

In sum, our data supported a net climate benefit of native tallgrass prairie restoration in the humid prairie-savanna ecotone of the Upper Midwest, driven primarily by strong CH 4 uptake in woodland soils. Future work should add further non-growing season and multi-site flux observations of native and restored ecosystems, including vegetation CO 2 flux components, in the Upper Midwest prairie savanna to assess more accurately the potential climate benefits of restoration in this region.

Modeling of prairie-savanna GHG fluxes

Decision tree models were effective at explaining variability in GHG fluxes and captured non-linear interactions between cover type and soil conditions. Soil CO 2 emissions arise from the aerobic redox reactions that drive microbial respiration and are therefore almost universally sensitive to temperature via chemical reaction kinetics (Lloyd and Taylor, 1994). Our data and model results supported a strong non-linear positive temperature sensitivity of microbial respiration (Figure 4; Table 2). The soil temperate-based decision tree branches were defined by the temperate ranges of 10 °C–15 °C, 15 °C–20 °C, and 20 °C–25 °C and 91% of fluxes were observed between 10 °C–25 °C. Therefore, given terminal node CO 2 emission rates for each 5 °C bin, we estimate a Q 10 -equivalent of 3.6 for 12.5 °C–17.5 °C and a 3.23 for 17.5 °C–22.5 °C. Observed temperature sensitivity of soil respiration was therefore higher than commonly reported values, but is likely due to the confounding effect of declining grass and tree autotrophic respiration contributions from August to October. Key departures from this general trend with warmer soil temperatures, was the lower soil respiration rates at low soil moisture (<11%) across all cover types for temperatures <20 °C, as well as a 32% reduction in soil respiration in exceptionally warm (mean 29.0 °C) and dry (4.7%) soil conditions at mid-and-low elevation transects within Wd-1. Overall, our soil respiration findings from decision tree modeling are congruent with the Dual-Arrhenius Michael-Menten framework outlined in Davidson et al. (2012) whereby both temperature and moisture interact to regulate rates, while phenological effects on carbon substrate availability and autotrophic respiration can confound apparent sensitivity to these key drivers.

Decision tree models identified that woodland cover results in CH 4 uptake rates approximately four-fold greater than wet or upland prairie cover types at the study site, and these rates were exceptionally high when compared against other grassland and woodland observations (Figure 5; Table 2). Woodland cover at the site covers most higher elevation areas, specifically occupying knolls and outcrops underlain by sandstone-derived soils with limited understory. In other work focused on GHG flux regulation along catena-scale elevation gradients, lower moisture and higher oxygen availability on ridge and slope (as opposed to valley) landscape positions supported higher uptake of CH 4 (O’Connell et al., 2018). In accordance with these findings, the well-drained sandy soil on slopes and ridges of prairie-savanna woodland knolls may facilitate the exceptionally high CH 4 uptake rates in this study, underscoring the relatively under-studied geo-edaphic (rather that climatic) controls on microbially mediated soil GHG fluxes.

Microbial community associations

Defining GHG-related functional guilds related to CO 2 production (respiration type), CH 4 production and consumption, and N-cycling (Supplementary Tables S5, S6) enabled some key ecosystem cover-related insights to be distilled from the 599 genus-level and 357-family level observed amplicon sequence variants (ASVs). Given observed ASVs at the genus level, assigned communities were dominated by taxa commonly associated with aerobic heterotrophy (Figure 6), whereas taxa commonly associated with anaerobic metabolisms were rare and only slightly more abundant in prairie and wet-prairie samples. This compositional pattern is notable in light of the emerging anoxic microsite framework, which has documented anerobic metabolism, and even methanogenesis, in warm, moist, and finely textured soils (He et al. 2015; Keiluweit et al., 2018; Lacroix et al., 2023), including conditions similar to those underlying prairie cover in our study (McNicol et al., 2014; Angle et al., 2017). Despite soil moisture reaching 47.3% in wet prairie at ∼23 °C, no ASVs commonly assigned to methanogen taxa were detected at the family or genus level. One possible explanation is a time lag in the assembly of slow-growing methanogenic communities after restoration or wetting (Jerman et al., 2009; Cui et al., 2024). Another is that deep-rooted prairie vegetation may increase soil aeration through radial oxygen loss while also supplying substrate (Le Mer and Roger, 2001; Vroom et al., 2022). However, these mechanistic interpretations are conjectures that could be formulated as hypotheses. Taxonomic non-detection is not evidence of absence because ASV analysis of 16S sequencing data can miss rare or poorly classified taxa. Without functional-gene or metagenomic data, we cannot conclude methanotrophs or methanogens were absent, nor infer process rates from 16S composition. The flux data themselves showed consistent net CH 4 uptake and no net CH 4 emission. Together, the compositional and flux results point to a need to identify where field-to-landscape scale redox boundaries fall between aerobic, anaerobic, and methanogenic soil conditions, a transition previously described mainly at microsite to profile scales (Li et al. 2026; Yang et al., 2017; Keiluweit et al., 2018; Lacroix et al., 2023).

Our study provides initial empirical evidence that conservation and restoration of tallgrass prairie and oak savanna in the Upper Midwest could contribute to a potential climate benefit, primarily through very high rates of atmospheric CH 4 uptake in woodland soils developed on well-drained sandy knolls. Decision tree modeling revealed strong non-linear controls of temperature, moisture, and vegetation cover on GHG fluxes, while 16S microbial community composition was dominated by taxa commonly associated with aerobic heterotrophy, and no ASVs commonly assigned to methanogen taxa were detected. Compositional patterns therefore are consistent with the flux measurements but do not themselves demonstrate metabolic function. Overall, our findings indicate that geo-edaphic factors, particularly coarse-textured soils on sandstone outcrops, can create conditions for strong CH 4 oxidation that exceed those reported for most other grassland and woodland systems. At regional scales, restoration of prairie-savanna ecosystems on similar parent materials could generate substantial CO 2 -equivalent sinks. Future work should integrate direct measurements or remote sensing estimates of gross primary production and avoided agricultural emissions to fully quantify net GHG balances and inform potential restoration strategies for climate mitigation in the Upper Midwest.