Come Rain, Come Shine: Peatland Carbon Dynamics Shift Under Extreme Precipitation

Precipitation patterns are becoming increasingly extreme, particularly at northern latitudes. Current climate models predict that this trend will continue in the future. While droughts have been repeatedly studied in many ecosystems over the last decades, the consequences of increasingly intense, but less frequent rainfall events, on carbon (C) cycling are not well understood. At northern latitudes, peatlands store one third of the terrestrial carbon and their functioning is highly dependent on water. Shifts in rainfall regimes could disrupt peatland C dynamics and speed-up the rates of C loss. How will these immense stocks of C be able to withstand and recover from extreme rainfall? We tested the resistance and resilience effects of extreme precipitation regimes on peatland carbon dioxide (CO2) and methane (CH4) fluxes, pore water dissolved organic carbon (DOC) and litter decomposition rates by exposing intact peat cores to extreme, spring-time rainfall patterns in a controlled mesocosm experiment. We find that more intense but less frequent rainfall destabilized water table dynamics, with cascading effects on peatland C fluxes. Decomposition and respiration rates increased with a deeper mean water table depth (WTD) and larger WTD fluctuations. We observed similar patterns for CO2 uptake, which were likely mediated by improved vascular plant performance. After a three-week recovery period, CO2 fluxes still displayed responses to the earlier WTD dynamics, suggesting lagged effects of precipitation regime shifts. Furthermore, we found that CH4 emissions decreased with deeper mean WTD, but this showed a high resilience once WTD dynamics stabilised. Not only do our results illustrate that shifting rainfall patterns translate in altered WTD dynamics and, consequentially, influence C fluxes, they also demonstrate that exposure to altered rainfall early in the growing season can have lasting effects on CO2 exchange. Even though the increased CO2 assimilation under extreme precipitation patterns signals peatland resistance under changing climatic conditions, it may instead mark the onset of vascular plant encroachment and the associated C loss.


INTRODUCTION
Climatic change includes shifting precipitation patterns along with increases in temperature and drought (Dai, 2013). Moreover, current observations of increased heavy precipitation confirm earlier predictions (Fischer and Knutti, 2016). In the past three decades, extreme precipitation events have become more common (Lehmann et al., 2015), particularly in North America and Europe (Ummenhofer and Meehl, 2017). As a consequence of more intense but less frequent rainfall events in northern temperate to arctic regions, local carbon (C) dynamics are expected to change (Blodau and Moore, 2003;Gerten et al., 2008;Knapp et al., 2008;Frank et al., 2015;Dai et al., 2018). Heavy rain might temporarily invoke waterlogged and anoxic soil conditions and thus limit processes that require oxygen such as respiration (by plants and soil biota) and decomposition (Knapp et al., 2008;Reichstein et al., 2013) and promote methanogenesis. Conversely, increasingly long intervals between heavy rain events could lead to drought stress, affecting both plant and soil microbial physiology (Zeppel et al., 2014;Schimel, 2018), as well as related ecosystem functions, including photosynthesis, respiration, and methanogenesis (Knapp et al., 2002;Reichstein et al., 2013;Zeppel et al., 2014;Helfter et al., 2015). While the effects of drought and flooding are reasonably well understood, the combined impact of more intense rain and increased duration of dry intervals on ecosystem C dynamics is little researched.
At northern latitudes, peatlands are one of the dominant ecosystems (Leifeld and Menichetti, 2018) and their C dynamics heavily depend on their hydrology (Limpens et al., 2008). Predominant shallow water tables lead to C sequestration because of an imbalance of C entering the system through photosynthesis, while little assimilated C is lost because of repressed decomposition and respiration (Dise 2009;Turetsky et al., 2012). As a result peatlands store the equivalent of 25-50% of the global soil organic C stock (i.e., 500-1000 GtC) (Yu, 2012;Wieder et al., 2013;Jackson et al., 2017;Nichols and Peteet, 2019). Peatland C stocks are vulnerable to ongoing climatic change (Dise 2009) and especially to hydrological shifts (Swindles et al., 2019). Hydrological shifts such as caused by extreme rainfall patterns or drought threaten peatland C dynamics (Radu and Duval, 2018b;Evans et al., 2021), pushing peatlands toward becoming a C source by increasing decomposition rates (Belyea, 1996;Moore et al., 2007) and CO 2 respiration (Lund et al., 2012;Helfter et al., 2015) and reducing CO 2 uptake (Kuiper et al., 2014;Helfter et al., 2015). As a consequence of extreme drought, methanogenesis may decrease while methanotrophy increases, thus reducing CH 4 emissions Turetsky et al., 2014;Evans et al., 2021).
Hydrological shifts do not always translate into peatland C loss Muhr et al., 2011;Nijp et al., 2014;Wang et al., 2015). Small but frequent rain events can be readily exploited by the crucial peat-forming Sphagnum mosses, maintaining the moisture level at the peat surface sufficiently for CO 2 uptake Strack and Price, 2009;Adkinson and Humphreys, 2011;Nijp et al., 2014). The influence of rain is particularly relevant when water tables are low (Strack and Price, 2009). Sphagnum mosses are prone to drying because they lack vascular tissue and roots to cope with dry conditions. Then again, Sphagnum photosynthesis can recover within 1 h when rewetted after a long period of drought (Jassey and Signarbieux, 2019). Yet, this quick recovery of Sphagnum C assimilation might not be sufficient to off-set autotrophic respiration costs of restarting photosynthesis (McNeil and Waddington, 2003;Nijp et al., 2014). At the same time, under dry conditions vascular plants may take-over the C sink function from Sphagnum (Huth et al., 2018;Schwieger et al., 2020;Radu and Duval, 2018b), making the effects of shifting rainfall regimes on net ecosystem exchange harder to predict. Influence of precipitation shifts on peatland C dynamics may also act through dissolved organic carbon (DOC). Flooding and heavy rains can result in a temporary increase in DOC concentration (Moore and Dalva, 2001;Blodau and Moore, 2003) and DOC loss by run-off (Pastor et al., 2003;Leach et al., 2016;Rosset et al., 2019). Furthermore, recent evidence showed that prolonged drought has little effect on the concentration of DOC, although lowering of the water table decreased DOC aromaticity and/or lower its molecular weight (Jassey et al., 2018). These observations highlight the uncertainty and the need for a better understanding of peatland C dynamics responses to more extreme precipitation patterns.
Besides its hydrological state, seasonal timing also controls an ecosystem's response to precipitation shifts. Indeed, disturbances to vegetation phenology can have cascading effects on peatland C uptake (Helfter et al., 2015;Koebsch et al., 2020). While repeated drought stress early in the growing season may lead to a loss of resilience regarding CO 2 uptake, its effect on other C fluxes, such as decomposition, DOC, respiration, and CH 4 efflux, is little studied. This underlines the need to investigate the influence of redistribution of rainfall early in the growing season on peatland C dynamics, which may be substantial given the hydrological control on C uptake, turn-over and loss. Moreover, effects of extreme rainfall patterns on peatland C dynamics in spring could have immediate as well as longer term consequences when the response mechanisms include lasting physiological and community changes among plants and microorganisms (Frank et al., 2015;Schimel 2018).
Here, we evaluated the capacity of peatland C dynamics to tolerate severe spring rainfall redistribution. We studied how peatland C dynamics continue to function during (i.e., resistance) and "bounce back" from (i.e., resilience) increasingly extreme precipitation regimes during spring (sensu Nimmo et al., 2015;Ingrisch and Bahn 2018), in order to understand how C storage in peatlands may be affected by prospective climatic change early in the growing season. For wet systems, such as peatlands, Knapp et al. (2008) postulated that low frequent but intense rainfall events will reduce the number of days during which the system experiences anoxic conditions. That is why we expected that 1) under more intense but less frequent rainfall C cycling rates speed up with increased assimilation, decomposition, and respiration but reduced CH 4 emissions as caused by larger fluctuations of the water table over time. Moreover, 2) peatland C dynamics are expected to be resilient and C dynamics recover after rewetting without lasting effects. We tested our hypotheses in a mesocosm experiment in which intact peat cores were exposed to increasingly intense, but less frequent rainfall for nearly 3 months, followed by rewetting and a three-week recovery period.

Study Site
To construct the experimental mesocosms, peat monoliths with homogeneous vegetation cover were extracted from a poor fen near Counozouls in the French Pyrenees (42 o 41′16 N; 2 o 14′18E, altitude 1,350 m, 7.9°C mean annual temperature, 1,027 mm mean annual precipitation). The site is situated within the Special Area of Conservation Natura 2,000 site "Massif du Madres Coronat." The vegetation has a moss layer dominated by Sphagnum warnstorfii Russow and S. palustre L., the vascular vegetation was characterized by Molinia caerulea (L.) Moench, Eriophorum angustifolium L., Calluna vulgaris (L.) Hull, Potentilla anglica Laich. and Drosera rotundifolia L. Intact peat monoliths (14 cm h, 27 cm l, 17 cm d) were extracted in November 2019 and climatized for 3 months in a greenhouse in watertight containers, during which time the water table was kept at approximately−4 cm below the peat surface by watering regularly (for timeline see Table 1).
Meteorological field observations at the Counozouls site served as a baseline for the precipitation regimes (Supplementary Figure S1). In spring 2019 (March-May) 248 mm of rain fell at the site, varying from 0 to 31 mm a day, while in summer (June-Aug) 177 mm of rain fell, varying from 0 to 33 mm a day. The longest consecutive period without rain was 12 days during spring and 11 days during the summer months. We further calculated the rain frequency (% of days with rain), rain intensity (average mm of rain on a rainy day), and drought intensity (average length of dry spell in days) to characterize the rainfall regime at the field site during spring and summer (Figure 1).

Experimental Design
A week before the start of the experiment, intact cores were taken from the acclimated peat monoliths to construct 30 mesocosms in watertight plant pots (13 cm d, 9.5 cm h). The depth of the mesocosms covered the actively growing Sphagnum capitula, recently senesced Sphagnum stems and approximately 5 cm bulk peat including vascular plant roots. For 11 weeks the mesocosms were exposed to experimental precipitation regimes. Each pot was randomly assigned to one of six regimes. The precipitation regimes entailed the same amount of precipitation over the 11 weeks, but administered in increasingly extreme patterns sensu Knapp et al. (2008), namely 133-187 mm except for the most extreme regime (regime 6) which received 106-160 mm (Supplementary Figure S2). Every regime included 8-14 drought intervals (6-7 for regime 6) that intermitted the rain events. The regimes became increasingly extreme in three aspects: increase in rainfall intensity, reduction in rainfall frequency and consequentially lengthening of dry intervals between rain events ( Figure 1). Dry intervals in the least extreme regime lasted maximum 5 days and increased to maximum 29 days in the most extreme regime. In all this, regime 2 most closely resembled precipitation patterns observed at the field site during spring. In order to mimic natural variation, the five replicates within a regime varied slightly in the amount of precipitation per rain event and the exact length of a drought spell, but the order of the drought spells was synchronized to prevent pseudo-replication. The mesocosms were watered according to the precipitation treatments with artificial rain water based on the protocol by Garrels and Christ (1965), with addition of micronutrients, see Supplementary Table S1.
The resistance and resilience responses to the experimental precipitation regimes were observed during two experimental phases ( Table 1). The precipitation regimes were applied during the resistance phase (I) which lasted 11 weeks. The resilience phase (II) immediately followed the resistance phase. This lasted 20 days during which the water table was restored to pre-experimental conditions in all pots (−4 cm) by watering the pots with approximately 50 ml every 2 days. We considered this duration to be sufficient to measure recovery as photosynthetic activity in Sphagnum mosses resumes within an hour after rewetting (Jassey and Signarbieux, 2019), and net CO 2 assimilation of intact peat cores was previously shown to recover in about 2 weeks , while photosynthetic potential of vascular plants is not or much less affected by reduced rainfall (Rastogi et al., 2019).
Peat C dynamics were measured as litter decomposition rates (during phase I), CO 2 and CH 4 gas fluxes on a weekly basis, during phase I, and biweekly basis during phase II, and as dissolved organic carbon (DOC) concentrations and quality in x 14 x x x Peat mesocosms were exposed to increasingly extreme precipitation regimes (see Figure 1) during the 11-week resistance phase (I) and allowed to recover during the 3-week resilience phase (II) in which water table was maintained at-4 cm below the Sphagnum capitulum. Water table depth (WTD) was measured every 2 days during phase I and twice a week during phase II, CO 2 and CH 4 fluxes were measured on a weekly basis, and vascular plant abundance (VP), DOC, and decomposition rate constant (k) were. measured as indicated.
Frontiers in Environmental Science | www.frontiersin.org June 2021 | Volume 9 | Article 659953 the pore water at five time points to cover the start, middle and end of phase I, and start and end of phase II ( Table 1). Mid-way the resistance phase (week 7), we counted the number of graminoid and ericoid shoots and their average length to verify whether the abundance of vascular plants was similar across the six precipitation regimes. The 30 mesocosms (5 reps, 6 regimes) were first placed in the greenhouse in a randomized design. Unfortunately, we were forced to move the experiment to a new location in week 5, because the greenhouse facilities closed due to the outbreak of the COVID-19 pandemic. The mesocosms were moved to a sheltered location with natural, albeit reduced light at ambient outdoor temperature. The precipitation regimes were maintained and continuous monitoring of the experimental conditions and inclusion of location in our statistical analysis (both see below) safeguarded the continuation and interpretation of our measurements.

Experimental Conditions
Throughout the experiment, peat temperature (at−5 cm), volumetric moisture content and light intensity were registered every 30 min. Soil temperature and moisture loggers (5TM sensor, Decagon devices Inc., United States) were placed in three pots from precipitation regimes 1, 3, and 6, randomly chosen from among the 5 replicates per regime. We monitored the light intensity by placing a photosynthetic active radiation (PAR) light sensor (QSO-S, Decagon Devices Inc., United States of America) in the middle between the mesocosms. Water table depth (WTD), as the distance (cm) from the top of the Sphagnum layer, was measured once every 2 days during the resistance phase, and twice a week during the resilience phase.

Decomposition Rates
Decomposition rates were measured during the resistance phase by incubating graminoid litter in litter bags (Cornelissen, 1996). The graminoid litter was collected at the field site in January 2020, dried at 60°C for 72 h, cut to fragments (1 cm) and put in nylon bags (5 × 6 cm, mesh size <1 mm), approx. 0.75 g per litterbag. The litterbags were inserted vertically into the peat at−5 cm at the start of the resistance phase and retrieved after 75 days, at the end of the resistance phase. The collected litterbags were cleaned with tap water and dried (60°C) before weighting. Decomposition rates were measured as constant k, following the exponential decay model by Olson (1963): where decomposition rate constant k is calculated based on the remaining mass (X t ) after t days of incubation as a fraction of the initial mass (X 0 ).

Dissolved Organic Carbon
We monitored the DOC concentration and quality by sampling 10 ml porewater from each mesocosm, and replacing the sampled amount by equal amounts of artificial rainwater to compensate for the water removal. Porewater was filtered at 1 µm to remove particles and micro-organisms (Whatman GmbH, Dassel, Germany) and stored at −20°C before determining DOC FIGURE 1 | Characteristics of rainfall patterns to which the mesocosms (circles) were exposed during the resistance phase, along with those observed at the reference site Counozouls (triangles) during spring and summer 2019. Experimental rainfall treatments fell into one of six precipitation regimes of reducing rainfall frequency but increasing intensity. Size of the bubbles is scaled to drought intensity, given as the average number of days between rainfall events.

Gas Fluxes (CO 2 , CH 4 )
Gas flux rates were measured using the chamber method (Riutta et al., 2007). Mesocosms were placed in airtight flux chambers (17.4 cm × 25.7 cm x 17.4 cm) to simultaneously measure CO 2 and CH 4 fluxes using a trace gas analyser (LI-COR 7810, Li-Cor Environmental-GmbH, Bad Homburg, Germany) set to record the CO 2 and CH 4 concentration every second during 5 min. Net ecosystem exchange of CO 2 (NEE) was measured using a transparent chamber while ecosystem CO 2 respiration (R eco ) was measured using a darkened chamber. Gross primary productivity (GPP), CO 2 uptake as a result of photosynthesis, was calculated as the difference between NEE and R eco . All gas fluxes were measured at ambient temperature and light conditions in the greenhouse (week 1-4) or outdoors (week 5-14). After visually confirming that for all measurements CO 2 and CH 4 concentrations in the chamber headspace changed linearly when plotted against time, gas flux rates were calculated with linear regression models using the R package gasfluxes (Fuss, 2020). Gas fluxes were corrected for the container volume, mesocosm surface area, light intensity and soil temperature following the package guidelines. Positive values for the CO 2 and CH 4 fluxes indicate C release from and negative values C uptake by the mesocosm.

Cumulative CO 2 and CH 4 fluxes
To assess the resistance and resilience responses of the C balance to extreme precipitation, we reconstructed the C fluxes on an hourly basis for each mesocosm. Often process-based models based on the physical processes that regulate energy, carbon and water fluxes are used to predict C fluxes (f.e. Lai, 2009). However, application of these models has some limitations due to the inherent complexity of the model structure and related parameters (Tramontana et al., 2015). Empirical models, and especially machine learning algorithms such as random forests, are increasingly employed to develop quantitative predictive models of C fluxes while giving more accurate predictions (Cutler et al., 2007). Therefore, we used random forests models trained with all the environmental variables and flux data recorded. Random forests (RF) algorithm for regression trees creates successive and independent regression trees using bootstrap samples of the data (i.e., a random forest) to classify data and make predictions (Liaw and Wiener, 2002;Breiman, 2001;Cutler et al., 2007). We trained a RF model with different combinations of environmental variables (soil temperature, PAR, water table depth, and volumetric water content) for each flux data type (i.e., NEE, R eco and CH 4 ; Supplementary Table S2). To get the best performing models, we iteratively tuned RF parameters in R using the randomForest package (Liaw and Wiener, 2002): number of trees to grow ("ntree") and the number of predictor variables randomly sampled as candidate at each split ("mtry") to determine the best predictor to split the response variable at each node. We further performed a 10-fold bootstrap resampling to test the performance and overfitting of the final model. For each fold, we created a bootstrap data set by randomly collecting observations from the original dataset with replacement. In other words, we repeatedly sampled data with replacement from the original training set in order to produce separate training sets. These sets were then used to allow the RF algorithm to reduce the variance of its prediction, thus greatly improving its predictive performance. We repeated this procedure 10 times and evaluated model performance using the root mean squared error (RMSE) and R 2 values calculated from the results of the bootstrapping. RMSE measures the ability of the model to predict new data, and the result are easily interpretable as they are on the same scale as the original data. The R 2 value describes the fit of the model to the data. The final RMSE and R 2 values for the ensembled C fluxes are given in Supplementary Figure S3.
The final RF models tended to underestimate NEE, R eco , and CH 4 fluxes at high flux rates (Supplementary Figure S4). Nevertheless, all models performed well in cross-validation and gave high predictive strengths (R 2 > 0.9 for all fluxes). The final RF models were used to predict NEE, R eco , and CH 4 fluxes per hour, during the resistance and resilience phases. Hourly flux rates for GPP were calculated as the difference between NEE and R eco , for day time hours, see Supplementary Figure S10 for daily averages. Finally, all hourly flux rates were summed to obtain cumulative fluxes for the resistance and resilience phases.

Calculation Resilience Response
We quantified the resilience responses to rewetting during the resilience phase for DOC concentration, S 250-465 , CO 2 and CH 4 flux rates. It is recommended to express resilience responses relatively to a baseline or a control (Ingrisch and Bahn 2018). Since regime 2 most closely resembled the rainfall regime we observed at our field site (Figure 1), we took the regime 2 as a reference treatment to control for seasonal change. We calculated the resilience response per treatment (T) as the difference in observed values for measurements in week 11 (W11, end resistance phase) and 14 (W14, end resilience phase) proportional to the median of the change in the reference treatments (R): Resulting response values that are > 1 express a larger increase (a greater capacity to "bounce back") in the treatment compared to the reference regime, and vice versa.

Statistical Analyses
All numerical analyses were performed using R 3.6.1 (R Core Team 2019).
Frontiers in Environmental Science | www.frontiersin.org June 2021 | Volume 9 | Article 659953 5 Differences in rainfall characteristics and vascular plant abundance between precipitation regimes were checked with separate generalized least squares (gls) models with regime as a fixed factor, using the nlme package (Pinheiro et al., 2016). Effects of the precipitation regimes on the WTD dynamics were tested with in separate linear-mixed effects (lme) models for the resistance and resilience phase, also using the nlme package. We used the square-root transformed WTD observations as a response variable, precipitation regime and sampling date and their interaction term as fixed factors, with mesocosm identity as a random factor to account for repeated measurements. Heterogeneity of variance between regimes was controlled by applying a varIdent variance structure (Zuur et al., 2010). We used the standard deviation (SD) of the WTD as an index for the WTD fluctuation and tested for differences between precipitation regimes with a gls model. Per mesocosm we calculated the mean and SD WTD during the resistance phase and used both as explanatory variables in subsequent analyses.
Prior to testing the responses, we verified whether the DOC concentration and quality (S 254-465 ) and C gas fluxes were equal upon the start of the experiment, by running a linear model on the measurements of week 0 with precipitation regime as an explanatory variable (see Supplementary Tables S3, S4, Supplementary Figure S9). Since R eco values differed between regimes (F 5,24 2.80, p 0.0398) and GPP values differed marginally (F 5,24 2.60, p 0.052), we normalized the data for DOC concentration, S 254-465 , NEE, R eco , GPP, and CH 4 . For each pot, the difference between the observation (X 0 ) and the overall mean (µ 0 ) at week 0 was subtracted from all observations (X n ).
The normalized data, excluding data from week 0, were used for subsequent analysis described below.
Decomposition rates were analyzed for effects of precipitation regimes with a gls model as above. We also used separate gls models to test the relationship between decomposition rates and the mean and SD WTD. For all normalized DOC concentration and S 254-465 observations we ran lme models with precipitation regime and sampling date and their interaction term as fixed factors, with mesocosm identity as a random factor. For both DOC concentration and S 254-645 the models included a varIdent on regime to account for heteroscedasticity.
The normalized gas flux rates were likewise tested for effects of precipitation regime and sampling date with lme models. Data for the resistance and resilience phase were tested separately. First, we tested whether experiment location was a significant random factor in addition to mesocosm identity for the resistance phase observations, which was the case for R eco , and GPP (Supplementary Table S5). We used the CH 4 flux rate observations from both the transparent and darkened chamber measurements and added date as an additional random factor to correct for repeated measurements. CH 4 flux rates were ln transformed before analysis. A varIdent variance structure on sampling date was applied to all gas flux rate models on the resistance phase data and for the R eco model on the resilience phase data.
Cumulative fluxes (NEE, R eco , GPP, and CH 4 ) were analyzed for effects of precipitation regime, mean, and SD WTD in separate gls models, for both the resistance and the resilience phase. Likewise, the resilience responses for all fluxes, DOC concentration and S 254-465 were tested for effects of precipitation regime, mean and SD WTD with gls models. Models testing the resilience responses of CH 4 and S 254-465 were fitted including varIdent on precipitation regimes.
Parameters of all gls and lme models were estimated with restricted maximum likelihood (REML). Normality and homogeneity of all model residuals were verified through Kolmorgorov-Smirnov and Levene's tests as well as diagnostic plots. Differences between precipitation regimes were inspected using Tukey post-hoc test. Variance explained by fixed factors (R m 2 ) and full mixed effects models (R c 2 ) were calculated according to Johnson (2014) using the MuMIn package (Barton, 2020). All graphs were made using the ggplot2 (Wickham, 2016), RColorBrewer (Neuwirth, 2014) and cowplot (Wilke, 2019) packages.

Experimental Conditions
Throughout the resistance phase the mesocosms received approximately the same amount of rainfall (145-159 mm), except for the most extreme regime, regime 6 (125 mm). Although this was part of the design to avoid overflowing of the pots, the difference was significant (Supplementary Figure  S2; F 5,30 4.24, p 0.007). As planned the regimes differed significantly in rain frequency (F 5,24 98.49, p < 0.0001), rain intensity (F 5,24 56.01, p < 0.0001) and drought intensity (F 5,24 62.36, p < 0.0001) (Figure 1). The vascular plant abundance that had developed during the resistance phase did not differ between regimes (Supplementary Table S6).
The experimental rainfall patterns during the resistance phase lead to a temporal differences in WTD (Time: F 1,1370 75.49, p < 0.0001). While the WTD differed between regimes over time (Regime*Time: F 5,1370 8.80, p < 0.0001), the overall differences between regimes were not significant (Regime: F 5,1370 0.92, p 0.48). We further inspected the WTD dynamics and found that the standard deviation (SD) of WTD per mesocosm (on average 2.2 cm), differed significantly between precipitation regimes (F 5,24 4.975, p 0.0029), indicating that the mesocosms underwent WTD fluctuations with different intensities. Therein, regimes 3 and 6 had highest SD WTD values ( Figure 2C).
During the resilience phase WTD was kept more or less constant, although a temporal trend was still significant (Time: F 1,114 8.65, p 0.004) and differed between regimes (Regime*Time: F 5,114 6.37, p < 0.0001). No overall differences between precipitation regimes were found (Regime: F 5,114 1.00, p 0.44, Supplementary Figure S5A). Notably, the WTD fluctuation during the resilience phase was small compared to the resistance phase (SD on average 0.94 cm), and the SD WTD did not differ between regimes (Supplementary Figure S5B).

Decomposition Rate
Graminoid litter lost between 2-22% of its mass during the 11 weeks of incubation in the mesocosms. This equated to values for the decomposition rate constant k varying between 0.00026-0.00328 days −1 (Figure 3, Supplementary Figure S6).  Although the decomposition rates did not significantly differ between precipitation regimes (F 5,24 1.14; p 0.366), it significantly decreased along with the mean WTD (F 1,28 7.01, p 0.0132, R adj 2 0.17), and increased with the SD WTD (F 1,28 12.05, p 0.0017, R adj 2 0.28), respectively. In other words, decomposition was faster in mesocosms with a deeper mean WTD and a high SD WTD (Figure 3).

Dissolved Organic Carbon
The concentration of DOC in the porewater of the mesocosms increased over time (F 1,53 38.34 p < 0.0001). The normalized DOC concentration initially increased, then leveled off during the resistance phase, and increased again during the resilience phase ( Figure 4A), with normalized values ranging from 32.78 ± 0 mg C L −1 at the start of the experiment, to a mean of 49.94 ± 2.78 mg C L −1 at the final sampling date. There were no overall differences in DOC concentration between precipitation regimes (F 5,53 1.02, p 0.426), nor did the temporal increase vary between precipitation regimes (Time * Regime: F 5,53 1.43, p 0.228). Values for S 254-465 became more negative as time progressed (F 1,54 13.03, p 0.0007, Figure 4B), ranging from −0.012 ± 0 to −0.014 ± 0.0004 nm −1 from the start to end of the experiment. We did not find significant differences between precipitation regimes (F 5,54 1.65, p 0.186) nor was the interaction term significant (F 5,54 Figure S8E) and related negatively to the SD WTD (F 1,23 10.16, p 0.0041, Supplementary Figure S8F).

CO 2 and CH 4 flux Rates
During the resistance phase, gas flux rates varied over time with (near) significant influence of precipitation regimes ( Table 2). Normalized NEE rates became more negative as time progressed, except for precipitation regimes 1 and 4 as indicated by the significant interaction term. R eco values displayed a significant temporal effect and overall differences between regimes ( Figure 5B, Supplementary Figure S11C). Specifically, R eco gradually increased over time, with higher respiration rates as precipitation regimes got more intense and displayed a significant increase of 55.5% between regimes 1 (mean R eco 198.94 mg CO 2 m −2 h −1 ) and 5 (mean R eco 309.45 mg CO 2 m −2 h −1 ). Because R eco peaked mid-march (week 4, Figure 5B), we ran the same analysis without the week 4 measurements. After removing the peak, the results were similar with an overall effect of time (F 1,201 25.33, p < 0.0001) and differences between precipitation regimes (F 5,201 2.42, p 0.0475) but no interaction between regime and time (F 5,201 0.87, p 0.504). GPP rates decreased with time with marginal significant differences between precipitation regimes, but the interaction term (Time * Regimes) was not significant. CO 2 uptake rates tended to be stronger the more extreme the precipitation regime became, with largest differences between regime 1 and 5 (Supplementary Figure S11E). CH 4 flux rates showed a significant interaction between time and precipitation regimes, along with a significant overall effect of time. CH 4 fluxes increased gradually during the resistance phase in most precipitation regimes except in regime 6, where the flux rates decreased by 49.7% from 144.40 to 72.61 µg CH 4 m −2 h −1 (Table 2; Figure 5D).
Gas flux rates measured during the resilience phase all increased over time (Table 2; Figure 5). Yet, none of the fluxes showed significant differences between regimes, neither as main effect nor in interaction with time. The resilience response of the CO 2 fluxes did not differ between precipitation regimes. The resilience response of CH 4 , however, tended to differ between regimes (F 4,43 2.24, p 0.080) with a smaller increase in CH 4 flux rates for regime 5 compared to regime 2 (e.g. resilience < 1, Supplementary Figure  S12D). Moreover, the resilience responses for NEE, R eco , and CH 4 did TABLE 2 | Outcome of linear mixed effects models testing the influence of precipitation regime (Regime, R) and sampling date (Time, T) on rates of net ecosystem exchange (NEE), ecosystem respiration (R eco ), gross primary productivity GPP), and methane (CH 4 ), during the resistance phase (I, week 1-10), and resilience phase (II, week 10-14). not relate to the earlier WTD dynamics during the resistance phase, whereas the resilience response of GPP tended to relate negatively to the mean WTD of the mesocosms during the resistance phase (F 1,23 3.25, p 0.085, Supplementary Figure S12).

Cumulative CO 2 and CH 4 fluxes
The total uptake and release of CO 2 and CH 4 during the resistance phase related to the experimental WTD dynamics but did not differ significantly between precipitation regimes ( Figure 6; Table 3). Most mesocosms released more CO 2 than they took up, as shown by the positive NEE values ( Figure 6A). Particularly, mesocosm which underwent little change in WTD (measured as SD WTD) emitted most CO 2 . The underlying fluxes, namely R eco and GPP, related to mean WTD (Table 3). Mesocosms with a deep mean WTD respired more CO 2 than mesocosms with a shallow mean WTD, yet CO 2 uptake was also Frontiers in Environmental Science | www.frontiersin.org June 2021 | Volume 9 | Article 659953 9 higher at deep mean WTD than at shallow mean WTD ( Figures  6C,E). Likewise, CH 4 emissions were positively related to mean WTD, with highest emissions at shallow mean WTD ( Figure 6G).
For the resilience phase, comparable relationships between CO 2 exchange and WTD dynamics were found. In other words, the WTD dynamics of mesocosms under the experimental precipitation regimes were predictive of the CO 2 dynamics when WTD dynamics were restored to pre-experimental conditions. Thus, mesocosms which underwent little variation in WTD (SD WTD) had a more positive overall CO 2 balance during resilience phase compared to mesocosms for which a high SD WTD was measured in the resistance phase. However, the relationship between mean WTD and CH 4 emission was not significant for the resilience phase; all mesocosm emitted 100-150 mg CH 4 during this three-week period.
FIGURE 6 | Relationships between cumulative gas fluxes during the resistance (I, panels A, C, E, G) and resilience (II, panels B, D, F, H) phases, and the WTD dynamics of the resistance phase. Panels (A), (B): net ecosystem exchange (NEE) is explained by the standard deviation water table depth (SD WTD). In (C), (D) ecosystem respiration (R eco ) is plotted against WTD mean, as are gross primary production (GPP, in E and F) and methane (CH 4 , in G and H). Drawn relationships are significant at p < 0.05, based on linear regressions given in Table 3. For details on precipitation regimes, please see Figure 1.
Frontiers in Environmental Science | www.frontiersin.org June 2021 | Volume 9 | Article 659953 10 DISCUSSION Precipitation regimes are becoming more extreme, as is currently observed (Knapp et al., 2008;Lehmann et al., 2015) and predicted by climate models (Dai, 2013;Fischer and Knutti, 2016). While the consequences of deeper water levels on peatland C dynamics are relatively well known (Dise, 2009;Page and Baird, 2016), the relationship between peatland water table dynamics and carbon cycling under shifting rainfall patterns is not. We aimed to explore how changes in water table depth (WTD) dynamics (overall mean and fluctuations) in response to more extreme precipitation regimes affects peatland C dynamics. By exposing intact peat cores to approximately the same amount of rain but according to increasingly extreme patterns, we found that C uptake (GPP), decomposition rate of graminoid litter and C loss (DOC, R eco , and CH 4 ) all relate to the altered WTD dynamics induced by the altered rainfall patterns. Most importantly, we showed that extreme rainfall early in the growing season had lasting effects on CO 2 exchange in our mesocosms, while the effects on CH 4 emission diminished quickly after restoration of the WTD dynamics. The limited depth of the mesocosms meant that hydrological conditions during the experiment could have been more extreme than under field conditions and the results represent amplified responses. Under natural conditions hydraulic lift and biological activity from deeper layers would contribute to the overall response. Nevertheless, WTD changes in the top 10 cm are considered relevant for peatland C dynamics (Evans et al., 2021). Our findings highlight the importance of WTD fluctuations under shifting rainfall patterns in addition to mean WTD in driving peatland C dynamics, and that short-term changes in conditions may have lasting effects.

Carbon Uptake Under Extreme Precipitation
Besides a general increase in C uptake as the growing season progressed (Figure 5), we found a near significant difference between precipitation regimes with higher uptake rates under the more extreme precipitation regime 5 than 1. In addition, cumulative C uptake (GPP) was found to be increased with deeper mean WTD (Figure 6), and tended to be higher with larger WTD fluctuations. Our findings contrast with those of Nijp et al. (2014) who showed that photosynthesis efficiency and CO 2 uptake in Sphagnum is most efficient at high moisture content and can decline steeply as the volumetric moisture content decreases (Nijp et al., 2014). While high volumetric moisture content benefits Sphagnum photosynthesis, vascular plants benefit from increased aeration of the peat profile (Breeuwer et al., 2009;Jassey et al., 2018). A decrease of rainfall frequency in a field experiment resulted in reduced C uptake by Sphagnum dominated peat, while C uptake in vascular plant dominated peat increased (Radu and Duval, 2018b). Therefore, the observed negative relationship between GPP and mean WTD is most likely caused by increased performance of vascular plants. This could be experimentally verified by monitoring the plant photosynthetic activity and biomass more closely in future studies.
Upon rewetting, C uptake rates slowed down as GPP became less negative during the resilience phase. This resilience response was similar for all precipitation regimes. As a consequence, the cumulative C uptake during the resilience phase could still be predicted by the mean WTD of the peat cores during the resistance phase. These observations fit with our supposition that vascular plants grew better under extreme precipitation. Indeed, altered vascular plant growth during spring affected carbon uptake later in the growing season under field conditions (Lund et al., 2012;Helfter et al., 2015).
Moreover, the reduction in C uptake rates upon rewetting could have been dampened by the capacity of Sphagnum to recover quickly from drying. Jassey and Signarbieux (2019) found that Sphagnum photosynthetic activity can recover from droughts within an hour after rewetting. Although this shows that peatland functioning can bounce back from extreme meteorological conditions, repeated restarting of the photosynthetic pathway after droughts when rainfall frequency is low could come at the cost of increased autotrophic respiration (McNeil and Waddington, 2003;Nijp et al., 2014Nijp et al., , 2015. We further note that under natural circumstances rain coincides with cloud cover and temporally limit photosynthesis. This could act as a secondary mechanism by which reduced rainfall frequency translates to reduced cloud cover and increased evapotranspiration thereby influencing peatland carbon dynamics (Nijp et al., 2015;Dai et al., 2018).

Litter Decomposition
We found that graminoid litter mass loss varied between 2-22% during the resistance phase, and most importantly, that both the mean and fluctuation of the WTD explained this variation (Figure 3). These results suggest that intense rain interspersed with long dry intervals promote decomposition in peatlands. We inserted our litterbags just below the overall WTD (−4.49 cm). This finding agrees with previous studies which showed that decomposition in peatlands is optimal in 3 | Outcome of separate linear regression models testing the influence of precipitation regime (Regime), WTD mean or SD (measured during phase I) on cumulative net ecosystem exchange (NEE), ecosystem respiration (R eco ), gross primary productivity (GPP), and methane (CH 4 ) during the resistance phase (I, week 0-10), and resilience phase (II,. See Figure 6 for supporting visualization. Bold and italic text indicate significant (p < 0.05) and marginally significant (p < 0.1) p values, respectively.
Frontiers in Environmental Science | www.frontiersin.org June 2021 | Volume 9 | Article 659953 the zone of water table fluctuation (Belyea, 1996;Belyea and Clymo, 2001). The graminoid litter we utilized represents recently senesced plant material entering the pool of organic matter. In general, the young organic C deposited in the acrotelm (surface layer of peat) has a higher potential to be lost by decomposition than older organic C in the deeper catotelm (Belyea 1996;Ise et al., 2008). Taken together, WTD fluctuations at the level of the acrotelm are likely to have more profound consequences than further draw down of the water table into the catotelm (Ise et al., 2008). Our results signify the sensitivity of C stored in the acrotelm to potential loss as a consequence of shift in rainfall patterns (Knapp et al., 2008;Leifeld and Menichetti, 2018). Altered WTD dynamics could also have influenced decomposition indirectly. Increasing decomposition rates was concomitant with increasing photosynthetic C uptake (GPP). As photosynthesis and root exudation are strongly linked (Kuzyakov and Cheng, 2001), plant exudates may have influenced decomposition of organic matter by catalyzing degradation of organic molecules, triggering nutrient competition or even preferential substrate utilization by decomposers (Kuzyakov, 2010;Dijkstra et al., 2013;Saar et al., 2016;Barel et al., 2019). Indeed, vascular plants have been shown to promote C loss through increase heterotrophic respiration and enhance decomposition of peat carbon due to rhizosphere priming (Gavazov et al., 2018). Further studies are however necessary to establish whether these indirect drivers play a role in decomposition on relatively short time-frames of 3 months in peatlands.

Dissolved Organic Carbon
C lost from peatlands as DOC can reach 1-50 g m −2 y −1 and represent up to 64% of organic C in streams (Blodau, 2002;Rosset et al., 2019). DOC concentrations observed in our mesocosms (19.47-88.23 mg C L −1 ) approximated earlier observations at the same field site in June 2018 (22.91-65.20 mg C L −1 ) and fall within in the range (20-60 mg C L −1 ) of concentrations reported for northern peatlands (Blodau, 2002). We observed a seasonal increase in DOC concentration. DOC concentrations are often found to increase with temperature as a consequence of increased microbial and vegetation activity (Blodau, 2002;Pastor et al., 2003;Rosset et al., 2019). Moreover, C loss by stream export varies seasonally which has been attributed to seasonal variation in hydrological and biological processes in addition to streamflow and run-off after heavy rain (Leach et al., 2016). However, we found no indication that shifts in precipitation patterns affected the DOC concentration. Yet, the limited size of our mesocosms could have excluded mechanisms from occurring that play a role in natural DOC flows, such as the release of DOC from deeper peat layers.
Additionally, we observed a seasonal decrease in S 254-645, indicating a decrease in low molecular weight compounds and/or increase in aromaticity (Hansen et al., 2016). Which is in line with long-term field observations that demonstrated an increase in DOC aromaticity during spring (Leach et al., 2016). The resilience response of S 254-645 values continued to decrease in the reference precipitation regime (regime 2), whereas mesocosms exposed to more extreme rainfall patterns did not display a decrease in S 254-465 values or to a lesser extent. As high molecular weight and aromatic compounds are presumably hard to break-down by microbial decomposers (Kalbitz et al., 2003;Cotrufo et al., 2013), these results could reflect the preferential break-down of low molecular weight/non-aromatic compounds under a normal compared to extreme precipitation regimes, regardless of differences in DOC concentration.

Ecosystem Respiration and Net Ecosystem Exchange
Besides a seasonal increase, CO 2 respiration rates differed significantly between precipitation regimes, with lowest R eco values for the least extreme regime (regime 1) and highest respiration under extreme precipitation (regime 5). Based on this observation we would have expected a relationship between cumulative R eco and WTD fluctuation, but there was not. Field observations demonstrated that R eco flux rates are less responsive to WTD dynamics than GPP and NEE rates, because R eco is the combined result of autotrophic and heterotrophic respiration (Olefeldt et al., 2017). However, we did find a significant relationship of increased cumulative R eco at deeper mean WTD during the resistance phase. This relationship follows the earlier relationships of increased heterotrophic (decomposition) and autotrophic (GPP) activity with deeper WTD.
Responses of overall CO 2 exchange (NEE) to the precipitation treatments reflected the observations in GPP and R eco . The seasonal increase in the CO 2 sink function (negative NEE rates) was stronger under extreme precipitation than the least extreme regimes. This significant interaction between time and precipitation regime echoes the differences between precipitation treatments observed in R eco and GPP flux rates. When inspecting the cumulative NEE responses, mesocosms with increasing WTD fluctuations became a weaker source of CO 2, even becoming a sink when WTD had a standard deviation of more than 3 cm, which represents the tendency for cumulative GPP to be controlled by SD WTD in a similar manner. Increased CO 2 uptake, as we found under extreme precipitation, could signal a long-term change in peatland ecology. Notably, repeated exposure of peatlands to extreme rainfall and drought is likely to trigger vascular plant encroachment (Buttler et al., 2015;Dieleman et al., 2015;Jassey et al., 2018;Lamentowicz et al., 2019) with cascading effects on overall ecosystem functioning (Bardgett et al., 2013;Dieleman et al., 2015;Jassey et al., 2018). The fact that we found no significant relationship between cumulative NEE and mean WTD stresses the importance to include a measure of WTD fluctuation as a potential controlling factor in peatland studies researching hydrological effects on peatland carbon dynamics, in addition to the average WTD.
During the resilience phase, CO 2 uptake rates lessened while respiration rates rose slightly and resulted in a rise in NEE rates. Moreover, the relationship between NEE and the WTD fluctuation observed during the resistance phase was still measurable during the resilience phase, as were the relationships between cumulative R eco and GPP values and the mean WTD. These lagged effects indicate that exposure to extreme precipitation early in the growing season may have prolonged effects on ecosystem functioning, even if meteorological conditions return to normal (Frank et al., 2015).

Methane Emissions
The evolution of CH 4 flux over time differed between precipitation regimes ( Table 2). In particular, CH 4 emissions were low in mesocosms exposed to the most extreme precipitation regime (regime 6). The limited CH 4 emission of these mesocosms in the second half of the resistance phase (April) coincided with an extended period without rainfall resulting in deeper WTD (Figures 2, 5). This link between WTD dynamics and CH 4 emission is further evidenced by the positive relationship between cumulative CH 4 flux and the mean WTD during the resistance phase ( Figure 6).
However, these results should be interpreted with some caution. Our overall measured average CH 4 flux rate of 162.2 μg m −2 hr −1 during the resistance phase, and 482.1 μg m −2 hr −1 during the resilience phase fall in the low end of the reported range (5-80 mg m −2 d −1 ) of rates in peatlands under field conditions (Blodau, 2002), probably because of the limited depth of our mesocosms. Under field conditions, CH 4 production in deeper peat layers would contribute to overall emissions, while hydraulic lift could maintain moisture content in the peat surface layer under dry conditions. Interestingly, our findings contrast those from Radu and Duval (2018a). They found that more extreme regimes increased CH 4 emission overall and stipulated that rainfall oscillations stimulated the build-up of labile substrate as result of aerobe respiration, which upon repeated rewetting boost CH 4 production. Our S 254-465 values for DOC indicate no such buildup, while the near-significant relationship of lower cumulative CH 4 emission with increased WTD fluctuation imply a reduction of CH 4 emission at more extreme precipitation. Instead, our findings are similar to patterns observed in other studies. WTD is a wellrecognized driver of CH 4 emissions Updegraff et al., 2001;Blodau et al., 2004;Turetsky et al., 2014;Riutta et al., 2020). Shallow WTD promote CH 4 emissions because methanogenesis requires anaerobic conditions, while deep WTD increases the portion of aerated surface peat accommodating consumption of CH 4 by methanotrophs as it diffuses toward the atmosphere Sundh et al., 1995).
Drying and rewetting of peat may lead to pulses of substrate and nutrients stimulating CH 4 emissions (Turetsky et al., 2014). Recently, Huth et al. (2018) observed a strong CH 4 increase after heavy precipitation in a German fen. We observed that rewetting during the resilience phase led to a steep increase in CH 4 flux rates compared to the resistance phase which tended to be restricted in the more extreme precipitation treatments. Although these results indicate that more extreme rainfall reduces the CH 4 emissions from peatlands, the contrasting results from Radu and Duval (2018a) and the variable rewetting responses between peatland types (Turetsky et al., 2014) make further investigation of gas flux responses to shifting rainfall across peatland types under field conditions advisable.

Synthesis
In summary, our results imply that shifts in precipitation patterns early in the growing season are likely to influence peatland C dynamics. In our mesocosm study, increased water table fluctuations induced by shifting rainfall patterns reduced overall C emission from peatlands, whereby a deeper overall water table improved CO 2 uptake more than it increased ecosystem respiration, and reduced CH 4 emissions. Despite the limitations of our mesocosm study, our findings are in line with observations under field conditions (Radu and Duval, 2018b). Although it may seem promising that precipitation shifts increase carbon uptake, this conclusion should be treated with caution as this may signal an increased dominance of vascular plants that could off-set peatland carbon dynamics on the long-term. Our findings further showed that extreme precipitation early in the growing season can invoke changes in CO 2 fluxes that last beyond the period of altered precipitation. Therefore, in order to fully understand the consequences of shifting precipitation patterns for the capacity of peatlands to sequester carbon, the fluctuations in the water table have to be considered as a driving factor and both the immediate and lagged effects on peatland carbon dynamics should be taken into account.

DATA AVAILABILITY STATEMENT
The data supporting the conclusions of this article are available at 10.6084/m9.figshare.c.5468205.

AUTHOR CONTRIBUTIONS
JB and VJ conceived the experiment. SH collected the intact peat cores in the field. VM and VJ executed the experiment with support of JB, AS, and SH. JB, VM, and VJ analyzed the data. JB wrote the manuscript together with VJ to which all other authors contributed.

FUNDING
This work was funded by the French National Research Agency (MIXOPEAT project, grant no. ANR-17-CE01-0007 to VJ) and Fondation pour la Recherche sur la Biodiversité which acted in cooperation with the French Agency for Biodiversity (AFB), and its partners (FRB-www. fondationbiodiversite.fr).