# Seasonal effects of long-term warming on ecosystem function and bacterial diversity

Melissa S Shinfuku, Luiz A Domeignoz-Horta, Mallory J Choudoir, Serita D Frey, Megan F Mitchell, Ravi Ranjan, Kristen M Deangelis

► To cite this version:

Melissa S Shinfuku, Luiz A Domeignoz-Horta, Mallory J Choudoir, Serita D Frey, Megan F Mitchell, et al.. Seasonal effects of long-term warming on ecosystem function and bacterial diversity. PLoS ONE, 2024, 19 (10), pp.e0311364. 10.1371/journal.pone.0311364. hal-05071277

HAL Id: hal-05071277

https://hal.inrae.fr/hal-05071277v1

Submitted on 16 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# Seasonal effects of long-term warming on ecosystem function and bacterial diversity

Melissa S. Shinfuku¹, Luiz A. Domeignoz-Horta¹,², Mallory J. Choudoir¹,³, Serita D. Frey⁴, Megan F. Mitchell⁵, Ravi Ranjan¹,⁶, Kristen M. DeAngelis¹*

1 Microbiology Department, University of Massachusetts, Amherst, MA, United States of America, 2 INRAE, AgroParisTech, UMR EcoSys, Université Paris-Saclay, Palaiseau, France, 3 Department of Plant and Microbial Biology, North Carolina State University, Raleigh, NC, United States of America, 4 Center for Soil Biogeochemistry and Microbial Ecology, Department of Natural Resources and the Environment, University of New Hampshire, Durham, NH, United States of America, 5 Graduate Program in Organismic and Evolutionary Biology, University of Massachusetts, Amherst, MA, United States of America, 6 Genomics Resource Laboratory, Institute for Applied Life Sciences, University of Massachusetts, Amherst, MA, United States of America

* deangelis@microbio.umass.edu

<figure><img src="image_2.png" /></figure>

## Abstract

Across biomes, soil biodiversity promotes ecosystem functions. However, whether this relationship will be maintained within ecosystems under climate change is uncertain. Here, using two long-term soil warming experiments, we investigated how warming affects the relationship between ecosystem functions and bacterial diversity across seasons, soil horizons, and warming duration. Soils were sampled from these warming experiments located at the Harvard Forest Long-Term Ecological Research (LTER) site, where soils had been heated +5°C above ambient for 13 or 28 years at the time of sampling. We assessed seven measurements representative of different ecosystem functions and nutrient pools. We also surveyed bacterial community diversity. We found that ecosystem function was significantly affected by season, with autumn samples having a higher intercept than summer samples in our model, suggesting a higher overall baseline of ecosystem function in the fall. The effect of warming on bacterial diversity was similarly affected by season, where warming in the summer was associated with decreased bacterial evenness in the organic horizon. Despite the decreased bacterial evenness in the warmed plots, we found that the relationship between ecosystem function and bacterial diversity was unaffected by warming or warming duration. Our findings highlight that season is a consistent driver of ecosystem function as well as a modulator of climate change effects on bacterial community evenness.

<figure><img src="image_3.png" /></figure>

### OPEN ACCESS

**Citation:** Shinfuku MS, Domeignoz-Horta LA, Choudoir MJ, Frey SD, Mitchell MF, Ranjan R, et al. (2024) Seasonal effects of long-term warming on ecosystem function and bacterial diversity. PLoS ONE 19(10): e0311364. https://doi.org/10.1371/journal.pone.0311364

**Editor:** Mette Vestergård, Aarhus University, DENMARK

**Received:** May 13, 2024

**Accepted:** September 9, 2024

**Published:** October 24, 2024

**Peer Review History:** PLOS recognizes the benefits of transparency in the peer review process; therefore, we enable the publication of all of the content of peer review and author responses alongside final, published articles. The editorial history of this article is available here: https://doi.org/10.1371/journal.pone.0311364

**Copyright:** © 2024 Shinfuku et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** The primer-trimmed bacterial sequencing data are available at the NCBI Sequence Read Archive under the BioProject accession number PRJNA957454. Ecosystem

## Introduction

Climate change is driving losses in biodiversity, and these losses threaten both ecosystem productivity [1] and stability [2, 3]. While there is evidence that bacterial diversity promotes ecosystem function, not enough is known about how increased temperatures associated with climate change will affect this relationship. Studies examining the effects of warming on
function data (beta-glucosidase activity, N-acetyl glucosaminidase activity, oxidative enzyme activity, microbial biomass carbon, and respiration) is available in the Harvard Forest Data Archive under the Dataset ID HF431. Total soil organic carbon and nitrogen, z-score tranformed ecosystem function measurements, and ecosystem multifunctionality indices are available in the Supporting information and in the Harvard Forest Data Archive under the Dataset ID HF449.

**Funding:** This work was supported by a National Science Foundation (NSF, https://www.nsf.gov/) Long-Term Research in Environmental Biology grant (DEB-1456528) to SDF and KMD, and an NSF CAREER award (DEB-1749206) to KMD. The soil warming experiments at Harvard Forest are maintained with support from the National Science Foundation (NSF) Long Term Ecological Research Program (DEB-1832110). This work was also supported in part by a Graduate School Dissertation Research Grant from the University of Massachusetts Amherst to MSS (https://www.umass.edu/graduate). This research was also supported [in part] by the intramural research program of the U.S. Department of Agriculture (https://www.usda.gov), National Institute of Food and Agriculture, Hatch Multistate, accession number 7004345. The findings and conclusions in this publication have not been formally disseminated by the U. S. Department of Agriculture and should not be construed to represent any agency determination or policy. None of the funders had any role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. There was no additional external funding received for this study.

**Competing interests:** The authors have declared that no competing interests exist.

microbial diversity and ecosystem function have revealed decreases in both microbial diversity and ecosystem function [4], a weaker relationship between microbial diversity and ecosystem function [5], or even abiotic modulators of the microbial diversity-ecosystem function relationship [6]. Increased temperatures also alter interactions between taxa in microbial communities [7, 8], which might explain why at higher temperatures decreases in diversity have disproportionate effects on ecosystem function [9]. Much of the basis for the current understanding of how warming affects the relationship between biodiversity and ecosystem function (BEF) has been gleaned from laboratory incubations or manipulated communities, which has limited our ability to accurately predict future ecosystem productivity. Clarifying whether the dual stresses of increased temperature and biodiversity loss affect ecosystem function in natural soils in the field is therefore critical to gaining a clearer understanding of the potential impacts of climate change.

Long-term global change experiments, like those at Harvard Forest Long-Term Ecological Research (LTER) site, provide a unique opportunity to examine long-term soil warming in a field setting. Given that the effects of warming on ecosystem function can change or compound over time [5, 10], maintaining and studying long-term climate change experiments is critical to the ability to predict ecosystem function. Additionally, microbial community compositional shifts resulting from climate change can take over a decade to manifest [11, 12]. How microbial communities respond to climate change is influenced by the duration of stress, with acclimation theorized to occur in the short-term and adaptation theorized to occur in the long-term [13]. Two soil warming experiments were established in 1991 and 2006 at the Harvard Forest, where soils have been continuously heated +5°C above ambient along with control soils that have not been heated. Both the soil carbon quantity and quality has declined in response to heating [10], resulting in an increase in lipid concentration in the heated mineral soils [14]. The effects of warming on soil respiration at these sites has been nonlinear over decadal time scales at Harvard Forest [10], where there are years when there is a large treatment effect of heating compared to controls and years where there are no differences between heated and control plots. These changes in ecosystem function and nutrient pools have occurred in conjunction with changes to the microbial community. The abundance of Alpha-proteobacteria as well as Acidobacteria have increased in the heated plots [11], and the overall diversity of the fungal community has decreased in the heated plots [7]. Altogether, these changes call into question how the relationship between ecosystem function and bacterial diversity will be altered by global warming.

Climate change alters both abiotic and biotic drivers of ecosystem function, which complicates predicting the response of ecosystem function to climate stresses, such as warming. One approach to addressing and condensing the dissimilar responses of individual ecosystem functions to climate change is the use of ecosystem multifunctionality (EMF) indices, which aggregate individual ecosystem functions into a single index. Compared to single function approaches, which examine ecosystem functions separately, EMF indices are more sensitive to biodiversity changes [15] and can better represent inherent trade-offs between ecosystem functions [16]. For instance in plant trait simulation study, ecosystem nitrogen use efficiency and litter quality had a negative relationship, indicating a trade-off between these two functions [17]. An EMF meta-analysis that incorporated results from 68 different warming studies observed that warming tends to have positive effects on EMF [18]. While kinetic effects are partially driving the positive effect of warming on EMF [19], warming-related changes in abiotic conditions, such as more limited substrate availability, also can increase certain ecosystem functions, like enzyme activity [14]. However, warming does not always have a positive effect on EMF, and when climate change alters abiotic conditions, warming can have a negative effect on EMF. For instance, warming decreases EMF under nitrogen limitation [20] or
precipitation declines [21]. Warming can also shift biotic drivers of EMF, like plant commu-
nity functional traits, which in turn will have negative effects on EMF [22]. Microbial diversity
and community composition, two of the most studied biotic drivers of EMF, also are influ-
enced by warming, where decreases in dominant taxa due to warming can result in decreases
in EMF [23]. The complex interactions between biotic and abiotic drivers of EMF have limited
the integration of climate change into the existing framework of BEF relationships.

Temperate deciduous forests experience marked seasonal differences, and further compli-
cating predictions, the effects of warming on ecosystem function and microbial diversity are
also season-dependent. Long-term warming increases soil CO2 efflux [10], and results in faster
upticks in soil respiration after the spring thaw [24]. Additionally, warmed soils exhibit higher
temperature sensitivity of respiration in the spring and winter [25], suggesting that with
increasing temperatures, soil CO2 efflux will increase non-linearly. Microbial biomass similarly
has a seasonal cycle that tends to peak in fall [26, 27], although decreases in microbial biomass
due to warming might alter this cycle [14]. Extracellular enzyme activity is likewise responsive
to the seasonal shifts in the availability of carbon. Following or coinciding with leaf fall, the
activity of peroxidase, phenol oxidase, and other oxidative enzymes involved in the decompo-
sition of lignin tends to increase [27, 28]. Microbial diversity likewise exhibits seasonal shifts;
however, the season with highest diversity varies from study to study. Summer, with its higher
soil temperatures and higher availability of photosynthates as root exudates, has been observed
as the season with the highest microbial diversity [29, 30]. Conversely winter has been noted as
having the highest microbial diversity, due to the availability of products from the decomposi-
tion of fresh leaf litter in the fall [31]. Moreover, season is consistently identified as a driver of
microbial community composition in forest soils [26, 32], where the community composition
shifts along with leaf fall. While during the spring and summer, the community is dominated
by organisms that can utilize root exudates, such as Proteobacteria, the fresh influx of leaf litter
that occurs during the fall selects for organisms that can decompose lignin, such as Actinobac-
teria and Mucilaginibacter [30, 31, 33]. Despite variability in microbial communities and eco-
system functions across seasons, most studies that examine the relationship between
ecosystem multifunctionality and microbial diversity utilize a single snapshot of the commu-
nity and ecosystem at the time of sampling.

Here, we examined how 13 or 28 years of long-term soil warming affected ecosystem multi-
functionality and bacterial diversity (Shannon's H), richness (Chao1), and evenness (Pielou's
J) in two different seasons. We utilized previously published ecosystem function data [27], as
well as additional ecosystem function data (soil total organic carbon and nitrogen) paired with
bacterial 16S rRNA gene amplicon sequencing to test the hypothesis that warmed soils would
have a stronger relationship between bacterial diversity and ecosystem multifunctionality than
the control soils.

**Materials and methods**

Site description and sample collection

Soil samples were collected from two long-term warming soil experiments located at the Har-
vard Forest in Petersham, MA (42°30′30″N, 72°12′28″W), as described in [27]. The Prospect
Hill Soil Warming Study was established in 1991 [34], and the Soil Warming x Nitrogen Addi-
tion (SWaN) Study was established in 2006 [24]. Warmed plots are heated continuously +5°C
above ambient using buried resistance cables placed 10 cm below the soil surface and spaced
20 cm apart. Soil samples were collected in 2019 on July 15th and October 19th, from Prospect
Hill and SWaN, which had been warmed for 28 and 13 years, respectively, at the time of sam-
pling. The experiments are located adjacent to one another and have the same dominant plant
overstory (Acer rubrum, Betula lenta, Betula papyrifera, Fagus grandifolia, Quercus velutina, Quercus rubra) and soil type (coarse-loamy incepitsols). The temperature ranges from a mean low of 3.3°C to mean high of 13.1°C [35]. The mean annual precipitation, including snow, is 1107 mm distributed evenly across seasons [35]. Both Prospect Hill and SWaN organic and mineral soils are acidic, with the organic horizon having pHs between 3.8 and 4.2 and the mineral soils having pHs between 3.9 and 4.4 [14, 36, 37]. The organic horizon in the control, unheated plots has a depth of approximately 5 cm, and the organic horizon in the heated +5°C plots has a depth of approximately 2–3 cm [14]. Duplicate cores were taken from each plot to 10 cm depth using a 5.7 cm diameter tulip bulb corer. Cores were separated into the organic and mineral horizons, duplicate soil cores were pooled by depth increment, roots and rocks were removed, and then the soil was sieved <2 mm. After sieving, samples were kept at ambient temperature, and within 4 hours of collection, samples were taken back to the lab for further analyses. The full sampling design was two sites (Prospect Hill, SWaN) x two treatments (control, heated) x two seasons (summer, fall) x two soil layers (organic horizon, mineral soil) x five replicate plots for a total of 78 samples (one sample from the SWaN experiment had a mineral horizon that was beyond the reach of the corer for the July sampling and this same plot was subsequently not sampled in October).

## Ecosystem functionality measurements

To assess ecosystem functionality, we utilized five soil properties or functions that were measured as a part of Domeignoz-Horta et al. (2023), including microbial biomass carbon, respiration, and the potential activities of four enzymes: phenol oxidase and peroxidase, β-glucosidase (BG), and N-acetyl-glucosaminidase (NAG). Additionally, we measured total soil organic carbon and total nitrogen. Altogether, there were seven different soil functions or properties measured for this study.

To measure total soil organic C and N and soil water content, soils were weighed into pre-weighed aluminum tins and dried in a 65°C oven until they reached a constant mass the same day the samples were brought back to the lab. Constant mass was verified by weighing the soils multiple times throughout the drying process. The drying temperature was selected so soils could be utilized for a separate analysis that was performed in Domeignoz-Horta et al. (2023) but not included in this study. After the soils achieved a constant mass, the tins were weighed again. The dry soil weight was subtracted from the wet soil weight, then divided by the dry soil weight to calculate water content. For total C and N, the dried soils were ground to a fine powder using a mortar and pestle. The soil was then weighed and packaged in duplicate into tins, which were run on a Perkin Elmer 2400 Series II CHN Elemental Analyzer with acetanilide as a standard at the University of New Hampshire Water Quality Analysis Lab. Total C and total N or duplicates were averaged.

Within three days of soil sampling, microbial biomass carbon was measured. Soils were stored at 15°C in the intervening time between sample collection and MBC measurement. The three day storage was to allow for similar treatment as samples used in carbon use efficiency calculations for Domeignoz-Horta et al. (2023). Four replicate soil samples were each split into two subsamples, with one group serving as a control and one group which was fumigated with chloroform vapors under vacuum pressure for 24 hours. Dissolved organic carbon (DOC) was then extracted from both the unfumigated and fumigated samples using 15 mls of 0.05 M K₂SO₄ and quantified on a Shimadzu TOC analyzer. Microbial biomass carbon was determined by subtracting the DOC concentration in the unfumigated subsample from the fumigated subsample.
Soils were stored at room temperature (20°C) overnight before being aliquoted for soil respiration measurements. Soil respiration was measured on triplicate subsamples (0.15 or 0.3 g for the organic horizon or mineral soil, respectively) that were placed into Hungate tubes. Tubes were then sealed and incubated at 15°C for 24 hours. A 30 ml headspace sample was then taken and injected into an infrared gas analyzer (Quantek 906) to measure CO2 concentrations.

Prior to extracellular enzyme activity assays, soils were stored for no more than 4 days at 15°C after soil sample collection. Potential extracellular enzyme activity was measured using fluorescent substrates. Soil slurries were prepared with 1.25 g wet weight soil and 175 mls of 50 mM pH 4.7 sodium acetate in a Waring blender. For the BG and NAG assays, 200 µls of soil slurry was pipetted into black 96 well plates, and for the oxidative enzyme assay, 500 µls of soil slurry was pipetted into deep well plates. Plates were then placed in a 15°C incubator for 25 minutes to allow for temperature acclimation. This temperature reflects the average air temperature between the summer sampling and the fall sampling. After temperature acclimation, either 50 µls of 4000 µM 4-methylumbelliferyl β-D-glucopyranoside, or 50 µls of 2000 µM 4-methylumbelliferyl N-acetyl-glucosaminidase were added to each well. For assessing phenol oxidase and peroxidase activity, 500 µls of 25 µM L-DOPA + 0.03% H2O2 were added to each well. Each plate contained a standard curve as well as a slurry-only control. All plates were read on a SpectraMax M2 plate reader. BG and NAG plates were measured at 360/450 nm excitation/emission wavelengths after substrate addition. Oxidative enzyme plates were incubated for 4 hours after substrate addition, and then 100 µls were removed, transferred to a clear 96 well plate, and read at 460 nm. Since phenol oxidase and peroxidase both act on L-DOPA and 0.03% H2O2, any enzyme activity measured using these substrates was called oxidative enzyme activity. All enzyme activity was normalized for each sample by the sample's microbial biomass carbon.

## Ecosystem multifunctionality calculation

Multifunctionality was calculated for the organic horizons and mineral soil in R [38] using the multifunc package [39]. We included microbial biomass carbon, soil respiration, total carbon, total nitrogen, N-acetyl glucosaminidase activity, β-glucosidase activity, and oxidative enzyme activity in the ecosystem multifunctionality index. These functions represent both process rates (respiration, potential enzyme activity) and nutrient pools (total C, total N, microbial biomass) [40]. The organic horizon and mineral soil layers were analyzed separately due to documented differences in soil parameters [14] and microbiome properties [11, 41]. A multifunctionality index was calculated for each sample by taking the average of the z-score transformed ecosystem functions. While there are different ways of assessing diversity-multipfunctionality relationships, we selected the averaging approach for ease of interpretation. Individual ecosystem functions' relationship with diversity was also investigated. For this approach, the raw untransformed function measurements were used.

## DNA extraction and library preparation

Bacterial diversity was measured using 16S ribosomal RNA (rRNA) gene amplicon sequencing. After sampling, soils were stored for three days at 15°C to keep treatments consistent for an analysis in Domeignoz-Horta et al. (2023) that was not included in this study. After this three day storage period, soils were stored at -80°C until extraction. DNA was extracted from soils using the Qiagen Powersoil kit following the manufacturer's protocol, and DNA concentration was measured using the Picogreen dsDNA kit (ThermoFisher Scientific).
We sequenced the 16S rRNA gene V4 region using the primers 515F (5' - GTG YCA GCM GCC GCG GTA A- 3') and 806R (5' - GGA CTA CNV GGG TWT CTA AT - 3') following the Earth Microbiome Project protocol [42], with minimal modifications. Forward primers contained the 5' Illumina adapter, 12 basepair unique Golay barcode, forward primer pad and linker, and 515F. Reverse primers contained the 3' reverse complement Illumina adapter, reverse primer pad and linker, and 806R. Templates were amplified in duplicate prior to sequencing using a 25 µl reaction containing 10 µM 2X Invitrogen Platinum Hot Start master mix, 0.5 µl 10 µM forward primer, 0.5 µl 10 µM reverse primer, 13 µl molecular grade H₂O, and 1 µl template. Samples were amplified in an Eppendorf Mastercycler Pro thermocycler with the program parameters of 94°C for 3 min, 35 cycles of 94°C for 45 s, 50°C for 60 s, 72°C for 90 s, and a final extension at 72°C for 10 min. Technical replicate amplicons were combined and visualized on a 1% agarose gel. We quantified the amplicons with the PicoGreen Assay for dsDNA (ThermoScientific) and pooled 300 ng of each sample. We used AMPure XP magnetic beads (Beckman Coulter) to clean the pooled library. Prior to sequencing, the quality of the amplicon library was checked using Qubit Fluorometer (Thermo Fisher Scientific Inc) and 2100 Bioanalyzer DNA 7500 assay (Agilent Technologies, Inc). The library quantification was done using NEB-Next Library Quant Kit for Illumina (New England Biolabs). The amplicon library was spiked with Illumina PhiX and sequenced on Illumina MiSeq platform using the v2-300 cycle kit with 156 bp paired end chemistry (Illumina Inc). Custom Read1, Read2, and Index Sequencing primers were used for sequencing. Negative controls with no template were PCR amplified and sequenced as well. Amplicon library quality assessment and sequencing was performed at the Genomics Resource Laboratory (RRID:SCR_017907), Institute for Applied Life Sciences, University of Massachusetts Amherst, MA, USA.

We were also interested in measuring fungal community diversity using ITS amplicon sequencing. However, we were unable to successfully PCR amplify enough samples to have triplicates for each combination of season, soil layer, warming treatment, and warming duration. Thus, the focus of this study considered only on the bacterial community.

## Sequencing processing and analysis

Raw FastQ files were demultiplexed and primers trimmed using cutadapt [43]. Bacterial 16S rRNA gene sequences were processed with the DADA2 pipeline (v.1.18.0) [44], which generates amplicon sequence variants (ASV), in R. Forward reads were trimmed to 150 bp, and reverse reads were trimmed to 140 bp. Any reads with expected error higher than 2 or with any ambiguous nucleotides were discarded, and all PhiX sequences were also removed. Error rates for the forward and reverse reads were estimated separately using the function 'learnErrors'. ASVs were called for the forward and reverse reads using their respective error models. Denoised forward and reverse ASVs were merged using the function 'mergePairs' with default parameters. Prior to assigning taxonomy, any chimera sequences were identified *de novo* and removed. Taxonomy was assigned using SILVA (v 138.1) [45]. All sequences that were classified at the family level as "mitochondria" or at the order level as "chloroplast" were filtered out from the final ASV table. To account for differences in sequencing depth between samples, raw ASV counts were normalized with library size estimation factors calculated in DESeq2 [46, 47]. Any ASVs that had 1 or 0 reads total across all samples after normalization were discarded.

Our final data set comprised 69 out of 79 possible samples due to difficulty with successfully PCR amplifying samples during library preparation. In the final data set, each combination of season, soil layer, warming treatment, and warming duration had at least three replicates out of the original five replicates, except the October organic horizon control treatment from the
13-year experiment, which had only two out of five replicates. At the end of sequence processing, the data set had 4,345,430 reads, with an average depth of 62,977 reads per sample, and 8,754 total ASVs. The no template negative controls had fewer than 500 reads. After correcting with a sample specific library size estimate calculated in DESeq2 [46], the total read count was 4,159,442, with an average read count of 60,281 per sample. Rarefaction curve analyses indicated that for the samples that were successfully sequenced, the sequencing efforts were adequate.

Since a variety of diversity metrics have been correlated with EMF, we calculated Simpson diversity, Shannon's H, Chao1 estimated richness, community structure (PCoA axes), and evenness (Pielou's J). Diversity metrics were calculated using the packages vegan [48] and phyloseq [49]. The distance matrix for the principal coordinate analysis was calculated using Bray-Curtis dissimilarity using the function 'vegdist' in the package vegan.

## Statistical analysis

All statistical analyses were done in R [38]. Prior to any statistical analysis, we evaluated whether warming treatment and warming duration could be combined into one variable, where a single “control” treatment (or “0 year warmed”) would derive from combined Prospect Hill and SWaN plot control samples, Prospect Hill warmed samples would be “28 year warmed”, and SWaN warmed samples would be “13 year warmed”. Both EMF and bacterial diversity metrics were tested for normality using a Shapiro-Wilk test. After checking for normality, we tested for significant differences between the control plot samples from Prospect Hill and SWaN using either a Wilcoxon rank sum test if the data were not normally distributed (as was the case for both the organic horizon and mineral soils’ EMF), or a t-test if the data were normally distributed (as was the case for the organic horizon and mineral soils’ bacterial diversity metrics). We found no significant differences between the control plots’ EMF in Prospect Hill and SWaN. However, we did find significant differences in bacterial diversity between the control plots from the two experiments (p = 0.0260). As a result, warming treatment and warming duration were combined into a single variable when EMF was the response variable, but not when bacterial diversity was the response variable.

We also investigated whether percent water content correlated with ecosystem multifunctionality or any of the bacterial diversity metrics in the organic horizon or mineral soils. Prior to correlation testing, ecosystem multifunctionality, Shannon H, Chao1 estimated richness, and Pielou J were tested for normality using a Shapiro-Wilk test. Of these four variables, the organic horizon bacterial diversity metrics and mineral soil Shannon's H were normally distributed. All other variables were not normally distributed. For variables that were normally distributed, we used a Pearson correlation test. For variables that were not normally distributed, we used a Spearman correlation test.

We built a set of general linear models using raw or log-transformed data to investigate the effects of warming, soil layer, season, and/or site on four response variables: individual ecosystem functions, nutrient pools, EMF, or bacterial diversity. When normality of model residuals was not satisfied or the response variable was not normally distributed, we constructed generalized linear models with non-Gaussian distributions. We considered both additive and interactive effects between warming, soil layer, season, and/or site.

A set of candidate generalized linear models were constructed with either EMF or bacterial diversity as the response variable. First, models were built using a Gaussian distribution and model residuals were tested for normality, except for models where Chao1 estimated richness was the response variable. When Chao1 estimated richness was the response variable, models utilized a Poisson or negative binomial distribution, which are both appropriate
for count data. Poisson models for Chao1 estimated richness were checked for overdispersion by looking at the ratio between residual deviance and degrees of freedom, and if overdispersion was observed (ratio > 1.1), a negative binomial distribution was used with the R package MASS [50]. For the mineral soils, any model where EMF was the response variable had a log-transformation applied to EMF to make the residuals normally distributed. The results from the log-transformed models are reported without back-transformation. We found that the different diversity metrics all performed similarly based on their ΔAICc, so we elected to utilized Shannon's H, Pielou's J, and Chao1 estimated richness as predictors since taxon richness and community evenness have both been found to be drivers of ecosystem function [51, 52].

We also tested whether warming treatment, warming duration, soil layer, and/or season affected the relationship between the individual ecosystem functions or nutrient pools and bacterial diversity by constructing linear models. Model residuals were tested for normality using a Shapiro-Wilk test. If residuals were not normally distributed, a log-transformation was applied. If the model residuals were still not normally distributed after a log-transformation, we constructed a set of generalized linear models using a gamma or inverse Gaussian distribution. Residuals for these models were validated using the DHARMa package [53]. Results for models where a log-transformation was applied or for models that did not use a normal distribution are reported without back-transformation.

Finally, we then investigated how warming treatment, warming duration, soil layer, and/or season affected the relationship between EMF and bacterial diversity. For all of these analyses, we constructed a set of candidate linear models. In these sets of candidate models, we considered all predictors and possible interactions between warming treatment and diversity, warming treatment and horizon, warming treatment and season, or warming treatment and warming duration. All model residuals were checked for normality, and if necessary a log-transformation was applied to the response variable. Estimates and errors from the log-transformed models are reported without back-transformation.

For all of the candidate sets of models, model fit was assessed using AICc in the package [54]. If models were within 2 ΔAICc units, we considered the models equally supported, and used a likelihood ratio test for the nested models to select a model. If the models with similar ΔAICc were not nested, we selected the model that explained more variance in the data. Model residual deviance was also compared to model null deviance to examine model fit. When significant interactions were found between warming treatment, warming duration, and/or season, a Tukey's HSD test was utilized to further investigate the interaction. Otherwise, p-values had a Benjamini-Hochberg correction applied to them to account for multiple comparisons. We utilized a p-value cutoff of 0.05 for statistical significance. If the p-values were between 0.05 and 0.1, we report it as a trend. All estimates are reported with a 95% confidence interval.

After analyzing the bacterial diversity and finding a significant difference in bacterial evenness between the warmed and control treatments in the summer organic horizon samples, we conducted a *post-hoc* analysis to see whether any of the dominant ASVs were differentially abundant between control and heated plots. We used three different methods: Analysis of Community Microbiomes with Bias Correction in the R package ANCOM-BC (v.2.0.3) [55], log2 fold change using the package DESeq2 [46], and a $\chi^2$ test between ASV counts in the heated and control plots. Dominance was determined based on a relative abundance threshold of ≥ 0.001 within each sample. Of the 8,629 ASVs in the organic horizon summer dataset, 256 ASVs were identified as dominant, which accounted for 2.967% of the total ASVs. For all three differential abundance methods, p-values were adjusted for multiple comparisons using a Benjamini-Hochberg correction.
<table><caption>Table 1. The effects of warming treatment, season, or site on diversity metrics in the organic horizon and mineral soil.</caption><thead><tr><th>soil layer</th><th>diversity metric</th><th>predictor</th><th>estimate</th><th>error</th><th>p value</th></tr></thead><tbody><tr><td rowspan="12">Organic horizon</td><td rowspan="3">Shannon diversity</td><td>season (15)</td><td>-0.076</td><td>0.157</td><td>n.s.</td></tr><tr><td>warming (16)</td><td>-0.134</td><td>0.151</td><td>n.s.</td></tr><tr><td>warming:season (9)</td><td>0.247</td><td>0.220</td><td>0.072</td></tr><tr><td rowspan="2">Chao1 estimated richness</td><td>site (16)</td><td>-0.110</td><td>0.086</td><td>0.019 *</td></tr><tr><td>warming (16)</td><td>0.038</td><td>0.086</td><td>n.s.</td></tr><tr><td rowspan="3">Pielou evenness</td><td>season (15)</td><td>-0.011</td><td>0.016</td><td>n.s.</td></tr><tr><td>warming (16)</td><td>-0.024</td><td>0.016</td><td>0.011 *</td></tr><tr><td>warming:season (9)</td><td>0.031</td><td>0.024</td><td>0.015 *</td></tr><tr><td rowspan="12">Mineral soil</td><td rowspan="2">Shannon diversity</td><td>season (19)</td><td>0.011</td><td>0.151</td><td>n.s.</td></tr><tr><td>warming (19)</td><td>0.001</td><td>0.015</td><td>n.s.</td></tr><tr><td rowspan="2">Chao1 estimated richness</td><td>site (19)</td><td>-0.047</td><td>0.122</td><td>n.s.</td></tr><tr><td>warming (19)</td><td>0.017</td><td>0.122</td><td>n.s.</td></tr><tr><td rowspan="2">Pielou evenness</td><td>site (19)</td><td>-0.006</td><td>0.022</td><td>n.s.</td></tr><tr><td>warming (19)</td><td>0.002</td><td>0.022</td><td>n.s.</td></tr></tbody></table>

Shannon diversity models and organic horizon Pielou evenness models used a Gaussian distribution, Chao1 estimated richness models used a negative binomial distribution, and mineral soil Pielou evenness models used a gamma distribution. Chao1 estimated richness and mineral soil Pielou evenness estimates are reported in the table without back-transformation. Estimates are reported where the reference for warming treatment is control non-heated plots, the reference for season is summer, and the reference for site is the 13 year experiment. Number of samples for the contrasts are reported in parentheses. Errors are reported for a 95% confidence interval. Any significant interactions were further investigated using a Tukey HSD test. Benjamini-Hochberg adjusted p-values are reported, and any adjusted p-values that were greater than 0.05 are reported as not significant (n.s.).

https://doi.org/10.1371/journal.pone.0311364.t001

# Results

## Bacterial diversity

The effect of warming on bacterial evenness varied by season in the organic horizon, where the best fitting model based on ΔAICc was the interaction model with warming and season (Table 1). Warming had a negative impact on Pielou's J in the summer but no effect in the fall (Tukey HSD summer: -0.024 ± 0.011, p = 0.026; Tukey HSD fall: p = 0.801). This pattern was not observed in the mineral soil, where warming did not affect bacterial evenness. None of the bacterial diversity metrics in the organic horizon had a significant correlation with soil water content (S1 Table). There was also a season-dependent warming trend on Shannon diversity, although this effect was not significant (p = 0.072). Since in the organic horizon in the summer warming decreased bacterial evenness but did not affect Chao1 richness, we carried out a *post-hoc* differential abundance analysis on dominant ASVs to determine if changes in the abundance of dominant ASVs could explain the decrease in evenness. None of the dominant ASVs were differentially abundant between the warmed and control plots.

Chao1 estimated richness significantly differed between the two different sites, Prospect Hill (28 year warmed) and SWaN (13 year warmed) in the organic horizon. The older experimental site had around 90% of the number of taxa compared to the younger site (Table 1). Soil water content did not have a significant correlation with any of the bacterial diversity metrics in the mineral soils (S1 Table). Unheated control plots from the 13-year warmed and 28-year warmed experiments had significantly different bacterial diversity metrics and so were not combined when bacterial diversity was the response variable.
Ecosystem multifunctionality

Ecosystem multifunctionality was affected only by season and not warming treatment or warming duration. In the organic horizon, EMF was significantly higher in fall compared to summer (0.089 ± 0.064, p = 0.022, Fig 1A). Similarly, in the mineral soil, EMF trended higher in fall compared to summer, although this difference was not statistically significant (p = 0.090, Fig 1B). Neither EMF in the organic horizon nor EMF in the mineral soils had a significant correlation with soil water content (S1 Table). Due to known differences between the organic horizon and mineral soil [14], ecosystem multifunctionality was calculated for the soils separately.

For both the organic horizon and mineral soils, neither 13 years nor 28 years of warming significantly altered EMF. In mineral soils, EMF tended to increase as warming duration

<figure><img src="image_4.png" /><figcaption>**Fig 1. Ecosystem multifunctionality across warming treatment and season.** Ecosystem multifunctionality (EMF) was calculated by averaging z-score standardized ecosystem function measurements for each sample, and EMF was calculated separately for the organic horizon (A) and mineral soils (B). In the mineral soils (B), EMF was log-transformed for general linear model construction, but original non-log-transformed values are displayed.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0311364.g001
increased, but these differences were not significant. In the organic horizon soils, we did not observe any trend in EMF with longer duration of heating.

## Ecosystem multifunctionality-bacterial diversity relationship

In the organic horizon, season had a significant effect on the relationship between EMF and bacterial Chao1 estimated richness, where fall had significantly higher intercept for EMF than summer, irrespective of warming treatment or duration (Fig 2A, 0.108 ± 0.074, p = 0.020). We

<figure><img src="image_5.png" /><figcaption>**Fig 2.** The relationship between ecosystem multifunctionality and bacterial Chao1 estimated richness. Ecosystem multifunctionality (EMF) was calculated by averaging z-score standardized function measurements for each sample. Shaded regions represent a 95% confidence interval. Mineral soil EMF was log-transformed and is displayed without back-transformation.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0311364.g002
observed a similar, seasonal trend in the relationship between EMF and Shannon's diversity (S2A Fig, 0.105 ± 0.075, p = 0.052) and Pielou's evenness (S2C Fig, 0.104 ± 0.074, p = 0.055). In the organic horizon soils, ecosystem multifunctionality did not correlate with bacterial Chao1 estimated richness, Shannon diversity, or Pielou evenness. The relationship between EMF and Chao1 estimated richness had a progressively lower intercept as warming duration increased. However, neither the 13-year warmed nor the 28-year warmed plots showed a relationship between EMF and Chao1 estimated richness that was significantly different from the control plots.

In the mineral horizon, Shannon diversity (S2B Fig), Chao1 estimated richness (Fig 2B), or Pielou evenness (S2D Fig) were not correlated with EMF, and the intercept for the relationship between EMF and Chao1 estimated richness did not change as warming duration increased. There also was no effect of season on the intercept for EMF for all three of the diversity metrics in the mineral soil.

## Individual ecosystem functions-bacterial diversity relationship

In the organic horizon, no individual ecosystem function or nutrient pool had a significant relationship with any diversity metric. MBC and oxidative enzyme activity both varied significantly by season, with MBC higher in fall compared to summer and oxidative enzyme activity higher in fall compared to summer for both the Shannon diversity and Chao1 estimated richness models (Table 2 and S2 Table).

In the mineral soil, β-glucosidase activity was the only ecosystem function that had a relationship with any bacterial diversity metric, where Shannon's H had a negative relationship with β-glucosidase activity (Table 3 and S3 Table). Similar to the organic horizon, significant seasonal differences were observed for MBC, oxidative enzyme activity, and β-glucosidase

<table><caption>Table 2. Single ecosystem function-diversity relationships in the organic horizon.</caption><thead><tr><th>ecosystem function</th><th>diversity metric</th><th>predictor</th><th>estimate</th><th>error</th><th>p value</th></tr></thead><tbody><tr><td rowspan="12">MBC</td><td rowspan="4">Shannon diversity</td><td>season (15)</td><td>2149.600</td><td>537.432</td><td>9.942e-08 ***</td></tr><tr><td>13 years heating (9)</td><td>-230.400</td><td>622.104</td><td>n.s.</td></tr><tr><td>28 years heating (7)</td><td>-470.800</td><td>675.612</td><td>n.s.</td></tr><tr><td>Shannon (32)</td><td>341.100</td><td>1654.436</td><td>n.s.</td></tr><tr><td rowspan="4">Chao1 estimated richness</td><td>season</td><td>2165.129</td><td>533.296</td><td>7.456e-08 ***</td></tr><tr><td>13 years heating (9)</td><td>-246.062</td><td>636.418</td><td>n.s.</td></tr><tr><td>28 years heating (7)</td><td>-477.112</td><td>678.174</td><td>n.s.</td></tr><tr><td>Chao1 (32)</td><td>0.1276</td><td>2.248</td><td>n.s.</td></tr><tr><td rowspan="12">OX activity</td><td rowspan="4">Shannon diversity</td><td>season</td><td>-11573.7</td><td>1724.212</td><td>1.468e-12 ***</td></tr><tr><td>13 years heating (9)</td><td>122.400</td><td>231.476</td><td>n.s.</td></tr><tr><td>28 years heating (7)</td><td>101.300</td><td>233.828</td><td>n.s.</td></tr><tr><td>Shannon (32)</td><td>245.900</td><td>545.664</td><td>n.s.</td></tr><tr><td rowspan="4">Chao1 estimated richness</td><td>season</td><td>-1.157e+04</td><td>1765.960</td><td>2.579e-12 ***</td></tr><tr><td>13 years heating (9)</td><td>1.419e+02</td><td>229.712</td><td>n.s.</td></tr><tr><td>28 years heating (7)</td><td>1.225e+02</td><td>229.908</td><td>n.s.</td></tr><tr><td>Chao1 (32)</td><td>0.186</td><td>0.706</td><td>n.s.</td></tr></tbody></table>

Microbial biomass (MBC) models used a Gaussian distribution and oxidative enzyme (OX) models used an inverse Gaussian distribution. Oxidative enzyme model estimates are reported without back-transformation. All model estimates are reported where the reference for warming treatment is control non-heated plots and the reference for season is summer. Number of samples are reported in parentheses. Errors are reported for a 95% confidence interval. Benjamini-Hochberg adjusted p-values are reported. Any adjusted p-values greater than 0.05 are reported as not significant (n.s).
<table><caption>Table 3. Single ecosystem function-diversity relationships in the mineral soil.</caption><thead><tr><th>ecosystem function</th><th>diversity metric</th><th>predictor</th><th>estimate</th><th>error</th><th>p value</th></tr></thead><tbody><tr><td rowspan="12">MBC</td><td rowspan="3">Shannon diversity</td><td>season (19)</td><td>0.944</td><td>0.278</td><td>8.355e-07 ***</td></tr><tr><td>13 years heating (10)</td><td>-0.294</td><td>0.333</td><td>n.s.</td></tr><tr><td>28 years heating (9)</td><td>-0.318</td><td>0.345</td><td>n.s.</td></tr><tr><td rowspan="3">Chao1 estimated richness</td><td>season (19)</td><td>0.942</td><td>0.276</td><td>3.833e-07 ***</td></tr><tr><td>13 years heating (10)</td><td>-0.309</td><td>0.333</td><td>n.s.</td></tr><tr><td>28 years heating (9)</td><td>-0.313</td><td>0.343</td><td>n.s.</td></tr><tr><td rowspan="3">Chao1 (37)</td><td>season (19)</td><td>0.942</td><td>0.276</td><td>3.833e-07 ***</td></tr><tr><td>13 years heating (10)</td><td>-0.309</td><td>0.333</td><td>n.s.</td></tr><tr><td>28 years heating (9)</td><td>-0.313</td><td>0.343</td><td>n.s.</td></tr><tr><td rowspan="12">OX activity</td><td rowspan="3">Shannon diversity</td><td>season (19)</td><td>1.927</td><td>0.247</td><td>1.479e-15 ***</td></tr><tr><td>13 years heating (10)</td><td>-0.329</td><td>0.296</td><td>n.s.</td></tr><tr><td>28 years heating (9)</td><td>-0.320</td><td>0.308</td><td>n.s.</td></tr><tr><td rowspan="3">Chao1 estimated richness</td><td>season (19)</td><td>1.926</td><td>0.247</td><td>7.401e-16 ***</td></tr><tr><td>13 years heating (10)</td><td>-0.334</td><td>0.298</td><td>0.058</td></tr><tr><td>28 years heating (9)</td><td>-0.318</td><td>0.298</td><td>0.063</td></tr><tr><td rowspan="3">Chao1 (37)</td><td>season (19)</td><td>1.926</td><td>0.247</td><td>7.401e-16 ***</td></tr><tr><td>13 years heating (10)</td><td>-0.334</td><td>0.298</td><td>0.058</td></tr><tr><td>28 years heating (9)</td><td>-0.318</td><td>0.298</td><td>0.063</td></tr><tr><td rowspan="12">BG activity</td><td rowspan="3">Shannon diversity</td><td>season (19)</td><td>-744.41</td><td>216.306</td><td>6.389e-07 ***</td></tr><tr><td>13 years heating (10)</td><td>-96.970</td><td>258.720</td><td>n.s.</td></tr><tr><td>28 years heating (9)</td><td>-157.880</td><td>268.598</td><td>n.s.</td></tr><tr><td rowspan="3">Chao1 estimated richness</td><td>season (19)</td><td>-747.101</td><td>232.613</td><td>2.316e-06 ***</td></tr><tr><td>13 years heating (10)</td><td>-91.320</td><td>280.160</td><td>n.s.</td></tr><tr><td>28 years heating (9)</td><td>-150.038</td><td>288.798</td><td>n.s.</td></tr><tr><td rowspan="3">Chao1 (37)</td><td>season (19)</td><td>-747.101</td><td>232.613</td><td>2.316e-06 ***</td></tr><tr><td>13 years heating (10)</td><td>-91.320</td><td>280.160</td><td>n.s.</td></tr><tr><td>28 years heating (9)</td><td>-150.038</td><td>288.798</td><td>n.s.</td></tr><tr><td rowspan="3">Chao1 (37)</td><td>season (19)</td><td>-747.101</td><td>232.613</td><td>2.316e-06 ***</td></tr><tr><td>13 years heating (10)</td><td>-91.320</td><td>280.160</td><td>n.s.</td></tr><tr><td>28 years heating (9)</td><td>-150.038</td><td>288.798</td><td>n.s.</td></tr></tbody></table>

A log-transformation was applied to microbial biomass carbon (MBC) and oxidative enzyme (OX) activity. MBC models, OX models, and β-glucosidase (BG) models all used a Gaussian distribution. MBC model estimates and OX model estimates are reported without back-transformation. All model estimates are reported where the reference for warming treatment is control non-heated plots and the reference for season is summer. Number of samples for contrasts/predictor are reported in parentheses. Errors are reported for a 95% confidence interval. Benjamini-Hochberg adjusted p-values are reported. Any adjusted p-values greater than 0.05 are reported as not significant (n.s.).

https://doi.org/10.1371/journal.pone.0311364.t003

activity. Both MBC and oxidative enzyme activity were higher in fall compared to summer, whereas the opposite was observed for β-glucosidase activity, irrespective of the diversity metric (Table 3).

Discussion

Season is a critical driver of ecosystem function, especially in deciduous forests. We observed that ecosystem multifunctionality was significantly higher in fall compared to summer. Indeed, season had a larger effect on EMF than warming treatment, which had no effect on EMF. The seasonal input of fresh litter could explain the difference in EMF between fall and summer. With leaf-fall, a priming effect could be initiated, resulting in higher microbial activity [27, 56, 57]. This priming could be substantial enough to mitigate the negative impact of warming on EMF. Further supporting a priming effect, we observed a significant increase in microbial biomass in fall, and an increased microbial biomass has also been noted following the addition of fresh carbon in priming experiments [58–60]. In the heated plots at Harvard Forest, soil organic matter is more depleted in simple sugars and is lower in quantity compared to soil
organic matter in the control plots [61], and the influx of litter in the fall might help alleviate this resource limitation [27]. Accordingly, seasonal influx of nutrients cancels out the warming-induced substrate limitation which can be observed by the response of the microbial biomass carbon [27]. This could subsequently increase enzyme production and activity and microbial biomass, both of which were included in our EMF index. While we only examined two seasons, summer and fall, our results highlight the importance of accounting for seasonal differences when measuring ecosystem function.

An additional contributor to the higher EMF observed in fall could be a two-day rainfall event that occurred two days prior to the fall sampling [35]. This rainfall could have initiated a Birch effect, which is an increase in respiration that occurs after rewetting a dried soil [62]. Furthermore, the rainfall could have resulted in a flush of dissolved organic matter, which was found in Domeignoz-Horta et al. (2023). Drying-rewetting cycles can increase microbial biomass [63, 64], which is in line with the significantly higher microbial biomass that we observed in fall (Tables 2 and 3). Warmed soils have significantly lower soil moisture in the spring and fall compared to the control plots in the [24], and the combination of warming and lower moisture can further increase CO₂ efflux after rewetting [65]. Yet, we did not observe a significant increase in respiration in fall [27], as would be expected with a Birch effect. This could be due to the fact that respiration was measured after the soil was taken back to the lab, up to six days after the rainfall. This falls outside of the 1.5 days when the respiration spike that is characteristic of the Birch effect is typically observed [66]. Enzyme activity can also increase after rewetting [67]. While we did observe increased phenol oxidase activity in the fall, we did not observe similar increases in activity for β-glucosidase and N-acetyl glucosaminidase [27]. Finally, we did not observe any significant correlation between soil water content and EMF (S1 Table). Even though we noted higher EMF in the fall, the ecosystem function responses we observed are not entirely consistent with a rewetting response.

Bacterial community evenness tends to be more sensitive to warming compared to bacterial community richness. We observed that bacterial community evenness was negatively affected by soil warming and that this effect was specific to the organic horizon in the summer (Table 1). This contrasts with previous work at our site, where warming increased bacterial community evenness in the organic horizon as a result of decreases in dominant taxa abundance [11]. Warming often drives shifts in community evenness, whether through increased dominance of a single taxon, decreases in dominant taxa, or losses of rare taxa [7, 68, 69]. We did not observe any significant effects of warming on bacterial community richness (Table 1), which is in alignment with other experimental warming studies or meta-analyses [18, 70–72]. Changes in abundance of dominant taxa were not the drivers of this decrease in evenness either.

Changes in substrate quality and quantity due to warming might also explain the observed changes in evenness. At the Harvard Forest, warming has negatively affected both substrate quality and quantity [27], and lower substrate availability likewise has negative effects on microbial community evenness [73]. Given that ecosystem function tends to be more resilient to stress in more even communities [52], the negative effects of warming on evenness ultimately could make ecosystem function less stable. Yet these negative effects might be masked or compensated for by changes to microbial community function as a result of heating, such as increased abundance of carbohydrate activated enzymes (i.e. CAZymes) involved in degrading complex carbohydrates [41] or an expanded number of bacteria able to degrade lignin [74].

In this study, we observed no significant relationship between bacterial diversity, richness, or evenness and EMF, irrespective of season, soil horizon, warming treatment, or warming duration. Typically, ecosystem multifunctionality positively correlates with bacterial diversity [75–77], but neutral, or marginally positive relationships are not uncommon [78, 79]. In fact,
abiotic or environmental factors, such as soil moisture or pH, are consistently identified as stronger drivers of EMF than diversity itself [51, 79–81]. Even the spatial scale of sampling can affect the relationship between EMF and diversity, where at larger scales, like continental or global scales, diversity and ecosystem function tend to be positively correlated, but at smaller scales, like plot level or regional scales, diversity and ecosystem function tend to be neutrally or negatively correlated [82]. Additionally, active microbial diversity [83] or functional trait diversity [84–86] may be a better predictor of ecosystem multifunctionality than total or taxonomic microbial diversity. This may explain why we did not see a relationship between bacterial diversity and EMF. Our DNA-based sequencing approach captured total taxonomic diversity, which comprises both the active and inactive community, and did not provide any functional trait information. Moreover, in soils many bacteria are dormant, with minimal metabolic activity [87], and these dormant cells are not likely contributing to EMF. Altogether, this suggests that when microbial diversity positively correlates with ecosystem function, this correlation tends to be weak, mainly driven by active members of the community or specific functional groups, and ultimately can be better explained by abiotic variables.

Functional redundancy within the bacterial community can also account for the lack of relationship between EMF and bacterial diversity. In theory, if an ecosystem function is functionally redundant within a community, the addition of more taxa that carry out that specific function will have negligible effects on the performance of that function at any given time. Forest soils are known to contain high levels of functional redundancy [88], therefore changes in diversity might not result in changes in function. Additionally, the functions included in the ecosystem multifunctionality index we calculated are known to be widespread within soil microbial communities. For instance, aerobic respiration is widespread throughout the phylogenetic tree, and the genes for β-glucosidase or N-acetyl glucosaminidase production are found or predicted in Proteobacteria, Chloroflexi, and Actinobacteria [89, 90], all of which were abundant in our samples. Conversely, ecosystem functions that are performed by a narrower group of taxa and therefore exhibit less functional redundancy, such as xenobiotic degradation or nitrogen fixation, are more sensitive to changes in diversity [80, 91]. Although, even within functions that are considered to be less widespread throughout the phylogenetic tree, such as nitrite oxidation, there does seem to be some level of functional redundancy or resilience in response to changes in diversity [92]. The lack of a relationship between bacterial diversity, richness, or evenness and EMF could partially be attributed to the functional redundancy of the ecosystem functions in this study.

In these temperate terrestrial forest soils, fungi could have a large role in ecosystem function, which would explain why we did not observe a relationship between any bacterial diversity metric and ecosystem multifunctionality. Both bacterial and fungal diversity are known drivers of ecosystem function [51, 77], and ecosystem multifunctionality-microbial diversity studies rely on the assumption that bacteria and fungi are equally involved in promoting ecosystem functions. While we selected our seven ecosystem functions as representatives of processes such as carbon cycling rates (β-glucosidase activity, respiration, oxidative enzyme activity), nitrogen cycling rates (N-acetyl glucosaminidase activity), or nutrient pools (total carbon, total nitrogen, microbial biomass), some of these functions might be more dependent on fungi. For instance, fungal biomass typically outnumbers bacterial biomass in temperate forests [93], although previous work at the sites in this study found that bacteria outnumbered fungi, based on 16S rRNA or ITS gene copy numbers [11]. Enzyme production and activity could also be skewed towards fungi since fungi are the main producers of oxidative enzymes, like phenol oxidase, in soils [94]. N-acetyl glucosaminidase targets chitin, which is a part of the fungal cell wall, and NAG activity is correlated with fungal biomass as well [95]. While we did initially set out to examine both the bacterial and fungal communities, we were unable to
successfully PCR amplify the fungal community during sequencing library preparation. As a result, we could not assess the role of the fungal community in driving ecosystem multifunctionality. Even though bacteria play a role in ecosystem function, fungi are still major contributors to EMF, possibly explaining why we did not observe a relationship between EMF and any bacterial diversity metric.

While we did not observe a relationship between EMF and bacterial diversity, richness, or evenness, we did observe that β-glucosidase activity had a significant relationship with bacterial Shannon diversity in the mineral soils. Unexpectedly, this relationship was negative (Table 3), meaning that soils with higher diversity tended to have lower β-glucosidase activity potential. We only observed this negative relationship for bacterial Shannon diversity and not bacterial Chao1 richness, suggesting that more even bacterial communities have lower BG activity. This could be explained by negative selection effects, which arise when the most abundant or dominant taxa do not perform or add to a specific ecosystem function, making the diversity-ecosystem function relationship negative or neutral [96]. In fact, β-glucosidase activity can be driven by individual taxa that are not the most abundant within a community [97], and the energetic costs of extracellular enzyme production might negatively impact growth of extracellular enzyme producers [98]. As a result, the presence or absence of individual taxa that are associated with high BG activity might not be reflected in diversity metrics that assess total community diversity. Further, metabolic overlap, which occurs when different taxa within a community can utilize the same substrate or compound, might also account for the negative relationship between diversity and function, since more diverse communities tend to have less metabolic overlap and more diverse metabolic pathways [99]. As it relates to β-glucosidase activity, higher diversity communities might have more diverse metabolisms that can utilize other substrates besides glucosyl, the substrate of BG, and therefore might have lower overall BG activity. Negative selection effects and metabolic overlap can both account for the negative relationship between bacterial diversity and β-glucosidase activity, which has also been observed in other studies [23, 100].

Controlled manipulative laboratory incubations and observational field studies are complementary methods to understand the drivers of the BEF relationship, and both of these methods have validated the relationship between ecosystem function and microbial diversity. Mycorrhizal and saprotrophic fungal diversity and bacterial diversity were all noted as having positive correlations with ecosystem multifunctionality in a study that sampled from 80 different sites across the globe [51]. Manipulative laboratory experiments, where diversity is altered via serial dilutions, have also observed a similar positive trend [101, 102]. However, the strength of the positive correlation tends to decrease when comparing manipulative studies to observational studies, due to factors such as lack of immigration and dispersal, reduced environmental complexity, or short experimental duration [84]. Additional variability in the strength of the relationship between microbial diversity and ecosystem function arises from abiotic factors, like season or soil type [79]. Indeed, we observed no relationship between bacterial diversity and ecosystem function in this observational field study. However, the removal of abiotic variables can affect how diversity confers benefits to ecosystem function. The benefits of diversity arise from mechanisms such as niche partitioning [5] or maintaining a pool of taxa that can thrive under different environmental conditions [52, 103, 104]. Niche partitioning, in particular, is dependent on the presence of environmental heterogeneity in order for organisms to colonize the niche they are best adapted to [105]. Further emphasizing the importance of abiotic variables, we found that season had the largest impact on ecosystem function. Thus, while controlled laboratory experiments have helped to elucidate the mechanisms driving the BEF relationship, observational field studies are still necessary to fully understand the role that abiotic factors, such as season, have on this relationship under more natural conditions.
## Conclusion

This study investigated if long-term soil warming influenced the relationship between ecosystem multifunctionality and bacterial diversity. We observed a strong seasonal effect on ecosystem function, where ecosystem multifunctionality was significantly higher in the fall than in the summer. The effects of soil warming on bacterial evenness were also dependent on season in the organic horizon, where in the summer warming negatively affected bacterial evenness, but in the fall warming did not affect bacterial evenness. However, we found that warming treatment or warming duration ultimately had no effect on the relationship between EMF and bacterial diversity. Indeed, we found that EMF had no relationship with bacterial diversity, which is not uncommon in natural or non-laboratory systems. The younger experimental site also had increased bacterial richness in the organic horizon. Overall, we found that season exerts a strong influence on ecosystem functionality, and season modulates the effects of warming in temperate deciduous forests. We emphasize that cross-season sampling is needed to best assess an ecosystem's total function and that diversity metrics other than taxonomic diversity may be better suited to capturing the role of diversity on ecosystem function.

## Supporting information

**S1 Fig. Relative abundance of the phyla of dominant bacteria based on 16S rRNA gene amplicon sequencing.** Taxonomy was assigned using SILVA (v 138.1). Sample amplicon sequence variant (ASV) counts were normalized using a sample-specific DESeq2 library size estimation correction factor. Dominant ASVs were determined as ASVs with a relative abundance of 0.001 or higher within a sample. Organic horizon samples are presented in panel A, mineral soil samples are presented in panel B.
(tiff)

**S2 Fig. Ecosystem multifunctionality bacterial diversity (Shannon H and Pielou J) relationship.** Ecosystem multifunctionality (EMF) was calculated by averaging a set of z-score transformed ecosystem functions or nutrient pools. Mineral soil EMF was log-transformed and is displayed without back-transformation. Bacterial diversity was measured using 16S rRNA gene amplicon sequencing. Shaded regions represent a 95% confidence interval. The reported p-values are Benjamini-Hochberg corrected to account for multiple comparisons. In the organic horizon, fall had a slightly higher intercept compared to summer for both Shannon H (A, p = 0.052) and Pielou J (C, p = 0.055). There was no seasonal trend in the mineral soils (B, D). In the organic horizon and the mineral soils, there was no significant relationship between EMF and bacterial Shannon H or between EMF and bacterial Pielou J, as well as no significant differences between the control, warmed for 13 years, and warmed for 28 years treatments.
(tiff)

**S1 Table. Correlation between soil water content and ecosystem multifunctionality and bacterial diversity metrics.** Ecosystem multifunctionality (EMF) and all three bacterial diversity metrics (Shannon H, Chao1 estimated richness, Pielou J) were tested for normality. Variables that were normally distributed (organic horizon Shannon H, Chao1 estimated richness, and Pielou J and mineral soils Shannon H) used a Pearson correlation test, and variables that were not normally distributed (organic horizon EMF and mineral soils EMF, Chao1 estimated richness, and Pielou J) used a Spearman correlation test.
(PDF)

**S2 Table. Non-significant single ecosystem function-diversity relationships in the organic horizon.** All models using a Gaussian distribution and model residuals were checked for
normality. To meet normality of model residuals, β-glucosidase (BG) activity, respiration, and total nitrogen in the Chao1 model were log-transformed. The model estimates and errors for the ecosystem functions that were log-transformed are reported without back-transformation. All model estimates are reported where the reference for warming treatment is control non-heated plots and the reference for season is summer. Number of samples are reported in parentheses. Errors are reported for a 95% confidence interval. Benjamini-Hochberg adjusted p-values are reported. Any adjusted p-values greater than 0.05 are reported as not significant (n.s.).

(PDF)

**S3 Table. Non-significant single ecosystem function-diversity relationships in the mineral soils.** All models used a Gaussian distribution, except total carbon and total nitrogen models, which used an inverse Gaussian distribution with an inverse link function and an inverse Gaussian distribution with an $1/μ^2$ link function, respectively. Model residuals were examined, and if normality of residuals for the Gaussian models was not met, a log-transformation was applied. Subsequently, N-acetyl glucosaminidase (NAG) was log-transformed to satisfy normality of model residuals. The model estimates and errors for the ecosystem functions that were log-transformed or that used an inverse Gaussian distribution are reported without back-transformation. All model estimates are reported where the reference for warming treatment is control non-heated plots and the reference for season is summer. Number of samples are reported in parentheses. Errors are reported for a 95% confidence interval. Benjamini-Hochberg adjusted p-values are reported. Any adjusted p-values greater than 0.05 are reported as not significant (n.s.).

(PDF)

**S1 Data. Soil organic carbon, soil organic nitrogen, z-score transformed ecosystem functions, and ecosystem multifunctionality indices metadata.** Metadata for the soil organic carbon, soil organic nitrogen, z-score transformed ecosystem functions, and ecosystem multifunctionality CSV file.

(CSV)

**S1 Dataset. Soil organic carbon, soil organic nitrogen, z-score transformed ecosystem functions, and ecosystem multifunctionality indices dataset.**

(CSV)

**Acknowledgments**

We would like to thank Jody Potter at the University of New Hampshire Water Quality Analysis Lab for assistance with the total organic carbon and nitrogen analysis.

MSS and KMD conceptualized the study. MSS, LADH, MJC, and RR conducted experiments and collected the data. SDF maintains the long-term soil warming experiments at Harvard Forest, facilitated sample collection, and provided input on data interpretation. MSS and MFM analyzed the data. MSS wrote the first draft of the manuscript. All authors contributed to editing and revising the manuscript.

**Author Contributions**

**Conceptualization:** Melissa S. Shinfuku, Kristen M. DeAngelis.

**Data curation:** Melissa S. Shinfuku, Luiz A. Domeignoz-Horta.
**Formal analysis:** Melissa S. Shinfuku, Luiz A. Domeignoz-Horta, Serita D. Frey, Megan F. Mitchell.

**Funding acquisition:** Melissa S. Shinfuku, Serita D. Frey, Kristen M. DeAngelis.

**Investigation:** Melissa S. Shinfuku, Luiz A. Domeignoz-Horta, Mallory J. Choudoir, Ravi Ranjan.

**Methodology:** Melissa S. Shinfuku, Megan F. Mitchell.

**Resources:** Melissa S. Shinfuku, Luiz A. Domeignoz-Horta, Serita D. Frey.

**Software:** Melissa S. Shinfuku, Mallory J. Choudoir, Megan F. Mitchell.

**Visualization:** Melissa S. Shinfuku.

**Writing – original draft:** Melissa S. Shinfuku.

**Writing – review & editing:** Melissa S. Shinfuku, Luiz A. Domeignoz-Horta, Mallory J. Choudoir, Serita D. Frey, Megan F. Mitchell, Ravi Ranjan, Kristen M. DeAngelis.

## References

1. Tilman D, Reich PB, Isbell F. Biodiversity impacts ecosystem productivity as much as resources, disturbance, or herbivory. Proc Natl Acad Sci USA. 2012; 109(26):10394–10397. https://doi.org/10.1073/pnas.1208240109 PMID: 22689971

2. Naeem S, Li S. Biodiversity enhances ecosystem reliability. Nature. 1997; 390(6659):507–509. https://doi.org/10.1038/37348

3. Rohner N, Lepori S, Loaiza V, Sebag D, Verrechia E, Nelson DB, et al. Erosion of community complexity increases temperature-dependency of microbial respiration, but not growth, in short-term incubations. Elementa. 2024; 12(1).

4. Wu L, Zhang Y, Guo X, Ning D, Zhou X, Feng J, et al. Reduction of microbial diversity in grassland soil is driven by long-term climate warming. Nat Microbiol. 2022; 7(7):1054–1062. https://doi.org/10.1038/s41564-022-01147-3 PMID: 35697795

5. Hong P, Schmid B, De Laender F, Eisenhauer N, Zhang X, Chen H, et al. Biodiversity promotes ecosystem functioning despite environmental change. Ecol Lett. 2022; 25(2):555–569. https://doi.org/10.1111/ele.13936 PMID: 34854529

6. Domeignoz-Horta LA, Pold G, Liu XJA, Frey SD, Melillo JM, DeAngelis KM. Microbial diversity drives carbon use efficiency in a model soil. Nat Commun. 2020; 11(1):1–10. https://doi.org/10.1038/s41467-020-17502-z

7. Anthony M, Knorr M, Moore J, Simpson M, Frey S. Fungal community and functional responses to soil warming are greater than for soil nitrogen enrichment. Elementa. 2021; 9(1):000059.

8. Van Nuland ME, Smith DP, Bhatnagar JM, Stefanski A, Hobbie SE, Reich PB, et al. Warming and disturbance alter soil microbiome diversity and function in a northern forest ecotone. FEMS Microbiol Ecol. 2020; 96(7):fiaa108. https://doi.org/10.1093/femsec/fiaa108 PMID: 32472932

9. García FC, Bestion E, Warfield R, Yvon-Durocher G. Changes in temperature alter the relationship between biodiversity and ecosystem functioning. Proc Natl Acad Sci USA. 2018; 115(43):10989–10994. https://doi.org/10.1073/pnas.1805518115 PMID: 30297403

10. Melillo JM, Frey SD, DeAngelis KM, Werner WJ, Bernard MJ, Bowles FP, et al. Long-term pattern and magnitude of soil carbon feedback to the climate system in a warming world. Science. 2017; 358 (6359):101–105. https://doi.org/10.1126/science.aan2874 PMID: 28983050

11. DeAngelis KM, Pold G, Topçuoğlu BD, van Diepen LT, Varney RM, Blanchard JL, et al. Long-term forest soil warming alters microbial communities in temperate forest soils. Front Microbiol. 2015; 6:104. https://doi.org/10.3389/fmicb.2015.00104 PMID: 25762989

12. Frey S, Drijber R, Smith H, Melillo J. Microbial biomass, functional capacity, and community structure after 12 years of soil warming. Soil Biol Biochem. 2008; 40(11):2904–2907. https://doi.org/10.1016/j.soilbio.2008.07.020

13. Martiny JB, Martiny AC, Brodie E, Chase AB, Rodríguez-Verdugo A, Treseder KK, et al. Investigating the eco-evolutionary response of microbiomes to environmental change. Ecol Lett. 2023; 26:S81–S90. https://doi.org/10.1111/ele.14209 PMID: 36965002
14. Pold G, Grandy AS, Melillo JM, DeAngelis KM. Changes in substrate availability drive carbon cycle response to chronic warming. Soil Biol Biochem. 2017; 110:68–78. https://doi.org/10.1016/j.soilbio. 2017.03.002

15. Gamfeldt L, Hillebrand H, Jonsson PR. Multiple functions increase the importance of biodiversity for overall ecosystem functioning. Ecology. 2008; 89(5):1223–1231. https://doi.org/10.1890/06-2091.1 PMID: 18543617

16. Giling DP, Beaumelle L, Phillips HR, Cesarz S, Eisenhauer N, Ferlian O, et al. A niche for ecosystem multifunctionality in global change research. Glob Chang Biol. 2019; 25(3):763–774. https://doi.org/10. 1111/gcb.14528 PMID: 30449061

17. Fiedler S, Monteiro JA, Hulvey KB, Standish RJ, Perring MP, Tietjen B. Global change shifts trade-offs among ecosystem functions in woodlands restored for multifunctionality. J Appl Ecol. 2021; 58 (8):1705–1717. https://doi.org/10.1111/1365-2664.13900

18. Zhou Z, Wang C, Luo Y. Meta-analysis of the impacts of global change factors on soil microbial diversity and functionality. Nat Commun. 2020; 11(1):1–10. https://doi.org/10.1038/s41467-020-16881-7 PMID: 32555185

19. García-Palacios P, Vandegehuchte ML, Shaw EA, Dam M, Post KH, Ramirez KS, et al. Are there links between responses of soil microbes and ecosystem functioning to elevated CO2, N deposition and warming? A global perspective. Glob Chang Biol. 2015; 21(4):1590–1600. https://doi.org/10.1111/gcb. 12788 PMID: 25363131

20. Dong H, Zhang S, Lin J, Zhu B. Responses of soil microbial biomass carbon and dissolved organic carbon to drying-rewetting cycles: A meta-analysis. Catena. 2021; 207:105610. https://doi.org/10.1016/j. catena.2021.105610

21. Ma Y, Tian L, Qu G, Li R, Wang W, Zhao J. Precipitation alters the effects of temperature on the ecosystem multifunctionality in alpine meadows. Front Plant Sci. 2022; 12:824296. https://doi.org/10. 3389/fpls.2021.824296 PMID: 35222457

22. Valencia E, Gross N, Quero JL, Carmona CP, Ochoa V, Gozalo B, et al. Cascading effects from plants to soil microorganisms explain how plant species richness and simulated climate change affect soil multifunctionality. Glob Chang Biol. 2018; 24(12):5642–5654. https://doi.org/10.1111/gcb.14440 PMID: 30239067

23. Chen W, Zhou H, Wu Y, Wang J, Zhao Z, Li Y, et al. Long-term warming impacts grassland ecosystem function: Role of diversity loss in conditionally rare bacterial taxa. Sci Total Environ. 2023; 892:164722. https://doi.org/10.1016/j.scitotenv.2023.164722 PMID: 37308011

24. Contosta AR, Frey SD, Cooper AB. Seasonal dynamics of soil respiration and N mineralization in chronically warmed and fertilized soils. Ecosphere. 2011; 2(3):1–21. https://doi.org/10.1890/ES10- 00133.1

25. Contosta AR, Frey SD, Ollinger SV, Cooper AB. Soil respiration does not acclimatize to warmer temperatures when modeled over seasonal timescales. Biogeochemistry. 2013; 112:555–570. https://doi. org/10.1007/s10533-012-9748-6

26. Contosta AR, Frey SD, Cooper AB. Soil microbial communities vary as much over time as with chronic warming and nitrogen additions. Soil Biol Biochem. 2015; 88:19–24. https://doi.org/10.1016/j.soilbio. 2015.04.013

27. Domeignoz-Horta LA, Pold G, Erb H, Sebag D, Verrecchia E, Northen T, et al. Substrate availability and not thermal-acclimation controls microbial temperature sensitivity response to long term warming. Glob Chang Biol. 2023; 29(6):1574–1590. https://doi.org/10.1111/gcb.16544 PMID: 36448874

28. Baldrian P, Šnajdr J, Merhautová V, Dobiášová P, Cajthaml T, Valášková V. Responses of the extra-cellular enzyme activities in hardwood forest to soil temperature and seasonality and the potential effects of climate change. Soil Biol Biochem. 2013; 56:60–68. https://doi.org/10.1016/j.soilbio.2012. 01.020

29. Du C, Xu CY, Jian JS, He WX, Hou L, Geng ZC. Seasonal dynamics of bacterial communities in a Betula albosinensis forest. Eur J Soil Sci. 2018; 69(4):666–674. https://doi.org/10.1111/ejs.12568

30. López-Mondéjar R, Voříšková J, Větrovský T, Baldrian P. The bacterial community inhabiting temperate deciduous forests is vertically stratified and undergoes seasonal dynamics. Soil Biol Biochem. 2015; 87:43–50. https://doi.org/10.1016/j.soilbio.2015.04.008

31. Landesman WJ, Freedman ZB, Nelson DM. Seasonal, sub-seasonal and diurnal variation of soil bacterial community composition in a temperate deciduous forest. FEMS Microbiol Ecol. 2019; 95(2): fiz002. https://doi.org/10.1093/femsec/fiz002 PMID: 30629168

32. Kaiser C, Koranda M, Kitzler B, Fuchslueger L, Schnecker J, Schweiger P, et al. Belowground carbon allocation by trees drives seasonal patterns of extracellular enzyme activities by altering microbial community composition in a beech forest soil. New Phytol. 2010; 187(3):843–858. https://doi.org/10. 1111/j.1469-8137.2010.03321.x PMID: 20553392
33. Kuffner M, Hai B, Rattei T, Melodelima C, Schloter M, Zechmeister-Boltenstern S, et al. Effects of season and experimental warming on the bacterial community in a temperate mountain forest soil assessed by 16S rRNA gene pyrosequencing. FEMS Microbiol Ecol. 2012; 82(3):551-562. https://doi.org/10.1111/j.1574-6941.2012.01420.x PMID: 22670891

34. Melillo J, Steudler P, Aber J, Newkirk K, Lux H, Bowles F, et al. Soil warming and carbon-cycle feedbacks to the climate system. Science. 2002; 298(5601):2173-2176. https://doi.org/10.1126/science.1074153 PMID: 12481133

35. Boose E, Gould E. Harvard Forest Climate Data since 2001; 2022. https://doi.org/10.6073/pasta/03dc1f107ca816675b4983daf7a97dbe.

36. Anthony MA, Stinson KA, Moore JA, Frey SD. Plant invasion impacts on fungal community structure and function depend on soil warming and nitrogen enrichment. Oecologia. 2020; 194(4):659-672. https://doi.org/10.1007/s00442-020-04797-4 PMID: 33141324

37. Pec GJ, van Diepen LT, Knorr M, Grandy AS, Melillo JM, DeAngelis KM, et al. Fungal community response to long-term soil warming with potential implications for soil carbon dynamics. Ecosphere. 2021; 12(5):e03460. https://doi.org/10.1002/ecs2.3460

38. R Core Team. R: A Language and Environment for Statistical Computing; 2023. Available from: https://www.R-project.org/.

39. Byrnes JE, Gamfeldt L, Isbell F, Lefcheck JS, Griffin JN, Hector A, et al. Investigating the relationship between biodiversity and ecosystem multifunctionality: challenges and solutions. Methods Ecol Evol. 2014; 5(2):111-124. https://doi.org/10.1111/2041-210X.12143

40. Garland G, Banerjee S, Edlinger A, Miranda Oliveira E, Herzog C, Wittwer R, et al. A closer look at the functions behind ecosystem multifunctionality: A review. J Ecol. 2021; 109(2):600-613. https://doi.org/10.1111/1365-2745.13511

41. Roy Chowdhury P, Golas SM, Alteio LV, Stevens JT, Billings AF, Blanchard JL, et al. The Transcriptional Response of Soil Bacteria to Long-Term Warming and Short-Term Seasonal Fluctuations in a Terrestrial Forest. Front Microbiol. 2021; 12:666558. https://doi.org/10.3389/fmicb.2021.666558 PMID: 34512564

42. Caporaso JG, Lauber CL, Walters WA, Berg-Lyons D, Huntley J, Fierer N, et al. Ultra-high-throughput microbial community analysis on the Illumina HiSeq and MiSeq platforms. ISME J. 2012; 6(8):1621-1624. https://doi.org/10.1038/ismej.2012.8 PMID: 22402401

43. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 2011; 17(1):10-12. https://doi.org/10.14806/ej.17.1.200

44. Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods. 2016; 13(7):581-583. https://doi.org/10.1038/nmeth.3869 PMID: 27214047

45. Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2012; 41(D1): D590-D596. https://doi.org/10.1093/nar/gks1219 PMID: 23193283

46. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014; 15(12):1-21. https://doi.org/10.1186/s13059-014-0550-8 PMID: 25516281

47. McMurdie PJ, Holmes S. Waste not, want not: why rarefying microbiome data is inadmissible. PLoS Comput Biol. 2014; 10(4):e1003531. https://doi.org/10.1371/journal.pcbi.1003531 PMID: 24699258

48. Oksanen J, Simpson GL, Blanchet FG, Kindt R, Legendre P, Minchin PR, et al. vegan: Community Ecology Package; 2022. Available from: https://CRAN.R-project.org/package=vegan

49. McMurdie PJ, Holmes S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One. 2013; 8(4):e61217. https://doi.org/10.1371/journal.pone.0061217 PMID: 23630581

50. Venables WN, Ripley BD. Modern Applied Statistics with S. 4th ed. New York: Springer; 2002. Available from: https://www.stats.ox.ac.uk/pub/MASS4/.

51. Delgado-Baquerizo M, Reich PB, Trivedi C, Eldridge DJ, Abades S, Alfaro FD, et al. Multiple elements of soil biodiversity drive ecosystem functions across biomes. Nat Ecol Evol. 2020; 4(2):210-220. https://doi.org/10.1038/s41559-019-1084-y PMID: 32015427

52. Wittebolle L, Marzorati M, Clement L, Balloi A, Daffonchio D, Heylen K, et al. Initial community evenness favours functionality under selective stress. Nature. 2009; 458(7238):623-626. https://doi.org/10.1038/nature07840 PMID: 19270679

53. Hartig F. DHARMA: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models; 2022. Available from: https://CRAN.R-project.org/package=DHARMA.
54. Mazerolle MJ. AICcmodavg: Model selection and multimodel inference based on (Q)AIC(c); 2020. Available from: https://cran.r-project.org/package=AICcmodavg.

55. Lin H, Peddada SD. Analysis of compositions of microbiomes with bias correction. Nat Commun. 2020; 11(1):3514. https://doi.org/10.1038/s41467-020-17041-7 PMID: 32665548

56. Lyu M, Xie J, Vadeboncoeur MA, Wang M, Qiu X, Ren Y, et al. Simulated leaf litter addition causes opposite priming effects on natural forest and plantation soils. Biol Fertil Soils. 2018; 54:925-934. https://doi.org/10.1007/s00374-018-1314-5

57. Zhang Y, Zou J, Meng D, Dang S, Zhou J, Osborne B, et al. Effect of soil microorganisms and labile C availability on soil respiration in response to litter inputs in forest ecosystems: A meta-analysis. Ecol Evol. 2020; 10(24):13602-13612. https://doi.org/10.1002/ece3.6965 PMID: 33391666

58. Jilling A, Keiluweit M, Gutknecht JL, Grandy AS. Priming mechanisms providing plants and microbes access to mineral-associated organic matter. Soil Biol Biochem. 2021; 158:108265. https://doi.org/10.1016/j.soilbio.2021.108265

59. Lyu M, Nie Y, Giardina CP, Vadeboncoeur MA, Ren Y, Fu Z, et al. Litter quality and site characteristics interact to affect the response of priming effect to temperature in subtropical forests. Funct Ecol. 2019; 33(11):2226-2238. https://doi.org/10.1111/1365-2435.13428

60. Zhou J, Wen Y, Shi L, Marshall MR, Kuzyakov Y, Blagodatskaya E, et al. Strong priming of soil organic matter induced by frequent input of labile carbon. Soil Biol Biochem. 2021; 152:108069. https://doi.org/10.1016/j.soilbio.2020.108069

61. VandenEnden L, Anthony MA, Frey SD, Simpson MJ. Biogeochemical evolution of soil organic matter composition after a decade of warming and nitrogen addition. Biogeochemistry. 2021; 156:161-175. https://doi.org/10.1007/s10533-021-00837-0

62. Birch H. The effect of soil drying on humus decomposition and nitrogen availability. Plant Soil. 1958; 10:9-31. https://doi.org/10.1007/BF01343734

63. Evans SE, Wallenstein MD. Soil microbial community response to drying and rewetting stress: does historical precipitation regime matter? Biogeochemistry. 2012; 109:101-116. https://doi.org/10.1007/s10533-011-9638-3

64. Fierer N, Schimel JP. A proposed mechanism for the pulse in carbon dioxide production commonly observed following the rapid rewetting of a dry soil. Soil Sci Soc Am J. 2003; 67(3):798-805. https://doi.org/10.2136/sssaj2003.7980

65. Schindlbacher A, Wunderlich S, Borken W, Kitzler B, Zechmeister-Boltenstern S, Jandl R. Soil respiration under climate change: prolonged summer drought offsets soil warming effects. Glob Chang Biol. 2012; 18(7):2270-2279. https://doi.org/10.1111/j.1365-2486.2012.02696.x

66. Canarini A, Kiær LP, Dijkstra FA. Soil carbon loss regulated by drought intensity and available substrate: A meta-analysis. Soil Biol Biochem. 2017; 112:90-99. https://doi.org/10.1016/j.soilbio.2017.04.020

67. Ouyang Y, Li X. Effect of repeated drying-rewetting cycles on soil extracellular enzyme activities and microbial community composition in arid and semi-arid ecosystems. Eur J Soil Biol. 2020; 98:103187. https://doi.org/10.1016/j.ejsobi.2020.103187

68. Deslippe JR, Hartmann M, Simard SW, Mohn WW. Long-term warming alters the composition of Arctic soil microbial communities. FEMS Microbiol Ecol. 2012; 82(2):303-315. https://doi.org/10.1111/j.1574-6941.2012.01350.x PMID: 22404643

69. Nottingham AT, Scott JJ, Saltonstall K, Broders K, Montero-Sanchez M, Püspök J, et al. Microbial diversity declines in warmed tropical soil and respiration rise exceed predictions as communities adapt. Nat Microbiol. 2022; 7(10):1650-1660. https://doi.org/10.1038/s41564-022-01200-1 PMID: 36065063

70. Wang Z, Lu G, Yuan M, Yu H, Wang S, Li X, et al. Elevated temperature overrides the effects of N amendment in Tibetan grassland on soil microbiome. Soil Biol Biochem. 2019; 136:107532. https://doi.org/10.1016/j.soilbio.2019.107532

71. Yu C, Han F, Fu G. Effects of 7 years experimental warming on soil bacterial and fungal community structure in the Northern Tibet alpine meadow at three elevations. Sci Total Environ. 2019; 655:814-822. https://doi.org/10.1016/j.scitotenv.2018.11.309 PMID: 30481708

72. Zi H, Hu L, Wang C, Wang G, Wu P, Lerdau M, et al. Responses of soil bacterial community and enzyme activity to experimental warming of an alpine meadow. Eur J Soil Sci. 2018; 69(3):429-438. https://doi.org/10.1111/ejss.12547

73. Dilly O, Bloem J, Vos A, Munch JC. Bacterial diversity in agricultural soils during litter decomposition. Appl Environ Microbiol. 2004; 70(1):468-474. https://doi.org/10.1128/AEM.70.1.468-474.2004 PMID: 14711676
74. Pold G, Melillo JM, DeAngelis KM. Two decades of warming increases diversity of a potentially lignolytic bacterial community. Front Microbiol. 2015; 6:480. https://doi.org/10.3389/fmicb.2015.00480 PMID: 26042112

75. Chen QL, Ding J, Zhu D, Hu HW, Delgado-Baquerizo M, Ma YB, et al. Rare microbial taxa as the major drivers of ecosystem multifunctionality in long-term fertilized soils. Soil Biol Biochem. 2020; 141:107686. https://doi.org/10.1016/j.soilbio.2019.107686

76. Mori AS, Isbell F, Fujii S, Makoto K, Matsuoka S, Osono T. Low multifunctional redundancy of soil functional diversity at multiple scales. Ecol Lett. 2016; 19(3):249-259. https://doi.org/10.1111/ele.12560 PMID: 26689733

77. Wagg C, Schlaeppi K, Banerjee S, Kuramae EE, van der Heijden MG. Fungal-bacterial diversity and microbiome complexity predict ecosystem functioning. Nat Commun. 2019; 10(1):1-10. https://doi.org/10.1038/s41467-019-12798-y PMID: 31649246

78. Jiang L, Pu Z, Nemergut DR. On the importance of the negative selection effect for the relationship between biodiversity and ecosystem functioning. Oikos. 2008; 117(4):488-493. https://doi.org/10.1111/j.0030-1299.2008.16401.x

79. van der Plas F. Biodiversity and ecosystem functioning in naturally assembled communities. Biol Rev Camb Philos Soc. 2019; 94(4):1220-1245. https://doi.org/10.1111/brv.12499 PMID: 30724447

80. Delgado-Baquerizo M, Giaramida L, Reich PB, Khachane AN, Hamonts K, Edwards C, et al. Lack of functional redundancy in the relationship between microbial diversity and ecosystem functioning. J Ecol. 2016; 104(4):936-946. https://doi.org/10.1111/1365-2745.12585

81. Jing X, Sanders NJ, Shi Y, Chu H, Classen AT, Zhao K, et al. The links between ecosystem multifunctionality and above-and belowground biodiversity are mediated by climate. Nat Commun. 2015; 6 (1):1-8. https://doi.org/10.1038/ncomms9159 PMID: 26328906

82. Hagan JG, Vanschoenwinkel B, Gamfeldt L. We should not necessarily expect positive relationships between biodiversity and ecosystem functioning in observational field data. Ecol Lett. 2021; 24 (12):2537-2548. https://doi.org/10.1111/ele.13874 PMID: 34532926

83. Bastida F, Torres IF, Moreno JL, Baldrian P, Ondoño S, Ruiz-Navarro A, et al. The active microbial diversity drives ecosystem multifunctionality and is physiologically related to carbon availability in Mediterranean semi-arid soils. Mol Ecol. 2016; 25(18):4660-4673. https://doi.org/10.1111/mec.13783 PMID: 27481114

84. Hillebrand H, Matthiessen B. Biodiversity in a complex world: consolidation and progress in functional biodiversity research. Ecol Lett. 2009; 12(12):1405-1419. https://doi.org/10.1111/j.1461-0248.2009.01388.x PMID: 19849711

85. Steudel B, Hallmann C, Lorenz M, Abrahamczyk S, Prinz K, Herrfurth C, et al. Contrasting biodiversity-ecosystem functioning relationships in phylogenetic and functional diversity. New Phytol. 2016; 212(2):409-420. https://doi.org/10.1111/nph.14054 PMID: 27301904

86. Trivedi C, Delgado-Baquerizo M, Hamonts K, Lai K, Reich PB, Singh BK. Losses in microbial functional diversity reduce the rate of key soil processes. Soil Biol Biochem. 2019; 135:267-274. https://doi.org/10.1016/j.soilbio.2019.05.008

87. Papp K, Mau RL, Hayer M, Koch BJ, Hungate BA, Schwartz E. Quantitative stable isotope probing with H218O reveals that most bacterial taxa in soil synthesize new ribosomal RNA. ISME J. 2018; 12 (12):3043-3045. https://doi.org/10.1038/s41396-018-0233-7 PMID: 30042501

88. Zhang Y, Cong J, Lu H, Deng Y, Liu X, Zhou J, et al. Soil bacterial endemism and potential functional redundancy in natural broadleaf forest along a latitudinal gradient. Sci Rep. 2016; 6(1):28819. https://doi.org/10.1038/srep28819 PMID: 27357005

89. Isobe K, Ise Y, Kato H, Oda T, Vincenot CE, Koba K, et al. Consequences of microbial diversity in forest nitrogen cycling: diverse ammonifiers and specialized ammonia oxidizers. ISME J. 2020; 14 (1):12-25. https://doi.org/10.1038/s41396-019-0500-2 PMID: 31481743

90. Moreno B, Cañizares R, Nuñez R, Benitez E. Genetic diversity of bacterial β-glucosidase-encoding genes as a function of soil management. Biol Fertil Soils. 2013; 49:735-745. https://doi.org/10.1007/s00374-012-0765-3

91. Li Y, Ge Y, Wang J, Shen C, Wang J, Liu YJ. Functional redundancy and specific taxa modulate the contribution of prokaryotic diversity and composition to multifunctionality. Mol Ecol. 2021; 30 (12):2915-2930. https://doi.org/10.1111/mec.15935 PMID: 33905157

92. Wertz S, Degrange V, Prosser JI, Poly F, Commeaux C, Guillaumaud N, et al. Decline of soil microbial diversity does not influence the resistance and resilience of key soil microbial functional groups following a model disturbance. Environ Microbiol. 2007; 9(9):2211-2219. https://doi.org/10.1111/j.1462-2920.2007.01335.x PMID: 17686019
93. He L, Rodrigues JLM, Soudzilovskaia NA, Barceló M, Olsson PA, Song C, et al. Global biogeography of fungal and bacterial biomass carbon in topsoil. Soil Biol Biochem. 2020; 151:108024. https://doi.org/10.1016/j.soilbio.2020.108024

94. Sinsabaugh RL. Phenol oxidase, peroxidase and organic matter dynamics of soil. Soil Biol Biochem. 2010; 42(3):391-404. https://doi.org/10.1016/j.soilbio.2009.10.014

95. Miller M, Palojärvi A, Rangger A, Reeslev M, Kjøller A. The use of fluorogenic substrates to measure fungal presence and activity in soil. Appl Environ Microbiol. 1998; 64(2):613-617. https://doi.org/10.1128/aem.64.2.613-617.1998 PMID: 9464399

96. Jiang L. Negative selection effects suppress relationships between bacterial diversity and ecosystem functioning. Ecology. 2007; 88(5):1075-1085. https://doi.org/10.1890/06-1556 PMID: 17536392

97. Bailey VL, Fansler SJ, Stegen JC, McCue LA. Linking microbial community structure to β-glucosidic function in soil aggregates. ISME J. 2013; 7(10):2044-2053. https://doi.org/10.1038/ismej.2013.87 PMID: 23719152

98. Ramin KI, Allison SD. Bacterial tradeoffs in growth rate and extracellular enzymes. Front Microbiol. 2019; 10:494333. https://doi.org/10.3389/fmicb.2019.02956 PMID: 31921094

99. Hester ER, Jetten MS, Welte CU, Lücker S. Metabolic overlap in environmentally diverse microbial communities. Front Genet. 2019; 10:476961. https://doi.org/10.3389/fgene.2019.00989 PMID: 31681424

100. Han S, Tan S, Wang A, Chen W, Huang Q. Bacterial rather than fungal diversity and community assembly drive soil multifunctionality in a subtropical forest ecosystem. Environ Microbiol Rep. 2022; 14(1):85-95. https://doi.org/10.1111/1758-2229.13033 PMID: 34962072

101. Osburn ED, Yang G, Rillig MC, Strickland MS. Evaluating the role of bacterial diversity in supporting soil ecosystem functions under anthropogenic stress. ISME Commun. 2023; 3(1):66. https://doi.org/10.1038/s43705-023-00273-1 PMID: 37400524

102. Wagg C, Bender SF, Widmer F, van der Heijden MG. Soil biodiversity and soil community composition determine ecosystem multifunctionality. Proc Natl Acad Sci USA. 2014; 111(14):5266-5270. https://doi.org/10.1073/pnas.1320054111 PMID: 24639507

103. De Mazancourt C, Isbell F, Larocque A, Berendse F, De Luca E, Grace JB, et al. Predicting ecosystem stability from community composition and biodiversity. Ecol Lett. 2013; 16(5):617-625. https://doi.org/10.1111/ele.12088 PMID: 23438189

104. Yachi S, Loreau M. Biodiversity and ecosystem productivity in a fluctuating environment: the insurance hypothesis. Proc Natl Acad Sci USA. 1999; 96(4):1463-1468. https://doi.org/10.1073/pnas.96.4.1463 PMID: 9990046

105. Cardinale BJ. Biodiversity improves water quality through niche partitioning. Nature. 2011; 472 (7341):86-89. https://doi.org/10.1038/nature09904 PMID: 21475199