Winter warming in Alaska accelerates lignin decomposition contributed by Proteobacteria

Background In a warmer world, microbial decomposition of previously frozen organic carbon (C) is one of the most likely positive climate feedbacks of permafrost regions to the atmosphere. However, mechanistic understanding of microbial mediation on chemically recalcitrant C instability is limited; thus, it is crucial to identify and evaluate active decomposers of chemically recalcitrant C, which is essential for predicting C-cycle feedbacks and their relative strength of influence on climate change. Using stable isotope probing of the active layer of Arctic tundra soils after depleting soil labile C through a 975-day laboratory incubation, the identity of microbial decomposers of lignin and, their responses to warming were revealed. Results The β-Proteobacteria genus Burkholderia accounted for 95.1% of total abundance of potential lignin decomposers. Consistently, Burkholderia isolated from our tundra soils could grow with lignin as the sole C source. A 2.2 °C increase of warming considerably increased total abundance and functional capacities of all potential lignin decomposers. In addition to Burkholderia, α-Proteobacteria capable of lignin decomposition (e.g. Bradyrhizobium and Methylobacterium genera) were stimulated by warming by 82-fold. Those community changes collectively doubled the priming effect, i.e., decomposition of existing C after fresh C input to soil. Consequently, warming aggravates soil C instability, as verified by microbially enabled climate-C modeling. Conclusions Our findings are alarming, which demonstrate that accelerated C decomposition under warming conditions will make tundra soils a larger biospheric C source than anticipated. Video Abstract


Background
Nearly half of global soil organic C is stored in the northern permafrost regions (1330-1580 Pg organic C) [1,2]. With rapid increase of temperature occurring in higher latitudes, this large C pool becomes vulnerable to microbial decomposition [3]. It was shown that an increase of 2°C would accelerate the decomposition of chemically recalcitrant C by 21%, compared with only a 10% rise for chemically labile C, suggesting that chemically recalcitrant C storage is more vulnerable to global warming [4]. Therefore, the chemically recalcitrant C pool in tundra soil may act as a source for further accumulation of atmospheric greenhouse gasses because warming significantly stimulates microbially driven degradation of vulnerable C in upland tundra ecosystem, forming a positive feedback to global climate change [3,5].
As a complex aromatic heteropolymer in plant cell walls and litters, lignin is an important component of chemically recalcitrant C [6]. In addition to fungi [7], bacteria, such as Streptomyces viridosporus T7A, Pseudomonas putida MT-2, Nocardia, and Rhodococcus jostii RHA1 [8][9][10][11], have recently been demonstrated to be potent lignin decomposers. Given that the fungal-tobacterial ratio was declined by warming [12] and diverse terminal electron acceptors can be utilized by bacteria during recalcitrant C decomposition, bacteria play a key role in affecting soil C stability in warmed soil [13,14]. However, the identity of bacterial lignin decomposers in tundra soils and their responses to warming remain elusive, which prevent accurate prediction of future C fate in tundra regions.
Warmer climate leads to permafrost soil thaw in tundra regions, which releases previously frozen organic C so as to be accessible for microbial decomposition [3]. In moist acidic tundra soils, microbial community subjected to warming treatment showed a higher C-decomposing capacity [15]. Abundances of chemically recalcitrant Cdecomposing genes were also increased by warming in tundra soils, which corroborated with higher ecosystem respiration [1]. Moreover, warming in tundra regions promotes plant root exudates such as organic acids, sugars, and amino acids [16], which could accelerate the chemically recalcitrant C vulnerability (termed the priming effect) [17]. Given that approximately 65% of total C in tundra soil is stored as lignin, chitin, and terpenes [1,2], it was hypothesized that warming would aggravate tundra C instability by shifting microbial community composition, increasing decomposer abundance and the decomposing capacity of chemically recalcitrant C.

Strong warming effects on lignin decomposers
To identify potential lignin decomposers and their responses to warming, the tundra soils, which were subjected to in situ winter warming for a 1.5-year period, in parallel with unwarmed/control soils, were collected. Since the turnover time of slow C pool of tundra soil is6 00 days [18], the soils were starved by laboratory incubation for 975 days to completely deplete chemically labile C reserves, a practice also recommended for helping accurate modeling of soil C kinetics [19]. Subsequently, stable isotope probing (SIP) experiments were performed to label active bacterial decomposers in warmed and control samples. Intuitively, lignin should be used. However, it is difficult to label all of C atoms of lignin owing to its complex structure. More importantly, lignin is resistant to biochemical breakdown [20]; thus, incubation time for the SIP experiment has to be extended, which inevitably causes cross-feeding. Therefore, vanillin, an intermediate of lignin decomposition widely used as a model aromatic substance to detect lignin depolymerization [20,21], was used as the SIP substrate to identify potential ligninolytic microorganisms.
The peak of 13 C-labeled DNA was detected around the density of 1.748 g/ml (heavy fractions, Supplementary Fig. 1), which was present only in samples incubated with 13 C-vanillin but not in samples incubated with 12 Cvanillin or no vanillin. The peak of 12 C-labeled DNA was detected around the density of 1.720 g/ml (light fractions). Total abundance of 13 C-labeled DNA was increased from (2.7 ± 0.5) × 10 5 copies/g soil in control samples to (1.0 ± 0.3) × 10 6 copies/g soil in warmed samples (P < 0.01), suggesting that potential lignin decomposers were stimulated (Fig. 1a). High-throughput sequencing of 13 C-labeled DNA showed that there were 63 operational taxonomic units (OTUs) in warmed samples, which were significantly (P < 0.05) higher than that of control samples (28 OTUs). In addition, warming enhanced total abundance of 12 C-labeled DNA from (4.7 ± 1.6) × 10 5 to (1.4 ± 0.4) × 10 6 copies/g soil (P < 0.05), suggesting that growth of the overall bacterial communities was also stimulated.
Two clusters comprised of the majority of lignin decomposers were identified in a phylogenetic tree constructed from potential decomposers (Fig. 1b). Cluster I contained 14 OTUs, of which 5 OTUs were shared by both warmed and control samples. Strikingly, cluster I was solely composed of β-Proteobacteria, with 4 OTUs belonging to the genus Burkholderia. Although warming slightly decreased relative abundance of Burkholderia in 13 C-labeled DNA from 99.4 to 90.8% (P < 0.01, Supplementary Fig. 2a), it increased absolute abundance of 13 C-labeled Burkholderia from 2.6 × 10 5 to 9.1 × 10 5 copies/g soil (Fig. 1c), reflecting strong stimulation. Burkholderiales dominated lignocellulosic decomposers in coniferous forest soils across North America [22]. It accounted for 64% of bacterial clone sequences in Canadian High Arctic soil [23] and was identified as keystone species in Arctic heathland soils using a network analysis [24]. Burkholderiales are also Fig. 1 The absolute abundances of various taxa, measured by 16S rRNA genes, after incubation with 13 C-vanillin, determined by quantitative polymerase chain reaction (qPCR) and amplicon sequencing. a The numbers of total 12 C-and 13 C-labeled 16S rRNA gene copies in 13 C-vanillin incubated samples (n = 3, biological replicates of warmed or control samples). b The circular maximum likelihood phylogenetic tree of the 13 Clabeled active decomposer operational taxonomic units (OTUs). Relative abundance is the modified sequence number in heavy fractions (MSNH; see the "Methods" section for details) of the OTU (n = 3, biological replicates of warmed or control samples), c The number of 13 C-labeled 16S rRNA gene copies of Burkholderia (n = 3, biological replicates of warmed or control samples). d The number of 13 C-labeled 16S rRNA gene copies of α-Proteobacteria (n = 3, biological replicates of warmed or control samples). *0.01 < P ≤ 0.05 and **0.001 < P ≤ 0.01, as determined by a twotailed t test. W, warmed soils; C, unwarmed/control soils. Data are shown as mean ± standard error abundant in bog ecosystems in the northern hemisphere [25,26]. Further, relative abundance of Burkholderia significantly increased by 31.5-folds after the 975-day incubation (P < 0.05) and then increased by another 25.4-folds after addition of vanillin (P < 0.001), demonstrating that Burkholderia is highly responsive to C substrate availability ( Supplementary Fig. 3). Clearly, Burkholderia is a keystone taxon in Arctic tundra soils ( Supplementary Fig. 4).
To directly prove Burkholderia as lignin decomposers, two Burkholderia strains were isolated from tundra soils. Both strains were able to grow with lignin as the sole C source ( Supplementary Fig. 5). Using whole-genome sequencing followed by genome annotation with Genome Taxonomy Database Toolkit (GTDB-tk) [27], the average nucleotide identity (ANI) between both isolate strain genomes was 99.34%, which belonged to the species of Burkholderia zhejiangensis [28]. A number of peroxidases were annotated in the genomes, including multiple copies of a gene encoding catalase-peroxidase well known in lignin decomposition [29]. Seven aromatic acid transporter genes and β-ketoadipate pathway genes associated with lignin decomposition were identified, which constitute a lignin metabolism pathway ( Supplementary Fig. 6).
Cluster II contained 35 OTUs. Only two of these OTUs were shared between warmed and control samples, with 26 OTUs exclusive to warmed samples. Thirty-two cluster II OTUs belonged to α-Proteobacteria, including genera Bradyrhizobium, Phenylobacterium, Sphingomonas, Novosphingobium, and Methylobacterium. To provide further evidence for lignin-decomposing capacity, 7 α-Proteobacterial genomes, including 3 Bradyrhizobium, 3 Phenylobacterium, and 1 Sphingomonas, were assembled from deep metagenomic sequencing data of starved soils despite their low abundance compared with cluster I OTUs (Supplementary Table 1). Genes encoding catalaseperoxidases are present across all the genomes. Several genes associated with the β-ketoadipate pathway are also present, which are responsible for metabolizing intermediates during lignin decomposition.
Warming considerably enhanced total abundance of α-Proteobacteria in 13 C-labeled DNA from 1.1 × 10 3 to 8.9 × 10 4 copies/g soil (Fig. 1d), and its relative abundance by 3.1folds ( Supplementary Fig. 2b). Among α-Proteobacteria, increases in representatives of genera Bradyrhizobium (from undetected to 4.0 × 10 4 copies/g soil) and Methylobacterium (from undetected to 3.2 × 10 4 copies/g soil) were most notable ( Supplementary Fig. 7). Warming did not significantly change relative abundance of active α-Proteobacteria before and after the 975-day incubation at 25°C (P > 0.05), but significantly stimulated it by 3.1-fold after incubation with vanillin at 25°C (P < 0.05, Supplementary Fig. 8), suggesting that α-Proteobacterial decomposers were considerably induced. Typically, α-Proteobacteria prefer nutrient-rich environments and exhibit fast growth rates [30]. Tundra soil thawing by warming exposes previously frozen organic C to microbial decomposition [31], which stimulates α-Proteobacteria. Similarly, warming has increased the abundance of α-Proteobacteria in Antarctic environments [30]. However, bearing in mind a caveat that fast-growing bacteria may be preferably stimulated in the SIP experiment, slow-growing bacteria may not incorporate sufficient 13 C-label within 6 days of incubation.
Members of Acidobacteria and Actinobacteria phyla were also identified as active lignin decomposers in coniferous forest soils across North America [22]. Consistently, two OTUs of Acidobacteria, which were detected only in warmed samples, and eight OTUs of Actinobacteria were identified as potential lignin decomposers (Fig. 1b).

The priming effect doubled by warming
Changes in α-Proteobacteria abundance were positively correlated with soil CO 2 production in the High Arctic [30]. In this study, higher abundances of α-Proteobacteria in the warmed samples ( Fig. 1d) also corroborated with significantly higher total CO 2 production measured in the SIP experiment (182.6 ± 6.2 μmol in warmed versus 174.9 ± 4.0 μmol in control samples) (Fig. 2). The 13 C content of 13 CO 2 production after a 6-day incubation (92.7 ± 1.1 μmol in warmed and 90.2 ± 0.5 μmol in control samples) was close to that of added vanillin (82.6 μmol 13 C in warmed or control samples), indicating that added vanillin was depleted by the end of the incubation period. A strong priming effect caused by supplementing with vanillin was detected, as the C molar of total CO 2 production was substantially higher than the amount calculated from theoretical oxidation of vanillin (110.1 μmol). The C primed by vanillin under warming (19.1 ± 2.2 μmol, the inset of Fig. 2) was significantly more than that under the control condition (9.9 ± 0.5 μmol). This result contrasts a previous study that detected no significant priming effect in the organic layer of tundra soil [32]. This is likely attributed to our prolonged 975-day incubation of soils to deplete the chemically labile C reserves [31] prior to supplementing with vanillin.
To examine whether higher CO 2 production in warmed samples arises from changes in functional genes associated with C decomposition, especially of aromatics and lignin, relative abundances of related functional genes in 13 C-labeled DNA were quantified by GeoChip 5.0. Strikingly, almost all detected lignin-decomposing genes, including the mnp gene encoding peroxidase, lcc genes encoding phenol oxidase, glx gene encoding glyoxal oxidase, vanA gene encoding vanillate demethylase A, and vdh gene encoding vanillin dehydrogenase, increased by 10.1%-26.3% under warming (Fig. 3a). This is consistent with higher 13 C ratio and greater total CO 2 flux under warming (Fig. 2). More broadly, 32 out of 43 aromatic-decomposing genes also significantly increased Fig. 2 The accumulated CO 2 flux and the priming effect during the 6-day incubation trial. Red lines/symbols represent in situ warmed soils (W) and blue lines/symbols represent unwarmed/control soils (C). Solid lines represent total CO 2 and dashed lines represent 13 C-CO 2 . The significance of the difference was determined by one-way ANOVA (n = 3, biological replicates of warmed or control samples) on CO 2 amounts and primed C amounts between warmed and control samples (shown in the inset of the figure). W, warmed samples; C, control samples; NS, not significant. Data are shown as mean ± standard error Aromatics-and lignin-decomposing genes belonging to Burkholderia in 13 C-labeled DNA. c Genes associated with other key C-decomposing pathways. Ppl, phospholipids; Ca, Camphor; Cu, Cutin; and Ta, Terpenes. Differences in relative abundances were determined using one-way ANOVA (n = 3, biological replicates of warmed or control samples). *0.01 < P ≤ 0.05, **0.001 < P ≤ 0.01, and ***P ≤ 0.001. Data are shown as mean ± standard error by 4.9-184.1% under warming (Fig. 3a). Among them, aromatic and lignin-decomposing genes derived from Burkholderia significantly increased by 11.9% (Fig. 3b), suggesting higher decomposition potentials of Burkholderia. Interestingly, 83.5% of the other C-decomposing genes irrelevant to lignin decomposition, including those associated with decomposing chitin (e.g., chitinase), terpenes (e.g., cdh encoding carveol dehydrogenase), pectin (e.g., pel encoding pectin lyase), cellulose (e.g., endoglucanase), hemicellulose (e.g., xylanase), and starch (e.g., amyA encoding α-amylase), also significantly increased under warming by 5.7-40.6% (Fig. 3c). This is likely to arise from a broad functional response of priming to warming. In contrast, no C-decomposing gene was significantly decreased in abundance by warming. Given that warming may promote plant root exudation, releasing more chemically labile C substrates to the soil [16], this strong priming effect implicates that there is an important positive feedback to global warming in tundra environments.

Model verification and data synthesis
To examine the scale of potential positive feedback of C cycling to warming, a Microbially ENabled Decomposition (MEND) model [33] was implemented to estimate the changes in C-decomposing rates. Simulated soil respiration for the 975-day laboratory incubation period agreed well with observed data (R 2 ≥ 0.85, Fig. 4a). Therefore, the average soil respiration rate was simulated over the long term. The C-decomposing rates were significantly higher under warming (Fig. 4b), which were verified by in situ respiration measurements in tundra regions [1,34]. Accordingly, tundra soil organic C is projected to decrease significantly (Fig. 4c), leading to greater net C loss under warming [1]. Higher soil respiration rates are likely caused by significantly more active microbial biomass (Fig. 4d) and higher decomposition rates of oxidative enzymes (Fig. 4e). Given that the oxidative enzymes mainly consist of ligninases [35], these modeling results are consistent with our SIP results (Fig. 1a), suggesting that warming would markedly increase the release of CO 2 from lignin or other chemically recalcitrant C sources to the atmosphere.

Conclusions
In this study, warming increased total abundance and functional capacities of all potential lignin decomposers in the Alaska tundra, resulting in a doubling of the priming effect on chemical recalcitrant C decomposition. Coupled with high-throughput sequencing and strain isolation, Burkholderia and several genera of α-Proteobacteria identified as major lignin decomposers were stimulated by warming. To our knowledge, this is the first study directly investigating active recalcitrant C-decomposing bacteria in response to warming, leading to establishment of explicit linkages between soil respiration and microbial community composition and functional capacity. The stronger priming effect under warming is alarming, which indicates a previously overlooked mechanism that accelerates climate warming in tundra regions. Given that the past 5 years are the five warmest years on record since 1880, our study provides important insights into tundra soil C stability under global warming.

Methods
Site description and soil sample preparation The warming experiment was carried out at the Carbon in Permafrost Experimental Heating Research (CiPEHR) site, which was established in 2008. As described previously [3], soils were warmed during winter months by increasing snow cover behind snow fences, which were perpendicular to the dominant south-easterly winter winds. Warmed plots were on the leeward side of snow fences, while the control plots were on the windward side. Snow fences trapped and accumulated an insulating snow layer on the warmed plots. As a result, the warmed soil was at an average temperature of 2.3°C higher than the control soil during winter. Soil samples were collected in May 2010. Intact soil core was collected from a depth of 15-25 cm in order to avoid litter and coarse root material that constituted most of the tundra soil from the 0-15 cm layer. Soil was dried in an ED-056 oven (Sigma-Aldrich, St. Louis, MO, USA) at 60°C until a constant weight was reached. Soil was briefly grinded before soil C and N contents were measured using an ECS 4010 Elemental Analyzer (Costech Analytical Technologies, Valencia, CA, USA). Since extended soil incubation of soil is important for identifying potential lignin decomposers and measuring the priming effect of recalcitrant C [19], intact soil samples were incubated in lightproof jars for 975 days to deplete chemically labile C. As described previously [31], soil was added to a perforated foil cup and placed over a bed of 3-mm glass beads inside the jar to allow drainage and maintain soil moisture. Because the temperature of surface soils at the CiPEHR site could reach 25°C during the growing season [31], jars were placed in a 25°C water bath to ensure complete depletion of chemically labile C within 975 days. Jars were covered with perforated lids to allow air exchange.
The stable isotope probing experiment and CO 2 flux measurement 13 C-vanillin (vanillin-(phenyl-13 C 6 ); 99 atom% 13 C) and 12 C-vanillin were used as stable isotope probe substrates (Sigma-Aldrich, St. Louis, MO, USA). To ensure even addition to the soil, vanillin was dissolved in water as 0.8% solution, not exceeding the solubility of 1%, and then injected evenly to the soil. Three incubation groups, i.e., (1) with 0.345 ml of 0.8% 13 C-vanillin in 2.76 g of soil (1 mg/g w/w) as isotopic treatment, (2) with 0.345 ml of 0.8% 12 C-vanillin in 2.76 g of soil as isotopic control, and (3) with 0.345 ml of water in 2.76 g of soil as the background, were set up for both warming and control samples. Each group was set up in three biological replicates. The final water content did not exceed 70% of the water-holding capacity (WHC) of the soil, and each replicate was sealed in a 25-ml lightproof bottle and incubated at 25°C for 6 days.
Five milliliters of headspace gas was collected daily into 12-ml evacuated vials (Labco Limited, Lampeter, UK), after which the bottle was opened and refreshed for 30 min. To generate positive pressure in relation to atmospheric pressure, sampled gas in vials were diluted by injecting 10 ml of N 2 gas into each vial. CO 2 and 13 CO 2 concentration were measured at the Stable Isotope Facility, University of California, Davis. The obtained CO 2 and 13 CO 2 concentrations in parts per million were transformed to obtain the number of moles (n) using the ideal gas law equation PV = nRT, in which P (pressure) was 101 kPa; V (volume) was 25 ml multiplying the obtained CO 2 or 13 CO 2 concentration in parts per million; R (the gas constant) was 8.314 J K −1 mol −1 , and T (temperature) was 298 K.
The percentage of the CO 2 -C derived from vanillin was calculated as: where δ C is the δ 13 C value of respired CO 2 from the soil with no added vanillin, δ T is the δ 13 C value of respired CO 2 from the soil with 13 C-vanillin, and δ L is the δ 13 C value of 13 C-vanillin. Because vanillin was added to the soil in the form of water solution, the amount of soil organic matter C primed by vanillin was calculated as total soil respiration after vanillin addition minus the amount of C respired from vanillin, and then minus the amount of C primed by water (C respired from the soil with no added vanillin).

Soil DNA extraction
Soil DNA was extracted using a freeze-grinding method as described previously [36] and purified by a MOBIO PowerSoil kit (MO BIO Laboratories Inc., Carlsbad, CA, USA) according to the manufacturer's protocol. DNA quality was assessed by a NanoDrop ND-1000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) based on spectrometry absorbance at wavelengths of 230 nm, 260 nm, and 280 nm. The absorbance ratios of 260/280 nm were larger than 1.8, and of 260/230 nm were around 1.7.

C-DNA separation
Density gradient ultracentrifugation of 13 C-labeled DNA was performed according to a previous protocol with minor modifications [37]. In brief, 5.1 ml of a solution was centrifuged, which was composed of 5 μg of soil DNA, 1.90 g ml −1 cesium chloride (CsCl) (MP Biomedicals, Santa Ana, CA, USA), and a gradient buffer of 1 mM EDTA, 0.1 M KCl, and 0.1 M Tris-HCl, reaching a final density of 1.725 g ml −1 . The solution was sealed in a polyallomer centrifuge tube (cat. no. 342412, Beckman Coulter, Brea, CA, USA) with a cordless tube topper, and centrifuged on a Vti 65.2 rotor of an Optima L-XP ultracentrifuge (Beckman Coulter, Brea, CA, USA) at 177,000 g and 20°C for 48 h. Solution from each centrifuged tube was then separated into twenty-four 220-ml fractions (14 drops per fraction) according to the density gradient. The buoyant density of each fraction was determined by an AR200 digital refractometer (Reichert, Depew, NY, USA). To this end, DNA in each fraction was precipitated with 20 μg of glycogen and 2 volumes of PEG solution (30% PEG 6000 and 1.6 M NaCl), washed with 70% ethanol, and re-suspended in 35 μl of ultrapure water.
Quantitative PCR (qPCR) was used to determine absolute abundances of 16S rRNA genes, thus identifying the fractions containing 13 C-DNA. Universal primers 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) were used for targeting the V4 region of the 16S rRNA genes. qPCR was performed in triplicate 20-μl reactions containing 10 μl SsoAdvanced™ Universal SYBR® Green Supermix (Bio-Rad, Hercules, CA, USA), 350 nM each primer and 1 μl of template, using a thermocycler program of 35 cycles of 95°C for 20 s, 53°C for 25 s, and 72°C for 30 s on an IQ5 Multicolor Real-time PCR Detection System (Bio-Rad, Hercules, CA, USA). Gene copy numbers were determined by a standard curve constructed with 16S rRNA gene segment of E. coli JM109 competent cells and TA cloning vector (Promega, Madison, WI, USA). According to the buoyant density and 16S rRNA gene copy number, a light DNA peak was identified as 12 C-labeled DNA, and a heavy DNA peak was identified as 13 C-labeled DNA.

16S rRNA gene amplicon sequencing
A two-step PCR was performed prior to 16S rRNA gene sequencing [38]. In the two-step PCR, the first step of the V4 region of 16S rRNA genes was amplified by the primer 515F and 806R in triplicate 25 μl reaction containing 2.5 μl of 10 × AccuPrime PCR buffer (containing dNTPs) (Invitrogen, Grand Island, NY, USA), 1 μl of 10 μM forward and reverse primer, 2 μl of template DNA and 0.2 μl of AccuPrime High-Fidelity Taq Polymerase. The thermocycler program was as follows: 94°C for 1 min, 10 cycles of 94°C for 20 s, 53°C for 25 s, and 68°C for 45 s, and a final extension at 68°C for 10 min. The second step of PCR also used a 25-μl reaction containing 2.5 μl of 10 × AccuPrime PCR buffer (including dNTPs), 1 μl of 10 μM 515F and 806R primer combined with the Illumina adaptor sequence, a pad and a linker of two bases, and a barcode sequences on the reverse primers, 15 μl of aliquot of the first step purified PCR product and 0.2 μl of AccuPrime High-Fidelity Taq Polymerase. The thermal cycling condition was the same as the first step except a cycle number of 20. PCR products from the second step were examined by agarose gel electrophoresis, and then triplicate PCR products were combined and quantified by Pico Green.
PCR products from each fraction were pooled at equal molarity and sequenced in the same MiSeq run. The pooled mixture was purified with a QIAquick Gel Extraction kit (Qiagen Sciences, Germantown, MD, USA) and re-quantified with Pico Green. The detailed protocol for MiSeq sequencing had been described previously [39].
The raw sequence reads were processed using an inhouse pipeline built on the Galaxy platform. First, the FastQC (http://www.bioinformatics.babraham.ac.uk/pro jects/fastqc/) was used to evaluate the quality of raw sequence data. Second, the spiked PhiX reads were removed by demultiplexing with E value < 10 −5 . Third, sequences were sorted to corresponding samples according to their barcodes on the primers, which allowed for 0 mismatch. Fourth, Btrim was performed for quality trimming [40], then forward and reverse reads of the same sequence with at least 20-bp overlap and < 5% mismatch were combined by FLASH v1.2.5 program [41]. Combined sequences were removed if they contained ambiguous bases or were less than 240 bp, and then Uchime was used to remove the chimeric sequences [42]. Finally, OTUs were clustered using UPARSE at the 97% similarity level [43]. Each fraction was randomly resampled according to the absolute abundance of 16S rRNA gene based on qPCR result, with a resampling size of 29,845. OTUs were annotated through Ribosomal Database Project (RDP) classifier 2.5 with minimal 50% confidence score [44].

Identification of potential lignin decomposers
The minor peak on the right side of the major peak was identified as 13 C-labeled DNA from the abundancedensity plot of 13 C-incubated samples. For each sample, four fractions contributing to the major peak were termed as the "heavy fractions", and four fractions contributing to the minor peak were termed as the "light fractions".
Not all OTUs detected in the heavy fraction can be labeled by 13 C due to interference of high GC content. In addition, OTUs are prone to sequencing errors. A strictly filtering method was developed to distinguish 13 C-labeled OTUs from false positives ascribed to high GC content or sequencing errors. This filtering method included three consecutive steps: (i) To exclude OTUs of high GC content, a non-parametric t test was performed to determine the significance of sequence number difference assigned to each OTU found in heavy fractions between 13 C-incubated samples and 12 C-incubated samples. OTUs with insignificant P values were removed (P > 0.10). (ii) The relative abundances of remaining OTUs were subtracted from the relative abundances in heavy fractions of the corresponding 12 C-incubated samples. The difference was termed as the modified sequence number in heavy fractions (MSNH). And (iii) to exclude OTUs from sequencing errors, OTUs were removed when its MSNH present in the heavy fractions was less than 20% of the total relative abundance. The remaining OTUs were regarded as potential lignin decomposers.
The BMM-defined medium [47] is comprised of 0.80 g L −1 NaCl, 1.0 g L −1 NH 4 Cl, 0. Burkholderia strains were inoculated into 50 ml of the BMM medium with 0.2% yeast and incubated at 28°C with constant shaking at 180 rpm to an OD 600 (optical density at 600 nm) of approximately 1.0. Then, 1 ml of culture was aseptically inoculated into three parallel culture flasks containing 50 ml of the BMM-defined medium. The flasks were incubated at 28°C with constant shaking at 180 rpm for 9 days. Uninoculated medium was used as a control.

Functional metagenomics analyses
Genomic annotation was performed using Automatic Genomic Analysis Pipeline (AGAP, version 1.2, an internal pipeline). Genomic annotation was conducted using PROKKA (version 1.11) [54]. First, finished genomes and draft genomes were submitted to gene calling using Prodigal (version 2.6) [54] with output of translated protein sequences, single mode and genetic code of bacteria and archaea. Then rRNA genes were predicted using Barrnap (version 0.7). Pseudogenes and coding sequences overlapping with tRNA and rRNA gene were removed by PROKKA. The 16S rRNA genes used for taxonomic classification were classified using RDP Classifier (version 2.12) [44]. Protein sequences were submitted to DIA-MOND [55] search (BLASTp) against NCBI NR database (version Jan 2016) with E value cutoff of 1e− 5, coverage cutoff of 0.5 and maximum target number of 50. The BLASTp results were imported into MEGAN6 (Ultimate Edition, version 6.6) [56] for functional profiling with output of SEED Subsystem, KEGG, and COG categories. Exported tables of functional profiles were integrated for comparison of genomes. Genome binning (assigning the non-overlapping contigs to genomes) was performed by tetra-nucleotide frequency and verified by differential coverage binning [57].

Experiments with GeoChip 5.0
Approximately 50 ng of DNA separated from heavy fractions in warming or control samples were amplified using a Templiphi kit (GE Healthcare, Little Chalfont, UK). The amplified DNA (2 μg) was labeled with fluorescent dye (Cy-3) dUTP using random primers and Klenow fragment of DNA polymerase I at 37°C for 6 h, followed by heating at 95°C for 3 min. Labeled DNA was then purified, dried in a SpeedVac at 45°C for 45 min, and re-suspended in 43.1 μl of hybridization buffer containing 27.5 μl of 2X HI-RPM hybridization buffer, 5.5 μl of 10X CGH blocking agent, 2.4 μl of cot-1DNA, 2.2 μl of universal standard, and 5.5 μl of formamide. DNA was hybridized with Geo-Chip 5.0 (60 K) in a SL incubator (Shel Lab, Cornelius, OR, USA) at 67°C for 24 h. Then, GeoChip arrays were washed and scanned by an MS 200 Microarray Scanner (Roche, Basel, Switzerland) at 532 nm and 635 nm. Raw signals from the scanning were processed by an online pipeline as previously described [58].

Statistical and phylogenetic analyses
Differences of relative abundances between isotopic treatment and control groups were determined by Wilcoxon rank sum and signed-rank tests [59]. The two-tailed t test was used to determine the difference of richness of active communities between warming and control samples. All analyses above were performed in R software (version 3.3.2). Unless otherwise stated, mean values are given ± standard error of the mean, significant differences are determined by one-way ANOVA, and values of P ≤ 0.05 were considered significant.
The maximum likelihood phylogenetic tree was constructed based on the representative sequence for each OTU, as determined by UPARSE. MEGA 6.05 [60] was used to construct the phylogenetic tree with MUSCLE alignment, maximum likelihood method and a bootstrap value of 1000. The visual of final tree was generated by iTOL [61].

Molecular ecological network analyses
Phylogenetic molecular ecological networks (pMENs) were constructed from the 16S rRNA gene sequencing data, using a random matrix theory (RMT)-based network approach [62]. To ensure reliability, only OTUs detected in at least 10 out of 12 samples were used for network construction. In brief, a matrix containing Pearson's rho correlation between any pair of OTUs was generated. The threshold for network construction was automatically determined when the nearest-neighbor spacing distribution of eigenvalues transitioned from GOE to Poisson distributions. Consequently, a threshold of 0.74 was used in both warming and control sample networks. Random networks corresponding to all pMENs were constructed using the Maslov-Sneppen procedure with the same network size and average number of links to verify the systemspecificity, sensitivity, and robustness of the empirical networks [63].

Microbially enabled decomposition modeling
A Microbially ENabled Decomposition model (MEND) [33,64] was used in this study. The MEND model is a sophisticated model with parameters representing microbial dormancy, resuscitation, and mortality. It simulates soil organic matter (SOM) decomposition in response to changes in environmental conditions such as temperature [33] by explicit microbe-mediated oxidative and hydrolytic processes. The temperature response is modeled by the Arrhenius equation characterized by the activation energy. Here, soil C pools were initialized by the measurements of soil organic C and microbial biomass C in the control and warmed samples at the beginning of the lab incubation experiments. To determine microbial parameters in soil samples, experimental data were combined at two experimental incubation temperatures (i.e., 15°C and 25°C) since differential temperature allows for more accurate estimate of model parameters. The stochastic shuffled complex evolution (SCE) algorithm was used [33] to determine model parameters by achieving the highest goodness-of-fit between modeled and measured CO 2 fluxes. The microbial traits or parameters (e.g., microbial growth, maintenance, mortality, C use efficiency, and active versus dormant fractions) represent the microbial community characteristics related to soil C mineralization. Then, the model simulated with those parameters was run for 10 years to predict the long-term warming effect on SOM decomposition.

Supplementary information
Supplementary information accompanies this paper at https://doi.org/10. 1186/s40168-020-00838-5. the genomes based on the metagenomic data from the soil the 975-day laboratory incubation period. R.T. conducted metagenomic assemblies and annotations. X.T., J.F., and F.F. carried out the 6-day incubation, and the related gas and 16S rRNA gene amplicon experiments and analyses. G.W. carried out the MEND model analysis. D.N. assisted with the statistical analyses. X.T., J.F., Y.Y., and J.Z. wrote the paper. C.T.B., LW.W., L.H., Q.G., LY.W., E.A.G.S., K.T.K., J.R.C., and J.M.T. edited the manuscript. All authors were given the opportunity to review the results and comment on the manuscript. The authors read and approved the final manuscript.