Introduction
Environmental stresses due to climate change and human activity are
threatening to degrade the biodiversity and services of terrestrial
ecosystems around the world (Singh et al. , 2010, De & Ashley,
2013, Jansson & Hofmockel, 2020). An ecosystem contains a myriad of
components that interact and affect each other in complex ways. Thus,
the degradation of a few, or even just one component, could potentially
change the state of others and lead to holistic ecosystem degradation
(Moreno-Mateos et al. , 2020, Cheng et al. , 2021). Among
ecosystem components, the microbial community is unique. On the one
hand, its members have extremely high richness and reproduction rate,
which allows the microbial community to be responsive to stress as its
composition can adapt rapidly when facing environment fluctuations
(Jansson & Hofmockel, 2020, Hall et al. , 2018). On the other
hand, the microbial community plays an irreplaceable role in
decomposition, nutrient cycling, and regulation of above-ground
community diversity and activity, therefore, degradation of the
microbial community is a potential threat to related ecosystem functions
and other ecosystem components (Berendsen et al. , 2012).With
rising awareness of the importance of the microbial community and
increasingly mature methodology, research on ecosystem degradation
related to microbial communities under climatic and anthropogenic stress
have begun to receive increasing attention (Firestone, 2006, Hallet al. , 2018).
However, most microbial community
degradation related research has been performed according to vegetation
or soil degradation gradients, rather than according to microbial
community degradation itself (Zhang et al. , 2010, Cheng et
al. , 2021). These studies have also mainly focused on the effects that
vegetation or soil degradation exert on the microbial community, either
directly or indirectly, rather than the effects that microbial community
degradation exerts on other ecosystem components (Zhang et al. ,
2010, Cheng et al. , 2021).
Yet, due to the much higher
reproduction rate of microbes than plants, the microbial community could
have already passed through numerous generations and established new
stages of succession before significant degradation of vegetation or
soil is detected (Curtis, 2006,
Konopka, 2006, Berendsen et al. , 2012). Thus, it is more
plausible that microbial community degradation occurs earlier than
vegetation degradation, and that the latter is a consequence of
microbial community degradation (Curtis, 2006, Konopka, 2006, Berendsen
et al., 2012). However, studies exploring the effects of microbial
community degradation on other ecosystem components, such as
ectomycorrhizal species, remain scarce.
Ectomycorrhizal fungi are fungal microbes that live in the interface
between plant roots and soils, and are possibly the first step that soil
microbial community degradation exert to plant community (Martin et al.,
2016). The cloze interactions between ectomycorrhizal fungi species and
host plants had been comprehensively studied from the molecular level to
ecological level (Martin et al., 2016, Hoeksema et al., 2018). However,
the interactions between ectomycorrhizal fungi species and soil
microbial community were seldom reported. Theoretically, since the
ectomycorrhizal fungi physically interact with the soil microbial
community, the ectomycorrhizal fungi species should also interact at the
biochemical or ecological level (Bowen & Theodorou, 1979). In recent
years, climate change has resulted in the decreasing production and the
regional extinction of ectomycorrhizal fungi species in several
ecosystems even if the vegetation has not yet shown significant
degradation (Collado et al. , 2019, Ohara & Hamada, 1967, Gillet al. , 2000). Taking the potential influence from the soil
microbial community on ectomycorrhizal fungi into consideration, the
extinction of ectomycorrhizal fungi is a potential sign of soil
microbial community degradation, thereby linking the degradation of the
soil microbial community to other ecosystem components.
In order to elucidate the relationship of the extinction of
ectomycorrhizal fungi species with the degradation of a soil microbial
community, we chose Hengduan mountain as the sampling site. It is
located in the eastern edge of Tibet plateau and is a hotspot of global
biodiversity. Several types ectomycorrhizal fungi species, especiallyT.matsutake , have shown signs of regional extinction due to
climatic and anthropogenic stresses (Ruth et al. , 2008, Yinet al. , 2020). We sampled the top-soil along a T.matsutakeextinction gradient and investigated the variation in
the associated microbial
communities.
Materials and Methods
Sample collection and basic soil physico-chemical properties
Soil samples were collected from the Baima Snow Mountain National Park
in Yunnan Province, near Nanren village (99°6´40.65” E, 28°34´59.06” N,
at 3,534 m). The soils were sampled from an old alpine forest dominated
by Cyclobalanopsis glauca (Thunb.) Oerst on August 9, 2019
(Fig. S1A). With help from native experienced mushroom collectors, soils
beneath T. matsutake were collected and considered as undegraded
areas (UD). Degrading areas (DI) were classified as location whereT. matsutake that did not fruit that year, and degraded areas
(DD) were where no T. matsutake growth had been observed for over
three years. Over three sites were selected for each soil type and
divided into 10 subsamples around the fruit, marked as 1 to 10, and each
subsample was about 10cm away from the fruit (Fig. S1B). Each subsample
was collected from a depth of 0–20 cm. For each soil type, the
subsamples from each site of the same mark were mixed to form 10 samples
(Zhou et al. , 2016). Using
this methodology, we were able to minimize the influences and errors
caused by heterogeneity and focus primarily on the state of the
microbial community under the same degree
of T. matsutake extinction.
Soil pH was measured using a pH meter (FE20-FiveEasyTM pH,
MettlerToledo, Germany) after shaking a soil water (1:5 w/v) suspension
for 30min. Soil moisture was measured gravimetrically. Total organic
carbon (TOC) and total nitrogen (TN) contents were measured using an
elemental analyzer (VarioMAX, Elementar, Germany). Ammonium
(NH4+-N) and nitrate
(NO3−-N) were extracted at a ratio of
10g fresh soil to 100mL 2M KCl. After shaking for 1h,
NH4+-N,
NO3−-N contents in the filtered
extracts were analyzed using a continuous flow analytical system (San++
System, Skalar, Holland). Total phosphorus content (TP) was measured
using the alkali fusion–Mo-Sb Anti spectrophotometric method.
DNA extraction and amplicon sequencing
Total DNA was extracted according to the protocol for the
FastDNATM SPIN kit (MP Biomedicals), followed by
measuring the DNA concentrations and quality using a Nano-100 NanoDrop
spectrophotometer. Primers for the 16S rRNA gene amplification targeted
the V3-V4 hypervariable region and included 338F
5´-ACTCCTACGGGAGGCAGCA-3´ and 806R 5´-GGACTACHVGGGTWTCTAAT-3´. The
primers for ITS region amplification, 3F 5´-GCATCGATGAAGAACGCAGC-3´ and
4R 5´-TCCTCCGCTTATTGATATGC-3´, targeted the ITS2 genetic region. PCR
products were mixed in equimolar ratios according to the GeneTools
Analysis Software (Version 4.03.05.0, SynGene), and the pooled PCR
products were purified using an E.Z.N.A. Gel Extraction Kit (Omega,
USA). Sequencing libraries were generated using the NEBNext® Ultra™ II
DNA Library Prep Kit for Illumina® (New England Biolabs, MA, USA)
following manufacturer’s recommendations, and indexing adapters were
added. Library quality was assessed with a Qubit@ 2.0 fluorometer
(Thermo Fisher Scientific, MA, USA). The library was then sequenced on
an Illumina Nova6000 platform to generate 250 bp paired-end reads
(Guangdong Magigene Biotechnology Co., Ltd., Guangzhou, China).
The raw sequencing data generated in this study are available in the
NCBI Sequence Read Archive under the accessions PRJNA687558 (bacterial
undegraded stage-BUD library), PRJNA690166 (bacterial degrading
stage-BDI library), PRJNA690139 (bacterial degraded stage-BDD library),
PRJNA690207 (fungal undegraded stage-FUD library), PRJNA690181 (fungal
degrading stage-FDI library), and PRJNA690169 (fungal degraded stage-FDD
library).
Sequence data preprocessing and bioinformatics analyses
Fastp (version 0.14.1) was used for quality control of the raw sequence
data using a sliding window (-W4-M20). Primers were also removed using
the cutadapt software program. Paired-end clean reads were then merged
using the usearch -fastq_mergepairs command based on overlap between
paired-end reads. Merging criteria included requiring at least 16 bp of
overlap from reads generated from the opposite end of the same DNA
fragment with a maximum of 5 mismatches being allowed in the overlap
region. The merged sequences were identified as raw sequences. Fastp
(version 0.14.1) was then used again to perform quality control on the
raw sequence data using a sliding window (-W4-M20) to obtain paired-end
clean sequences. Sequences with more than 97% nucleotide similarity
were clustered into operational taxonomic units (OTUs).
Bioinformatic analyses
Network analyses were conducted according to the molecular ecological
network analysis pipeline (MENA,
http://ieg2.ou.edu/MENA/) protocols.
A cutoff threshold of 0.88–0.89 was used, and threshold values ranging
from 0.01 to 0.99 over 0.01 intervals were applied to the Spearman
correlation matrix (Deng & Zhou, 2015). Module hubs and connectors were
determined based on the criteria of Zi > 2.5 and Pi
> 0.625, respectively. The constructed network and
sub-networks were visualized using Cytoscape 3.3.0 (Xie et al. ,
2020). Statistically significant differences in basic soil properties
and diversity indices among the three soil types were measured using
one-way ANOVA tests. In addition, differences in the relative abundances
of OTUs and functional gene modules were evaluated using the
Kruskal-Wallis method. Calculations of richness, Shannon index, and
equitability, in addition to Linear discriminant analysis Effect Size
(LEfSe) tests, principal coordinates analysis (PCoA) based on
Bray-Curtis distances, and 16S rRNA gene-based functional predictions
(by comparison to the database COG-cluster of orthologous groups), were
all conducted with the online platform
http://www.magichand.online/h5-BioCloud-site/#/. Factor
extractions were conducted using a correlation matrix in the SPSS
software suite.
Results
Variation in basic soil abiotic and biotic properties during
degradation
The degradation of T. matsutake habitat was accompanied by
significant variation in soil physico-chemical and microbial community
properties. The contents of TC, TN, TP, ammonia, nitrate, and pH
significantly increased with degradation status (Table S1). The most
prominent increase was for ammonia content that was nearly 42 times
higher in DD soils than in UD soils, while TC, TN, and TP increased by
only approximately 1.0-, 2.0-, and 0.5-fold, respectively. TN and
nitrate contents did not significantly differ between UD and DI soils.
The Shannon indices of the bacterial and fungal communities from UD
soils were 2.28 and 0.97, respectively, but increased to 2.77 and 1.60
in DI soils. The Shannon index values were lowest in DD soils, at 2.03
and 0.72 for bacterial and fungal communities, respectively. Similar
trends were also observed for richness and equitability metrics (Fig.
S2). Variations in the relative abundances (RA) of fungal and bacterial
populations were observed among communities (Fig. S3A and B). The RA of
α-Proteobacteria (the bacterial phylum with the highest RA)
increased from 21.9% (UD) to 31.6% (DI), and then decreased to 13.6%
(DD) in the respective soil communities. Acidobacteria RAs
significantly increased to as high as 62.7% in DD soils. The RAs of
β-Proteobacteria decreased from 17.7% in UD soils to 1.1% in DD
soils, while the RAs of Actinobacteria also decreased from 16.9%
in UD soils to 4.4% in DD soils.
The RAs of “other” fungi decreased from 28.1% in UD soils to 6.7% in
DD soils. Furthermore, the RAs of Eurotiales ,Thelephorales, and Hypocreales exhibited similar declining
trends. Helotiales RAs slightly increased from 4.2% in UD soils
to 7.5% in DD soils. The dominant Agaricales (41.5% RA) of UD
soils was replaced by Atheliales (73.3% RA in DD soils).Sebacinales and Russulales exhibited the highest RAs in DI
soils (18.8% and 33.6%, respectively) which were significantly higher
(p < 0.05) than in UD and DD soils. The RA ofT.matsutake in the fungal community decreased from 41% in UD
soils to 1.6% in DD soils, highlighting the decreasing dominance ofT.matsutake during the degradation process (Fig. S3C). LEfSe
analysis indicated that T.matsutake was the most enriched
representative species in UD soils, with decreasing prominence in
degraded soils (Fig. S3D). Only the RAs of Piloderma olivaceumwere significantly enriched in DD soils compared to the other two
degradation gradients.
The microbial community in DD exhibited a greater number of predicted
gene modules related to stress resistance and basic metabolism than in
DI and UD, but also exhibited fewer predicted modules related to
interaction, reproduction, and cofactor metabolism. In the T.
matsutake degraded habitats, the RAs of gene modules related to basic
metabolites (i.e., carbohydrates, lipids, and amino acids) were
enhanced, while carbohydrate metabolite-related gene modules exhibited
the largest increased abundances, increasing to approximately 25%
higher in the DD soils compared to UD soils. The predicted RAs of
defense and motility related gene modules were also higher in DD (Fig.
1). In contrast, the predicted RAs of replication-related gene modules
decreased from UD to DD, the largest decrement was chromatin structure
and dynamics-related gene modules, with the largest decreases associated
with chromatin structure and dynamics-related gene modules the RA in DD
soils were only 77.5% of levels in UD soils. The RAs of community
interaction and cofactor metabolite-related gene modules were also
predicted as lower in DD (Fig.1).
Microbial community network variation during degradation and
interactions with T.matsutake
Microbial networks become less interactive along theT.matsutake degradation
gradient, although the largest network was observed for DI communities.
The network topological parameters (Table S2) demonstrated that DI soil
networks exhibited the most nodes (1,280) and interactions (1,414),
although it is not the most interactive network among the three soil
types. Across the degradation gradient, the average degree (avgK)
decreased from 3.54 (UD) to 1.63 (DD), indicating a lower interaction
density (Deng & Zhou, 2015). Furthermore, the average clustering
coefficient (avgCC) decreased from 0.192 (UD) to 0.112 (DD), while the
centralization of eigenvector centrality (CE) increased from 0.231 (UD)
to 0.572 (DD). These observations indicate that the size proportion of
unit modules in the network decreased across the degradation gradient,
with an increase in the number of smaller modules (Deng & Zhou, 2015).
The transitivity (trans) decreased from 0.506 (UD) to 0.241 (DD) across
the gradient, reflecting a potentially decreased ability to communicate
and cooperate among network members (Deng & Zhou, 2015). In addition, a
mantel test indicated the presence of significant differences between
the empirical network and random networks (p < 0.05),
suggesting that the empirical networks were reasonable.
The form and state of interactions in the network varied during
degradation (Fig. 2). The UD network contained four module hubs (BOTU79,
BOTU257, BOTU99 and FOTU7) and two bacterial connectors (BOTU175 and
BOTU2394), while the DI network contained twelve module hubs and two
fungal connectors (FOTU1608 and FOTU1160), and the DD network had no
connectors and only a single module hub (BOTU37). The ratio of module
hubs/connectors increased from UD to DD, indicating that the interaction
types shifted from an inter-module focused state to an intra-module
focused state (Fig. 2A and Fig. 2B).
The ratios of total positive to negative interactions (p/n) in the UD,
DI, and DD networks were nearly equivalent, while the DI network had a
slightly higher ratio of negative interactions than the UD and DD
networks (Fig. 3A). However, the p/n of bacterial-fungal interactions
decreased from 1.18 (UD) to 0.82 (DD), while the p/n of
bacterial-bacterial interactions increased from 0.89 (UD) to 1.08 (DD).
The p/n of fungal-fungal interactions was lowest in the DI network
(1.07) compared to 1.72 and 1.85 in the UD and DD networks, respectively
(Fig. 3A). This indicated interactions between bacterial and fungi
transitioned from a more cooperative to a more competitive state during
degradation. The proportion of bacterial-fungal interactions to total
interactions was similar among the three networks, but the proportion of
fungal-fungal interactions decreased from 28.18% in the UD network to
24.31% in the DD network, while the proportion of bacterial-bacterial
interactions increased from 26.01% in the UD network to 31.33% in the
DD network (Fig. 3B). This observation indicates the potential
decreasing relative importance of fungal interactions compared to
bacterial interactions in the network along the degradation gradient.
Furthermore, the RAs of OTUs that were directly associated withT.matsutake encompassed 13 negatively interacting OTUs that had
significantly lower RAs in the DD communities compared to the UD
communities (Fig. 3C). However, half of the positively associated OTUs
(80% were bacterial) exhibited increased RAs along the degradation
gradient, indicating that they improved their fitness
once they no longer had
interactions with T.matsutake . Only BOTU43 and FOTU304 exhibited
decreased RAs with the loss of T.matsutake dominance, indicating
that they may have intimate and irreplaceable interactions withT.matsutake (Fig. 3C). Thus, the location of T.matsutakein the networks became marginalized as degradation increased (Fig. 2B),
moving from a large central module in the UD network that had tight
interactions, to a small peripheral module that had few interactions
with other nodes and modules in the DD network. This change indicated
the loss of functional interactions with T.matsutake during
habitat degradation.
The contribution of microbial community interactivity to
T.matsutake extinction
A strong linear relationship was observed between T.matsutake RAs
and network structures, with R2 values as high as
0.909 (p < 0.05), further implicating the influence of
network interactions on T.matsutake RAs (Fig. 4). However,
correlations of soil properties in addition to the diversity and
composition metrics for the bacterial and fungal communities exhibited
weak linear relationships to T.matsutake RAs. The lowest
R2 was observed for the relationship between community
structure and T.matsutake RAs (0.006). Furthermore, a weak
positive relationship was observed between community diversity andT.matsutake RAs (R2=0.302), while a weak
negative relationship was observed between soil properties andT.matsutake RAs (R2=0.417).
Discussion