Plant community composition explains spatial variation in year-round methane fluxes in a boreal rich fen

Climate change is expected to impact the methane (CH₄) budget of boreal peatlands, highlighting the need to understand the factors that influence methane cycling, including plant community structure. In northern peatlands, the majority of methane is transported through plants, and the magnitude of this process is strongly linked to plant community composition. Therefore, detailed information about the role of plants in regulating year-round methane fluxes is highly valuable. This paper explores the causes of spatial variability in plot-scale methane fluxes in a northern boreal rich fen. Methane fluxes were measured using the manual chamber technique in the context of fine-scale biomass variations in plant community composition from 36 study plots over 232 d throughout a full year. The mean methane flux rates for snow-free and snow seasons were 2.55 and 0.21 mg CH₄ m⁻² h⁻¹, respectively. We found a significant correlation between methane fluxes and a plant community type associated with the presence of sedge Carex rostrata during three studied periods: year-round, snow-free, and snow season. C. rostrata was present at 13 out of 36 plots, and these 13 plots contributed 44 %–49 % of the total measured methane fluxes. The biomass of vascular plants, sedges, and C. rostrata, as well as the ratio of vascular plant to bryophyte biomass, was also significantly correlated with methane fluxes in year-round and snow-free seasons. By identifying vegetation-driven emission hotspots, these results can enhance efforts to upscale emission predictions and improve ecosystem-scale methane modeling. Thus, our findings provide valuable insights for predicting realistic future changes in peatland methane emissions throughout the year.

1 Introduction

Northern peatlands are an intrinsic part of the global carbon cycle and currently, these peatlands store more than a third of all terrestrial carbon, act as strong sinks of carbon dioxide (CO₂), and are among the main natural terrestrial sources of methane (CH₄) (Ramage et al., 2024; Schuur et al., 2022). Indeed, wetlands produce about 25 % of the total methane emissions, which are still considered an important source of uncertainty in the global methane budget (Saunois et al., 2025). The uncertainty arises from several factors: the relative contributions of methane emissions from tropical and northern wetlands, how these regions respond to rising temperatures, and the spatiotemporal patterns of the emissions (Christensen, 2024; Yuan et al., 2024).

Climate change is predicted to affect the hydrology of peatlands by increasing the water table depth (WTD) (Evans et al., 2021; Helbig et al., 2020; Swindles et al., 2019), which is one of the most well-known regulators of methane fluxes along with temperature and vegetation (Turetsky et al., 2014). According to several studies, an increased WTD would decrease methane fluxes (e.g., Pearson et al., 2015; Riutta et al., 2020) and increase the rate of decomposition and soil CO₂ emissions (Ma et al., 2022). These ecosystem-level processes are complex, though, as the increasing level of atmospheric CO₂ is predicted to enhance plant productivity (Forkel et al., 2016) and thereby the rate of root exudation (Nielsen et al., 2017) that, via priming, leads to greater methane emissions (Turner et al., 2020; Waldo et al., 2019). The expected rise in methane production could be balanced by increased methane oxidation in the topsoil, which is a probable response to enhanced microbial activity and oxygen availability caused by rising temperatures (Zhang et al., 2021). A warming climate may, however, also increase the areal cover of wet fens in the Arctic region due to permafrost thaw, which could potentially create new sources of methane release (e.g., Christensen et al., 2023; Grimes et al., 2024).

Vegetation type and responses to environmental changes are highly relevant for methane flux dynamics, as up to 90 % of ecosystem-level methane is transported through plants in northern peatlands (Ge et al., 2023; Korrensalo et al., 2022). The aerenchymatous tissues of certain vascular plant species allow methane to move from anoxic soil through the plant, thus avoiding oxidation in oxic soil layers (Ge et al., 2023; Joabsson et al., 1999), and, indeed, plant species and their specific traits have been found to reliably predict methane flux rates (Korrensalo et al., 2022). Sphagnum mosses of wet environments can also host methanotrophic microbes and thus have the potential to oxidize methane and affect the magnitude of the total emissions (Larmola et al., 2010). Climate change is predicted to accelerate the natural vegetational succession in boreal rich fens toward Sphagnum-dominated plant communities even in stable hydrological conditions (see Kolari et al., 2021, and references therein), which could have major impacts on methane production and release. To improve our understanding and predictions of future methane emissions, it is important to have a stronger focus on the vegetation composition and the specific plant species controlling the magnitude of the fluxes (Riutta et al., 2020).

The relationship between plant community composition and methane fluxes remains an important topic of study (e.g., Lai et al., 2014; Riutta et al., 2007; Ström et al., 2015), and research on individual plant species has shown significant variation in the magnitude of flux rates and transport efficiencies (e.g., Bhullar et al., 2013; Koelbener et al., 2010; Korrensalo et al., 2022). However, year-round methane flux data from northern boreal rich fens at the plot scale are limited. This study aims to better understand the causes behind local spatial variability in methane fluxes year-round and to provide a new perspective on assessing plant-mediated methane emissions. To this end, we focus on fine-scale variations in plant community composition based on species' biomass, using non-destructive in situ methods. We intend to answer the following questions: (1) Do methane flux variations correlate with plant community type at a study plot scale? (2) Does plant community composition correlate with methane flux variability alone or in combination with other environmental factors? We hypothesize that (1) the plant community composition affects the methane flux and that (2) the flux is highest on study plots with the largest biomass of vascular plants in absolute terms or in proportion to the biomass of bryophytes.

2 Materials and methods

2.1 Study site

This study was implemented in Puukkosuo, an open and slightly sloping calcareous fen located in the northern boreal zone at Oulanka National Park in Kuusamo, Northeast Finland (66.377299° N, 29.308062° E) (Fig. 1). The long-term (1992–2022) mean annual, January, and July temperatures were 0.6, −13.0, and 15.3 °C, respectively, and the mean annual precipitation was 557.4 mm. The study period (19 October 2021–31 October 2022) was slightly warmer and drier than the long-term mean, with temperatures 0.9, −11.4, and 16.5 °C and with total annual precipitation of 528.8 mm. During the snow-free season of 2022, the mean pH (n=179) measured from peat porewater at approximately 10 cm depth was 7.0, ranging from 6.74 to 7.38. These values are relatively high and place our site at the upper end of the pH range typically observed in northern boreal rich fens, which often exhibit pH values below 7 (e.g., Hájek et al., 2021; Olefeldt et al., 2017; Laitinen et al., 2021). The deepest measured water table was 9 cm below the peat surface, while the highest was 7 cm above the surface (Fig. A1). The plot-scale variation in WTD during the snow-free season of 2022 was 4–9 cm with an average of 6 cm. The vegetation at the site is dominated by vascular plants typical of rich fens (Carex spp., Trichophorum spp., Molinia caerulea, Potentilla erecta, and Menyanthes trifoliata), as well as brown mosses (Scorpidium cossonii, Campylium stellatum, and Cinclidium stygium) and peat mosses (Sphagnum spp., mostly S. warnstorfii).

Figure 1 A map of Puukkosuo rich fen showing the location of boardwalks, the experimental area, and the fence excluding reindeer. The small inset shows the location of the study area, Oulanka, in NE Finland. Orthomosaic © Petra Korhonen 2024.

2.2 Experimental design

The study area was approximately one hectare in size and included 12 spatial blocks. Each block had three study plots, resulting in 36 study plots. The plots were established in summer 2018, and the size of a plot was 2 m × 3.5 m, which included a 0.5 m wide buffer zone. Wooden boardwalks, built to minimize stepping on the surface of the peat, led to the plots. Half of the study plots (n=18) were located inside a fence, built in spring 2019, to exclude grazing by reindeer (Fig. 1). At the time of this study, the exclusion had lasted for 2–3 years. The location of the plots and the fence followed the hydrological gradient of the fen (Fig. 1) and the differences in vegetation and methane fluxes between inside and outside the fence were likely related to the hydrological conditions, rather than the effects of exclusion (Väisänen et al., unpublished data). The study plots were also assigned to snow level manipulations that began in January 2019. Within each block, one plot was an untreated control with ambient snow level, one was a snow removal plot where the snow depth was maintained at 0.25 m throughout the snow season, and one plot was a snow addition plot, where the snow from the removal plot was placed. The snow treatments had no statistically significant effects on the methane fluxes (Fig. B5).

2.3 Methane flux measurements

For methane flux measurements, a round PVC collar (inner diameter 29.5 cm) was inserted approximately 5 cm into the ground in the rear end of each plot in September 2020. During 19 October 2021–31 October 2022, we measured methane fluxes (mg CH₄ m⁻² h⁻¹) over 232 individual days between 08:00 and 18:00 LT, conducting one measurement per plot per day. On most days, only half of the plots were measured, using randomized plot selection. We used a manual, closed chamber technique (e.g., Christensen et al., 2000) with a portable LI-COR CH₄/CO₂/H₂O Trace Gas Analyzer (LI-7810) and a transparent polycarbonate chamber (height 38 cm, diameter 29 cm) equipped with a small fan to circulate the air inside the chamber (Fig. 2a). For each measurement during the snow-free season, the chamber was placed on the PVC collar for an airtight seal (Fig. 2b). Each measurement lasted for 5 min. During winter, when snow covered the collars and it was not possible to place the chamber directly on them, we measured the fluxes on top of the snowpack, also known as the floating chamber technique, to avoid disturbance to the snow and methane diffusion (Björkman et al., 2010) and consequently to the experimental plots. In comparison to the snow-free season, these measurements were taken at a slightly different spot closer to the boardwalks. Due to these limitations, snow season flux measurements should be examined with caution. Altogether, 4121 individual measurements were used in this study.

Figure 2 A flux measurement carried out with (a) a portable LI-COR Trace Gas Analyzer and (b) a clear polycarbonate chamber placed on a round collar installed at each plot in 2020. Photos by Eeva Järvi-Laturi.

To calculate methane flux rates for each plot, we used a Python script that calculated the slope of methane concentration change during the three centermost minutes of a five-minute measurement period and computed the flux in mg CH₄ m⁻² h⁻¹ using ambient air pressure and air temperature at the time of each measurement (linear regression model, e.g., Pirk et al., 2016a). We accepted the measurements with an R² value ≥ 0.95 (n=3589) and inspected all the rest (n=691) individually, leaving out measurements showing very strong non-linearity or any other sign of failed measurement (n = 159). We examined the fluxes in three periods: (1) year-round (19 October 2021–31 October 2022), (2) snow-free (13 May–26 October 2022), and (3) snow season (19 October 2021–12 May 2022 and 27–31 October 2022). The snow-free season presented a period when we were able to measure fluxes on the collar (part of the experimental plot or the experimental site could be covered with snow). The snow covered season presented a season when we measured the fluxes on top of the snow, not on the collar. Annual accumulated flux (1 November 2021–31 October 2022) was estimated by calculating a 24 h accumulated flux for each available data point by multiplying the hourly mean flux by 24. These daily flux values were then summed to obtain the annual total. The days that were missing a measurement were given the value from a previous measurement, assuming the fluxes did not vary remarkably over the diel cycle or between the days. These assumptions were based on data from two automatic chambers located at the same site near our study plots, which showed limited diel variation in fluxes (Mastepanov et al., unpublished data).

2.4 Plant community data

We studied the plant communities by identifying species and their abundance within each collar, collecting separate biomass samples from approximately 50–150 m distance from the main research area, and using the mean biomass (BM) value of each species' samples to determine the community composition within each collar. Plant surveys and sample collection were conducted from 25 July to 12 August 2022. To survey the vascular plant (VP) species within the collars, we counted each aboveground shoot individually and separated them into flowering and non-flowering categories, as the flowering shoots may hold higher biomass than the non-flowering ones. There were altogether 31 vascular plant species identified at the study site. Then we collected separate samples of each identified species. The biomass sampling was randomized so that we first selected an area where vegetation heights resembled the heights of the vegetation within the collars. Then, we randomly threw a marker and selected the first 10 non-flowering individuals of the target vascular plant species close to the marker. We additionally collected 10 flowering individuals for those species that were found flowering within the collars. Shoots were cut at the peat surface. With uncommon species, randomization could not be put into practice (Angelica sylvestris, Carex dioica, Carex panicea, Dactylorhiza sp., Drosera sp., Eriophorum angustifolium, Festuca ovina, Pinguicula sp., Saussurea alpina, and Viola epipsila). For these, samples were collected from where they could be found. For bryophytes, we visually estimated the percentage cover of each species within each collar. There were altogether 10 bryophyte species identified at the study site. We took either three samples representing 5 % (diam. 6.6 cm) or one sample representing 1 % (diam. 2.95 cm) of the collar area with small plastic jars. The bigger sample size was collected for species that could easily be found in pure patches (Aulacomnium palustre, Campylium stellatum, Scorpidium cossonii, Sphagnum spp., and Tomentypnum nitens) and the smaller sample size for the remaining, more scattered and sparse species (Aneura pinguis, Cinclidium stygium, Fissidens adianthoides, Mesoptychia rutheana, and Paludella squarrosa). The sampling locations were selected so that the target bryophyte species could be found as “pure monoculture” as possible. We then cleaned the samples of other species and litter, removed the non-living parts, and included only the colorful or leafy parts to represent the aboveground biomass.

For BM estimation, we dried the samples (n=454) in a hot air circulation oven at 40 °C for a minimum of two days and weighed them with a four decimal precision (Denver Instrument SI-234, Table A1). The dry weights were normalized either by shoot (vascular plant, g per shoot) or by cover percent (bryophyte, g per 1 %). For each vascular plant species, we used the total number of shoots, separated by flowering and non-flowering status, within a collar, multiplied by the mean dry mass per shoot to obtain the total species BM for each collar. For each bryophyte species, we used the total percentage coverage within a collar multiplied by the mean dry mass per 1 % area to obtain the total species BM for each collar. We also calculated a VP to bryophyte BM ratio for each collar using these BM values (Table A2) to study whether this parameter could potentially be used as a predictive parameter for methane emissions. We acknowledge that studying the plant communities using non-destructive methods and a mean BM of the separate samples, instead of the actual BM of the plant communities within the collars, introduces a margin of error in the calculations.

2.5 Environmental variables

We used weekly measures of WTD from June to October from regular measuring points located approximately 1 m from the methane flux measurement points. As the microtopographic variation and WTD fluctuation at the site are minimal (see Sect. 2.1), WTD measured 1 m from the plot likely provides a reasonable estimate of WTD at the plot. At every plot, soil temperature was recorded at 10 min intervals at 5 cm depth (model p107 CS CR1000, Campbell Scientific Inc., Logan, UT, USA). In October 2023, we measured peat layer thickness at the rear edge of each study plot with a thin, metallic 300 cm long auger. To determine peat chemistry, we collected pore water during frost-free periods. At the beginning of the snow-free season in 2022, Rhizon samplers (Rhizosphere Research Products, the Netherlands) were installed (at 10 cm depth) for pore water sampling in the middle of the collar in each experimental plot. Pore water was sampled five times during the frost-free periods (31 May, 29 June, 22 July, 28 August, and 29 September). Samples were collected into evacuated opaque syringes over a period of 24–48 h, filtered (0.45 µm, sterile nylon, Sarstedt, Germany), and frozen (−18 °C). Thawed samples from all sampling campaigns were analyzed for pH (913 pH/DO Meter, Metrohm), dissolved organic carbon (DOC) and inorganic carbon (DIC; Shimadzu DOC-VCX, Trios), dissolved organic nitrogen (DON), ammonium (NH${}_{\mathrm{4}}^{+}$), nitrite + nitrate (NO${}_{\mathrm{2}}^{-}$ + NO${}_{\mathrm{3}}^{-}$; AA500 Seal Analytical), and for DOC and DON reported as mg L⁻¹ and for NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{2}}^{-}$ + NO${}_{\mathrm{3}}^{-}$ as µg L⁻¹. In analyses, we used the mean value of each measured variable. Additionally, we estimated the percentage cover of litter inside the collar while identifying the vegetation.

2.6 Data analyses

To analyze the vegetation, we divided the plant communities into three species combinations: (1) all species, (2) VPs, and (3) bryophytes. We analyzed these combinations separately using plant BM estimates (see Sect. 2.4) with hierarchical cluster analysis using Sorensen (Bray–Curtis) distance measure and Flexible Beta group linkage method (McCune et al., 2002) with a beta value of −0.25. Clusters with fewer than six samples (study plots) were discarded. To evaluate which species were statistically most connected to the different clusters, we carried out an indicator species analysis separately for all three species combinations. We tested the differences between the clusters with Tukey's Honestly Significant Difference (HSD) test following an analysis of variance (aov function, R Core Team, 2024) to assess the significance of clusters in relation to methane fluxes in different time periods (snow-free, year-round, and snow season). With clustered data, we performed detrended correspondence analysis (DCA) to identify patterns in species composition and canonical correspondence analysis (CCA) to relate the species composition to environmental variables (McCune et al., 2002). We conducted the analyses separately for VPs and bryophytes to determine the main environmental factors characterizing the composition of these plant communities. Pairwise correlations between multiple BM and environmental variables were examined with a separate CCA, using snow-free season averages. The resulting Pearson correlation matrix (Table B2) provides insights into these relationships. Finally, we conducted local regression models (LOESS, locally estimated scatterplot smoothing) to explore nonlinear trends between various variables and methane fluxes in relation to VP clusters.

Cluster, indicator species, and correspondence analyses were implemented with PC-ORD version 7.09 (McCune and Mefford, 2018). Regression models and Tukey's HSD tests were performed using RStudio version 2024.4.2.764 (Posit team, 2024). Both programs were used for data visualization: PC-ORD for cluster dendrograms, boxplots, and ordination graphs and RStudio for line graphs and scatterplots. Significance was defined as p<0.05. We used the following R packages: readxl v1.4.3 (Wickham and Bryan, 2023), dplyr v1.1.4 (Wickham et al., 2023a), tidyr v1.3.1 (Wickham et al., 2024), ggplot2 v3.5.1 (Wickham, 2024), forcats v1.0.0 (Wickham, 2023), scales v1.3.0 (Wickham et al., 2023b), paletteer v1.3.0 (Hvitfeldt, 2021), ggnewscale v0.5.0 (Campitelli, 2024), and viridis v0.6.5 (Garnier et al., 2024). R scripts were created with the assistance of Microsoft 365 Copilot, an AI-powered productivity tool.

Figure 3 Year-round methane fluxes (mg CH₄ m⁻² h⁻¹) of Puukkosuo with average soil temperature at 5 cm depth. The gray flux range shows the overall variation of the fluxes. Smoothed conditional means (LOESS method) and data points, which represent individual plot-scale flux measurements, are segregated by vascular plant clusters defined in the text in Sect. 3.2. The vertical dashed lines mark the start and end of the snow-free season (13 May–26 October 2022).

3 Results

3.1 Methane fluxes

During the snow-free season, methane fluxes were the highest, ranging from 0.02 to 9.17 mg CH₄ m⁻² h⁻¹ (i.e., 0.48–220 mg CH₄ m⁻² d⁻¹), with an overall mean of 2.55 mg CH₄ m⁻² h⁻¹ (i.e., 61 mg CH₄ m⁻² d⁻¹). During this season, the mean of plot-scale methane fluxes ranged from 0.51 to 4.67 mg CH₄ m⁻² h⁻¹ (Fig. B7), and the daily coefficient of variation spanned from 38.9 % to 85.4 %. The highest individual fluxes per plot varied between 1.51 and 9.17 mg CH₄ m⁻² h⁻¹. The lowest flux values of the snow-free period (0.02–1.18 mg CH₄ m⁻² h⁻¹) were measured in May after a spring burst (1 April–12 May 2022) and at the end of the season in late October (Fig. 3). The magnitude of the spring burst differed among the plots, with maximum individual fluxes ranging from 0.15 to 6.65 mg CH₄ m⁻² h⁻¹. During the snow season, flux measurements fluctuated between 0 and 6.65 mg CH₄ m⁻² h⁻¹ (i.e., 0–160 mg CH₄ m⁻² d⁻¹), with an overall mean of 0.21 mg CH₄ m⁻² h⁻¹ (i.e., 5 mg CH₄ m⁻² d⁻¹). The mean of plot-scale fluxes ranged from 0.07 to 0.56 mg CH₄ m⁻² h⁻¹ (Fig. B8), and the daily coefficient of variation spanned from 39.3 % to 300.4 %. Further, the snow season fluxes accounted, on average, for 8.2 % of the estimated annual accumulated flux, with values ranging from 2.3 % up to 21.3 % across the study plots (Table B3). On a year-round scale, the overall mean flux was 1.37 mg CH₄ m⁻² h⁻¹, with mean fluxes varying between 0.29 and 2.52 mg CH₄ m⁻² h⁻¹, and the daily coefficient of variation ranged from 38.9 % to 300.4 %.

3.2 Plant community structures

The plots differed in their plant community structures, with wide variation in plant BM. VP BM per plot ranged from 62.1 to 486 g m⁻², while the bryophyte BM ranged from 65.1 to 269 g m⁻². The lowest total BM for an individual plot was 167 g m⁻², while the highest total BM was 674 g m⁻². The percentage of VPs in the total BM varied between 21 % and 72 % (Table A2). The species producing the most BM, based on combined dry mass estimates from all plots, were sedges Carex lasiocarpa (122 g) and C. rostrata (84 g) and bryophytes Scorpidium cossonii (193 g), Campylium stellatum (113 g), and Sphagnum spp. (92 g) (Table A3).

The cluster and indicator species analyses for all species, VPs, and bryophytes yielded three clusters (groups of plots), each with significant (p<0.05) indicator species statistically connected to the clusters (Figs. B1–B3). The cluster analyses indicated that bryophytes (Sphagnum spp., S. cossonii, and C. stellatum) showed the strongest connection to the clusters when analyses were conducted with all plant species or bryophytes (Figs. B1 and B3). On the other hand, when analyzing VPs alone, the clusters connected mostly with different sedges (C. rostrata with the first, C. chordorrhiza with the second, and Trichophorum cespitosum, C. lasiocarpa, and Potentilla erecta (a forb) with the third cluster). From here on, the clusters are referred to as the C. rostrata cluster, the C. chordorrhiza cluster, and the T. cespitosum cluster, respectively. The characteristics of each VP cluster, including their community structure and indicator species, were studied by comparing the indicator values of the species (Table B1). The results demonstrated that the plant communities differed between the clusters, as most species were abundant in only one or two clusters, and only three species (Campylium stellatum, Vaccinium oxycoccos, and Scorpidium cossonii) were common in all clusters (Table B1). The spatial division of the vascular plant communities was clear – the C. rostrata community plots were all located in the upper part of the sloping fen, while the other two clusters were more widely distributed (Fig. B6).

When all plant species were considered, bryophytes defined the vegetation clusters (Fig. B1). Bryophyte communities and the biomass of Sphagnum mosses correlated with WTD and pH, but these did not, however, explain the distribution of the vascular plant communities (Table B2). Additionally, there was an association between Sphagnum mosses and C. rostrata in the plant communities: the C. rostrata cluster had the highest abundance of Sphagnum mosses among the three vascular plant clusters (Table B1).

Figure 4 Methane fluxes in the snow-free season (left panel), the whole year (middle panel), and the snow season (right panel), as divided by cluster analyses of (a) vascular plants (VPs), (b) all species (ALL), and (c) bryophytes (BRYO). Abbreviations for the clusters and the species with the highest fidelity according to indicator species analysis: C.Ros = C. rostrata, C.Cho = C. chordorrhiza, T.Ces = C. lasiocarpa, P. erecta, and T. cespitosum, Sco = S. cossonii, Sph = Sphagnum spp., and Cam = C. stellatum. Asterisks above the bars denote significant differences between the clusters at the level $p<{\mathrm{0.01}}^{**}$ and $p<{\mathrm{0.001}}^{***}$.

3.3 Clusters' relation to methane fluxes

An analysis of variance (ANOVA) revealed that VP clusters significantly explained variation in methane fluxes during the snow-free (F=10.71, p<0.001) and year-round seasons (F=10.92, p<0.001) but not during the snow season (F=2.14, p>0.05). Tukey's HSD test indicated that the C. rostrata cluster, which had the highest fluxes, especially during the snow-free season (Fig. 3), differed significantly from the C. chordorrhiza and T. cespitosum clusters in the snow-free (p<0.01 and p<0.001, respectively) and year-round periods (p<0.01 and p<0.001, respectively) but not in the snow season (Fig. 4a). There were no significant differences between the C. chordorrhiza and T. cespitosum clusters in any of the periods (Fig. 4a–c). The cluster with all species considered did not significantly explain variation in methane fluxes in any period (F=1.51, p>0.2 for snow-free, F=1.45, p>0.2 for year-round, and F=0.57, p>0.5 for snow season, Fig. 4b). Similarly, bryophyte clusters did not show significant effects on methane fluxes (F=1.23, p>0.3 for snow-free, F=1.26, p>0.2 for year-round, and F=0.90, p>0.4 for snow season, Fig. 4c).

Figure 5 Ordination graphs of vascular plant data. (a) DCA with rare species downweighted (eigenvalue axis 1 = 0.463, axis 2 = 0.154) and (b) CCA (eigenvalue axis 1 = 0.423, axis 2 = 0.163). Plant community clusters: C.ros = Carex rostrata, C.cho = Carex chordorrhiza and T.ces = Carex lasiocarpa, Potentilla erecta, and Trichophorum cespitosum. Time spans for methane fluxes: CH₄ year-round (19 October 2021–31 October 2022), CH₄ snow-free (13 May–26 October 2022), and CH₄ snow (19 October 2021–12 May 2022, 27–31 October 2022). NO₃+NO₂ is the combined concentration of nitrate (NO${}_{\mathrm{3}}^{-}$) and nitrite (NO${}_{\mathrm{2}}^{-}$), BMratio is the ratio of vascular plant and bryophyte biomasses, and Peat (m) is the peat layer thickness.

3.4 Ordination analyses

Within VP communities, the C. rostrata cluster diverged from the C. chordorrhiza and T. cespitosum clusters in both detrended correspondence analysis (DCA) and canonical correspondence analysis (CCA), while the C. chordorrhiza and T. cespitosum clusters diverged in DCA but not in CCA (Fig. 5a and b). In bryophyte communities, the Sphagnum cluster diverged from the S. cossonii and C. stellatum clusters in both DCA and CCA, while the S. cossonii and C. stellatum clusters diverged in DCA but not in CCA (Fig. B4a and b). The main compositional gradients of VP and bryophyte communities displayed different correlation patterns with environmental variables. The first ordination axis of the VP communities correlated with snow-free and year-round methane fluxes (r=0.775 and 0.782 in DCA, $r=-\mathrm{0.856}$ and −0.866 in CCA, respectively, Fig. 5a and b). Additionally, the first ordination axis correlated with snow season fluxes (r=0.445 in DCA, −0.402 in CCA). The correlation was also significant for the ratio of VP to bryophyte BM (r=0.562 in DCA, −0.730 in CCA), peat layer depth (r=0.342 in DCA and −0.580 in CCA), and the combined nitrate (NO${}_{\mathrm{3}}^{-}$) and nitrite (NO${}_{\mathrm{2}}^{-}$) concentration (r=0.529 in DCA, −0.453 in CCA). In contrast, the strongest compositional gradient of the bryophyte data correlated with WTD and pH in both DCA ($r=-\mathrm{0.606}$ and 0.473, respectively) and CCA ($r=-\mathrm{0.614}$ and 0.408, respectively) (Fig. B4a and b). The correlation between the bryophyte communities and methane was not significant in any of the periods (r=0.267, 0.262, and 0.132 in DCA, r=0.316, 0.326, and 0.130 in CCA in snow-free, year-round, and snow season, respectively).

3.5 BM, environmental variables, and methane

In the correlation matrix (Table B2), the total BM of VPs correlated with the total BM of sedges (r=0.98) and the total BM of C. rostrata (r=0.93), indicating that sedges were the main functional group of VPs, and C. rostrata was the main VP species producing BM at our study site. The BM variables that had a significant pairwise correlation with methane fluxes in year-round and snow-free seasons were the total BM of VPs (p<0.001), the total BM of sedges (p<0.001), the total BM of C. rostrata (p<0.001), and the ratio of VP to bryophyte BM (p<0.01) (Fig. 6a–d). Significant correlations during the snow season were not found. Environmental variables that had a significant pairwise correlation with methane fluxes in snow-free and year-round periods were pH (p<0.05 and p≤0.01, respectively, Fig. 6e) and combined concentration of NO${}_{\mathrm{3}}^{-}$ and NO${}_{\mathrm{2}}^{-}$ in peat pore water (p<0.05 for both periods, Fig. 6f). There was no significant correlation between methane fluxes and WTD or soil temperature in any period. All correlation coefficients are listed in a correlation matrix in Table B2.

Figure 6 Scatterplots showing the relationships between methane (CH₄) flux (mg m⁻² h⁻¹) and environmental variables: (a) total biomass (BM) of vascular plants (p<0.001, r=0.54), (b) total BM of sedges (p<0.001, r=0.57), (c) BM of C. rostrata (p<0.001, r=0.60), (d) BM ratio of vascular plants and bryophytes (p<0.01, r=0.45), (e) pH (p<0.05, $r=-\mathrm{0.38}$), and (f) NO₃ + NO₂ concentration (p<0.05, r=0.34) across different vascular plant clusters. Vascular plant clusters are defined in the text in Sect. 3.2. The smoothing lines were fitted using local polynomial regression (LOESS) with a span value of 1 for moderate smoothing. Biomasses are in grams of dry weight.

We found that when C. rostrata BM was 0 g, there was a large variation in the magnitude of methane flux (Fig. 6c). Therefore, we additionally examined the clusters' relation to fluxes in plots without any C. rostrata shoots (23 plots out of 36), as well as the relationship between methane fluxes and environmental variables using linear regression models (lm function, R Core Team, 2024). This revealed that, in the absence of C. rostrata, pH and litter coverage were the primary environmental drivers of snow-free season methane fluxes, showing significant relationships (p<0.001 and p=0.01, respectively) and Pearson correlation values of −0.66 and −0.29, respectively.

4 Discussion

4.1 Methane flux variation over time and space

During the snow-free season, methane fluxes ranged from 0.48 to 220 mg CH₄ m⁻² d⁻¹, with an overall mean of 61 mg CH₄ m⁻² d⁻¹. Similar fluxes have been measured from a northern boreal rich fen (Jammet et al., 2017). We also observed that the snow season fluxes accounted for 2.3 %–21.3 % of the estimated annual accumulated flux (Table B3). Even though we measured fluxes on top of an undisturbed snowpack, our results are in line with a study from fens in central Finland, where 6 %–17 % of the annual methane release was observed during wintertime, when measuring the fluxes on the peat surface after clearing the spots from snow (Alm et al., 1999). We acknowledge that direct measurements from the ground surface would have improved our ability to connect winter fluxes with the studied plant communities. However, we do not expect significant differences in flux values obtained using a closed chamber on top of the snowpack, as methane diffusion through the snowpack may be relatively unrestricted (Alm et al., 1999; Pirk et al., 2016b). Since winter may account for up to 20 % of the annual methane flux in boreal fens (this study, Alm et al., 1999), any changes in wintertime processes may impact future methane emissions from these regions. Our results highlight the importance of including winter in methane flux studies, as this information may help reduce the current large uncertainties in the net carbon balance. Further, the reduced variability of fluxes and the lack of significant variation among vegetation clusters during mid-winter suggest that winter fluxes could be effectively captured even with lower but more optimized sampling efforts, aiding in the estimation of annual emissions (Vargas and Le, 2023).

Flux rates during the snow-free season reached up to 9.17 mg CH₄ m⁻² h⁻¹ at certain study plots, while they never exceeded 2 mg CH₄ m⁻² h⁻¹ in some others (Fig. 3). The high BM of vascular plants, and especially the presence of C. rostrata, contributed to the spatial variability in fluxes: all plots with more than five flux measurements exceeding 6 mg CH₄ m⁻² h⁻¹ (n=9) contained C. rostrata shoots, but only three out of 19 plots that never recorded fluxes above this threshold had C. rostrata present. However, C. rostrata biomass cannot be the sole factor in high methane fluxes, as one plot without any C. rostrata shoots still measured fluxes as high as 9.05 mg CH₄ m⁻² h⁻¹. The possible explanations, other than the effect of aboveground BM of vascular plants or C. rostrata, are discussed further in Sect. 4.2 and 4.3. Unlike the snow-free season, the high flux values during the spring burst at the end of the snow season could not be explained by plant biomass or any other studied environmental variable. Other reasons for the spatial variation in flux rates may include species-specific plant traits such as rooting characteristics (Ge et al., 2023), ecohydrological aspects such as peat water holding capacity (Zhang et al., 2020), or microbial metabolic interactions such as nutrient cycling (Kujala et al., 2024; Yavitt et al., 2012), which contribute to soil conditions, substrate availability, and microbial activity. Interestingly, the methane fluxes from C. rostrata community plots were the highest from late July to late August (Fig. 3), when vegetation at our site remained predominantly green, even though fluxes from C. rostrata shoots have been reported to be the highest when leaves are senescing (Ge et al., 2024). The period from late July to late August may coincide with the peak development of permeable root surface area (Reid et al., 2015), whose extent, in turn, is a key factor influencing methane transport in plants (Henneberg et al., 2012). Therefore, the seasonal changes in methane flux rates associated with C. rostrata may be controlled by the belowground parts of the plant (Ge et al., 2024).

4.2 Plant communities, biomass, and methane

During the growing season, vegetation is known to correlate substantially with methane fluxes (e.g., Lai et al., 2014; Riutta et al., 2007; Ström et al., 2015), but we found significant correlations year-round, supporting our first hypothesis. Thus, our results indicate that the influence of plant communities on methane flux dynamics is not limited to the active growing season or plant senescence, and that some species are likely to be more efficient at supporting methane production and transport, even under the snowpack. Indeed, continuous production and spatial variation in methane flux rates have been observed throughout the cold season, and the variation has been associated with differences in plant community composition and the adjoining differences in substrate quantity and quality (Pirk et al., 2016b). However, the significant correlation between vegetation and methane fluxes during the snow season, identified by CCA but not by linear regression models, suggests that the relationship between vegetation and snow season methane fluxes is complex and nonlinear.

Both the total vascular plant biomass and the biomass ratio of vascular plants and bryophytes correlated significantly with methane fluxes during snow-free and year-round periods, supporting our second hypothesis (Fig. 6a and d; Table B2). There was a strong correlation between the biomasses of vascular plants, sedges, and C. rostrata (Table B2), and the biomasses of vascular plants, sedges, and C. rostrata also correlated with methane fluxes (Fig. 6a–c; Table B2). These findings suggest, similar to earlier studies (e.g., Ge et al., 2023; Korrensalo et al., 2022), that plant functional type and species largely determine the magnitude of the fluxes. High amounts of vascular plant, sedge, and C. rostrata biomass likely enhance methane production and release by supplying labile organic carbon substrates for methanogenesis through deep root systems throughout the year (Alm et al., 1999; Joabsson et al., 1999; Saarinen, 1996). High flux rates from C. rostrata dominated plots (Fig. 6c) may be due to the species' high methane transport rate (Ge et al., 2023), the high porosity of its roots (Ge et al., 2023), and a low capacity to oxidize methane into CO₂ in the rhizosphere (Ström et al., 2005). Additionally, the perennial nature and deep rooting traits of C. rostrata (Saarinen, 1996, 1998) could support methane production and transport during the cold season by providing substrates for microbial processes in deeper peat layers and a potential pathway from belowground to the atmosphere. Approximately 40 % of C. rostrata shoots at our study site overwinter green (Cunow et al., unpublished data), indicating the potential to transport gases during winter as well. Moreover, the significant correlation between the ratio of vascular plant to bryophyte biomasses and methane fluxes brings a new perspective to the discussion, as previous studies have mainly focused on finding differences between single plant species (e.g., Bhullar et al., 2013; Ge et al., 2023; Koelbener et al., 2010; Korrensalo et al., 2022) or their role in a plant community (e.g., Lai et al., 2014; Riutta et al., 2007; Ström et al., 2015). This ratio could be obtained remotely (e.g., multispectral imaging, Wolff et al., 2023) and potentially used as a predictive parameter for modeling peatland methane fluxes.

A majority (75 %) of the highest methane fluxes observed in this study originated from C. rostrata dominated plots. However, the magnitude of the fluxes was not solely dictated by C. rostrata biomass; the fluxes increased rapidly in a non-linear manner and only slightly with higher biomass, remaining relatively high and stable in most C. rostrata cluster plots (Fig. 6c). As a species with high gas transport efficiency, even a single C. rostrata shoot can transport a substantial volume of gases, potentially the same amount as a community with multiple shoots, where gas transport is distributed among many individuals (Koelbener et al., 2010; Korrensalo et al., 2022). This division of gas transport among multiple shoots could explain the observed saturation in methane flux rates: as the biomass of C. rostrata increases, the total flux remains consistent regardless of the number of shoots. In our analysis, an increase in biomass corresponds to a higher number of shoots, and therefore we cannot argue that the size of individual C. rostrata shoots did not impact the plant's methane transport efficiency. This study focused only on the aboveground parts of the plant communities, and therefore the multiple effects that the belowground parts may have on methane production, consumption, transport, and emission (e.g., Määttä and Malhotra, 2024) were not considered in the analyses. However, we acknowledge that belowground plant characteristics likely play a significant role in explaining the variation in our methane flux data that could not be explained by the other studied variables.

4.3 Environmental factors and methane

We also asked if the spatial variability of methane fluxes could be explained by the combination of plant communities and environmental variables and found that vegetation composition was the primary driver of this variation. However, in some cases, the role of an environmental variable seemed to depend on the presence of plants. Firstly, we analyzed snow-free season methane fluxes from plots without any C. rostrata shoots and found that pH and litter cover emerged as the main drivers of methane fluxes. However, as the significance of pH decreased when including the plots with C. rostrata, it is likely that the significant relationship between methane fluxes and pH in snow-free and year-round periods is explained through vegetation. Secondly, when C. rostrata was present, NO${}_{\mathrm{3}}^{-}$ and NO${}_{\mathrm{2}}^{-}$ levels correlated positively with methane fluxes during snow-free and year-round periods. In plots without C. rostrata shoots, this significant relationship disappeared, suggesting that the positive correlation is related to the presence of the species. The positive correlation was unexpected, as these compounds usually inhibit methanogenesis (e.g., Knorr and Blodau, 2009). More precisely, higher levels of NO${}_{\mathrm{3}}^{-}$ may promote methane oxidation, leading to lower fluxes (Song et al., 2022), but the impact of these electron acceptors remains debated (Zhang et al., 2021). Thus, our findings underscore the need to further investigate the interaction between C. rostrata and NO${}_{\mathrm{3}}^{-}$ and NO${}_{\mathrm{2}}^{-}$ in driving methane fluxes. Water table depth is commonly thought to regulate methane fluxes in peatlands (Lai, 2009; Turetsky et al., 2014) but did not correlate with methane fluxes at our site in any period. Indeed, in fens with a stable water table, such as our study site (standard deviation of WTD 1–3 cm), the depth of the water table may not control the variability of the fluxes (e.g., Ge et al., 2023). The thickness of the peat layer did not show a significant pairwise correlation with methane fluxes, which supports the findings that most of the methane is produced from fresh root litter and root exudates instead of old, recalcitrant peat (Ström et al., 2012). However, peat layer thickness correlated significantly with the first ordination axis in CCA and had a significant pairwise correlation with the biomass of vascular plants, sedges, and C. rostrata, as well as with the biomass ratio of vascular plants and bryophytes (Table B2), all of which were proxies for higher methane fluxes. This suggests that peat layer composition and depth can also impact methane production through an indirect relationship. Lastly, methane fluxes did not correlate with peat temperature at 5 cm depth. Indeed, methane fluxes in boreal rich fens are associated with deeper soil temperatures, which are connected to water table position, rather than with surface temperatures influenced by air temperature (Olefeldt et al., 2017). All these findings highlight that vegetation, rather than environmental factors, was the main driver of methane fluxes at our site.

5 Conclusions

Our year-round methane flux data showed significant spatial and temporal variations in flux magnitude in a northern boreal rich fen. Plant community composition, particularly the biomass of C. rostrata, explained the plot-scale spatial variation in methane fluxes well during snow-free and year-round periods. Multivariate analysis also revealed that a vascular plant cluster depicted by C. rostrata correlated with snow season methane fluxes. In addition, the total biomass of vascular plants and the ratio of vascular plant to bryophyte biomasses also showed a significant positive relationship with methane fluxes in both year-round and snow-free seasons. On the other hand, the role of other environmental factors, such as pH and nitrogen, in driving methane fluxes was subtle and seemed to depend on vegetation. Our findings suggest that the biomass ratio of vascular plants and bryophytes could be used as a parameter for predicting methane emissions at study sites with similar vegetation to ours. Importantly, these findings help to predict future changes in peatland methane emissions more accurately and realistically, which are essential for estimating the potential impacts of ongoing climate change.

Appendix A

Figure A1 The depth of the water table (WTD) at 36 study plots in Puukkosuo fen during June–October 2022. The average minimum value for WTD was −4.8 cm, and the average maximum was +1.8 cm from the peat surface.

Table A1 Distribution (Dist), mean (M), median (Mdn), and standard deviation (SD) of the heights (h) and dry biomass (BM) of vascular plant (VP) and bryophyte (Bryo) species identified from the 36 study plots (P) in July 2022. The number of VPs within the studied area was counted by the number of shoots (Shoots_P). BMs are reported only for samples (S) collected for BM estimations. Separate samples were collected for plants in different growth stages: flowering (FL) and non-flowering (NF) VP species, except for four species for which these categories were combined (Comb). The sample size (Sample) of bryophytes depended on the commonness of the species and was either 1 % or 5 % of the studied area. For all Bryo species, a BM for 1 % coverage of the studied area was calculated (BM 1 %). Abbreviations for the plant species: Andromeda polifolia (AndPol), Angelica sylvestris (AngSyl), Betula nana (BetNan), Carex chordorrhiza (CarCho), Carex dioica (CarDio), Carex flava (CarFla), Carex lasiocarpa (CarLas), Carex limosa (CarLim), Carex panicea (CarPan), Carex rostrata (CarRos), Dactylorhiza species (DactSp.), Drosera anglica (DroAng), Drosera rotundifolia (DroRot), Equisetum fluviatile (EquFlu), Equisetum variegatum (EquVar), Eriophorum angustifolium (EriAng), Eriophorum gracile (EriGra), Eriophorum latifolium (EriLat), Festuca ovina (FesOvi), Menyanthes trifoliata (MenTri), Molinia caerulea (MolCae), Pedicularis palustris (PedPal), Pinguicula species (PingSp.), Potentilla erecta (PotEre), Saussurea alpina (SauAlp), Selaginella selaginoides (SelSel), Tofieldia pusilla (TofPus), Trichophorum alpinum (TriAlp), Trichophorum cespitosum (TriCes), Vaccinium oxycoccos (VacOxy), Viola epipsila (VioEpi), Aneura pinguis (AnePin), Aulacomnium palustre (AulPal), Campylium stellatum (CamSte), Cinclidium stygium (CinSty), Fissidens adianthoides (FisAdi), Mesoptychia rutheana (MesRut), Paludella squarrosa (PalSqu), Scorpidium cossonii (ScoCos), Sphagnum warnstorfii (SphWar), and Tomentypnum nitens (TomNit). Plot-scale vegetation data can be obtained from the corresponding author.

Table A2 Total plot-scale dry mass weight biomass (BM) of vascular plants (VP) and bryophytes (Bryo) in grams inside the collar (A=660.5 cm²) at each experimental plot (37–72), a VP to Bryo BM ratio (VP : Bryo ratio), and the portion of VP BM in the total BM (VP of total BM).

Table A3 Combined dry mass weight estimation (g) of all the study plots for each species. Vascular plants are listed at the top and bryophytes at the bottom of the list in alphabetical order. For abbreviations of the species, see Table A1.

Appendix B

Figure B1 Cluster dendrogram of all species data. Indicator species by clusters: Sco = Scorpidium cossonii, Sph = Sphagnum spp. and Paludella squarrosa, and Cam = Campylium stellatum.

Figure B2 Cluster dendrogram of vascular plant data. Indicator species by clusters: C.ros = Carex rostrata, C.cho = Carex chordorrhiza, and T.ces = Carex lasiocarpa, Potentilla erecta, and Trichophorum cespitosum.

Figure B3 Cluster dendrogram of bryophyte data. Indicator species by clusters: Sco = Scorpidium cossonii, Sph = Sphagnum spp. and Paludella squarrosa, and Cam = Campylium stellatum.

Figure B4 Ordination graphs of bryophyte data. (a) DCA with rare species downweighted (eigenvalue axis 1 = 0.950, axis 2 = 0.364) and (b) CCA (eigenvalue axis 1 = 0.711, axis 2 = 0.510). Abbreviations for plant community clusters: Sco = Scorpidium cossonii, Sph = Sphagnum spp. and Paludella squarrosa, and Cam = Campylium stellatum. WTD stands for water table depth and CarRos for the biomass of C. rostrata.

Figure B5 The relation of the snow level manipulation experiment's treatments and methane flux in different time periods. There were no significant (p<0.05) differences between the treatments in relation to the fluxes. Control plots represent the natural snow level, while in removal plots, the snow depth was maintained at 0.25 m by shoveling the snow and distributing it evenly onto the addition plots.

Figure B6 Location of the study plots at Puukkosuo fen segregated by vascular plant clusters. Indicator species most strongly connected to the vascular plant clusters: C.ros = C. rostrata, C.cho = C. chordorrhiza, and T.ces = C. lasiocarpa, P. erecta, and T. cespitosum. Orthomosaic © Petra Korhonen 2024.

Figure B7 Spatial variation of snow-free season (13 May–26 October 2022) methane fluxes (mg CH₄ m⁻² h⁻¹). Orthomosaic © Petra Korhonen 2024.

Figure B8 Spatial variation of snow season (19 October 2021–12 May 2022 and 27–31 October 2022) methane fluxes (mg CH₄ m⁻² h⁻¹). Orthomosaic © Petra Korhonen 2024.

Table B1 Indicator values (% of perfect indication, based on combining the values of relative abundance and relative frequency) of all identified species (n=40) of Puukkosuo segregated by the vascular plant clusters. The values for the clusters' significant indicator species are bolded. Species with an indicator value lower than 20 % were excluded from the table. Abbreviations for the species most strongly connected to the clusters: C.ros = C. rostrata, C.cho = C. chordorrhiza, and T.ces = C. lasiocarpa, P. erecta, and T. cespitosum. Abbreviations for the species: CarRos = Carex rostrata, MenTri = Menyanthes trifoliata, TriAlp = Trichophorum alpinum, CarLim = Carex limosa, CinSty = Cinclidium stygium, CarCho = Carex chordorrhiza, VacOxy = Vaccinium oxycoccos, BetNan = Betula nana, AndPol = Andromeda polifolia, EriLat = Eriophorum latifolium, TriCes = Trichophorum cespitosum, PotEre = Potentilla erecta, CarLas = Carex lasiocarpa, EquFlu = Equisetum fluviatile, MolCae = Molinia caerulea, TofPus = Tofieldia pusilla, PingSp. = Pinguicula sp., CarDio = Carex dioica, CamSte = Campylium stellatum, ScoCos = Scorpidium cossonii, and SphWar = Sphagnum warnstorfii.

Table B2 Pearson's correlation matrix among the environmental and biomass variables. Significant correlations are bolded, p-value ≤ 0.05 (r≥0.33) marked with one asterisk (^*), p-value ≤ 0.01 (r≥0.42) with two asterisks (${}^{**}$) and p-value ≤ 0.001 (r≥0.53) with three asterisks (${}^{***}$). CH4_y = year-round methane flux, CH4_sf = snow-free season methane flux, CH4_s = snow season methane flux, WTD = water table depth, NO₃ + NO₂ = nitrate and nitrite, NH₄ = ammonium, TOTN = total nitrogen, TOC = total organic carbon, IC = inorganic carbon, BMratio = vascular plant to bryophyte biomass ratio, BMvp = vascular plant biomass, BMbryo = bryophyte biomass, BMsed = biomass of all sedges, BMsph = biomass of Sphagnum mosses, Peat = peat layer thickness, Litter = litter percentage, CarRos = biomass of Carex rostrata, SoilT_y = average year-round soil temperature at −5 cm, SoilT_sf = average snow-free season soil temperature at −5 cm, and SoilT_s = average snow season soil temperature at −5 cm. WTD and porewater chemistry concentrations are values of the snow-free season.

Table B3 Accumulated annual (1 November 2021–31 October 2022, AccAnnual) and snow season (1 November 2021–12 May 2022 and 27–31 October 2022, AccSnow) methane flux rates in mg CH₄ m⁻² and the portion of snow season fluxes in the annual flux (%SnowOfAnnual) calculated for each plot from the manual static chamber measurements. Ordered by the highest portion of the annual flux.