Atmospheric Methane Consumption and Methanotroph Communities in West Siberian Boreal Upland Forest Ecosystems

: Upland forest ecosystems are recognized as net sinks for atmospheric methane (CH 4 ), one of the most impactful greenhouse gases. Biological methane uptake in these ecosystems occurs due to the activity of aerobic methanotrophic bacteria. Russia hosts one-ﬁfth of the global forest area, with the most extensive forest landscapes located in West Siberia. Here, we report seasonal CH 4 ﬂux measurements conducted in 2018 in three types of stands in West Siberian middle taiga–Siberian pine, Aspen, and mixed forests. High rates of methane uptake of up to − 0.184 mg CH 4 m − 2 h − 1 were measured by a static chamber method, with an estimated total growing season consumption of 4.5 ± 0.5 kg CH 4 ha − 1 . Forest type had little to no effect on methane ﬂuxes within each season. Soil methane oxidation rate ranged from 0 to 8.1 ng CH 4 g DW − 1 h − 1 and was negatively related to water-ﬁlled pore space. The microbial soil communities were dominated by the Alpha - and Gammaproteobacteria , Acidobacteriota and Actinobacteriota . The major group of 16S rRNA gene reads from methanotrophs belonged to uncultivated Beijerinckiaceae bacteria. Molecular identiﬁcation of methanotrophs based on retrieval of the pmoA gene conﬁrmed that Upland Soil Cluster Alpha was the major bacterial group responsible for CH 4 oxidation.


Introduction
Atmospheric level of methane, the second most important greenhouse gas after CO 2 , started to rise actively after a period of no growth in 2000-2006 [1][2][3]. The attribution of this trend to particular sources and sinks is still an unresolved issue for the scientific community. The most likely explanation is a combination of processes related to different components of the global methane budget [1]. Consumption in upland soils is one of these components and the only biological methane sink estimated in a range of 20-40 Tg yr −1 [2]. Recent top-down inventories suggested that a decreasing trend in soil methane uptake may partly explain the post-2006 renewed growth of atmospheric CH 4 and simultaneous decrease in the ratio of stable carbon isotopes of CH 4 ( 13 C/ 12 C) [4]. Thus, reducing uncertainty in soil methane consumption is important for understanding atmospheric global methane trends and future climate prediction. methane oxidation and composition of methanotroph communities in soils of the West Siberian middle taiga. This study was aimed to fill in the gap in our knowledge of the seasonal CH 4 fluxes from three typical middle taiga forest ecosystems of West Siberia, the controls staying behind observed seasonal (not inter-annual) and ecosystem scale methane consumption variability, as well as of the composition of methanotroph communities in taiga soils.

Site Description
Field measurements and soil sampling were conducted in 2018 in three typical middle taiga forest ecosystems of West Siberia: Siberian pine, mixed and small-leaved aspen forests. Study plots were situated several kilometers from each other on the second terrace of Ob' river near the city of Khanty-Mansiysk ( Figure 1). atmospheric methane oxidation and composition of methanotroph communities in soils of the West Siberian middle taiga. This study was aimed to fill in the gap in our knowledge of the seasonal CH4 fluxes from three typical middle taiga forest ecosystems of West Siberia, the controls staying behind observed seasonal (not inter-annual) and ecosystem scale methane consumption variability, as well as of the composition of methanotroph communities in taiga soils.

Site Description
Field measurements and soil sampling were conducted in 2018 in three typical middle taiga forest ecosystems of West Siberia: Siberian pine, mixed and small-leaved aspen forests. Study plots were situated several kilometers from each other on the second terrace of Ob' river near the city of Khanty-Mansiysk ( Figure 1). The climate of the region is subarctic Dfc according to Köppen climate classification, with long cold winters (average winter air temperature is −17.5 °C), short warm summers (average summer air temperature is 16.0 °C), and annual average air temperature of −0.4 °C for the 1981-2010 period (www.pogodaiklimat.ru, accessed on 5 November 2021). The annual average precipitation is 549 mm, concentrated in the period from June to October. Snow cover lasts for 187 days on average, from October to May (www.pogodaiklimat.ru, accessed on 5 November 2021). Ground water depth is more than 5 m. According to WRB classification, all plots have the same soil type-Albic Podzol with sandy loam texture. Soil profile consisted of a thin O layer (mean depths are 0-2 cm) of plants litter, a topsoil The climate of the region is subarctic Dfc according to Köppen climate classification, with long cold winters (average winter air temperature is −17.5 • C), short warm summers (average summer air temperature is 16.0 • C), and annual average air temperature of −0.4 • C for the 1981-2010 period (www.pogodaiklimat.ru, accessed on 5 November 2021). The annual average precipitation is 549 mm, concentrated in the period from June to October. Snow cover lasts for 187 days on average, from October to May (www.pogodaiklimat.ru, accessed on 5 November 2021). Ground water depth is more than 5 m. According to WRB classification, all plots have the same soil type-Albic Podzol with sandy loam texture. Soil profile consisted of a thin O layer (mean depths are 0-2 cm) of plants litter, a topsoil A horizon (2-7 cm), an eluviated E horizon (7-12 cm), an illuviated B horizon (12-20 cm) and a C horizon (parent material deeper than 20 cm). In Siberian pine forest, an E horizon is the most pronounced, while in mixed and aspen forest, it is intermittent. Soil properties of the sampling plots are given in Table 1. Siberian pine forest. A mature dense stand (crown cover of 90%-95%) is dominated by Siberian pine (Pinus sibirica Du Tour) interspersed with Siberian fir (Abies sibirica Ledeb.) and Siberian spruce (Picea obovata Ledeb.). The overstory trees have an average diameter of 40 cm at breast height and an average height of 27 m. The grass layer is sparse (projective cover is less than 10%) with Equisetum sylvaticum L. and Oxalis acetosella L. as main components. The moss layer is fragmentary in windthrow gaps and formed by Polytrichum commune and Pleurozium schreberi.
Aspen forest. A dense stand (crown cover of 60%-70%) with dominant aspen (Populus tremula L.) appeared after a uniform clear-cut that was conducted about 30 years ago [25]. Silver birch (Betula pendula Roth) is also common in overstory while Pinus sibirica forms sparse (1%-5%) understory. The overstory trees have an average diameter of 7 cm at breast height and an average height of 12 m. The dominant species in a grass-shrub layer (projective cover of 10%-15%) are Vaccinium vitis-idaea L., Equisetum sylvaticum and Calamagrostis canescens.
Mixed forest. A dense conifer-deciduous mixed stand (crown cover of 70%-80%) is formed by Pinus sibirica, Abies sibirica, Populus tremula. The overstory trees have an average diameter of 30 cm at breast height and an average height of 22 m. The understory is presented by the same species and Sorbus aucuparia L. The grass layer is sparse (projective cover of 10%-15%) and consists solely of Vaccinium myrtillus L. The moss layer (30%-50%) is formed by Hylocomium splendens and Pleurozium schreberi. Grass and moss layers are the most pronounced under deciduous species.

CH 4 Flux Measurements
CH 4 fluxes were measured using the static chamber method [49]. Three static chambers were randomly installed in each forest type. Field flux measurements were conducted three times per year in 2018: 25-28 of May (just after melting of seasonal soil frost), 10-15 of July (warmest week of the year), and 9-12 of September (beginning of the abscission). Methane fluxes were measured in three consecutive replicates for each chamber in each forest type between 12 and 4 p.m. with a total of 81 measurements (3 forest types × 3 times per year × 3 chambers × 3 replicates).
The chamber consisted of a permanently installed square stainless steel collar (37 cm × 37 cm, embedded 10-15 cm deep into the soil) and a removable plexiglas box (30 cm height). A groove on the collar rim-a hydro lock against leaks-was filled with water to a depth of 5 cm before the measurement. To minimize changes in ambient temperature, the box was covered with reflecting aluminum fabric. The air inside the chamber was mixed by a battery-operated internal fan. Initial pressure shock during the chamber setting was  Four gas samples were collected  0, 20, 40 and 60 min after closure by flushing gas-tight 20 mL polypropylene syringes  (KD-JECT III, KDM, Germany) 10 times with headspace air through a tube in a rubber stopper inserted tightly into the hole on the chamber top. After sampling, the syringes were immediately sealed with rubber stoppers and stored in the dark at +4 • C.
Methane concentration in the samples was analyzed in the laboratory within 48 h after sampling. Flux density in mg CH 4 m −2 h −1 was calculated using first order kinetics model as described in [17]. A methane concentration of 1.33 mg CH 4 m −3 was used to calculate a flux value as a mean value of initial chamber headspace CH 4 concentration for all measurements. Since methane fluxes usually have a non-normal distribution [5], a median of three consecutive flux measurements were used for statistical analysis.

Determination of the Soil Methane Oxidation Rate
The soil was sampled at each forest plot in each season at four depths (3-5, 10-12, 20-22 and 30-32 cm) in two spatial replicates between installed chambers (total of 3 forest types × 3 times per year × 4 depths × 2 replicates = 72 samples); for sampling, we used Edelman Soil Auger (Eijkelkamp Soil & Water, Netherlands). Samples were taken randomly from the soil core at the given depth, pooled in the 100 mL high-density polyethylene jars and stored in the dark at +4 • C for 2-3 days. Before incubation, the jars with soil were kept in the climate chamber MK-53 (BINDER, Germany) for 1 h to equilibrate to experimental temperatures. Then, 3-5 g of a sieved (mesh size 2 mm) soil was placed from the jars to 200 mL preliminary autoclaved glass bottles sealed with butyl stoppers. Initial methane concentration in bottles was ambient and varied from 0.92 to 1.21 mg CH 4 m −3 . These bottles were incubated at in situ soil temperature during the field sampling (±0.5 • C) in the same climate chamber. Four gas samples were taken from the bottle headspace by removing 2 mL of air at 3-5 h intervals. Soil CH 4 oxidation rate in ng CH 4 g DW −1 h −1 was calculated using first order kinetics model as described in [17] for methane concentration of 1.33 mg CH 4 m −3 . The soil from glass bottles was dried at 70 • C and the results were normalized to dry weight. For soil in each of the jars, incubation was conducted in 2-3 replicates; the average value of replicates was used for further calculations.

Methane Concentration
The methane concentration in gas samples from both chamber measurements and incubation experiments was determined using a gas chromatograph Kristall-5000 (Khromatek, Yoshkar-Ola, Russia) equipped with a flame ionization detector. Nitrogen served as a carrier gas with a flow rate of 35 mL min −1 . A stainless steel column with a length of 1 m and an internal diameter of 1 mm was filled with HayeSep Q (80-100 mesh) and held at 80 • C. Three external standards (2.28, 14.6 and 93 ppm, Ugra-PGS, Russia) were used for calibration each day after analysis. Precision (standard deviation) at 2.28 ppm was ±0.01 ppm for ten replicates. Each gas sample was analyzed twice; an averaged value was used for calculations of methane fluxes and oxidation rates.

Soil Chemical and Physical Properties
Soil temperature and volumetric moisture were measured in the field by Hydra Probe II (Stevens, CIIIA) in each soil core from 0 to 50 cm at 5 cm step. Total pore space was estimated by measuring volumetric water content in saturated samples. First, intact soil at 0 (surface), 10, 20 and 30 cm depths were sampled into 200 cm 3 rings in 2 spatial replicates at each study plot (after all measurements in September). Then, soil rings were submerged in a wide pan with foam rubber on the bottom for water infiltrating from both sides of the ring. The pan was accurately shaken for 10 h using an oscillating agitator to achieve full water saturation. Volumetric water content, measured by Hydra Probe II in saturated samples, was considered as the total pore space. Water-filled pore space was estimated as a ratio of in situ volumetric water content to the total pore space at the given depth. To determine pH and soluble NH 4 + and NO 3 − content of incubated soil samples, 1:6 slurries of soil and deionized H 2 O (w/v) were mixed using a vortex shaker for 5-10 s. pH was measured in the supernatant using a pH glass electrode (Mettler Toledo, Switzerland). After 2 h of extraction by the oscillating agitator, slurries were centrifuged for 15 min at 7000× g. Concentrations of NH 4 + and NO 3 − in the supernatant were assessed by the Nessler method using reagent set #2458200 and by cadmium reduction method using reagent set NitraVer ® 5, respectively (both sets-Hach-Lange, Düsseldorf, Germany), following the manufacturer guide. Quantification was performed using spectrophotometer DR5000 (Hach-Lange, Germany). The analytical accuracy of the method is ±10%.
Total organic carbon content in soil was measured using cuvette test system LCK 381 and spectrophotometer DR5000 (both Hach-Lange, Germany); it was estimated as a difference between total carbon and total inorganic carbon according to the manufacturer guide. Before measurements, 0.5 g of soil was homogenized with 5 mL of distilled water in a 50 mL propylene flask for 1 min using the oscillating agitator. The obtained mixture was immediately added to the test system in a volume of 50 µL. The analytical accuracy of the method is ±20%.

Statistical Analysis
We used 3-way ANOVA (function anovan in Matlab) with type III sum of squares [50] to evaluate the factors of depth, season and forest type, affecting the variability of methane flux (without depth as a factor) and soil CH 4 oxidation rate. Before the analysis, data were checked for normality by Anderson-Darling test (adtest, p > 0.05) and were Box-Cox transformed when applicable (boxcox). Omega-squared (corrected explained-to-total dispersion ratio) was used as an effect size measure [50]. Reported comparisons between groups were obtained using function multcompare with a significance level of 0.05 corrected by Tukey's honestly significant difference procedure [50].
We conceptualized the effect of different environmental controls with a help of regression modeling. To address non-linearity in controls (e.g., soil moisture [37]) we tested both linear and non-linear functions to explain soil CH 4 oxidation rate variability on a depth-, season-and ecosystem-specific basis. Linear, log-normal, power, exponential and 2nd order polynomial functions were fitted using function fitnlm. Reported models were checked for statistical significance (p < 0.05), parameters significance (p < 0.05), parameters non-zero condition (coefTest, p < 0.05) and normal distribution of residuals, estimated by Anderson-Darling test (adtest, p > 0.05). Both adjusted Pearson's correlation coefficient (R adj 2 ) and squared Spearman's rank correlation coefficient (Spearman 2 ) were used to estimate model performance (function corr).
Seasonal and ecosystem effects on models linking water-filled pore space and soil CH 4 oxidation rate were examined by non-linear mixed effect model; three seasons (May, July, September) and three forest types (Aspen, Siberian pine and mixed) were treated as random effect variables. Water-filled pore space was used as a predictor; it explained more variance in soil CH 4 oxidation rate than other tested controls. The calculation was made by means of function nlmefit. Obtained models were compared using the Akaike information criterion (AIC) and the likelihood ratio test (lratiotest).
Significant differences between groups were checked with a Mann-Whitney test (ranksum). All calculations were made in Matlab R2016b (The MathWorks, Natick, MA, USA). Soil methane consumption was assigned to negative flux values and positive values of oxidation rate.

Soil DNA Extraction
Extracts of total DNA used for molecular diversity studies were obtained from the samples collected from topsoil A horizon layers (3-5 cm depth). Four individual soil samples (of 0.5 g wet weight) were taken for the analysis from each of the studied forest sites and processed separately. Total DNA from samples was extracted using FastDNA SpinKit (MPBio, Santa Ana, CA, USA) according to the manufacturer's instructions.

Illumina Sequencing and Analysis of 16S rRNA Gene Fragments
To assess the microbial community composition in the forest soils including methanotrophic bacteria, fragments of 16S rRNA genes corresponding to the V4 region of were amplified from total DNA extracts. Libraries of the 16S rRNA gene fragments for highthroughput sequencing were prepared according to the earlier published protocol [51] and were sequenced on Illumina MiSeq platform (Illumina, San Diego, CA, USA) at Biospark (Moscow, Russia).

Illumina Sequencing and Analysis of pmoA Gene Fragments
Two different PCR assays were used to assess methanotroph diversity in the forest soils based on the retrieval of the pmoA gene fragments. The primer set A189/A682r [57] offers a broader coverage which covers USCα methanotrophs, but also targets amoA genes of ammonia-oxidizing bacteria. The primer set A189f/A650 was designed to target USCα methanotrophs [58]. The combination of these primer sets allowed in depth evaluation of methanotroph community composition in the forest soils. The obtained mixtures of amplicons were purified by agarose gel electrophoresis using Cleanup Standard Kit and sequenced on Illumina MiSeq platform (Illumina, San Diego, CA, USA) at Evrogen (Moscow, Russia). Demultiplexing and further processing of the resulting data set was carried out using QIIME 2 v.2020.8 package. DADA2 plugin was used for sequence quality control, merging of paired-end reads and chimera filtering. DADA2 outputs were translated to proteins using Framebot (http://fungene.cme.msu.edu/FunGenePipeline/framebot/form.spr, accessed on 4 October 2021) to detect and correct frameshifts in the reads. After correction step, the sequences were analyzed with QIIME 2 v.2020.8 Operational taxonomic units (OTUs) were clustered applying VSEARCH with dedicated reference database of 7809 unaligned pmoA nucleotide sequences with 86% identity [59,60].
Soil temperatures and moisture in May, July and September are similar to those in October, June and August, respectively, for studied plots [25]. Thus, we assumed the total methane sink in May, July and September to match the one in October, June and August, respectively. Thus, we estimated CH 4 consumption by studied soils at 4.5 ± 0.5 kg CH 4 ha −1 for the whole growing season. It is higher than 86% of estimates for boreal forest soils and 75% of estimates for all forest soils according to the database from [22]. Therefore, the taiga forest might be a strong methane sink at high latitudes.
Forest type affected CH 4 flux mostly in May, after the soil frost melt: consumption in Aspen forest was lower than in Siberian pine (p = 0.002) and mixed (p = 0.010) forests ( Figure 2). In July and September, soil sink showed no significant variability. Prior studies found the same seasonal emission patterns in temperate [62][63][64] and boreal [65,66] forest soils; forest type was reported both as significant [6,65,67] and non-significant [64,66,68] driver of methane consumption. Interaction of soil physical, chemical and biological factors trigger contrasting emission patterns. In studied soils, most variability in flux drivers was attributed to seasonal and depth scales; hence, forest type explained only a minor part of the total flux variance ( Table 2). Among ecosystem-scale controls, litter thickness tends to decrease soil methane consumption among different forest types [30,65], but it varied only slightly among study plots (Table 1). Soil pH and nitrogen content negatively correlated with methanotroph abundance in upland soil [69]. Both of these parameters were significantly lower in Siberian pine compared to Aspen forests (p = 0.048 and p = 0.009, respectively), but not to mixed forest (p = 0.161 and p = 0.156, respectively). It could potentially explain differences in consumption rates at studied forests in May.  0.04 * NS 0.53 *** * p < 0.05, ** p < 0.01, *** p < 0.001, NS, not significant; NM, not measured. a Data on methanotroph abundance are taken from the earlier published study [61].
Soil temperatures and moisture in May, July and September are similar to those in October, June and August, respectively, for studied plots [25]. Thus, we assumed the total methane sink in May, July and September to match the one in October, June and August, respectively. Thus, we estimated CH4 consumption by studied soils at 4.5 ± 0.5 kg CH4 ha −1 for the whole growing season. It is higher than 86% of estimates for boreal forest soils and 75% of estimates for all forest soils according to the database from [22]. Therefore, the taiga forest might be a strong methane sink at high latitudes.
Forest type affected CH4 flux mostly in May, after the soil frost melt: consumption in Aspen forest was lower than in Siberian pine (p = 0.002) and mixed (p = 0.010) forests (Figure 2). In July and September, soil sink showed no significant variability. Prior studies found the same seasonal emission patterns in temperate [62][63][64] and boreal [65,66] forest soils; forest type was reported both as significant [6,65,67] and non-significant [64,66,68] driver of methane consumption. Interaction of soil physical, chemical and biological factors trigger contrasting emission patterns. In studied soils, most variability in flux drivers was attributed to seasonal and depth scales; hence, forest type explained only a minor part of the total flux variance ( Table 2). Among ecosystem-scale controls, litter thickness tends
Water-filled pore space was the most powerful predictor for soil methane oxidation rate (dotted line in Figure 3, Table S1) with the power law as the best fit model. It substantiates the belief that soil methane consumption is diffusion limited during the growing season [37,62]. A power function of air-filled porosity is a commonly accepted approach to model soil gas diffusion [71,72]. Since air-filled porosity is a simple difference between total porosity and water-filled pore space, the latter works as a proxy for soil gas diffusion rate. Soil temperature and methanotroph abundance are significant predictors of the soil CH 4 oxidation rate as well, but they explain less variance (Table S1).
Water-filled pore space was the most powerful predictor for soil methane oxidation rate (dotted line in Figure 3, Table S1) with the power law as the best fit model. It substantiates the belief that soil methane consumption is diffusion limited during the growing season [37,62]. A power function of air-filled porosity is a commonly accepted approach to model soil gas diffusion [71,72]. Since air-filled porosity is a simple difference between total porosity and water-filled pore space, the latter works as a proxy for soil gas diffusion rate. Soil temperature and methanotroph abundance are significant predictors of the soil CH4 oxidation rate as well, but they explain less variance (Table S1). Potential cross-correlation of controls mask their individual effect on soil CH4 sink; its sensitivity to controls can also be ecosystem-specific [21,66]. To assess the effect, a nonlinear mixed effect model with a water-filled pore space as a fixed factor was compared with the power model (Table S2): all models performed similarly (for both, p = 0.36, df = 2 under likelihood ratio test), but the simple power model had lower (i.e., better) value of AIC. Thus, models had similar predictive power regardless of the data used: for each forest type (Figure 3), for each season ( Figure S1), for the whole dataset. Thus, the most important soil oxidation driver was the water-filled pore space; significant variations in temperatures and soil chemistry, throughout the season and among forest types, did not affect the soil CH4 oxidation rate and its sensitivity to soil moisture.
Being a biological process, methane consumption is controlled by temperature [37], but in forest soils, its temperature sensitivity could be negligible. Crill [62] reported Potential cross-correlation of controls mask their individual effect on soil CH 4 sink; its sensitivity to controls can also be ecosystem-specific [21,66]. To assess the effect, a nonlinear mixed effect model with a water-filled pore space as a fixed factor was compared with the power model (Table S2): all models performed similarly (for both, p = 0.36, df = 2 under likelihood ratio test), but the simple power model had lower (i.e., better) value of AIC. Thus, models had similar predictive power regardless of the data used: for each forest type (Figure 3), for each season ( Figure S1), for the whole dataset. Thus, the most important soil oxidation driver was the water-filled pore space; significant variations in temperatures and soil chemistry, throughout the season and among forest types, did not affect the soil CH 4 oxidation rate and its sensitivity to soil moisture.
Being a biological process, methane consumption is controlled by temperature [37], but in forest soils, its temperature sensitivity could be negligible. Crill [62] reported growing season Q 10 ranged from 0.9 to 1.4 with a mean value of 1.2 both for CH 4 flux and oxidation rate. Lind et al. [73] found no CH 4 flux dependence on soil and air temperatures at the seasonal scale in forests. Mean annual air temperature also did not correlate with CH 4 sink [22] or even correlated negatively [28] in boreal and temperate forest soils. We demonstrated that the gas diffusion controls methane consumption during the growing season; even in the early spring, when evapotranspiration is negligible and the soil retains snowmelt water, the effect is evident. Found correlation between temperature and soil CH 4 oxidation rate could be explained by high cross-correlation between temperature and water-filled pore space (Spearman 2 = 0.53).
Residuals for power fit model were not distributed normally (p = 0.008, N = 70 and p = 0.002, N = 22, respectively) across Siberian pine forest and all data. Outliers within medium water-filled pore spaces (denoted by red ellipse on Figure 3) could indicate the effect of other controls. We assessed them on a depth-specific basis as the depth explained most of the variance in the soil methane oxidation rate ( Table 2). Further regression modeling revealed its key drivers at a certain depth (Figure 4, Table S3).
Water-filled pore space was the only significant predictor in topsoil (3-5 cm) samples but had no effect at other depths ( Figure 4A). Methanotroph abundance strongly correlated with the soil methane oxidation rate in subsoil samples (20-22 and 30-32 cm) but had no effect in upper soil horizons ( Figure 4B). Methane concentration affected CH 4 oxidation rate only in samples from 10-12 cm ( Figure 4C); soluble NH 4 + and NO 3 − had no effect at any depth ( Figure 4D). Methane concentration and methanotroph abundance explained observed outliers (denoted by red ellipse on Figure 3). Another outlier caused strong correlation (R adj 2 = 0.82, Table S3) presented in Figure 4B, but Spearman 2 that is nonsensitive to them was also high (0.40) for this model. Divergent processes drive methane consumption within the soil profile. In the topsoil, methane concentration is high throughout the growing season because of the developed macropore system and close surface. Due to substrate availability, a population of methanotrophs reaches maximal abundance occupying a broad range of niches with different growth favorability [48,74]. Diversity of niches may result in high variability of methanotrophs abundance in the topsoil ( Figure 4B). High-affinity methanotrophs can also utilize other substrates [48,75] derived from root exudates. Thus, methane oxidation is limited by gas diffusion through the free pore space in soil aggregates.
In the subsoil (samples from 20-22 and 30-32 cm), methane concentration is low; it is close to the threshold for high-affinity methanotrophs [48,74]. Other substrates are less abundant than in the topsoil as well. Hence, their growth is limited for most of the season except periods of high methane availability, e.g., (i) temperature-induced lowered consumption in upper layers in early spring and late fall [62] and (ii) active methane production after heavy rain events [37]. Total methane consumption in the subsoil is limited by methanotroph abundance, which is significantly lower compared to 3 and 10 cm depths (p < 0.001), while water-filled pore space has medium values ( Figure 4A).
Methane consumption at the 10 cm depth occurs under intermediate conditions. Methane concentration is significantly less than in the topsoil (p < 0.001), but still higher than the threshold for high-affinity methanotrophs [48,74]. Methanotroph abundance drops comparing to the 3 cm depth (p < 0.001), but still exceeds the subsoil (p < 0.001). Soil methane oxidation rate negatively correlates with the methane concentration at this depth. If the oxidation were limited by the methanotrophic community, the correlation would be positive: higher CH 4 mixing ratios promote an extensive growth of methanotrophs. The observed negative correlation could be explained by preferential paths of gas transport through roots and macropores in the eluvial horizon (10 cm depth). Cm-scale heterogeneity of eluvial horizon in podzols is caused by the vertical migration of both organic and mineral substances from the topsoil [76]. Roots and macropores induce gas and water transport at the 10 cm depth, where a pore network is less developed compared to the topsoil [77]. Such preferential paths produce hot spots of microbial activity in soils [78]. In our study, macropores and hot spots might have been overlooked by measuring soil moisture and methanotroph abundance in small patches. In contrast, CH 4 concentrations integrate methane oxidation and production over the larger soil volume. We suggest that their small values correspond to active spots of microbial activity or preferential paths in the soil layer of 10 cm depth, where active methane oxidation takes place.

Prokaryote Diversity Patterns in Forest Soils
A total of 105 820 partial 16S rRNA gene sequences (mean amplicon length 250 bp) were retrieved from the examined forest soils. Of these, 89,533 reads were retained after quality filtering, denoising and removing chimeras. According to the alpha-rarefaction analysis, the richness of the examined samples has been fully sequenced and observed. The microbial diversity was highest in the mixed forest soil (average Shannon index 6.85 and Pielou evenness 0.90), followed by the aspen forest soil (average Shannon index 6.40 and Pielou evenness 0.89) and Siberian pine forest soil (average Shannon index 6.24 and Pielou evenness 0.89). As revealed by the UniFrac analysis and a further Permanova test, the microbial assemblages within the particular forest type were highly similar to each other but were significantly (p ≤ 0.05) different between the various forest types.
The pools of reads obtained from the examined forest soils were dominated by 16S rRNA gene sequences of bacterial origin. The relative abundance of archaeal 16S rRNA gene reads ranged from 0.086% to 0.12% of all sequences. Archaeal populations were represented by members of the Crenarchaeota (uncultured representatives of the class Nitrososphaeria).
The pools of reads obtained from the examined forest soils were dominated by 16S rRNA gene sequences of bacterial origin. The relative abundance of archaeal 16S rRNA gene reads ranged from 0.086% to 0.12% of all sequences. Archaeal populations were represented by members of the Crenarchaeota (uncultured representatives of the class Nitrososphaeria).
The bacterial communities in the three forest sites were dominated by members of the Alphaproteobacteria (22-35% of the total number of 16S rRNA gene fragments), Acidobacteriota (15-35%), Actinobacteriota (12-28%), and Gammaproteobacteria (10-21%). These four major bacterial groups comprised 76-88% of the total prokaryote diversity revealed in these soils ( Figure 5). This pattern of bacterial diversity is characteristic for acidic forest soils [79]. Minor bacterial groups included Verrucomicrobiota (1.5-4.2%), Planctomycetota (1.5-3.6%), Bacteroidota (1.2-4.6%), Myxococcota (0.2-3.4%), Chloroflexi (0.4-4.8%), Gemmatimonadota (0.2-1.7%) and others. The spectrum of minor bacterial groups also included the Candidate group RCP2-54, which was especially well represented (5% of the total number of 16S rRNA gene fragments) in the soil of Siberian pine forest. This as-yet-uncultivated group of bacteria was named after the environmental 16S rRNA gene sequence (GenBank accession number AF523886), which was retrieved in 2002 from a forested wetland [80]. The number of species-level OTUs determined at 97% sequence identity ranged between 295 in the Siberian pine forest soil and 488 in the mixed forest soil. The most abundant OTUs identified in each of the forest sites are listed in Table 3. Notably, one particular OTU that exhibited 100% sequence similarity to Bradyrhizobium sp. 175LB2PYPT obtained from the root nodules of Acacia pycnantha (GenBank accession number HQ698291), was identified in all forest sites in high relative abundances (5.7, 6.2 and 9.1% of the total number of 16S rRNA gene fragments retrieved from Siberian pine, mixed and aspen forest soils, respectively). Members of the genus Bradyrhizobium are dinitrogen-fixing soil bacteria, which commonly occur in association with various plants. Several other major OTUs  Table 3. Notably, one particular OTU that exhibited 100% sequence similarity to Bradyrhizobium sp. 175LB2PYPT obtained from the root nodules of Acacia pycnantha (GenBank accession number HQ698291), was identified in all forest sites in high relative abundances (5.7%, 6.2% and 9.1% of the total number of 16S rRNA gene fragments retrieved from Siberian pine, mixed and aspen forest soils, respectively). Members of the genus Bradyrhizobium are dinitrogen-fixing soil bacteria, which commonly occur in association with various plants. Several other major OTUs affiliated with the Acidobacteriia (Subdivisions 1 and 2), which are characteristic inhabitants of soils rich in plant-derived organic matter [81]. All major gammaproteobacterial OTUs belonged to the order-level group WD260, which does not contain cultivated representatives. This group was named after the environmental clone sequence WD260 (GenBank accession number AJ292673) retrieved two decades ago from an acidic polychlorinated biphenylpolluted soil near Wittenberg, Germany [82]. Members of WD260 group are common inhabitants of soils and peatlands. Despite the recent success in culturing various groups of elusive soil bacteria, members of WD260 soil group resist cultivation efforts till now. No information is available about the physiology of these bacteria. Acid-tolerant actinobacteria of the genus Mycobacterium and uncultured members of the family Solirubrobacteraceae were also among the most abundant OTUs identified in all three forest types. Thus, with exception of Bradyrhizobium and Mycobacterium species, all most abundant OTUs in studied forest soils represented as-yet-uncultivated groups of bacteria.

Identification of Methanotrophs Based on 16S rRNA Gene Analysis
The search for OTUs that represent well-studied methanotrophic bacteria revealed a single OTU, which was affiliated with the order Methylococcales. This OTU was detected in a low relative abundance in soils of mixed and aspen forests and was absent from Siberian pine forest soil ( Figure 6A,B). The only relatively abundant group of reads (~2.4% of all bacterial 16S rRNA gene fragments) that could potentially belong to methanotrophs was classified as uncultivated Beijerinckiaceae bacteria ( Figure 6A,B). The latter were represented by seven species-level operational taxonomic units (OTUs), with three OTUs shared between the three forest types ( Figure 6C). These Beijerinckiaceae-affiliated phylotypes displayed 93.2% and 100% sequence similarity to 16S rRNA gene sequence of Candidatus Methyloaffinis lahnbergensis [47] and 93.2% and 98.2% sequence similarity to 16S rRNA gene sequence of 'Methylocapsa gorgona' MG08 [48], respectively.

PmoA-Based Identification of Methanotrophic Bacteria
A total of 212,869 and 44,964 partial pmoA gene sequences (mean amplicon length 300 bp) were retrieved from the examined forest soils using A189/A682r and A189f/A650 primer sets, respectively. Of these, 3170 and 4379 pmoA gene sequences (mean amplicon length 530 bp) were retained after quality filtering, denoising, merging paired-end sequences, removing chimeras, removing amoA gene sequences and frameshift corrections via Framebot. All reads were clustered at 86% nucleotide similarity [59,60]. The obtained pool of pmoA fragments was represented by 23 species-level OTUs, three of which were common for all forest types ( Figure 7A). The highest number of the pmoA-based OTUs was detected in the mixed forest soil, followed by the Siberian pine forest soil and aspen forest soil. All of these fragments belonged to the alphaproteobacterial methanotrophs of the family Beijerinckiaceae and displayed 86.7%-99.4% sequence identity to the pmoA of Candidatus Methyloaffinis lahnbergensis, which was identified in a deciduous forest soil in Germany [47]. Notably, OTUs 22 and 23, which displayed highest similarity to the pmoA of Candidatus Methyloaffinis lahnbergensis, were detected only in the soils of Siberian pine and mixed forests. Of the 23 OTUs identified in the total pool of pmoA sequences, 20 OTUs were obtained using the primer set A189f/A650r, while application of the primer set A189f/A682r allowed retrieval of 15 OTUs only ( Figure 7B). all bacterial 16S rRNA gene fragments) that could potentially belong to methanotrophs was classified as uncultivated Beijerinckiaceae bacteria ( Figure 6A,B). The latter were represented by seven species-level operational taxonomic units (OTUs), with three OTUs shared between the three forest types ( Figure 6C). These Beijerinckiaceae-affiliated phylotypes displayed 93.2% and 100% sequence similarity to 16S rRNA gene sequence of Candidatus Methyloaffinis lahnbergensis [47] and 93.2% and 98.2% sequence similarity to 16S rRNA gene sequence of 'Methylocapsa gorgona' MG08 [48], respectively.  As revealed by molecular analysis, microbial communities in the studied Siberian boreal forest soils have a large proportion of bacteria belonging to as-yet-uncultivated groups, such as the Candidate phylum RCP2-54 or the order-level group WD260 of the Gammaproteobacteria. No information is currently available regarding the biology and environmental functions of these bacteria; the latter are attractive objects for further cultivation efforts or metagenome analyses. Methanotroph populations in these soils are represented by uncultivated Beijerinckiaceae bacteria of the USCα clade. These types of methanotroph communities, which are composed exclusively of high-affinity USCα methanotrophs, were earlier reported for soils of a beech forest in Denmark, a rainforest in Brazil, a mixed hardwood forest in the United States [39], a sub-boreal pine forest in Canada [83], temperate forests with European beech and Norway spruce in Germany [44], a lichen-dominated pine forest of Russian tundra [14] and a wide range of forest soils in China [46]. Given the large areas occupied by the upland forests in the world, USCα methanotrophs appear to represent one of the most environmentally relevant methanotroph populations in terrestrial ecosystems.

PmoA-Based Identification of Methanotrophic Bacteria
Overall, our study characterized West Siberian boreal upland forest ecosystems as a strong, earlier underestimated sink for atmospheric methane, which is oxidized by USCα methanotrophs.
the family Beijerinckiaceae and displayed 86.7-99.4% sequence identity to the pmoA of Candidatus Methyloaffinis lahnbergensis, which was identified in a deciduous forest soil in Germany [47]. Notably, OTUs 22 and 23, which displayed highest similarity to the pmoA of Candidatus Methyloaffinis lahnbergensis, were detected only in the soils of Siberian pine and mixed forests. Of the 23 OTUs identified in the total pool of pmoA sequences, 20 OTUs were obtained using the primer set A189f/A650r, while application of the primer set A189f/A682r allowed retrieval of 15 OTUs only ( Figure 7B). As revealed by molecular analysis, microbial communities in the studied Siberian boreal forest soils have a large proportion of bacteria belonging to as-yet-uncultivated groups, such as the Candidate phylum RCP2-54 or the order-level group WD260 of the Gammaproteobacteria. No information is currently available regarding the biology and environmental functions of these bacteria; the latter are attractive objects for further cultivation efforts or metagenome analyses. Methanotroph populations in these soils are represented by uncultivated Beijerinckiaceae bacteria of the USCα clade. These types of Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/f12121738/s1, Figure S1: Power models fitted for soil CH 4 oxidation rate versus water-filled pore space for all data and for each season separately, Table S1: Regression models for soil CH 4 oxidation rate (N = 72), Table S2: Comparison of simple regression model and mixed effect models with season and forest type as tested random variables, Table S3: Regression models at the specified depth, depicted on a Funding: This study was supported by the Russian Science Foundation (RSF), with field studies and laboratory analyses of methane oxidation rates supported by the RSF project no. 19-77-10074, and bioinformatic analyses of microbial and methanotroph diversity supported by the RSF project no. 21-14-00034.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: The 16S rRNA and pmoA gene reads retrieved using Illumina sequencing from the studied forests soil have been deposited under the Bioproject number PRJNA779228 in the NCBI Sequence Read Archive, with the accession numbers SAMN23019536, SAMN23019537, SAMN23019538 for rRNA dataset and SAMN23019615, SAMN23019616, SAMN23019617 for pmoA dataset.