Pathways for Methane Emissions and Oxidation that Influence the Net Carbon Balance of a Subtropical Cypress Swamp

We evaluated the major pathways for methane emissions from wetlands to the atmosphere at four wetland sites in the Big Cypress National Preserve in southwest Florida. Methane oxidation was estimated based on the δ13C-CH4 of surface water, porewater, and bubbles to evaluate mechanisms that limit surface water emissions. Spatially-scaled methane fluxes were then compared to organic carbon burial rates. The pathway with the lowest methane flux rate was diffusion from surface waters (3.50 ± 0.22 mmol m−2 d−1). Microbial activity in the surface water environment and/or shallow oxic sediment layer oxidized 26 ± 3% of the methane delivered from anerobic sediments to the surface waters. The highest rates of diffusion were observed at the site with the lowest extent of oxidation. Ebullition flux rates were 2.2 times greater than diffusion and more variable (7.79 ± 1.37 mmol m−2 d−1). Methane fluxes from non-inundated soils were 1.6 times greater (18.4 ± 5.14 mmol m−2 d−1) than combined surface water fluxes. Methane flux rates from cypress knees (emergent cypress tree root structures) were 3.7 and 2.3 times higher (42.0 ± 6.33 mmol m−2 d−1) than from surface water and soils, respectively. Cypress knee flux rates were highest at the wetland site with the highest porewater methane partial pressure, suggesting that the emergent root structures allow methane produced in anaerobic sediment layers to bypass oxidation in aerobic surface waters or shallow sediments. Scaled across the four wetlands, emissions from surface water diffusion, ebullition, non-inundated soils, and knees contributed to 14 ± 2%, 25 ± 6%, 34 ± 10%, and 26 ± 5% of total methane emissions, respectively. When considering only the three wetlands with cypress knees present, knee emissions contributed to 39 ± 5% of the total scaled methane emissions. Finally, the molar ratio of CH4 emissions to OC burial ranged from 0.03 to 0.14 in the wetland centers indicating that all four wetland sites are net sources of atmospheric warming potential on 20–100 yr timescales, but net sinks over longer time scales (500 yr) with the exception of one wetland site that was a net source even over 500 yr time scales.


INTRODUCTION
Inland aquatic ecosystems actively cycle carbon from the terrestrial biosphere (Battin et al., 2009), which results in large emissions of greenhouse gases (GHGs) such as carbon dioxide (CO 2 ) (Raymond et al., 2013;Sawakuchi et al., 2017) and methane (CH 4 ) (Baker-Blocker et al., 1977;Matthews and Fung, 1987) from inland waters to the atmosphere. Wetlands play a prominent role in the inland water carbon cycle due to their close connection between water, vegetation, and land. High rates of primary production by terrestrial and aquatic vegetation and low redox potential conditions in wetlands favor both organic carbon (OC) burial and methanogenesis, the relative balance of which regulate wetland impacts on the atmospheric GHG balance (Mitsch and Gosselink, 2007). Wetlands are the largest natural source of methane to the atmosphere with a global flux of 55-231 Tg C yr −1 (Houweling et al., 2000;Wuebbles and Hayhoe, 2002;Neef et al., 2010). The balance between carbon burial and GHG emissions in wetlands and aquatic ecosystems along the landocean continuum remain poorly constrained in the global carbon cycle and poorly represented in Earth system models Ward et al., 2020).
The role of CH 4 in this balance is particularly difficult and important to constrain because of high spatiotemporal variability in field observations and high global warming potential compared to CO 2 ; CH 4 has 96 and 11 times the warming potential of CO 2 when CH 4 emissions are sustained for 20 and 500 years, respectively (Neubauer and Megonigal, 2015). An evaluation of three wetlands ecosystems along a latitudinal gradient showed that the molar ratio of CH 4 emitted to net atmospheric CO 2 fixation by vegetation on an annual basis ranged from 0.05 to 0.20 (Whiting and Chanton, 2001). Considering the relative global warming potential of CH 4 , these wetlands exacerbate atmospheric greenhouse effect over short (20 years) time scales, but serve as net sinks of warming potential (i.e., produce a global cooling effect) over 100-500 years timescales since the atmospheric residence time of methane is an order of magnitude shorter than CO 2 (Whiting and Chanton, 2001).
Global estimates of CH 4 emissions from wetlands are poorly constrained because complex interactions between hydrology, vegetation, and soil/sediment properties result in emissions that can be orders of magnitude different even at short separation distances (Bridgman et al., 2013). Gaps in mechanistic understanding of how methane is produced, consumed, and emitted from wetlands further complicate predictions for how these ecosystems will respond to future climate and restoration scenarios (Riley et al., 2011). Fluctuating water table dynamics in wetlands means that saturated, inundated, and aerated zones exist and shift over time, with implications for where and when CH 4 is produced. However, general predictions about these dynamic controls on CH 4 emissions remains unclear because the functional relationship between water table height and CH 4 emissions is inconsistent across studies and wetland types. For example, CH 4 emissions respond differently to water table height in bogs, fens, and swamps while antecedent conditions (e.g., dry periods) appear to exert an important influence (Turetsky et al., 2014).
Wetland vegetation plays an important role in transporting CH 4 from belowground to the atmosphere because the soil and porewater that are the venue for plant roots contain high CH 4 levels (Schutz, 1991), and because vascular plant tissues are known to offer reduced resistance to atmospheric gas exchange. This phenomenon has been well-studied in herbaceous plants (Nouchi et al., 1990;Nisbet et al., 2009), but woody plants have recently been recognized as also playing an important role in transporting CH 4 to the atmosphere both in wetlands and upland forests (Barba et al., 2019a). In tropical wetlands, for example, emissions from tree stems accounted for up to 65% of total ecosystem CH 4 emissions (Pangala et al., 2017).
Despite the clear role of plant vascular systems in gas transport, few studies have considered some of the unique morphological features of highly flood tolerant trees. For example, flood tolerant trees of the genus Taxodium have "cypress knees" that visually appear similar to pneumatophores found in mangroves. However, the function of cypress knees remains unknown. Despite the ability for Taxodium to thrive in low oxygen soil environments compared to other genera (Pezeshki et al., 1996;Anderson and Pezeshki, 2000) there is little evidence that their emergent root structures provide oxygen to the root system. Another proposed function is that the knees provide stability in soft muddy soil (Briand, 2000). Considering these emergent structures provide a direct link between the anaerobic subsurface and the atmosphere, they may be an active, but largely unknown pathway for CH 4 emissions. Several studies have measured CH 4 flux rates from cypress knees, registering significantly higher evasion rates than nearby surface waters and soils (Pulliam, 1992;Bianchi et al., 1996). Likewise, pneumatophores have been shown to be conduits for methane emissions in mangrove forests (Purvaja et al., 2004). An experimental study showed that methane emissions from Taxodium distichum trees increased under both elevated ambient CO 2 levels and water table depth, suggesting that wetland methane emissions will increase under future climate scenarios (Vann and Megonigal, 2003). Despite their ubiquity in many coastal plain swamps, it remains unclear how important knees are for overall CH 4 emissions, particularly since this pathway would bypass microbial processes that might otherwise reduce methane export.
The microbial oxidation of methane (MOX) has been shown to dramatically reduce emissions from surface water environments such as rivers, in some cases by ∼99% . The magnitude and mechanisms for MOX, however, remain poorly constrained and are not frequently measured alongside emissions (Segers, 1998). Aerated surface waters are the most likely venue for MOX to occur, but MOX has also been attributed to processes occurring within herbaceous plants (Frenzel and Rudolph, 1998;Frenzel, 2000). If wetland surface waters play an important role in reducing CH 4 emissions similar to rivers and lakes (Bastviken et al., 2002), then emissions from vegetation may represent an important pathway to "bypass" this oxidative layer. A similar phenomenon has been observed in brackish coastal floodplains, where transpiration by trees allowed CH 4 produced deep in the soil profile (below shallow layers where sulfate inhibits methanogenesis) to bypass oxidation in the shallow soil layers .
Small aquatic ecosystems such as ponds and wetlands are increasingly recognized as exerting outsized influence on global aquatic GHG emissions (Holgerson and Raymond, 2016;DelSontro et al., 2018;Peacock et al., 2019), but settings with extensive organic matter turnover, such as cypress dome wetlands, remain poorly studied (Villa and Mitsch, 2014;Pereyra and Mitsch, 2018). Here, we examine methane cycling and carbon burial dynamics in the Big Cypress National Preserve (BICY) in southwest Florida, a subtropical wetland characterized by isolated forest wetland depressions (cypress domes) that are evenly distributed within short-hydroperiod marshes and pine uplands (Watts et al., 2014). The cypress domes in BICY are thought to form via feedback between water storage and carbonate rock dissolution, mediated and accelerated by organic matter cycling (Cohen et al., 2011;Dong et al., 2019). The resulting regular pattern of wetland sizes and arrangement  illustrates the importance of these feedbacks in creating geomorphic structure, a process thought to have occurred over the last 10,000 years (Chamberlin et al., 2018;Zhang et al., 2019).
Given the central importance of carbon cycling to the evolution of landscape depressions, our overarching goal was to enumerate carbon fluxes in these wetlands. The central relevance of methane fluxes in BICY emerges from a global understanding that wetlands are important control points for land-atmosphere methane exchanges, and from the local need to understand the hydrologic controls on organic matter cycling. Our objective was to quantify the components of CH 4 losses from small seasonal wetlands and evaluate how carbon burial balances with these CH4 flux pathways over 20-500 yr timescales. This required us to quantify the mass fluxes of methane in saturated and inundated settings, with varying water table conditions, to enumerate the mass fluxes through emergent root structures, and also to quantify mass retention due to water column oxidation. Measured methane emissions and existing carbon burial data (Zhang et al., 2019) was then scaled across the wetlands to evaluate the net radiative balance of the system. We hypothesized that the presence of oxic surface waters reduces overall wetland emissions, while emissions from cypress knees and ebullition bypass MOX and contribute significantly to wetland-scale emissions. We further hypothesized that the combination of these methane flux pathways will result in the system being a net source of atmospheric warming over decadal timescales.

Study Site
The Big Cypress National Preserve is a flat (mean surface slope 0.03 m km −1 ) karst landscape in southwestern Florida, and part of the Greater Everglades ecosystem ( Figures 1A-C). The combination of low relief and a humid subtropical climate with mean annual precipitation of 1.33 m (McPherson, 1974;Shoemaker et al., 2011) results in significant water storage on the landscape. Thousands of small (< 2 ha) karst depressions dot the landscape in a remarkable regularly patterned landscape of distributed water storage . These depressions have far thicker soils (Watts et al., 2014) with far higher OC content (Zhang et al., 2019) than the surrounding uplands because they hold water for much of the year, thereby supporting hydrophytic vegetation such as pond cypress (Taxodium ascendens) and a variety of wetland taxa (e.g., Salix caroliniana, Cladium jamaicense). These cypress "domes," so called because the cypress trees at the margins are typically shorter than those in the center due to growth rate differences (Katherine, 1988), are crucial elements in the complex mosaic that makes BICY an important biodiversity hotspot. The mosaic of uplands within which these depressions occur are dominated by cabbage palm (Sabal palmetto) and south Florida slash pine (Pinus elliottii var. densa) and, in stark contrast to the wetlands, are rarely inundated and lack any soil development.
Water levels in the wetlands follow pronounced seasonal patterns in rainfall and evapotranspiration; high water levels typically occur in the summer and fall, with drydown over the winter and the lowest water levels in late spring ( Figure 1D). Hydrologic connectivity in this landscape occurs when the wetlands fill sufficiently to exceed a critical stage threshold where water spills out and generates landscape sheetflow ; below this spill threshold, wetlands hold water but do not appear to exchange water regionally. While construction of elevated roads and canals for logging, oil drilling, and recreation over the last century has resulted in localized hydrologic modification (McPherson, 1974), the hydrologic and biological conditions in BICY are far less impacted than in the adjacent Everglades system.
We studied four cypress domes across two contrasting lithologies within BICY. Two were in the region known as Raccoon Point, near the end of 80.9279 W and RP3: 25.9144 N,80.9392 W) on the Pleistocene Fort Thompson Formation, a more recent carbonate lithology with a high proportion of insoluble material. The RP domes were shallower, and lacked surface water during the dry season. The other two domes were located near Turner River Road (TR2: 26.1197 N, 81.2691 W, andTR3: 26.0208 N, 81.2661 W), located on the Miocene Tamiami Formation. The TR domes were deeper, and each wetland center was inundated during the entire study period. All sampling described below was performed during daylight hours, generally between 10:00 and 15:00.

Above and Belowground CH 4 Concentrations and Stable Isotopic Composition
The partial pressure of CH 4 (pCH 4 ) and stable isotopic composition (δ 13 C-CH 4 ) was measured in surface waters, bubbles collected from the sediment-water interface, and two porewater depths at three locations per wetland site-the edge of the wetland near the upland-wetland boundary, midway into the wetland, and the center of the wetland ( Figure 1C and Figure 2). For each of these locations, nested porewater samplers were installed to the soil-bedrock interface and a depth halfway between the bedrock and sediment-water interface ( Table 1). The samplers were constructed of 5.1 cm diameter schedule 80 PVC and had 10 cm long sections with 125 µm well screening at their base. Water was pumped from each sampler through an overflow cup which contained calibrated (daily) YSI ProPlus pH, DO, specific conductivity, and temperature probes. Water was pumped until all parameters stabilized, to ensure pristine porewaters were sampled.
Sampling for measurements of pCH 4 and δ 13 C-CH 4 values included overfilling 635 ml polycarbonate bottles for about three bottle volumes ensuring they contained no bubbles. The bottles were closed with butyl stoppers equipped with two stop-cocks and short/long straws to allow creation of a headspace by removing 60 ml of water and injecting 60 ml of N 2 to fill the headspace without contact with the atmosphere . The bottle was shaken vigorously for 2 min and the headspace was extracted into a 60 ml syringe.
After collecting the above samples and measuring CH 4 fluxes (described below), bubbles were collected from the sedimentwater interface to measure non-microbially degraded δ 13 C-CH 4 values. Bubbles were collected by manually disturbing the sediment and capturing the released gas in an inverted funnel equipped with a 3-way luer valve . When the funnel filled with sediment bubbles, the gas was collected in 60 ml syringes. Prior to collecting the bubbles, care was taken to not disturb sediments as porewater were collected and flux measurements were made. δ 13 C-CH 4 data from the bubbles was used in our estimations of MOX (Estimations of CH 4 Oxidation).
Gas samples were injected from the syringe into a Picarro G2201-i Cavity Ring-Down Spectrometer (CRDS) within ∼12 h of sampling to measure pCH 4 and δ 13 C-CH 4 . Porewater and bubble samples were diluted 10-1,000 times with pure N 2 prior to analysis to achieve values within the CRDS's detection range. pCH 4 values were corrected for dilution with N 2 based on the  Table 3). The flux of CH 4 was measured using various opaque chambers depending on the source of the CH 4 flux ( Figure 2). For surface waters we used a "floating dome" design, which was an inverted high density polyethelene planter pot with polyethelene foam glued to its rim for flotation (Sawakuchi et al., 2017). The chamber was wrapped in tin foil to avoid heating of the headspace, which was sampled via an air-tight 2-way luer valve. The surface water chambers had a volume of 9.9 L, diameter of 0.33 m where the chamber meets the water, and surface area of 0.086 m 2 where the chamber meets the water. Five replicate chambers were deployed simultaneously for 5-6 min, after which the headspace was sampled using a 60 ml syringe. Ambient air was also sampled in the same manner. When surface waters were not present, e.g., at the edge of the wetland during certain sampling periods, we utilized triplicate soil flux chambers constructed from white PVC pipe with a sealed endcap on one end and air-tight 2-way luer valve. Chambers were gently placed on the soil surface for 4-5 min and samples were collected in 60 ml syringes. The soil chambers had a volume of 2.0 L, diameter of 0.10 m where the chamber meets the soil, and a surface area of 0.008 m 2 where the chamber meets the soil.

Frontiers in Earth
Another chamber, made out of the same PVC materials used for the soil chambers, was designed to fit over cypress knees protruding from the water surface ( Figure 2). At inundated sampling stations the chamber's opening was placed just below the water surface so as to not touch the sediment-water interface for 5 min prior to sampling with a 60 ml syringe. One chamber was deployed in triplicate over one knee per sampling location when present; no cypress trees or knees were present at TR3. In order to calculate how much of the chamber's volume was submerged under water and filled with the knee above the water-air interface, we measured how deep the chamber was submerged underwater, the knee height above the water surface, FIGURE 2 | The concentration and stable isotopic composition of methane was measured in surface waters, porewaters sampled near the bedrock interface and midway between the bedrock and sediment surface, and bubbles extracted from sediments. Methane fluxes were measured from surface waters, non-inundated soils when present, and cypress tree knees. These measurements were made in the center of the domes, the edge of the dome, and a midpoint between these locations. and the circumference at the bottom and top of the knee assuming the geometry of a tapered cylinder. For each chamber design, a Licor-820 with external diaphragm pump was interfaced to one chamber to monitor CO 2 fluxes (not reported here) and to ensure that the chamber's headspace was initially in equilibrium with the atmosphere and did not reach equilibrium with the surface waters/soils/knees during the deployment (i.e. gas fluxes remained linear). Water and air temperatures were measured using a YSI Pro Plus multiparameter sonde. Samples collected from the flux chambers were analyzed by direct injection on the CRDS as previously described. The rate of change ( dCH4 dt ; atm L −1 min −1 ) was calculated as follows: where pCH f 4 is the partial pressure of CH 4 inside the flux chamber at the end of the deployment (atm L −1 ), pCH i 4 is the partial pressure of CH 4 in ambient air collected with a syringe (atm L −1 ), and t is the length of time between placing the chamber on the water, soil, or knee and sampling the chamber's headspace (minutes). CH 4 fluxes (F CH4 ; μmol m −2 min −1 ) were then calculated based on the ideal gas law as follows: where, V is the volume of the chamber in L units, R is the universal gas constant (0.082057 L atm°K −1 mol −1 ), T is temperature in°K, and A is the surface area of the chamber opening in m 2 units. We then converted F CH4 values to mmol m 2 d −1 for reporting. The above flux calculations were applied to our three chamber types to represent: (Eq. 1) the total methane flux from surface waters (F T ), which include ebullitive and diffusive fluxes, (Eq. 2) the total methane flux from soils (F S ), and (Eq. 3) the total methane flux from cypress knees (F K ) to the atmosphere. We then calculated the relative contributions of diffusive vs. ebullitive fluxes (F D and F E , respectively) in the case of surface waters using gas transfer velocities similar to other river and lake studies (Bastviken et al., 2004Sawakuchi et al., 2014). To do this, we first converted surface water and atmospheric pCH 4 from μatm L −1 units to μmol L −1 concentrations ([CH 4 ] SW and [CH 4 ] ATM , respectively) based on the temperature dependent Henry's Law constant (K H ) (Wiesenburg and Guinasso, 1979): where T is temperature in°C units. Gas transfer velocity (k CH4 ) was then calculated in m d −1 units for each deployed chamber as: k CH4 was then converted to k 600 based on the following equation (Wanninkhof, 1992;Wanninkhof, 2014): where Sc CH4 is the temperature-dependent Schmidt number for methane calculated as follows: where T is temperature in°C. We then calculated the ratio of k 600 in each individual chamber to the minimum k 600 observed per deployment of five chambers (k 600 /k 600-min ) to determine which chambers captured bubbles from ebullition (Bastviken et al., 2004;Bastviken et al., 2010;Sawakuchi et al., 2014). For example, a high k 600 /k 600-min value indicates a large contribution of ebullitive vs. diffusive methane fluxes. The frequency distribution of k 600 /k 600-min for all chamber deployments indicates two distinct groups of ratios, whereby a ratio of 1.75 is the point of distinction (Figure 3). We chose the threshold value of 1.75 based on change point analysis (described in Data Analysis). Thus, we assumed that chambers with a k 600 / k 600-min greater than 1.75 experienced ebullition. We first calculated the average rate of CH 4 diffusion for chambers with a k 600 /k 600-min below 1.75. We then calculated ebullition as the difference between total and diffusive fluxes. For comparison across the study domain, we calculated the average ebullition rate for all five simultaneously deployed chambers, using a value of zero for chambers with a k 600 /k 600-min below 1.75.

Estimations of CH 4 Oxidation
Microbial oxidation of CH 4 in either oxic shallow sediment layers or surface waters results in an increase in δ 13 C-CH 4 values relative to its origin (Bastviken et al., 2002). Thus, we calculated the extent to which MOX limits diffusive surface water methane emissions by comparing surface water δ 13 C-CH 4 to both bubble and porewater δ 13 C-CH 4 . We applied two different models to assess the largest range of uncertainty. We used a model that represents an open system at steady state (Eq. 7) (Happell et al., 1994) and a second Rayleigh model for closed systems (Eq. 8)  to calculate the fraction of CH 4 oxidized in shallow sediments and surface waters (f): where δ SW is the δ 13 C-CH 4 of surface waters, δ b is the δ 13 C-CH 4 of bubbles and/or porewater, and a is the isotopic fractionation factor. Each calculation was made using two fractionation factor values-1.025 and 1.033-representative of the range of literature values available for similar systems (Tyler et al., 1997;Zhang et al., 2013). We also made the calculations using δ 13 C-CH 4 from bubbles and both porewater depths considering it was not always logistically possible to collect bubbles in shallow waters. We present MOX estimations as the average ± 1 SE of all calculations (i.e., two equations, two a values, and one or more belowground δ 13 C-CH 4 values).

Spatial Scaling of CH 4 Fluxes
We estimated total annual CH 4 emissions from each wetland in order to evaluate 1) the contribution of each flux pathway to wetland-scale CH 4 emissions and 2) the sustained global warming potential impact of CH 4 emissions compared to OC burial in the wetland sediments. The hydrology and topographic setting of the studied wetlands are extensively described by McLaughlin et al. (2019). In short, we measured water levels continuously for the entire study period using pressure transducers (Solinst Gold Levelogger) deployed in shallow groundwater wells situated in the deepest part of each wetland. Barometric pressure correction was performed using data collected from barometric pressure transducers (Solinst Barologger) deployed in a dry well adjacent to each water level recorder. LIDAR digital elevation models were provided by the National Center for Airborne Laser Mapping, processing of the LIDAR digital elevation models are described in detail in Quintero and Cohen (2019). Elevations are reported with an estimated mean accuracy of 5 cm across each domain, and with a spatial resolution of 5 m. The surface area of the center, intermediate, and edge sections of each wetland were determined based on elevation thresholds that aligned with shifts in vegetation coverage and inundation extent. Edge delineations capture the entirety of the wetland depression within the delineated extent while making sure to exclude any upland area. Intermediate delineations capture the bald/pond cypress extent of each wetland while making sure to exclude dwarf cypress and pine at higher elevations. Center delineations capture the center communities in our domes, these include open water areas, which typically differ from the dominant cypress community. Surface water CH 4 fluxes were multiplied by the percentage of the year that each section was inundated and soil CH 4 fluxes were multiplied by the percentage of the year that was dry. To scale the knee surface area we used an available literature value for the density of knees (0.315 knees m −2 ) in a similar cypress swamp setting since knee density was not measured in this study (Pulliam, 1992). Pulliam (1992) determined knee density across two 1,250 m transects in the frequently inundated lower floodplain of the Ogeechee River with 1 m 2 plots established every 8m (i.e., 156 plots per transect).We then multiplied this knee density by the average surface area of knees that we made flux measurements on (0.57 m 2 knee −1 ), and finally multiplied by the landscape surface area of each compartment. TR3 had no cypress trees or knees, whereas the other wetlands had knees in qualitatively similar abundance from the center to edge. Due to logistical constraints (i.e., time constraints), knee CH 4 flux measurements were not able to be made in each wetland section, thus, in those cases the average flux rate for a given wetland was applied to each of its sections.
Carbon burial rates were estimated for the center of each wetland based on previously reported 210Pb sediment accumulation rates, which were not measured for the intermediate and edge wetland sections (Zhang et al., 2019). 210Pb dating was performed down to 42 cm at the RP3 wetland and we applied this same rate to the other three wetlands. OC abundance in sediments (%OC) and bulk density were measured every 1 cm (Zhang et al., 2019). We used the average %OC and bulk density down to 42 cm to match the 210Pb time frame for carbon burial estimates. Carbon burial rates were determined as follows: OC Burial %OC x Bulk Density x Sediment Accumulation Rate x Surface Area The ratio of methane emissions to OC burial was calculated to compare their relative rates. This ratio was then multiplied by the sustained global warming potential of CH 4 for 20, 100, and 500 yr time frames (Neubauer and Megonigal, 2015) to determine the net balance in atmospheric radiative forcing from CH 4 emissions FIGURE 3 | The frequency distribution of the ratio of the observed gas transfer coefficient in a given chamber to the minimum gas transfer coefficient observed during a given set of chamber deployments (k 600 /k 600 Minimum) was used to differentiate between diffusive and ebullitive surface water methane emissions. The threshold of 1.75 was determined by statistical change point detection.
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 573357 and CO 2 uptake by vegetation and subsequent burial as OC in sediments.

Data Analysis
All statistical analyses were performed in the statistical computing language R using R Studio version 1.3.1093 (RStudio Team, 2020). Calculations of the k 600 /k 600-min threshold for ebullition to occur described in CH 4 fluxes From Surface Water, Soils, and Cypress Knees were performed using the "changepoint" R package. The changes in mean (cpt.mean) function was used with the At Most One Change (AMOC) method and the Modified Bayesian information criterion (MBIC) penalty. The change point of k 600 /k 600-min was determined to be 1.75 with a confidence value of 1.00. Field data was determined to not have a normal distribution based on visual inspection of density and QQ plots. Thus, we chose the non-parametric unpaired Wilcoxon test to evaluate the statistical significance of differences in measured parameters (e.g., pCH 4 , δ 13 C-CH 4 , CH 4 fluxes, and MOX) between wetland sites, sampling stations (i.e., center, intermediate, and edge), measurement types (e.g., soil vs. surface water fluxes), and sampling season. For example, we used the Wilcoxon test to determine if a given flux pathway was significantly higher at the site with the highest average flux rate compared to each of the other three sites individually. Thus, reported p values represent comparisons of one group to another (e.g., TR2 vs. RP1, TR2 vs. RP3, etc.). Significant differences were considered to fall within a 95% confidence interval (i.e., p < 0.05). In cases where all comparisons yielded insignificant differences (e.g., no wetland sites were statistically distinct from one another) we report the range of computed p values as (maximum p value < p > minimum p value) for the sake of brevity.
We performed linear regressions to determine correlations between measured parameters (i.e., pCH 4 vs. δ 13 C-CH 4 , fluxes vs. temperature, and fluxes vs. wetland stage) and report R 2 values.
All data is reported as the average ± standard error of the mean. This standard error was propagated through our estimations of total wetland methane emissions. The distribution and variability of field data is shown in box plots (Figures 4-6).
Porewater CH 4 concentrations ranged from 0.748 to 75.2 ppt with an average of 16.1 ± 1.85 ppt. The δ 13 C-CH 4 of porewaters ranged from −69.1‰ to −43.3‰ with an average of −57.0 ± 0.6‰ throughout the study period. Average pCH 4 in pore waters was 20 times higher than surface waters (p 0.00) and porewater δ 13 C-CH 4 was more 6.7‰ more depleted than surface waters (p 0.00) (Figure 4). Porewater collected at mid-depth had an average pCH 4 that was 1.1 times higher than at the bedrock interface, but this difference was not significant (p 0.35). Likewise, average porewater δ 13 C-CH 4 at mid-depth (-57.7 ± 0.7‰) was slightly more depleted than at the bedrock interface (-56.5 ± 0.9‰; p 0.28), but this difference was not significant.
The pCH 4 of bubbles collected from the sediment-surface water interface ranged from 55.0 to 1,000 ppt (i.e., pure CH 4 ) for all sites with an average value of 376 ± 71.1 ppt. The average δ 13 C-CH 4 of bubbles was −61.3 ± 0.68‰, which was 4.3% more depleted than porewaters (p 0.00). Because bubbles were not always able to be collected (e.g., when surface water level was too low for the sampling funnel), we did not attempt statistical comparisons between wetland sites, seasons, or sampling stations.
The average flux rate of CH 4 from soils (F S ) that were not inundated during the time of sampling (18.4 ± 5.14 mmol m −2 d −1 ) was 1.6 times higher than the total flux from surface waters, but this difference was not statistically significant (p 0.77). However, there was a more significant difference between F S and F D and F E when considered separately (p 0.06 and 0.00, respectively). F S was substantially higher in May 2016 (34.8 ± 8.31 mmol m −2 d −1 ) compared to both May 2015 (3.56 ± 1.01 mmol m −2 d −1 ; p 0.02) when water levels were substantially lower and October 2015 when only non-inundated soils were present at RP1 (2.85 ± 1.69 mmol m −2 d −1 ; p 0.29) ( Figure 5C) Figure 5A).

Microbial Oxidation of CH 4
On average, we estimate that 26 ± 3% of the CH 4 diffused from soil porewaters into surface waters is oxidized to CO 2 before evasion to atmosphere, resulting in an observed surface water diffusive CH 4 flux (i.e., F D ) that is 26% lower than the rate of CH 4 input to surface waters from sediments. The extent of MOX was higher in the TR wetlands than the RP wetlands ( Figure 6A); TR2 had the highest extent of MOX on average (40 ± 5%), which was 1.4 times higher than TR3 (28 ± 7%; p 0.46), 2.0 times higher than RP1 (20 ± 5%; p 0.01), and 2.7 times higher than RP3 (15 ± 7%; p 0.04). There was no significant variability between sampling period for the wetlands when all wetlands were considered together (1.00 < p > 0.37) ( Figure 6C). MOX was highest in the center sampling stations (35 ± 6%) compared to intermediate stations (25 ± 5%; p 0.17) and edge stations (18 ± 7%; p 0.11) ( Figure 6B), but these differences were not significant.

Spatially-Scaled CH 4 Fluxes and Carbon Burial
The low elevation center of each dome comprised 3-12% of the surface area of the studied wetlands. The intermediate and edge (i.e., mid to high elevation) sections comprised of 25-47% and 49-88% of each wetland, respectively ( Table 2). The center, intermediate, and edge components were non-inundated (i.e., no surface water) 9-17%, 15-27%, and 16-46% of the study period, respectively. The total flux of CH 4 from each wetland ranged from 239 ± 70 kg C yr −1 for the small RP3 wetland to 1741 ± 220 kg C yr −1 for the larger RP1 wetland. Diffusion and ebullition of methane from surface waters to the atmosphere contributed to 14 ± 2% (range 8-17%) and 25 ± 6% (range 12-34%) of the total spatially-scaled flux estimates from all four wetlands, respectively ( Table 2). Methane emissions from  non-inundated soils contributed to 34 ± 10% of the total spatiallyscaled flux estimates from all four wetlands, but had a much larger range (1-71%) for each individual wetland compared to surface water fluxes ( Table 2). Cypress knees contributed to 26 ± 5% of the total spatially-scaled flux estimates from all four wetlands. Knees were not present at TR3. When only considering the other three wetlands, knees accounted for 39 ± 5% of the total scaled methane fluxes with a range of 24-75% for each individual wetland. These totals do not include fluxes from the stems and leaves of the cypress trees, which is likely to increase the total. OC burial rates over the last century were estimated for the dome centers based on a sediment accretion rate of 0.45 cm yr −1 (Zhang et al., 2019). %OC averaged across the top 42 cm ranged from 13.5% to 19.1% for each wetland center and bulk density ranged from 1.20 to 1.37 g cm −3 ( Table 3). Carbon burial rates ranged from 160 kg C yr −1 for the small RP3 wetland to 1,263 kg C yr −1 for TR3. The molar ratio of CH 4 emissions to OC burial ranged from 0.03 to 0.14, or, in other words, spatially-scaled OC burial rates were 7-36 times greater than total spatially-scaled CH 4 emissions from the wetland centers. When ecosystem CH 4 emissions are sustained (as opposed to a one-time pulse), CH 4 has a sustained global warming potential of 96, 45, and 11 times that of CO 2 on 20 yr, 100 yr, and 500 yr timescales (Neubauer and Megonigal, 2015). Over the 20 and 100 yr time frames, the wetland centers of all four wetland sites were a net source of warming potential to the atmosphere; i.e., the methane they emitted to the atmosphere had a larger greenhouse gas impact than the CO 2 removed from the atmosphere by vegetation that gets buried in sediments as OC ( Table 3). The ratio of methane warming potential to carbon burial ranged from 2.7-13.3 and 1.3-6.3 for 20 yr and 100 yr time frames, respectively ( Table 3). Over the 500 yr time frame the wetland centers at three wetland sites were net sinks of warming potential (i.e., OC burial outweighed the impact of CH 4 emissions) with the ratio of methane warming potential to carbon burial ranging from 0.3-0.9. TR2 was a net source of warming potential even over 500 yr time scales with a warming potential ration of 1.5 (Table 3).

Pathways for Methane Emissions and Oxidation
Our results show that the total flux of CH 4 from cypress dome surface waters was not significantly different from non-inundated soil CH 4 fluxes on the wetland edge or during periods of low water levels ( Figure 5). Likewise, our observed methane flux rates from soils and surface waters (diffusion + ebullition) were nearly equivalent to other chamber-based measurements made in similar subtropical wetlands in Florida such as Corkscrew Swamp, FL (11.2 mmol m 2 d −1 compared to 11.3 mmol m 2 d −1 in this study) (Villa and Mitsch, 2014;Pereyra and Mitsch, 2018). 2 | Summary of surface area (SA), elevation delineations (Elev), the percentage of the year without inundation (Dry Pd), and spatially scaled total surface water methane fluxes (F T ), diffusive surface water fluxes (F D ), ebullitive surface water fluxes (F E ), methane fluxes from cypress knees (F K ), and methane fluxes from soils (F S ). The deep slough of the Corkscrew Swamp riverine cypress strand, which accumulates peat, is most similar to the cypress domes we studied, whereas the other four topographic settings studied by Villa and Mitsch (2014) were more similar to the unstudied upland settings surrounding our study domain and had an order of magnitude lower flux rates. Our average flux rates for surface water diffusion were also comparable to average rates reported for subtropical wetlands in a synthesis of global wetland studies, whereas our observed flux rates for total surface water emissions (diffusion + ebullition) was 2.8 times higher, soil emissions were 4.6 times higher, and knee emissions were 10.5 times higher (Turetsky et al., 2014). Only two subtropical wetlands were included in this synthesis, highlighting the paucity of comprehensive GHG data needed to constrain global wetland CH 4 budgets, which represent nearly half of all natural sources of CH 4 to the atmosphere (Saunois et al., 2020). Cypress knees were the most rapid pathway for methane release in our study with flux rates 3.7 and 2.3 times higher than total surface water and soil flux rates, respectively ( Figure 5). The highest knee flux rates were observed at the wetland site with the highest porewater pCH 4 (TR2) suggesting that subsurface methane content is an important driver of emissions through emergent root structures and perhaps also tree stems (Barba et al., 2019a), which were not studied here. Cypress knees represent a relatively small feature of the landscape in a two-dimensional plane, however, unlike the soil and water surface, gas exchange occurs from its entire three-dimensional structure. This attribute, along with the high measured flux rates, scales to represent a substantial percentage (39 ± 5%) of wetland CH 4 emissions at the three studied wetland sites with cypress trees present. This potentially important emission pathway has only been addressed by several previous studies. Observations made in the floodplain of a southeastern, USA blackwater river concluded that cypress knees contributed to less than 1% of the total methane emissions from the floodplain (Pulliam, 1992). The small contribution of cypress knees compared to floodplain surface CH 4 emissions in this case was likely due to substantially lower porewater CH 4 concentrations that resulted in average knee CH 4 flux rates that were ∼45 times lower than our measured rates. Knee CH 4 flux rates measured in the Sabine River floodplain in S.E. Texas were similar in magnitude to our measurements and it was reported that vegetation contributed to 64-96% of total ecosystem methane release (Bianchi et al., 1996). In that case the measured porewater CH 4 concentrations at depths greater than ∼10 cm were the same order of magnitude as our observations, whereas surface soils had nearly undetectable CH 4 levels, and surface waters had measurable rates of CH 4 oxidation.
Understanding the mechanisms that drive methane cycling through each pathway is a prerequisite for determining how effectively wetlands act as sources or sinks of GHGs and how they feedback and interact with changes to the Earth system. Water level is considered to be the primary factor that controls soil methane emissions in wetlands, particularly in tropical and subtropical wetlands that have wet and dry seasons (Whalen, 2005;Mitsch et al., 2010). Mitsch et al. (2010) observed that tropical wetlands with seasonal pulses in water level (e.g., monsoon cycles) had higher total annual methane emissions than permanently flooded wetlands, and that the highest seasonal emissions occurred when water levels were between 15-30 cm above the ground surface. This relationship to depth has also been observed in ombrotrophic peatlands and in both cases was attributed to enhanced MOX in shallow soils when water levels are low enough to allow soil aeration and enhanced MOX in surface waters when surface water depths are greater (Jauhiainen et al., 2005), which aligns with our observations of higher MOX in the wetland centers compared to edge ( Figure 6B).
The open water environments in BICY are small and shallow, which favor CH 4 ebullition due to low hydrostatic pressure and plant mediated fluxes from the predominantly anoxic porewater environments. MOX decreased the diffusive flux from open water by 26 ± 3% in the BICY wetlands ( Figure 6). However, the extent of MOX observed in our study area was considerably lower than reported in a large tropical floodplain lake in the Amazon (Barbosa et al., 2018), which found that in some cases up to 100% of the methane delivered from sediments to surface waters was oxidized. Future studies are needed to resolve the drivers of the large variability in MOX across different wetlands, and thereby better predict the conditions under which MOX significantly limits wetland methane emissions. Studies in rivers and lakes have identified some factors that influence the extent of MOX. For example, a lab incubation study showed that methane formation rates in lake sediments increased with increasing temperatures, whereas MOX potential was more closely related to CH 4 concentrations rather than temperature or sediment type (Duc et al., 2010). In rivers, NH 4 + concentrations above 20 μmol L −1 and light exposure have 3 | The percent organic carbon (%OC) and bulk density of sediments, OC burial rates, the molar ratio of total methane fluxes to OC burial, and net warming potential of the wetland centers from 20-500 yr timescales. Net warming potential values above 1 indicate the wetland is a net source of atmospheric warming and values below 1 indicate that the wetland is a net sink (i.e., has a net cooling effect). both been shown to inhibit MOX (de Angelis and Scranton, 1993;Dumestre et al., 1999), whereas MOX generally increases with higher suspended sediment concentrations (Weaver and Dugan, 1972). The high light exposure to the shallow waters at BICY, low suspended sediment concentrations (mean 5-7 mg L −1 ; , and historic mean NH 4 + concentrations of 13-20 μmol L −1 (Lietz, 2000) are all factors that may result in low MOX compared to other wetland systems.
Methane emissions from woody vegetation, which has only recently received significant attention (Barba et al., 2019a), are linked to several mechanisms relevant to our measurements. One mechanism for methane emissions from woody plants appears to be similar to herbaceous plants, whereby they take up soil-derived CH 4 by the root system (Covey and Megonigal, 2019). Stem methane emissions are not related to transpiration (Barba et al., 2019b), suggesting that diffusion through the plant's cells is the physical mechanism for methane exchange. For example, though plant adaptations such as aerenchyma and lenticels for increased O 2 diffusion in anoxic or water-logged environments can facilitate root-mediated transport of soil-produced methane (Pangala et al., 2013), transport of soil-derived CH 4 through tree stems has also been found in plants without these structures (Maier et al., 2018). Thus, regardless of whether or not cypress knees actively provide O 2 to roots in anoxic soil layers, diffusion of methane from soils to the atmosphere through knees is likely to occur, resulting in the high flux rates we observed. Methane can also be produced inside the biomass of woody vegetation by methanogenic archaea (Zeikus and Ward, 1974;Yip et al., 2019). This phenomenon is typically attributed to trees in upland forests with high CH 4 flux rates even when soils exhibit net uptake of methane (Barba et al., 2019b), though it is certainly possible that trees in wetland environments both transport methane from soils and produce it internally depending on the redox state and microbiome of woody biomass.
Woody vegetation plays an important role in CH 4 emissions in the case of environments with elevated porewater concentrations and high amounts of oxidation in shallow soils and/or surface waters. For example, average porewater pCH 4 was 20 times greater than surface water pCH 4 in our study (Figure 4). The roots of trees, which sample water and solutes from deeper soils than herbaceous plants, act as a direct conduit for methane to vent from deep soil porewaters with limited exposure to settings where oxidation occurs. On the other hand, O 2 may enter sediments from tree roots, oxidizing the rhizosphere and enabling MOX to occur deeper in the sediment than one might expect. This may explain several of the relatively high porewater δ 13 C-CH 4 values observed (i.e., maximum of −43.3‰), while the average porewater δ 13 C-CH 4 values (−57.0 ± 0.6‰) fall in the range of acetoclastic rather than hydrogenotrophic methanogenesis. Tree stems may contribute even more than knees to whole ecosystem methane fluxes because of their greater surface area and height, though few studies have shown how CH 4 fluxes change along the radial and vertical axes of a stem (Jeffrey et al., 2020). Tree-mediated methane emissions contribute 9-27% and 44-65% of the total ecosystem CH 4 flux in forested temperate and tropical wetlands, respectively (Pangala et al., 2017). Similar to methane emissions, the soil depths where methane oxidation can occur are linked to water level (Whalen and Reeburgh, 2000). Thus, the extent of wetland tree CH 4 emissions are likely linked to the interplay between rooting depth, water table depth, and redox zonation along the soil profile.

Comparison of Wetland Methane Emissions and Organic Carbon Burial
Are freshwater wetlands net sources or sinks of atmospheric GHGs and over what timescales? This critical question has been addressed in previous studies using the molar ratio between wetland CH 4 emissions and CO 2 uptake. These values ranged from 0.05-0.06 for two subtropical wetlands, 0.09-0.11 for two temperate wetlands, and 0.13-0.20 for three high-latitude wetlands (Whiting and Chanton, 2001). Our estimates were in the same range as subtropical wetlands studied by Whiting and Chanton (2001) (0.04-0.08) with the exception of TR2 (0.14) ( Table 3), which had the highest knee flux rates and highest rates of ebullition from sediments ( Figure 5A). Whiting and Chanton (2001) used global warming potential values of CH 4 based on the lifetime of a single pulse of CH 4 in the atmosphere (20 yr 21.8, 100 yr 7.6, 500 yr 2.5) (Lelieveld et al., 1993) to reach the conclusion that all of their studied wetlands were net sinks of warming potential over a 500 yr time frame and the subtropical wetlands were also net sinks over a 100 yr time frame. We used updated global warming potentials based on a sustained flux of methane from an ecosystem, which are considerably higher (20 yr 96, 100 yr 45, 500 yr 11; (Neubauer and Megonigal, 2015) to reach the conclusion that, over 100 yr timescales, the subtropical wetlands we studied are net sources of warming potential in contrast to prior studies.
There is not a clear consensus on what timescales are most important to consider with respect to global warming potential (Myhre et al., 2013;Kirschbaum, 2014). The 100 yr time frame represents the lifetime of CO 2 in the atmosphere-an order of magnitude longer than the residence time of CH 4 -and the time it generally takes terrestrial ecosystems to reach successional steady state (Odum, 2014). The 500 years time frame is similar to the time period over which ocean circulation exerts an important influence on climate change (Lelieveld et al., 1998). For context, the BICY patterned landscape appears to have begun forming in the middle to late Holocene (Chamberlin et al., 2018) and the dominant vegetation shifted from herbaceous to woody vegetation ∼3,000 years before present as precipitation and hydroperiods increased (Zhang et al., 2019). Radiocarbon dating of long-chain fatty acids showed that sediment accretion rates in the BICY wetlands were an order of magnitude slower from 1,600 to 3,500 years before present (0.03 cm yr −1 vs. 0.45 cm yr −1 present-day) (Zhang et al., 2019). Although we have no way of measuring past methane emissions, it is possible that carbon burial and methane emission rates have maintained a similar relative balance Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 573357 throughout the various stages of the wetland system's development, i.e., as carbon burial increases, so do CH 4 emissions. On the other hand, if woody vegetation is in fact a dominant pathway for deep soil methane to bypass oxidation, the wetlands may have shifted from a net sink to net source of greenhouse potential after the vegetation shift towards woody species around 3,000 years ago. Results from this study build upon a growing body of work demonstrating the importance of wetlands for methane emissions and thus their contributions to climate change (Zhang et al., 2017). Model simulations suggest that natural wetland methane emissions may be a positive feedback for warming. This feedback may be amplifying anthropogenic radiative forcing by up to 5% by 2100, similar to increases in anthropogenic methane emissions over the same time scales (Gedney, 2004). Our results further demonstrate that emergent vegetation is an important flow path for releasing methane from wetland soils as it limits the amount of methane oxidized in surface waters and shallow soils. While the role of herbaceous plants on wetland methane emissions is represented in current global scale Earth system models, the role of woody vegetation is crudely parameterized (Riley et al., 2011). Adding mechanistic detail to predictive models is paramount for fully constraining the role of wetlands as a positive feedback for warming and, likewise, predicting the impact of natural and anthropogenic changes to hydrological regimes on wetland carbon cycling.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession numbers can be found below: https://doi.org/10. 6084/m9.figshare.c.4791546.

AUTHOR CONTRIBUTIONS
TB, JM, and MC led the Big Cypress project and oversaw field and lab activities. CQ organized field logistics and NW performed methane measurements and sample collections. NW and HS performed calculations and data analyses. All authors discussed the results and commented on the manuscript.

FUNDING
This research was funded by National Science Foundation (DEB#1354783). The Jon and Beverly Thompson Endowed Chair of Geological Sciences to TB also provided support.