Changes in Soil Microbial Community and Carbon Flux Regime across a Subtropical Montane Peatland-to-Forest Successional Series in Taiwan

Abstract

Subtropical montane peatland is among several rare ecosystems that continue to receive insufficient scientific exploration. We analyzed the vegetation types and soil bacterial composition, as well as surface carbon dioxide and methane fluxes along a successional peatland-to-upland-forest series in one such ecosystem in Taiwan. The Yuanyang Lake (YYL) study site is characterized by low temperature, high precipitation, prevailing fog, and acidic soil, which are typical conditions for the surrounding dominant Chamaecyparis obtusa var. formosana forest. Bacterial communities were dominated by Acidobacteriota and Proteobacteria. Along the bog-to-forest gradient, Proteobacteria decreased and Acidobacteriota increased while CO₂ fluxes increased and CH₄ fluxes decreased. Principal coordinate analysis allowed separating samples into four clusters, which correspond to samples from the bog, marsh, forest, and forest outside of the watershed. The majority of bacterial genera were found in all plots, suggesting that these communities can easily switch to other types. Variation among samples from the same vegetation type suggests influence of habitat heterogeneity on bacterial community composition. Variations of soil water content and season caused the variations of carbon fluxes. While CO₂ flux decreased exponentially with increasing soil water content, the CH₄ fluxes exhibited an exponential increase together with soil water content. Because YYL is in a process of gradual terrestrialization, especially under the warming climate, we expect changes in microbial composition and the greenhouse gas budget at the landscape scale within the next decades.

1. Introduction

Soil is a highly complex ecosystem with numerous macrobial and microbial taxa. Microbes are the main driver of biogeochemical cycling in soils. Bacteria, archaea, and fungi are the main components in the diverse biochemical processes of organic matter turnover. Soils can differ substantially in their physicochemical characteristics, which may affect major microbial physiological properties such as survival and growth, as well as ecological properties such as distribution and dispersal [1]. In turn, this can affect the microbial composition and consequently ecological functions.

Forests store large amounts of carbon in the form of living biomass, as well as labile and recalcitrant plant-derived polymers accumulated in soils. Forest soil accounts for a total of 1699 Gt carbon and ranks as the largest of all terrestrial carbon storages [2]. Forest soil bacterial communities are usually dominated by Acidobacteriota (formerly Acidobacteria), Actinobacteriota (formerly Actinobacteria), Proteobacteria, Bacteroidota (formerly Bacteroidetes), and Firmicutes [3]. Acidobacteriota (Acidobacteria), Actinobacteriota (Actinobacteria), and Proteobacteria have been shown to comprise over 80% of the actual active community in the soil of temperate forests, on the basis of rRNA-based analysis [4,5]. Proteobacteria and Bacteroidota (Bacteroidetes) have been shown to exhibit a preference for using labile carbon substrates [6]. Soil pH is usually the dominant factor influencing the structure of the soil bacterial community, with other important factors including organic matter content, nutrient availability, climate conditions, and biotic interactions with plants [7]. Acidobacteriota usually dominate at lower pH, and Firmicutes dominate at higher pH [8].

Microbial taxa lack a distinct global distribution pattern. However, the compositions and functions of dominant phyla in peatland play a crucial role in current and changing global carbon cycles and therefor require scientific investigation [9]. Peatlands, defined as exhibiting a peat layer of >30 cm organic matter [10], occupy only 2.84% of the total land area [11] but store approximately 612 Gt of carbon [12]. Given the total terrestrial soil carbon of 2273 Gt [2], the relatively small land area of peatland stores 27% of it and is therefore a vital component of the global carbon pool. The microbial communities in peatland may therefore be crucial factors in global climate change. Proteobacteria, Acidobacteriota (Acidobacteria), and Chloroflexi were reported to be the dominant phyla in peatland [13]. In tropical peat swamp, Acidobacteriota (Acidobacteria), Proteobacteria, Actinobacteriota (Actinobacteria), and Firmicutes have been reported as the main phyla [14].

In peatlands, the organic matter remains mostly in anoxic condition, which on the one hand suppresses the aerobic decomposition of organic matters, and on the other hand, facilitates the production of methane by archaeal methanogens [15]. In a recent study of the global methane budget, Saunois et al. [16] estimated the atmospheric source of methane to account for 550 to 593 Tg CH₄ yr⁻¹, 31.4% of which was emitted by natural wetlands (peatlands and floodplains). The remainder of the atmospheric methane originated predominantly from anthropogenic sources such as agriculture, waste, and fossil fuels. On the drier part of the terrestrial lands, the soils acted oppositely as a net sink of atmospheric methane. Although the spatial extent of the land source and sink of atmospheric methane exhibits a seemingly clear boundary, most of the soils, regardless of whether they are situated in peatlands or in upland forests, exhibit a complex three-dimensional structure that constitutes a large spatiotemporal heterogeneity of redox conditions [17]. The methane produced in the soil can therefore be released to the atmosphere, or it can be consumed by diverse methylotrophic bacteria along the diffusion pathway [18].

Among upland forests, those that are frequently immersed in fog may exhibit wetland-like soil conditions, which are the result of high atmospheric water input from precipitation and fog deposition [19] and low evaporation [20]. In tropical mountains, forest soils at elevations with higher precipitation were demonstrated to have a weaker methane sink strength compared with soils at other altitudes [21]. Peatlands located within cloud forests may therefore exhibit a less distinct transition from peatland to the adjacent upland forests. However, microbial communities and the greenhouse gas fluxes in tropical montane cloud forests containing peatlands are insufficiently studied and therefore poorly understood.

For this study, we investigated the microbial communities and carbon fluxes along a peatland–forest continuum in the subalpine Yuanyang Lake (YYL) Nature Reserve in northern Taiwan. YYL is an upstream small and shallow lake surrounded by a subtropical montane cloud forest. The succession from open water through marsh, bog, and swamp to the upland forest within a short distance allowed an investigation under similar micro meteorological conditions. The adjacent Chi-Lan Mountain long-term ecosystem research site provides background information about nutrient cycling, fog meteorology, fog deposition, and flux tower measurements of the forest ecosystem. Being classified as montane cloud forest [22], the site receives 300 mm yr⁻¹ fog deposition in addition to abundant precipitation [23,24]. This high soil water content in combination with a relatively low air temperature (annual mean approximately 13 °C) are hypothesized to lead to a low decomposition rate of aboveground litter [25], and an extremely low soil respiration rate [26].

Soils of the YYL upland forest exhibited lower diversity of bacterial communities than those of the nearby disturbed forest and were dominated by Proteobacteria, Acidobacteriota (Acidobacteria), and Actinobacteriota (Actinobacteria) [27,28]. Chang et al. [29] studied the soil microbiota following a disturbance gradient from a natural Chamaecyparis forest, disturbed Chamaecyparis forest, and a secondary Chamaecyparis forest, to a reforested secondary Cryptomeria forest. Bacterial and fungal compositions changed along this gradient, as determined on the basis of denatured gradient gel electrophoresis and phospholipid fatty acid analysis. However, no information exists on microbial communities along the vegetation succession series from the peatland to the upland forest. Similarly, whether the greenhouse gas fluxes vary because of the change of the soil water regime is unclear regarding this subtropical montane landscape.

The objectives of this study were thus to investigate the microbial composition and the greenhouse gas (carbon dioxide and methane) flux in a series of vegetation types from peatland to upland forest in a little-explored subtropical montane area. Methane production can be influenced by microbial activities because methane can be produced only by archaeal methanogens, some of which are consumed by methanotrophic bacteria. We hypothesized that methane release would be accompanied by increases in methanotrophic bacteria. Carbon dioxide is produced mainly through heterotrophic respiration. Copiotrophic heterotrophic bacteria, which are represented by Proteobacteria, tend to grow faster when nutrient availability is high, which is expected in bog and swamp. However, CO₂ production is controlled by oxygen availability; therefore, higher CO₂ production is generally expected in forest.

2. Materials and Methods

2.1. Site Description

This research was conducted on the peatland and its surrounding upland forest of the subalpine YYL watershed. The lake is located at an elevation of 1650 m above sea level in northern Taiwan (121°24′ E, 24°35′ N) (Figure 1). The warm temperate climate of this site is characterized by a low annual mean air temperature of 12.7 °C, and a high annual precipitation of up to 3396 mm. The prevailing foggy condition of up to 40% of the year [22,30] adds approximately 300 mm yr⁻¹ water to the ecosystem [24]. YYL is a humic mesotrophic shallow lake with a surface area of 3.6 ha and maximum water depth of 4.5 m [31]. The surrounding peatland area is approximately 2.2 ha and the upland forest of the watershed is approximately 374 ha.

Sparganium fallax and Schoenoplectus mucronatus are the two dominant water plants in the lake [32]. The peatland is covered by Sphagnum mosses and herbaceous plants and gradually transitioning to the pristine Taiwan yellow false cypress (Chamaecyparis obtusa Sieb. & Zucc. var. formosana (Hayata) Rehder, hereafter referred to as “yellow cypress”) forest. Yellow cypress is an endemic, giant, and long-lived tree species and was one the most important timber sources in Taiwan [33]. The soil of the peatland is water-saturated Histosol. The soils of the upland forest are Typic Hapludult (Ultisols), Typic Dystrochrept (Inceptisols), and Lithic Medihemist (Histosols) on the summit, footslope, and toeslope, respectively [34].

2.2. Plant Community Survey

The plant communities were surveyed along a transect line from the open water region to swamp forest in September 2017. The physiognomic group of the vegetation types on the peatland was identified following the classification of peatland ecosystem [10]. All soil at the lakeshore was saturated or flooded with water. Marsh was the place permanently flooded by lake water; bog was the place not flooded by lake water and that exhibited stunted and sparsely located trees; the swamp forest was dominated by trees (

Table 1. Vegetation composition variables of the 10 vegetation types in the Yuanyang Lake (YYL) Nature Reserve and nearby forest.

Vegetation Type	A	B	C	D	E	F	G	H	I *	J
Ecosystem	Peatland								Forest
Physiognomic group	Marsh		Bog				Thicket swamp forest	Conifer swamp forest	Upland in watershed	Outside watershed
Water level (cm)	46.0	36.8	2.5	2.8	7.1	2.0	0	0	0	0
Plant species number	1	1	5	4	5	4	14	6	-	15
Wetland plant species	1	1	3	3	4	3	3	0	0	0
Percent wetland species	100	100	60	75	80	75	21.4	0	0	0
Herbaceous species number	1	1	5	4	5	4	7	4	-	5
Percent herbaceous species number	100	100	100	100	100	100	50	67	-	33
Coverage of wetland species (%)	84	82	78	28	90	74	25	0	0	0
Coverage of herbaceous species (%)	84	82	100	100	100	100	100	3	-	23
Mosses and lichen (% ground area)	0	0	19.7	19.3	65.7	22.7	41.0	30.0	-	83.0
Woody species number	0	0	0	0	0	0	7	2	42	10
Woody species density (individuals ha−1)							9333	1200	6688	2367
Woody species basal area (m2 ha−1)	0	0	0	0	0	0	18.4	86.3	58.6	54.5
Yellow cypress density (individuals ha−1)							0	300	848	1367
Individual yellow cypress basal area (cm2)							0	1667	480	360
Microbial sampling plot						F	G	Ha,b,c	I	Ja, b
CO2 and CH4 flux sampling plot	A				E	F	G	H		J

* Values were cited from [35].

). The vegetation types A to H were identified on the basis of patches of apparently distinct dominant plant species. Two vegetation types of upland forest ecosystem were included in this study: the vegetation type I within the YYL watershed and the vegetation type J located 1.5 km away from the watershed, where the southeastern aspect of the wider slope allows high water vapor exchange through advection. Both forest sites are dominated by yellow cypress trees. Three separated 10 × 10 m² quadrats for vegetation type H and J, three separated 5 × 5 m² quadrats for vegetation type G, and three 1 × 1 m² separated quadrats for each vegetation types in the marsh and bog were established for detailed plant species and abundance surveys. In each quadrat, the diameter at breast height of the tree trunk was measured for all the woody species and the coverage (the percentage of plant projected area to the ground area) was recorded for all the herbaceous species. Total mosses and lichen coverage on the ground surface was surveyed in each quadrat without distinguishing their species. Regarding the mosses and lichens, distinct herbaceous species were distributed at different height layers in a quadrat; thus, they could shade each other and render the total coverage of all species in a sampling quadrat over 100%. We considered that recording the coverage of each species by the ratios of their projected area to the ground area presented the real situation more accurately than minimizing the portion of all species to fit the total coverage of 100%. This method is also more precise when comparing the coverage of one species across several quadrats. The water level of each quadrat was averaged by measuring the water depth at two opposite corners. Because detailed vegetation surveys have already been reported for the upland forest in the YYL watershed [32,35], we did not survey the vegetation type I.

2.3. Microbial Characterization

Samples for microbial characterization were collected in September 2017. We collected soil samples from eight sampling plots corresponding to vegetation type F, Ga, Gb, Gc, H, I, Ja, and Jb (Table 1). Ha, Hb, and Hc are three nearby plots along a water content gradient (Ha > Hb > Hc). Ja and Jb are two nearby plots of Ja (mature forest plot) located downslope of Jb (young forest plot). At each plot, three nearby soil samples (within 1 m) were collected as replicate samples. After removing covering leaves and twigs, a core of the soil was sampled with a stainless-steel corer (7 cm in diameter), and the top 2–5 cm soil were stored in a plastic jar. In some sampling plots, the depth of the soil layer was <10 cm to the rock bottom; therefore, we decided to use the 2–5 cm surface soil (organic layer) as our samples for analysis. F samples contained organic matter with plant litter taken from the top layer of peat, approximately 50 cm below water, because it was prohibitively difficult to take samples from a proper depth with a soil corer. These samples were transferred to a laboratory at Tzu-Chi University one day after collection, aliquoted, and stored at −20 °C until further use.

Before DNA extraction, each soil sample was mixed and sieved with a metal net (pore size 1 mm) to remove plant debris. A 0.25 g fraction was taken from each sieved soil sample for DNA extraction using PowerSoil kit (MoBio, Carlsbad, CA, USA), following the manufacturer’s instructions. The 16S rDNA V3-V4 fragment was amplified using primer pair 341F-805R. The first polymerase chain reaction (PCR) was performed for 25 cycles using 341F and 805R primers, with DNA extracted from soil as the template. We added the Nextera style tags P5 and P7 to the 5-end of amplicons with five cycles of second PCR runs using P5-341F and P7-805R primers. Eight μL of the first run product was used as the template for the second run in a 40 μL PCR mixture. Ex Taq DNA polymerase (Takara Bio Inc., Shiga, Japan) with proofreading activity was used in all PCR reactions. The PCR thermal procedure was one cycle at 95 °C for 3 min, 25 cycles at 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s, and then 72 °C for 5 min. Amplicons of multiple PCR reactions from identical DNA samples were pooled and purified using a GenepHlow Gel/PCR Kit (Geneaid, New Taipei City, Taiwan) and then sent to NGS High Throughput Genomics Core Laboratory, Biodiversity Research Center, Academia Sinica, for library construction and pair-end sequencing using an Illumina MiSeq platform. A total of 24 libraries from each of the 24 soil samples were barcoded and pooled using a Nextera DNA preparation kit (Illumina, San Diego, CA, USA). The sequences generated in this project are available through NCBI under BioProject number PRJNA846692.

Basic sequence processing was performed using Mothur version 1.41 [36,37], according to MiSeq SOP provided on the Mothur website (https://mothur.org/wiki/miseq_sop/) (accessed on 1 May 2022). Usearch version 11 [38] was used to define 97% of operational taxonomic units (OTUs) and construct the OTU table. We used Mothur to classify OTUs on the basis of the SILVA database [39] version 138.1 (SILVA 138). Sequences assigned to Archaea, Eukaryote, mitochondria, and chloroplast, and unclassified sequences were removed from further analysis. SILVA 138 uses the new naming system, following the decision made by International Committee on Systematics of Prokaryotes [40]. In the remainder of this paper, we use the new phylum names and put the old ones in parenthesis when citing previous studies.

2.4. Analysis of Microbial Community

Statistical analyses of bacterial community data were conducted on the MicrobiomeAnalyst website (https://www.microbiomeanalyst.com) (accessed on 1 May 2022) [41,42]. We used the default setting for data filtering: for the low-count filter we set the minimum count to four and prevalence at 20%; for the low-variance filter set the percentage for removable at 10% on the basis of the inter-quantile range. No rarefaction and transformation were applied, and total sum scaling was used to adjust for distinct sequencing depth differences. MicrobiomeAnalyst was used for the following analyses: alpha diversity (Chao1, Shannon index, OTU richness), beta diversity (PCoA), and LEfSe [43].

Methanotrophs were determined on the basis of their phylogenetic affiliation. We classified taxa in Methylomirabilota (formerly NC10), Methylobacterium-Methylorubrum, Methylocella and Methylocystis (Beijerinckiaceae, Alphaproteobacteria), Methylocystaceae (Alphaproteobacteria), Methylococcaceae (Gammaproteobacteria), Methylomonadaceae (Gammaproteobacteria), Methylothermaceae (Gammaproteobacteria), and Methylacidiphilaceae (Verrucomicrobia) as putative methanotrophs for further analysis, according to previous studies [18,44]. Known methylotrophic taxa Methyloligellaceae (Alphaproteobacteria) and Methyloversatilis (Gammaproteobacteria) were used for comparison.

2.5. Surface CO₂ and CH₄ Fluxes

To understand the spatial variation of surface carbon fluxes along the gradient of water regime from the marsh, through swamp, to the upland forest, the CO₂ and CH₄ fluxes were measured at selected vegetation types along the gradient. Measurements at five vegetation types were conducted between the summer of 2017 and 2019, on the marsh (vegetation type A, (Table 1), in the bog (two vegetation types; E and F, Table 1), and in the swamps (G and H, Table 1). For each vegetation type, three randomly chosen points were used for measuring carbon flux. Twelve measurement campaigns were conducted during this period (7 September 2017, 27 November 2017, 25 January 2018, 30 May 2018, 30 March 2019, 31 July 2019, 23 September 2019, 5 November 2019, 1 January 2020, 7 March 2020, 9 May 2020, and 11 July 2020). By contrast, the measurements for vegetation types E, F, and G were performed only between 7 September 2017 and 31 July 2019. Because equipment had to be repaired, a large measurement gap exists between June 2018 and February 2019. For the young yellow cypress forest (vegetation type J, Table 1), the soil CO₂ and CH₄ fluxes were measured on points along two 100 m transects with consecutive points separated by 10 m. Four measurement campaigns were conducted (8 September 2017, 26 November 2017, 26 January 2018, and 25 May 2018).

2.5.1. Flux Measurement

The exchange of CO₂ and CH₄ between the soil and atmosphere was measured using a static closed chamber system. The system comprises a laser-based spectroscopy (Ultraportable Greenhouse Gas Analyzer, Los Gatos Research, Inc., San Jose, CA, USA) (hereafter UGGA) and a self-made transparent chamber. The cylinder-shaped chamber has two parts: the upper part is made of acrylic and measures 15 cm in height; the lower part is a stainless-steel cylinder and has an inner diameter of 15.2 cm and height of 20 cm. For measurements, the lower part was first carefully inserted into the soil. However, in the marsh sites where Trichophorum subcapitatum or S. fallax grew densely on the soil surface, the plant bodies along the outer rim of the cylinder were cut to allow the proper insertion of the cylinder into the soil. The upper part was then set on top, ensuring the cylinder was thoroughly inserted into the acrylic chamber. A rubber gasket fixed on the bottom of the inner wall of the chamber ensures a gas-tight connection between the two parts. For each measurement, the height of the headspace was measured for the calculation of the total chamber volume. Two PVC tubes (4 mm internal diameter) connecting UGGA and the chamber enable the inflow of chamber headspace air to UGGA and the outflow of measured air back to the chamber. The UGGA gas analyzer has an internal gas pump that circulates the air between the UGGA and the chamber headspace. The flow rate of the air was fixed at 1 L min⁻¹.

Soil volumetric water content and soil temperature of the top 10 cm of soil were recorded using portable equipment (HH2 moisture meter and WET-2 sensor, Delta-T Devices, UK). During the measurement of soil carbon flux, the WET-2 sensor was inserted into the soil proximal to the chamber.

2.5.2. Analysis of Carbon Flux Data

For each measurement, the changing CO₂ and CH₄ concentrations in the chamber head space were calculated to CO₂ flux in μmol m⁻² s⁻¹ and CH₄ flux in nmol m⁻² s⁻¹. The first 30 s of the 120 s measurement period were excluded from the flux calculation to minimize the influence of the initial unstable gas concentrations.

To explore the seasonal dynamics of the soil CO₂ and CH₄ fluxes, the measurements were integrated into four seasons, resulting in a data frame consisting of 24 groups (six vegetation types × four seasons). For each group, the outliers of the CO₂ and CH₄ fluxes were separately determined using the median and median absolute deviation (MAD) method with the upper (T_max) and lower (T_min) thresholds [45]:

$${\mathrm{T}}_{\max },{ \mathrm{T}}_{\min }= \mathrm{median}\left(\mathrm{x}\right) \pm 3\times \mathrm{MAD}$$

(1)

$$\mathrm{MAD} =1.4826\times \mathrm{median}\left(\left|{\mathrm{x}}_{\mathrm{i}}- \mathrm{median}\left(\mathrm{x}\right)\right|\right)$$

(2)

where x_i is the individual measurement within one group. A total of 13.3% of the CH₄ and 5.2% of the CO₂ fluxes were defined as outliers and excluded in the following analysis.

The influence of seasons and vegetation types on the soil CO₂ and CH₄ fluxes was explored using two-way analysis of variance (ANOVA) via the aov() function of R [46].

3. Results

3.1. Plant Community

The vegetation survey revealed a total 34 vascular plant species from 25 families, among which 17 were woody and 17 were herbaceous species. A thorough vegetation inventory in the upland forest of the YYL watershed (vegetation type I in this study), revealed 42 tree species and the basal area values of the dominant 11 species [35] (

Table 2. Dominant plant species in the 10 vegetation types in this study. Dominance of the woody species was determined on the basis of the basal area of the ground area (m² ha⁻¹) and the herbaceous species by their coverage of ground area (%). * Data for vegetation type I were cited from [35].

Life Form	Species	A	B	C	D	E	F	G	H	I *	J
Woody species
	Chamaecyparis obtusa var. formosana								50.0	40.7	49.0
	Rhododendron formosanum								36.3	9.8
	Illicium anisatum									3.5	0.13
	Adinandra formosana var. formosana									1.1
	Schefflera taiwaniana									0.3
	Barthea barthei									0.2
	Viburnum sympodiale									0.1
	Dendropanax dentiger									0.7
	Skimmia arisanensis									0.04
	Ternstroemia gymnanthera									0.8
	Neolitsea acuminatissima									0.1
	Hydrangea paniculata							14.02
	Rhododendron chilanshanense							2.51
	Prunus campanulata										1.97
	Acer serrulatum							1.38
	Lindera communis										1.28
	Dendropanax dentiger							0.1			1.0
	Vaccinium wrightii var. wrightii										0.47
	Ilex hayatana										0.3
	Rhamnus nakaharae							0.2
	Eurya crenatifolia							0.09			0.06
	Machilus thunbergia										0.14
	Cleyera japonica var. japonica										0.1
	Lyonia ovalifolia var. ovalifolia							0.04
	Viburnum luzonicum							0.01
Herbaceous species
Aquatic	Schoenoplectus mucronatus ssp. robustus	83.7		88.3	4.3	21.8
	Trichophorum subcapitatum				28.0	3.2	78.0	21.3
	Rhynchospora alba					81.8	0.7
	Sparganium fallax		82.3
	Juncus effusus var. decipiens				6.0	2.3	7.0	4.3
	Persicaria sagittate			15.3
	Persicaria thunbergii			4.3				0.3
Terrestrial	Miscanthus transmorrisonensis			21.3	100	11.5	29.3	76.7
	Plagiogyria glauca								2		12
	Viola adenothrix var. adenothrix			9.0
	Elatostema trilobulatum										6
	Ophiopogon intermedius										3
	Hymenophyllum paniculiflorum								0.33		2
	Sarcopyramis napalensis var. delicata								0.33		1
	Tripterospermum lanceolatum							0.33			<0.01
	Rubus corchorifolius							0.33			<0.01
	Nertera nigricarpa								0.33

). The vegetation in the marsh, bog, swamp, and forest exhibited highly distinctive compositions. Among the herbaceous species, seven species were known wetland aquatic species. In the marsh with deeper water (vegetation type A and B), the plant community was dominated by only one wetland species (Table 1). S. fallax appeared only in deep water, but Schoenoplectus mucronatus ssp. robustus was distributed in both the deep and shallow areas (Table 2). The reduced water level could be due to the invasion of a large terrestrial grass species, Miscanthus transmorrisonensis. This species was abundant in the bog zone (vegetation types C to F) and distributed into the thicket swamp forest where the broad-leaved species Hydrangea paniculata was distributed sparsely (vegetation type G). The number of plant species increased in the swamp and forest with the reduction in water level. The thicket swamp forest (vegetation type G) presented the highest density of woody species. These were all slim and small trees with only 18.4 m² ha⁻¹ basal area. The other three forest vegetation types were all dominated by yellow cypress. However, the conifer swamp forest (vegetation type H) had the lowest number and density of the woody species. The yellow cypress in this vegetation type was obviously larger and older than that in the other two forests (much higher basal area of individual yellow cypress trees). Although the sizes of yellow cypress trees differed among the three forest vegetation types (H to J), the total basal area of yellow cypress among them were similar. Apparent higher total basal area value of all tree species in vegetation type H was contributed to by Rhododendron formosanum (36.3 m² ha⁻¹) because this species thrives exclusively in conifer swamp forests. The yellow cypress trees in the watershed, both in the swamp and upland forest, exhibited a lower density but higher individual basal area than the trees outside of the watershed in vegetation type J. This indicates that the yellow cypress trees were larger and older in the YYL watershed. In the yellow cypress forest, the coverage of ground herbaceous species was reduced. The understory of the yellow cypress forest was dominated by Plagiogyria glauca. Although no herbaceous species survey data was reported by Liao et al. [35], the dominance of Plagiogyria glauca is mentioned by Chou et al. [32] with four other species, Coptis quinquefolia, Ardisia japonica, Asarum spp., and Pellionia trilobulata in the yellow cypress forest community.

3.2. Phylogenetic Characterization of Microbiota

Samples representing soil from vegetation type F, G, H, I, and J were processed for 16S rDNA-based prokaryotic community composition analysis (Table 1). Sequencing of these samples yielded 1,767,440 high-quality sequences, among which 1,746,324 (99%) were bacterial. We removed archaeal sequences before further analysis because their abundance was low and may not be representative of the archaeal community. We defined 3787 bacterial OTUs using 97% similarity as the threshold. On average, 72,763 (49,276–89,205) sequences and 1567 (1189–2079) OTUs were determined per sample. Only 688 (18.2%) of the OTUs could be properly assigned to known genera. The taxonomic assignment for each OTU was based on SILVA 138 taxonomy.

Alpha diversity figures, as indicated by OTU richness, the Shannon index, and the Chao 1 index, are presented in Figure 2A. Alpha diversity differed significantly among the sampling plots (ANOVA: richness, p = 0.008; Shannon index, p = 0.010; Chao 1 index, p = 0.000). The vegetation type significantly affected alpha diversity (ANOVA: richness, p = 0.000; Shannon index, p = 0.009; Chao 1 index, p = 0.000).

Figure 2B depicts the phylum-level composition and indicates the major bacterial groups constituting these communities. Across these samples, bacterial communities were composed of similar bacterial phyla. These communities were dominated by Acidobacteriota (mean relative abundance 42.6 ± 8.5%) and Proteobacteria (26.0 ± 4.8%, with 17.9% as Alphaproteobacteria and 8% as Gammaproteobacteria). Dominance of Acidobacteriota over Proteobacteria indicates that these plots were relatively nutrient-poor habitats. Other major bacterial phyla detected included Verrucomicrobiota (6.2 ± 1.8%), Planctomycetota (3.9 ± 1.1%), Actinobacteriota (3.0 ± 2.1%), Bacteroidota (2.7 ± 2.7%), Desulfobacterota (2.1 ± 3.4%), Myxococcota (2.0 ± 0.9%), and Nitrospirota (2.0 ± 2.6%).

3.3. Conservation of Taxa within the YYL Watershed

We constructed a rank-abundance plot to elucidate the consistency in the phylogenetic composition exhibited in each vegetation type with the watershed mean (Figure 3). We ranked the OTU, genera, and families with their respective mean abundances within the watershed (averages among the F, G, H, and I samples), and plotted the cumulative abundance of these taxa in the F, G, H, and I samples against this rank. The presence/absence and the cumulative relative abundance of bacterial genera was almost identical among the F, G, H, and I plots for abundance taxa. For example, almost all of the top 200 most-abundant genera can be found in the F, G, and I plots, and these genera constituted practically the whole communities and accounted for over 99% of sequences in all plots. This result indicates the presence of common taxa in every plot along the gradient, which already constituted a major proportion of the respective communities.

3.4. Microbial Community Types

We performed beta diversity analysis on these bacterial communities using PCoA, based on the OTU compositions of each sample (Figure 4A). Samples collected from the same vegetation type were located closer to each other on the PCoA plot, suggesting a major influence of vegetation types on soil bacterial composition.

PCoA ordination and statistical analysis clearly reveal that these samples form four clusters (PERMANOVA, R-squared: 0.64004; p < 0.001). Cluster 1 included bog (vegetation type F, three out of three) and wetter swamp (G, thicket swamp forest, two out of three) samples; Cluster 2 included dryer swamp (H, conifer swamp forest, seven out of nine). Clusters 3 and 4 contained most mature forest (I and Jb) samples. The bacterial community composition to some extent correlated with the vegetation types. The alpha diversity among clusters differed significantly (ANOVA: richness, p = 0.000; Shannon index, p = 0.000; Chao 1 index, p = 0.000). A heatmap of the genus-level taxa turnover is presented in Figure 4B. Samples in each cluster exhibited distinct phylogenetic composition. We considered it more reasonable to conduct further community analysis on the basis of PCoA clustering because this is based on the similarity in bacterial composition.

3.5. Microbiota along the Bog-to-Forest Gradient

Our samples allowed us to compare the bacterial community composition along the bog-to-forest gradient. The four PCoA clusters roughly represent samples along this gradient. Clusters 1 to 4 comprised samples from the bog, swamp, forest, and forest outside of the watershed, respectively (Figure 4). LEfSe analyses were performed to identify differential class-level taxa among the vegetation types and among the four clusters (Figure 5A). Major soil phyla including Proteobacteria, Nitrospirota, Bacteroidota, and Verrucomicrobiota are good indicators in both vegetation type-based and PCoA-based groups. By contrast Acidobacteriota and Actinobacteriota were selected by LEfSe analysis only in the PCoA-based grouping.

Evidence of habitat preference is evident from both vegetation type-based and PCoA cluster-based grouping (Figure 5B). Relative abundance of Proteobacteria (especially Gammaproteobacteria) and Myxococcota tend to be higher in wetter habitats, and Bacteroidota and Desulfobacterota tend to be enriched in drier habitats. Verrucomicrobiota, Planctomycetota, and Actinobacteria seem to prefer swamps.

The bog-to-forest gradient exhibits a gradual decrease in diversity and F (bog) samples exhibited the highest diversity (Figure 2). Samples obtained from Ja and Jb were collected from the forest outside of the watershed. The I, Ja, and Jb samples have distinct bacterial community composition, although they were all collected from forests (Figure 2 and Figure 4). Ja and Jb samples did not form a separate cluster, indicating that geography was not a major factor in soil bacterial community assembly.

3.6. Methanotrophic Community

Methanotrophs capture methane and use it as both an energy and carbon source, and therefore impact the net methane release to the atmosphere. All major methanotrophic bacterial groups, including Type I (Gammaproteobacteria) and Type II (Alphaproteobacteria) methanotrophs, Verrucomicrobial methanotroph, and bacteria conducting anaerobic methane oxidation (AOM-bacteria, Methylomirabilota), were detected in our samples. No information on the archaeal methanotrophic group is described here because we did not obtain enough sequences for reliable analysis. The mean relative abundance for AOM, alphaproteobacterial, gammaproteobacterial, and verrucomicrobial methanotrophs are 0.10%, 0.52%, 0.11%, and 0.07%, respectively, of bacterial sequences. Methane production requires an anaerobic environment. Therefore, we expect methanotroph abundance to be higher in bog samples and lower in forest samples. We did not observe this common trend in either the vegetation type-based nor PCoA-based grouping (Figure 5B). Each taxon seems to have a different habitat preference, and therefore niche differentiation is more likely.

3.7. Surface CO₂ and CH₄ Fluxes

Along the lake–forest transect of the YYL peatland, the soil CO₂ fluxes exhibited an increasing trend as the vegetation types transition gradually from peatland-species-dominated plant communities to upland forests (Figure 6A). Two-way ANOVA revealed that both vegetation types and seasons exerted significant control over the CO₂ fluxes (both factors; p < 0.01). The Tukey honestly significant difference (Tukey-HSD) test further reveals that the vegetation types in the swamp and marsh emitted less CO₂ than both forest soils. The swamp forest (vegetation type H) exhibited the highest average CO₂ emission rate especially in the wintertime, which might lead to the observed significantly higher average CO₂ flux among the four seasons (Figure 6A).

While the CO₂ emission through aerobic respiration increases with the distance to the lake surface, the soil along the transect was concurrently changing from CH₄ sources in the peatland to CH₄ sinks in the upland forests (Figure 6B). The two factors vegetation type and season significantly influenced the CH₄ fluxes as revealed by two-way ANOVA (p < 0.01). However, the CH₄ fluxes varied substantially in the swamp forest as revealed by the springtime measurements. Among the six vegetation types, the submerged type A exhibited the highest average CH₄ emission rate and the two forests (Type I and J) had the lowest average values (Tukey-HSD, p < 0.05).

The soil water content was obviously the most varying physical environment for the vegetation types along the transect. Over the measurement period, a significant difference in average soil volumetric water content by vegetation types was determined by a one-way ANOVA test (F = 171.2, p < 0.001). The Tukey-HSD test revealed significant differences in the soil water content, and the six vegetation types were separated into four different groups of soil water contents (legend of Figure 6). The yellow cypress forest outside of the YYL watershed exhibited the lowest soil water content among all vegetation types.

The CO₂ flux as a function of the soil water content was expressed using an exponential decay function, y = a × exp(b × x). Interestingly, the results of the equation fitted with all data points (y = 1.082 × exp(−1.938 × x), R² = 0.242) were different when the yellow cypress forest outside of the watershed was excluded (y = 71.801 × exp(−9.313 × x), R² = 0.619) (Figure 7A). The relation between the CH₄ flux and soil water content revealed an opposite pattern, namely an exponential growth equation, y = 0.226 × exp(5.062 × x) (R² = 0.163) (Figure 7B).

4. Discussion

Categorizing ecosystems by vegetation types is simple. In this study, we demonstrated that a vegetation-type-based grouping of samples is a good proxy for microbiota-based PCoA clustering (Figure 5). Plant–microbe interaction is known to be a major driver of soil microbial composition in the rhizosphere [47,48]. This effect may extend to bulk soil and contribute to the changes we observed.

However, disagreement exists regarding the relationship between vegetation type and PCoA clustering in our study site (Figure 4A). Physiochemical factors that shape the vegetation type comprise its microclimate and soil properties. Microbes, because of their substantially smaller sizes, are more likely to be affected by environmental conditions at a much smaller spatial scale and are therefore prone to being influenced by habitat heterogeneity. In this study, although the three sampling points of the same plot were very close to each other (within 1 m), they did exhibit an apparent variation in bacterial composition (Figure 2). This trend was more apparent in the mature forest (I, Jb) soil samples. This indicates the crucial influence of the local environment on bacterial composition. We conclude that although vegetation type can be used as a predictor of bacterial composition and functions, the influence of other factors representing habitat heterogeneity should be considered, too.

Soil bacterial communities in the YYL region have previously been studied using a clone library-based approach [27,28], phospholipid fatty acid and denaturing gradient gel electrophoresis analysis [29]. The lakeshore communities described in Lin et al. [28] had higher Proteobacteria (59.2% vs. 25.3% in our H samples, Figure 2) and Actinobacteriota (15.1% vs. 4%), but considerably fewer Acidobacteriota (11.2% vs. 48.4%). Their native Chamaecyparis forest samples had much more Proteobacteria (55% vs. 21.4% in our I samples) and Actinobacteriota (Actinobacteria) (15.5% vs. 1.2%), but fewer Acidobacteriota (19.1% vs. 43%) and Verrucomicrobiota (Verrucomicrobia) (1.5% vs. 7.7%). Their secondary Chamaecyparis forest samples had similar Proteobacteria (31.4% vs. 29.8% in our J samples), but significantly less Acidobacteriota (Acidobacteria) (35.8% vs. 46.4%). They reported that Actinobacteria-to-bacteria ratio did not significantly differ between native and secondary Chamaecyparis forests based on the phospholipid fatty acid analysis result [29]. By contrast, we determined a considerably higher Actinobacteriota abundance in H samples (Figure 5). Their data consistently revealed more Proteobacteria and fewer Acidobacteria. This discrepancy likely results from differences in methodology such as the choice of PCR primers or the database used for identification.

Along the bog-to-forest gradient, we observed more Proteobacteria in wetter plots and more Acidobacteriota in drier plots (Figure 5B). Kitson and Bell [49] reviewed wetland microbial changes during drainage and re-wetting, and determined that drainage of fens led to an increase in Acidobacteriota (Acidobacteria) and decrease in Proteobacteria and Bacteroidota (Bacteroidetes). Compared to forest samples in the YYL region, higher microbial respiration and biomass have been reported in lakeshores and marshes [50]. Based on their laboratory incubation experiments, Imberger and Chiu [50] demonstrated that the lakeshore (swamp forest, vegetation type H in our study) soil CO₂ flux was higher than those in upland original or secondary forests. Our in situ soil CO₂ flux measurement confirms their estimations. The soil CO₂ flux in swamp forest was significantly higher than that in the yellow cypress forest outside of the watershed (Figure 6A). In addition, we determined that the microbial diversity was higher in the samples from wetter plots. The alpha diversity of the microbial community was highest in the bog (vegetation type F) and decreased gradually along the bog-to-forest gradient (Figure 2A). However, the plant alpha diversity in the bog was low, with fewer plant species present and a higher coverage of single species (Table 1 and Table 2). More plant species were distributed in both the thicket swamp forest and the yellow cypress forest outside of the watershed. It seems that the alpha diversity of the microbial community was not correlated to the plant species diversity. Nevertheless, plant roots can provide some oxygen and create oxic microhabitats in more anoxic sediment [51], providing more diverse microhabitats and thereby increasing microbial diversity. We propose that higher nutrient availability and habitat diversity in wetland contribute to this.

Although we determined the distinct characteristic bacterial composition along the bog-to-forest gradient, these communities were in fact relatively similar (Figure 3). As shown in the rank abundance plot, lines of F, G, H, and I are very close to each other, apart from J, suggesting these abundant taxa are commonly distributed within the watershed. This conservation of bacterial taxa is evident even at the OTU level. Because the majority of bacterial taxa are highly conserved (similar among plots), we propose that these distinct types of bacterial communities are likely interchangeable when changes in environmental conditions occur. Our observation that sometimes one of the three samples collected at the same plot exhibited a substantially different bacterial composition (G, Ha, Hb, I, and Jb in Figure 2B) corroborates this proposition.

The microbial composition within a peatland is affected by the vegetation and water level. Plant functional groups and water level considerably affected peatland bacterial and fungal community structures in a mesocosm study [52]. Similarly, manipulation of water levels has a strong effect on bacterial and fungal communities in a boreal peatland [53]. Our samples were obtained along the marsh-to-forest gradient that exhibited changes in both vegetation and soil water content. Both CO₂ and CH₄ productions are affected by soil water content (Figure 7), which in turn can affect microbial composition and activities.

All major bacterial methanotrophic taxa, including alphaproteobacterial, gammaproteobacterial, and verrucomicrobial taxa, were detected in our samples (Figure 5B). We expected that their relative abundance would vary depending on methane availability, which is higher in F and lower in I (Figure 6). However, we did not observe a corresponding trend in the relative abundance of total methanotroph or each specific taxon (Figure 5B). Ivanova et al. [54] indicated that the release of methanol from pectin degradation in acidic peatlands may lead to an increase in proteobacterial methanotroph because of methanol consumption. We observed an increase in two methylotrophic taxa, Methyloligellaceae (Alphaproteobacteria) and Methyloversatilis (Gammaproteobacteria) in the bog and swamp samples (Figure 5B), suggesting the possible involvement of methylated compounds. This may be one reason why we did not observe a correlation between methane production and the relative abundance of methanotroph. On the basis of the in situ CH₄ flux measurements, we did not determine significant differences among the vegetation types from the bog to the thicket swamp forest to the conifer swamp forest in the peatland area (Figure 6B). Another possibility is the influence of local non-methanotrophs, which have been reported to have a greater impact than methane availability on methanotrophic communities in various habitats [55]. In addition, we cannot exclude the potential influence of archaeal methanotrophs, for which we do not have sufficient data in this study.

The YYL region is extremely wet, with high annual precipitation and frequent fog. This results in higher soil water content and low pH, which was 3.5–4.6 in the forest soil [50,56] and 3.4–3.5 in the lakeshore [27,50]. We hypothesize that a low pH greatly influences the composition of bacterial communities. Acidobacteriota are considered oligotrophs and are well-adapted to low-nutrient environments [57]. Therefore, they are likely outcompeted by copiotrophic taxa such as Proteobacteria. Although we did observe an increase in Acidobacteriota in the G and H samples, no apparent increase in Proteobacteria was observed in the marsh (Figure 5B). The low pH in YYL region likely restricts the expansion of Proteobacteria and enables the more acidotolerant or even acidophilic Acidobacteriota to remain dominant. Kalam et al. indicated that in grasslands, as well as in forest or woodlands, Subdivision 4 (now class Blastocatellia) are the dominant Acidobacteriota [58]. However, in our samples, the major Acidobacteriota taxa were unclassified Acidobacteriales (mean abundance 14.17%) and Subgroup 2 (10.98%), and Subdivision 4 constituted only 0.005% of the community. Foesel et al. [59] indicated that the abundance of Subdivision 4 tends to increase with increasing soil pH. The low pH in our sampling sites may explain the low abundance of Subdivision 4. Planctomycetota are likely active degraders of complex polysaccharides, and their sequences comprise 5% to 22% of the total 16S rRNA gene sequences from peat samples in boreal and subarctic wetlands (reviewed in [60]). In our samples, they usually constitute no more than 5% of total reads. These bacteria are known to be more abundant in mildly acidic (pH 5.3–6.0) habitats and less abundant in more acidic (pH 3.7–4.5) bogs [60]. Low Proteobacteria and low Planctomycetota communities were also observed in a mixed peat swamp forest in Malaysia, which was also characterized by a fairly low pH of 3.9 [14].

Potential predatory bacteria were detected in our samples. Sequences of phylum Bdellovibrionia were enriched in bog samples (0.48% in F and 0.27% in others). A representative genus of this taxon, Bdellovibrio, is known for its ability to prey on prokaryotic cells (reviewed in [61]). Sequences of Vampirivibrionia (in the phylum Cyanobacteria) were also enriched in the bog samples (0.73% in F and 0.24% in others). The only cultured species of this taxon, Vampirovibrio chlorellavorus, is a predator of single-celled green algae Chlorella [62]. This suggests that a predatory lifestyle may be favored in aquatic habitats compared to terrestrial ones.

On the basis of an analysis of a series of aerial photographs spanning from 1980 to 2017, we determined that the areas of YYL and the surrounding peatland have changed over time. The high shoreline-to-area ratio of the morphometric nature of YYL together with the high watershed-to-lake area ratio, high average slope, and finally the high precipitation, especially the precipitation brought by typhoons, cause the high rate of accumulation of terrestrial material up to 15 t yr⁻¹ in the lake and on the lakeshore [63]. As Figure 8 indicates, the lake surface area has decreased from 40.6% in 1980 to 32.1% in 2017. This was determined on the basis of a reference area comprising the lake, peatland, and a surrounding 10 m buffering forest area in 1980. The lost lake area became peatland and forest, with a respective average annual increase rate of 0.14% yr⁻¹ and 0.084% yr⁻¹.

The observed terrestrialization of YYL may have accompanied a changing local water regime and resulted in the gradual change of greenhouse gas budgets. Although landscape-scale annual CH₄ and CO₂ fluxes were not exclusively estimated in this study, a gradual reduction in CH₄ emission from the watershed can be inferred on the basis of our flux measurements along the marsh–bog–swamp–forest transect. The terrestrialization of the lake (e.g., to forest) may on the one hand decrease the CH₄ emission and on the other hand increase the respiratory CO₂ emission. The net change in global warming potential depends on the magnitude of these two gas fluxes. The anthropogenic drainage of peatland provides an analogy of lake terrestrialization in view of a sinking water table and CO₂ versus CH₄ flux tradeoff. Huang et al. [64] estimated that the current degraded (i.e., drier) peatland area under the representative concentration pathway (RCP) 8.5 scenario of increasing drought stress, may emit 0.86 GtCO₂-eq yr⁻¹ more greenhouse gases at the end of the 21st century. How the landscape-scale greenhouse gas budget will change under the threatening global climate change, especially the change of hydrological regime by the lifting of the cloud forest zone, remains an important topic for further research.

5. Conclusions

In this study, we characterized the bacterial composition and the greenhouse gases (CO₂ and CH₄) production dynamics in soils along a peatland-to-forest gradient in a subtropical montane forested landscape. Because the carbon budgets of peatland may have substantial impacts on the global climate, considerable effort has been invested in the understanding the greenhouse gas budget of northern and tropical peatlands and the underlying microbial communities. To our knowledge, this study is one among recent pioneering works on subtropical montane peatland. The main conclusions of this study are presented in the following:

Bacterial communities changed along the bog-to-forest gradient, with more Acidobacteriota in drier samples and more Proteobacteria in wetter samples. The vegetation type commonly used for categorizing ecosystems can be used to predict bacterial communities. We demonstrated that genus-level bacterial taxa were highly conserved among the plots. The high humidity in YYL may help to sustain the more aquatic bacterial species in a terrestrial environment. This indicates that soil microbial communities are likely to be more susceptible to environmental changes because they already have the needed microbes to switch, for example, from a swamp-type to forest-type microbiota. Further studies are needed to determine the resilience of these microbial communities against disturbance and how general this phenomenon is in other ecosystems.

The increasing and decreasing trend in CO₂ and CH₄ emissions, respectively, along the peatland-to-upland-forest gradient were to be expected. However, a lifting cloud-forming belt over subtropical mountainous landscape caused by a warming climate could change the hydrologic regime of the current cloud forest and the embedded lakes and peatlands. The overall effect might mimic the anthropogenic draining of northern peatlands. Therefore, the feedback to the global climate system deserves closer investigation.