\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{Winter warming rapidly increases carbon degradation capacities of fungal
communities in tundra soil: potential consequences on carbon stability}



\author[1]{Jingmin Cheng}%
\author[1]{Yunfeng Yang}%
\author[2]{Mengting Yuan}%
\author[3]{Ziyan Qin}%
\author[1]{Qun Gao}%
\author[2]{Liyou Wu}%
\author[2]{Zhou Shi}%
\author[4]{Edward Schuur}%
\author[5]{James R Cole}%
\author[5]{James Tiedje}%
\author[1]{Jizhong Zhou}%
\affil[1]{Tsinghua University}%
\affil[2]{University of Oklahoma}%
\affil[3]{Affiliation not available}%
\affil[4]{Northern Arizona University}%
\affil[5]{Michigan State University}%


\vspace{-1em}



  
  \date{May 26, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
High-latitude tundra ecosystems are increasingly affected by climate
warming. As an important fraction of soil microorganisms, fungi play
essential roles in carbon (C) degradation, especially the old,
chemically recalcitrant C. However, it remains obscure how fungi respond
to climate warming and whether fungi, in turn, affect C stability of
tundra. In a two-year winter soil warming experiment of 2 \selectlanguage{ngerman}°C by snow
fences, we investigated responses of fungal communities to warming in
the active layer of the Alaskan tundra. Although fungal community
composition, revealed by 28S rRNA gene amplicon sequencing, remained
unchanged (P \textgreater{} 0.05), fungal functional gene composition,
revealed by a microarray named GeoChip, was altered (P \textless{}
0.05). Changes in functional gene composition were linked to winter soil
temperature, thaw depth, soil moisture, and gross primary productivity
(Canonical Correlation Analysis, P \textless{} 0.05). Specifically,
relative abundances of fungal genes encoding invertase, xylose
reductase, and vanillin dehydrogenase significantly increased (P
\textless{} 0.05), indicating higher C degradation capacities of fungal
communities under warming. Accordingly, we detected changes of fungal
gene networks under warming, including higher average path distance,
lower average clustering coefficient, and lower percentage of negative
links, indicating that warming potentially changed fungal interactions.
Together, our study revealed higher C degradation capacities of fungal
communities under short-term warming and highlights the potential
impacts of fungal communities on mediating tundra ecosystem respiration,
and consequently future C stability of high-latitude tundra.%
\end{abstract}\selectlanguage{ngerman}%



\sloppy


\textbf{Introduction}

High-latitude tundra is one of the most vulnerable ecosystems to climate
warming (Myers-Smith et al., 2015), wherein permafrost thawing and
changes of vegetation patterns have been observed during past decades
(Schuur, Crummer, Vogel, \& Mack, 2007; Sturm, Racine, \& Tape, 2001;
Xue et al., 2016). Given the fact that permafrost soil of high-latitude
tundra contains up to half of the global soil carbon (C) storage
(Tarnocai et al., 2009), warming-induced permafrost C loss to atmosphere
could be one of the most severe concerns for accelerating climate
warming by greenhouse gas effects (Natali et al. 2011; Schuur et al.
2009). In a warmer world, vast C may be released to the atmosphere as a
result of stimulated microbial activities to accelerate permafrost C
degradation (Mann et al., 2015; Zhou et al., 2012).

As the major soil C decomposers, fungi may affect the C stability of
high-latitude tundra considerably (Christiansen et al., 2016). For
instance, saprotrophic fungi are important regulators of soil nutrient
cycling via organic C decomposition, especially lignin (Buckeridge,
Banerjee, Siciliano, \& Grogan, 2013; Yuste et al., 2011). In winter,
fungi are the largest proportion of soil microbial biomass and
predominantly affect wintertime C degradation (Buckeridge et al., 2013;
Schadt, Martin, Lipson, \& Schmidt, 2003). Additionally, mycorrhizal
fungi form mutually beneficial associations to facilitate the growth of
arctic plants, especially those warming-responsive plant species such
as\emph{Salix} spp., and \emph{Betula nana} (Clemmensen, Michelsen,
Jonasson, \& Shaver, 2006), and thus can regulate C stability by
affecting primary productivity. Therefore, responses of fungal
communities to climate warming would be determinative for the C fate in
tundra ecosystems.

Using PCR-based marker gene sequencing, several \emph{in situ} field
studies have demonstrated that experimental warming changes taxonomic
composition and environmental group abundances of fungal communities in
tundra soil, with increased diversity of saprotrophic fungi but
decreased diversity of ectomycorrhizal fungi (Morgado et al., 2016;
Mundra et al., 2016; Semenova et al., 2016). However, such taxonomic
information provides little insight into the relationship between fungal
functions and ecological processes. Functional gene analysis has been
applied to quantify functional genes responsible for key ecosystem
processes. The main findings of previous studies (Levy-Booth, Prescott,
\& Grayston, 2014; Yergeau, Kang, He, Zhou, \& Kowalchuk, 2007), i.e.,
close relationships between functional gene abundances and ecosystem
processes, have provided important insights on the roles of microbial
communities in ecosystem functioning.

In addition to taxonomic composition, it is possible that fungal
functional genes, serving as markers of potential functional capacities,
are sensitive to climate warming. The justification includes 1)
Temperature affects the syntheses of cryoprotective compounds essential
for maintaining fungal growth and functions in cold environments
(Weinstein, Montiel, \& Johnstone, 2000); 2) Temperature sensitivity of
soil heterotrophic respiration (measured as Q\textsubscript{10}) is
widespread (Malcolm, Lopez-Gutierrez, Koide, \& Eissenstat, 2008),
implying that C degradation capacities of soil microorganisms including
fungi are variable; and 3) Plant litter input would be higher under
warming (Rinnan, Michelsen, Baath, \& Jonasson, 2007), which may
stimulate fungi-mediated C degradation.

As the change of functional gene composition can inevitably affect
interactions among microbial members, it is important to elucidate
complex interactions among microbial functional genes or species
containing these genes in order to understand ecosystem functioning (Wu
et al., 2016). However, it remains challenging due to myriad members in
microbial communities and difficulties in cultivating individual
species. In recent years, network analysis has been widely used to
reveal potential interactions among co-existing microbial members
(Banerjee et al., 2016; Banerjee, Schlaeppi, \& van der Heijden, 2018;
de Vries et al., 2018). Accumulating evidence shows that network
properties can reflect responses of microbial communities to
environmental disturbances (de Vries et al., 2018; De Vries et al.,
2012), as drought treatment in grassland mesocosms led to higher
connectivity and lower modularity of soil bacterial networks than fungal
networks, which suggested lower stability of bacteria than fungi under
drought (de Vries et al., 2018).

In the present work, we analyzed the effects of winter warming on fungal
communities in the active layer of the Alaskan tundra, via 28S gene
amplicon sequencing and GeoChip technologies. Through installing snow
fences to increase snow cover thus winter soil temperature, we tested
the following hypotheses: 1) Warming would affect fungal community
composition and potential interactions; 2) Fungal C degradation
capacities, represented by the abundances of C degradation genes, would
be increased by warming; 3) Increased fungal C degradation capacities
could explain the increase of ecosystem respiration (\emph{R}
\textsubscript{eco}), as recently witnessed in several studies (Liu et
al., 2015; Trivedi et al., 2016; Zhao et al., 2014).

\textbf{Methods \& Materials}

\emph{Site description and experimental design}

This study was carried out in Carbon in Permafrost Experimental Heating
Research (CiPEHR) site, established in September 2008. This site is a
typical moist acidic tundra located in the northern foothills of
Interior Alaska (63\selectlanguage{ngerman}°52'59''N, 149°13'32''W), near the Denali National
Park and Preserve of Eight Mile Lake. The annual mean temperature of the
year 1976-2009 was -1.0 °C, with monthly mean temperature ranging from
-16 °C (December) to 15 °C (July). The annual mean precipitation of
1976-2009 was 378 mm. The soil type is gelisol, with a 45-65 cm thick
organic horizon covered by a mixture of glacial till and loess. The
active layer thickness is 50-60 cm. Dominant plant species included the
tussock-forming sedge \emph{Eriophorum vaginatum} (Cyperaceae) and the
small deciduous shrub \emph{Vaccinium uliginosum} (Ericaceae). Other
plants included \emph{Carex bigelowii}(Cyperaceae), bryophytes
\emph{Pleurozium schreberi} (Hylocomiaceae),\emph{Sphagnum} spp.
(Sphagnaceae), and lichens \emph{Cladonia} spp. (Cladoniaceae).

The experimental soil warming study was initiated in October 2008. A
total of 6 snow fences (1.5 m tall × 8 m long) were established and
arranged into 3 blocks, which was \textasciitilde{} 100 m apart to
mitigate small-scale spatial variability. For each snow fence, the
warmed plot was located at the leeward side with 5 m distance from the
fence, and the paired control plot was located at the windward side with
8-14 m distance from the fence. Snow fences increased soil temperature
via trapping insulating layers of snow in the leeward side. For example,
snow fences led to a threefold increase in snow depth of warmed plots in
2009 (1.35 m in warmed plots versus 0.38 m in control plots) (Natali et
al., 2011). Different from traditional snow fence experiments, both snow
and snow fences were removed before snowmelt to ensure that soil
hydrology, light conditions, and plant phenology were not affected by
extra snow. On March 8-15 of 2010, extra snow in warmed plots was
shoveled to ensure comparable snow cover depth between warmed and
control plots. Snow fences were removed at the time of snow removal to
avoid shading of plots during the growing season (Natali, Eag, \& Rubin,
2012; Natali et al., 2011). Snow removal allowed this experiment to
focus on the effects of soil warming.

\emph{Plant survey and ecosystem C flow measurements}

Plants were identified at the species level and classified as shrubs,
graminoids, bryophytes, and lichens. Plant biomass and net primary
productivity (NPP) of each plant group were measured by a point-frame
method (Walker, 1996), based on previously developed allometric
equations (Schuur et al., 2007).

Soil temperature at 15-25 cm depth in winter and soil moisture in the
growing season were measured \emph{in situ} every 30 minutes using
constantan-copper thermocouples with CR1000 data loggers (Campbell
Scientific, Logan, UT, USA). Thaw depth, which is the thickness of
unfrozen soil in the growing season, was measured weekly using a metal
depth probe. In the growing season, net ecosystem exchange (NEE) and
ecosystem respiration (\emph{R} \textsubscript{eco}) were measured using
an automated CO\textsubscript{2} flux system. Gross primary productivity
(GPP) was calculated as the sum of NEE and \emph{R}
\textsubscript{eco}of the growing season. Winter\emph{R}
\textsubscript{eco} was modeled from the measured CO\textsubscript{2}
flux data of March-April 2009 using an infrared gas analyzer, following
detailed description of the modeling methodology shown in the Appendix
of a previous paper (Natali, Schuur, Webb, Pries, \& Crummer, 2014). In
brief, an exponential equation between winter\emph{R}
\textsubscript{eco} and soil temperature was developed, i.e., winter
\emph{R} \textsubscript{eco} =\emph{R} ×exp\selectlanguage{greek}\textsuperscript{(β×Ts)}\selectlanguage{ngerman},
where \emph{R} is basal respiration, \selectlanguage{greek}β \selectlanguage{english}is the rate of respiration
change, and Ts is soil temperature at 20 cm depth. For CiPEHR site,
parameter values estimated for this model were \emph{R} = 0.47 (0.03)
and \selectlanguage{greek}β \selectlanguage{english}= 0.18 (0.06). This model was validated by comparing the
simulation results with all observed \emph{R} \textsubscript{eco} values
available from flux sites (\emph{R} \textsuperscript{2} = 0.70, \emph{P}
\textless{} 0.001). Since no photosynthesis occurs in winter, winter
\emph{R} \textsubscript{eco}mainly represents heterotrophic respiration.
\emph{R} \textsubscript{eco}in this study is annual \emph{R}
\textsubscript{eco} that is a sum of winter \emph{R} \textsubscript{eco}
and growing season\emph{R} \textsubscript{eco}. The period ranges of GPP
and NEE data were May 2010 - September 2010. The period ranges of other
environmental data, including\emph{R} \textsubscript{eco}, soil
temperature, soil moisture, and thaw depth, were October 2009-September
2010.

\emph{Soil sampling and biogeochemical measurements}

Soil samples were collected in early May of 2010. To dismiss within-plot
heterogeneity, we collected a relatively large soil core (\textgreater{}
2 Kg) from each plot. A total of 12 soil cores at 15-25 cm depth from
warmed plots and their paired control plots were collected using
electric drills and were transported to the laboratory on ice. In the
lab, we divided each sample into two parts using a sterile saw: one for
soil DNA extraction (stored at -80 \selectlanguage{ngerman}°C); the other for soil
biogeochemical measurements (stored at 4 °C), including bulk density,
labile carbon pool (LCP, mainly polysaccharides and cellulose),
recalcitrant carbon pool (RCP, mainly lignin), total organic carbon
(TOC), soil nitrogen, NH\textsubscript{4}\textsuperscript{+}-N, and
NO\textsubscript{3}\textsuperscript{-}-N. Soil bulk density was
calculated as its dry mass (oven-drying at 60\selectlanguage{english} degC for three\selectlanguage{english} days) divided
by its frozen volume. A two-step hydrolysis procedure separated LCP and
RCP. LCP was quantified by a Shimadzu TOC-V CPH PC-Controlled TOC
analyzer (Shimadzu Corporation, Tokyo, Japan). RCP was quantified by a
PerkinElmer Optima 2000DV ICP-OES spectrometer (PerkinElmer Inc.,
Waltham, MA, USA) (Rovira \& Vallejo, 2002).
NH\textsubscript{4}\textsuperscript{+}-N and
NO\textsubscript{3}\textsuperscript{-}-N were measured using anion (AG
1-X8, 20-50 mesh) and cation (AG 50W-X8, 20-50 mesh) exchange bags
(Biorad, Hercules, CA, USA). In each sampling plot, a cation bag and an
anion bag were placed at a depth of 10 cm during May-September 2010.
Then samples in the bags were collected, and inorganic nitrogen
components were extracted in the 0.1M HCl/2.0M NaCl solution for an hour
in the lab and analyzed on an Astoria Analyzer (Astoria-Pacific Int.,
Clackamas, OR, USA).

\emph{Soil DNA extraction}

DNA was extracted from 5 g of soil samples stored at -80 degC, following
the instruction of PowerMax Soil DNA Isolation Kit (MoBio Laboratories,
Carlsbad, CA, USA). DNA quality was assessed via a NanoDrop ND-1000
Spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE, USA),
based on UV spectrometry absorbance ratios at wavelengths of 260 nm/280
nm and 260 nm/230 nm. DNA concentration was measured by PicoGreen using
a FLUOstar OPTIMA fluorescence plate reader (BMG Labtech, Jena,
Germany).

\emph{Pyrosequencing of 28S rRNA gene amplicons}

28S rRNA gene was applied to investigate fungal taxonomy in this study.
After combining with a single sample-specific 8-mer tag sequence
separated by a 2-bp spacer sequence, the primers LR0R
(5'-ACCCGCTGAACTTAAGC-3') and LR3 (5'-CCGTGTTTCAAGACGGG-3') were applied
to amplify the D1-D2 region of 28S rRNA gene, which generated a gene
amplicon of about 650 bp. The 100 \selectlanguage{greek}µ\selectlanguage{english}L PCR system was composed of 50 mM
KCl, 10 mM Tris-HCl at pH of 8.3, 0.1\% Triton X-100, 2.5 mM
MgCl\textsubscript{2}, 200 \selectlanguage{greek}µ\selectlanguage{english}M dNTPs (Invitrogen Life Technologies,
Carlsbad, CA, USA), 5 units of \emph{Taq}polymerase, 0.1 \selectlanguage{greek}µ\selectlanguage{english}g BSA (New
England Biolabs Inc., Ipswich, MA, USA), 0.4 \selectlanguage{greek}µ\selectlanguage{english}M each primer and 12.5 ng
DNA template. PCR thermocycling conditions were as follows: 3 min at 94
\selectlanguage{ngerman}°C for initialization, followed by 35 cycles of 1 min at 94 °C for
denaturation, 45 secs at 53 °C for annealing, 1 min at 72 °C for
extension, and finally 15 min at 72 °C before maintained at 4 °C. For
each sample, five to ten PCR reactions were performed. PCR products were
combined and purified by agarose gel electrophoresis, and then collected
with QIAquick Gel Extraction Kit (Qiagen, Venlo, Netherlands). PCR
products were quantified by PicoGreen using a FLUOstar Optima (BMG
Labtech, Jena, Germany). Equimolar uniquely-tagged amplicons from each
sample were then combined into one library. Tags separated raw sequences
via the function ``Split by tag'' of the IEG Galaxy pipeline
(http://zhoulab5.rccc.ou.edu:8080). They were further cleaned via LUCY
based on the following thresholds: maximum average error at 0.01,
maximum error at the ends at 0.01, minimum score at 21, window size at
50, minimum length at 200, and allowable ambiguous base number (`N') at
10. Only sequences with no `N' were accepted and were clustered into
operational taxonomic units (OTUs) using UCLUST at the 97\% identity.
OTU singletons were removed and OTUs annotation was performed by the RDP
Classifier with the latest Fungal LSU training set 11
(http://rdp.cme.msu.edu/classifier/classifier.jsp) at the 50\%
confidence interval. We used raw OTU abundance in the following data
analyses.

\emph{GeoChip 4.2 hybridization}

Fungal functional genes were accessed by GeoChip 4.2 (Roche NimbleGen
Inc., Madison, WI, USA), which is a functional microarray containing
more than 107,950 oligonucleotide probes of 792 bacterial, archaeal and
fungal genes associated with C, nitrogen (N), phosphorus (P) and sulfur
(S) cycling, stress, virulence, energy process, to name a few\emph{.}
These oligonucleotide probes are designed from functional gene sequences
of microbial organisms with known taxonomy in the NCBI database (He et
al., 2007). Before hybridization, 1 \selectlanguage{greek}μ\selectlanguage{english}g of purified DNA was labeled with
the fluorescent dye Cy3. After overnight hybridization, unbound DNA was
removed. Then GeoChip was scanned by a NimbleGen MS 200 Microarray
Scanner (Roche, Basel, Switzerland) at 532 nm and 635 nm. Probes with
signal-to-noise ratio lower than 2 were considered to be unreliable and
hence removed. The normalized signal intensity of each detected probe
was calculated by dividing total signal intensity by each detected probe
and log-transformed.

\emph{Network construction}

To detect warming effects on potential fungal interactions, we used
functional gene data to construct functional molecule ecological
networks (fMENs) for control and warmed samples separately, via the
pipeline of Molecular Ecological Network Analyses (MENA)
(http://ieg4.rccc.ou.edu/mena/) (Deng et al., 2012). Only functional
genes detected in all six biological replicates were included to ensure
network reliability. To compare topological characteristics between
control and warmed networks, we chose the same threshold of 0.98 during
network construction, which refers to the Spearman's correlation
coefficients between pairwise OTUs. Network modules, which are groups of
nodes that highly connected within the group but less connected outside
the group, were separated by the fast-greedy modularity optimization.
Network topological characteristics were calculated via MENA, including
total nodes, total edges, average connectivity (node degree), average
clustering coefficient, average path distance, and modularity. Network
randomization (a total of 100 random networks) was performed via MENA.
To test whether fMENs of control and warmed samples have typical
properties of biological networks (e.g., small-world, hierarchy,
modularity), network topological characteristics were compared with
their corresponding 100 random networks. To test statistical differences
of fMENs between control and warmed samples, we performed two-tailed
Student's \emph{t} -tests by taking their network topological
characteristics as averages and standard deviations of their
corresponding 100 random networks as standard deviations (Wu et al.,
2016). Cytoscape (version 3.5.1) was used for network visualization.
Some network scores (e.g., connectivity) were used to identify putative
key genes/species that have large influences in maintaining community
structure and function, despite that species identities can be
inaccurate owing to horizontal gene transfer. To select network key
genes, we set thresholds based on within-module connectivity (Zi) and
among-module connectivity (Pi) of nodes: Zi [?] 2.5 and Pi [?] 0.62 for
network hubs, Zi [?] 2.5 and Pi \textless{} 0.62 for module hubs, Zi
\textless{} 2.5 and Pi [?] 0.62 for connectors, Zi \textless{} 2.5 and Pi
\textless{} 0.62 for peripherals (Guimera \& Amaral, 2005). Previous
studies have shown that network key genes/species were inclined to be
more responsive to environmental changes than others (Wu et al., 2016).
Therefore, we checked whether the connectivity of network key genes was
linked to warming-induced changes in environmental factors by performing
partial Mantel tests between the matrix of connectivity of network key
genes and the matrix of environmental factors.

\emph{Statistical analyses}

Statistical analyses were carried out via R software (version 3.4.2).
Diversity indices (i.e., Shannon-Wiener index and Richness) and Plelou's
evenness of functional genes and OTUs were calculated. Two-tailed
Student's\emph{t} -tests on diversity indices, relative abundances of
each functional gene and gene category, and relative abundance of each
OTU and higher taxonomic levels were performed between warmed and
control samples. \emph{P} -values were adjusted by the Benjamini \&
Hochberg (BH) correction with the false discovery rate of 0.10. The
multi-response permutation procedure (MRPP) and analysis of similarity
(ANOSIM) were performed via \emph{mrpp} and \emph{anosim} functions in
\emph{vegan}package in R to examine overall differences in functional
gene composition or OTU composition between warmed and control samples.
If significant differences were found, then the similarity percentage
(SIMPER) analysis, via the \emph{simper} function in the
\emph{vegan}package in R, was applied to reveal which functional genes
or OTUs caused the differences as well as their contributions. The
principal coordinate analysis (PCoA) was performed to visualize the
overall composition of functional genes and OTUs. The canonical
correlation analysis (CCA) and the variation partitioning analysis
(VPA), via the \emph{cca} function in \emph{vegan} package in R, were
performed to reveal which environmental factors were closely linked to
differences between control and warmed samples and their contributions.
We only selected environmental factors with the variance inflation
factor (VIF) \textless{} 20 in the CCA model and tested the overall
significance of the CCA model by ANOVA. Relationships between\emph{R}
\textsubscript{eco} and relative abundances of fungal functional genes
for C degradation were examined by linear regression.

\textbf{Results}

\emph{Environmental factors}

Snow fences and resulting thicker snow cover significantly increased
winter soil temperature from -5.6\selectlanguage{ngerman}±0.6 °C to -3.7±0.9 °C (Table S1,
\emph{P}\textless{} 0.01), and exerted legacy effects on soil and plants
in the subsequent growing season. For instance, soil moisture
significantly increased from 47.8±4.9\% to 52.7±2.0\% (\emph{P} = 0.03),
and thaw depth significantly increased from 35.3±0.8 cm to 38.4±2.6 cm
(\emph{P}= 0.03). Gross primary productivity (GPP) significantly
increased from 394.5±7.0 g m\textsuperscript{-2} to 514.7±7.4 g
m\textsuperscript{-2}(\emph{P} = 0.01), mainly attributed to the
increase in graminoid biomass (146.7±38.7 g m\textsuperscript{-2} in
control plots vs. 231.0±86.5 g m\textsuperscript{-2} in warmed plots,
\emph{P} = 0.05). As high plant biomass usually leads to more plant
respiration, \emph{R} \textsubscript{eco}, the sum of aboveground plant
and soil respiration, was significantly enhanced from 363.3±3.7 g
m\textsuperscript{-2} to 449.2±4.4 g m\textsuperscript{-2}(\emph{P}
\textless{} 0.01). There were no significant changes of soil C and N
pools, including LCP (\emph{P} = 0.94), RCP (\emph{P} = 0.94),
NH\textsubscript{4}\textsuperscript{+}-N (\emph{P} = 0.47), and
NO\textsubscript{3}\textsuperscript{-}-N (\emph{P} = 0.59).

\emph{Fungal taxonomic composition}

The C1 and C2 samples of control plots were not saturated in sequence
depths (less than 10,000 sequences, based on the rarefaction curves in
Fig. S1), probably due to technical problems of DNA sequencing.
Therefore, they were discarded for data analyses. After resampling of
raw 28S rRNA gene reads, a total of 126,050 high-quality sequence reads
were generated and were clustered into 991 OTUs at the 97\%
identity.\emph{Ascomycota} (70.66±11.80\%) was the dominant phylum,
followed by \emph{Basidiomycota}(14.21±10.19\%), \emph{Chytridiomycota}
(0.51±0.64\%) and\emph{Blastocladiomycota} (\textless{} 0.01\%). Those
OTUs were assigned to 135 genera and the most abundant genera included
\emph{Kriegeria}(9.14±9.64\%), \emph{Hyaloscypha} (6.63±5.41\%),
\emph{Collophora}(3.43±6.50\%), \emph{Mollisia} (2.20±3.68\%),
\emph{Russula}(0.80±3.40\%), \emph{Sorocybe} (0.57±1.77\%),
\emph{Rigidoporus}(0.46±1.53\%), \emph{Chloroscypha} (0.45±0.78\%),
\emph{Hymenoscyphus}(0.41±0.61\%) and \emph{Crinula} (0.40±1.00\%). The
fungal genera were further classified into 6 ecological groups based on
their roles in ecosystems (Tedersoo et al., 2014) - saprotrophs (57\%),
plant pathogens (24\%), ectomycorrhizal fungi (6\%), animal parasites
(5\%), mycoparasites (2\%) and mycobionts (6\%), among a total of the
135 genera. However, soil warming did not change the overall fungal
community composition, reputing part of our first hypothesis (Table 1 \&
Fig. S2a). Additionally, no differences in diversity and environmental
group abundances were found between warmed and control samples (Table
S3).

\emph{Fungal functional gene composition}

A total of 2,590 probes of 82 fungal functional genes responsible for
nutrient cycling, stress, energy processes, and virulence were detected.
In line with our first hypothesis, we found that warming altered the
overall functional gene composition significantly (Table 1,
\emph{P}\textless{} 0.05), which was further verified by the PCoA plot
(Fig. S2b). Therefore, we performed the SIMPER analysis to identify
major functional genes responsible for the differences. A total of 12
functional genes had significant contributions (Table 2,
\emph{P}\textless{} 0.05), which were associated with C degradation
(i.e., genes encoding cellobiase, xylanase, glyoxal oxidase,
metalloprotease, phospholipase A2, vanillin dehydrogenase, xylose
reductase, invertase), S cycling (i.e., genes encoding sulfate
transporter and PAPS reductase), P utilization (i.e., the gene encoding
phytase) and stress tolerance (i.e., the gene encoding superoxide
dismutase).

In line with our second hypothesis, 3 out of a total of 53 fungal C
degradation genes showed significantly (\emph{P} \textless{} 0.05)
higher relative abundances under warming , which encode invertase for
sucrose degradation (increased by 16.7\%), xylose reductase for
hemicellulose degradation (increased by 231.62\%), and vanillin
dehydrogenase for lignin-derived aromatics degradation (increased by
47.64\%) (Fig. 1). The gene encoding phytase for P utilization was also
significantly increased by 11.37\% (\emph{P} \textless{} 0.05).
Additionally, a total of 30 oligonucleotide probes of C degradation
genes were only detected in warmed samples (Table S2)\emph{.}

\emph{Network topological characteristics}

Using functional gene data, we generated fMENs for control and warmed
samples. Both networks were scale-free (power-law
R\textsuperscript{2}\textgreater{} 0.85), and possessed other typical
properties of biological networks such as small-world and modularity
(Table 3). In line with Hypothesis 1, the fMEN of warmed samples had
significantly different topological characteristics with the fMEN of
control (\emph{P} \textless{} 0.001). The average path distance of
warmed samples (18.718) was significantly higher than that of control
(10.325) (\emph{P} \textless{} 0.001), whereas the average clustering
coefficient of warmed samples (0.298) was significantly lower than that
of control (0.329) (\emph{P} \textless{} 0.001). A lower percentage of
negative links was found in the fMEN of warmed samples (19.78\%) than in
control (36.32\%) (Fig. 2 \& Table 3).

Based on Zi and Pi values of both networks, 12 genes were identified as
key genes (9 module hubs and 3 connectors) of warmed samples, while 14
genes were identified as key genes (11 module hubs and 3 connectors) of
control samples (Fig. 2 \& Table S4). A much higher percentage of key
genes were associated with C degradation in warmed samples (91.7\%, 11
out of 12 genes) than in control (64.3\%, 9 out of 14 genes) (Table S4).
Partial Mantel tests revealed that thaw depth was significantly linked
to network connectivity of key genes (\emph{P} \textless{} 0.01).

\emph{Linkages between environmental factors and fungal functional gene
composition}

The increase of winter soil temperature and associated changes of soil
hydrologic conditions (including soil moisture and thaw depth) and GPP
in the growing season significantly contributed to the differences of
functional gene composition (CCA model: F = 1.57, \emph{P} = 0.01) (Fig.
3a). Those factors explained 47.25\% of the differences, among which
winter soil temperature explained 10.35\%, soil hydrologic conditions
explained 20.10\%, and GPP explained 12.48\% (VPA, Fig. 3b).

\emph{Linkages between fungal C degradation genes and R}
\textsubscript{eco}

In line with the third hypothesis,\emph{R} \textsubscript{eco} had
positively linear regressions with the relative abundances of fungal C
degradation genes encoding xylose reductase (\emph{R}
\textsuperscript{2} = 0.56, \emph{P} \textless{} 0.01), vanillin
dehydrogenase (\emph{R} \textsuperscript{2} = 0.44,\emph{P} = 0.02),
glyoxal oxidase (\emph{R} \textsuperscript{2} = 0.32,\emph{P} = 0.03),
and aspartate protease (\emph{R} \textsuperscript{2} = 0.33, \emph{P} =
0.03) (Fig. 4).

\textbf{Discussion}

Understanding fungal response to climate warming is important for
revealing their roles in nutrient cycling processes and evaluating the
ecological consequences of climate warming on high-latitude tundra. In
this study, we hypothesized that climate warming would affect fungal
community composition and potential interactions. Our hypotheses were
partially upheld: fungal functional gene composition and potential
interactions among genes/species were significantly altered, whereas
fungal community composition was unchanged (Table 1 \& Fig. 2).
Different responses of functional gene composition and community
composition were similar to the findings of a previous warming study in
a tall-grass prairie ecosystem (Cheng et al., 2017), suggesting that
GeoChip is more sensitive for detecting subtle changes in functional
genes. Due to horizontal gene transfer among microbial species, species
identities detected in GeoChip do not reliably represent phylogenetic
species in environmental samples. However, if both sequencing and
GeoChip data point to the same species, the chance of disconnection
between sequencing and GeoChip is slim. For example, we detected both
functional genes and OTUs belonging to\emph{Hypocrea} in GeoChip data
and amplicon sequencing data. Similar to the overall changes of fungal
composition change, we found significant increases in gene abundances
from \emph{Hypocrea} (e.g., the genes encoding endochitinase,
exopolygalacturonase, and xylanase) (\emph{P}\textless{} 0.05), whereas
no changes of OTUs from this genus. Nonetheless, the lack of detectable
changes in fungal community composition should be interpreted with
caution because 28S rRNA gene has a relatively low taxonomic resolution
(Brown, Rigdon-Huss, \& Jumpponen, 2014), which might fail to detect
subtle changes in fine taxonomic levels (Hultman et al., 2015).

\emph{Ascomycota} was the dominant fungal phylum in our study and other
studies conducted in sedge-occupied tussock tundra (Christiansen et al.,
2016; Deslippe, Hartmann, Simard, \& Mohn, 2012; Timling, Walker,
Nusbaum, Lennon, \& Taylor, 2014). However, this phylum was less
dominant in shrub-occupied tundra (Wallenstein, McMahon, \& Schimel,
2007), indicating that vegetation pattern is probably closely related to
fungal community composition. Climate warming has caused vegetation
shifts in the tundra, with diminishing sedges and bryophytes but
expanding shrubs and trees (Fraser, Lantz, Olthof, Kokelj, \& Sims,
2014; Sturm et al., 2001). Considering the potential relationship
between fungi and plants, it is likely that dominant fungi in
high-latitude tundra will shift accompanying the shrub expansion.

The increase in fungal C degradation capacities, based on higher
relative abundances of functional genes encoding invertase, xylose
reductase, and vanillin dehydrogenase (Fig. 1), might imply that
fungi-mediated C degradation was accelerated under warming. Those genes
are associated with the degradation of complex plant-derived saccharides
(Culleton, Mckie, \& de Vries, 2013), probably resulting from higher
plant productivity (Table S1). Moreover, higher vanillin dehydrogenase
gene abundance might result from the ``priming effect,'' i.e., labile C
input can stimulate recalcitrant C degradation (de Graaff, Classen,
Castro, \& Schadt, 2010; Mau, Dijkstra, Schwartz, Koch, \& Hungate,
2018). Warming-induced increases of fungal C degradation capacities were
also found in other ecosystems including grasslands (Cheng et al.,
2017), subtropical freshwater wetlands (Wang et al., 2012), deserts
(Weber et al., 2011), showing a consistent response pattern of fungal
communities across ecosystems. For example, in the biological soil
crusts of a desert, soil warming increased the abundance of fungal
functional gene \emph{cbhI} , which encodes cellobiohydrolase for
cellulose degradation (Weber et al., 2011).

Higher average path distance but lower average clustering coefficient of
the warmed fMEN (Table 3) suggested that the complexity of potential
interactions decreased under warming. As the percentage of negative
links could represent competition within community members (Fuhrman,
2009), a lower percentage of negative links in warmed fMEN concurred
with higher soil nutrient availability that alleviates resource
competition (Banerjee et al., 2016). Higher nutrient availability
appears to stimulate C degradation, as almost all key genes, the most
important members in maintaining network structure and strongly
influence community stability and functions (Banerjee et al., 2018; Shi
et al., 2016), in warmed fMEN were associated with C degradation (Table
S4). However, it is important to bear in mind a caveat that most
networks cannot distinguish true ecological interactions based on
positive or negative correlations between nodes since it remains largely
intractable to analyze \emph{in situ} microbial interaction
experimentally in most communities of natural environments (Faust \&
Raes, 2012). When interpreting networks in ecological terms, topological
properties (e.g., average node connectivity, average path length,
clustering coefficient, and modularity), which reflect whole-network
changes, could be more reliable.

Legacy effects of winter warming included increases in thaw depth, soil
moisture, and GPP in the growing season (Table S1). The increase of thaw
depth, a sign of permafrost degradation, was also observed in previous
experimental warming studies using snow fences (Hinkel \& Hurd, 2006;
Nowinski, Taneva, Trumbore, \& Welker, 2010). Since we excluded snow
addition effects on soil hydrological conditions by removing snow before
melting, the most likely cause of soil moisture increase was ice-wedge
melting in permafrost (Liljedahl et al., 2016). Soil temperature, soil
moisture, thaw depth, and GPP imposed the strongest influence on fungal
functional gene composition (Fig. 3), suggesting that they were very
important in affecting fungal functional capacities. Oxygen availability
could be affected by soil hydrology and temperature, which is shown to
change species diversity, ecological functions, and survival of most
aerobic fungi (Wang et al., 2012; Zak \& Kling, 2006).

Permafrost thawing increases the potential of old C degradation, leading
to stronger heterotrophic respiration and CO\textsubscript{2} emission
(Nowinski et al., 2010; Schuur et al., 2009). Consistently, we found a
significant warming-induced increase in \emph{R}
\textsubscript{eco}(Table S1 \& Fig. 4), primarily resulting from the
increase in winter\emph{R} \textsubscript{eco} that mainly represents
heterotrophic respiration (by 103.2\%, \emph{P}\textless{} 0.05, Table
S1). Our winter \emph{R} \textsubscript{eco} data were derived from a
site-specific model that assumes winter \emph{R}
\textsubscript{eco}increases with soil temperature exponentially, which
is commonly used when simulating the temperature dependence of
heterotrophic soil respiration. The model fitted \emph{in situ} measured
data well (\emph{R} \textsuperscript{2}\selectlanguage{english} =\selectlanguage{english} 0.70, \emph{P}
\textless{}\selectlanguage{english} 0.001), which was consistent with observations elsewhere
(Tuomi, Vanhala, Karhu, Fritze, \& Liski, 2008).

Increased \emph{R} \textsubscript{eco} shared positive regressions with
fungal functional genes for C degradation (Fig. 4), suggesting that
fungal functional capacities in C degradation might be very important in
mediating \emph{R} \textsubscript{eco,} and thus the C stability of
tundra ecosystems. According to the central dogma, DNA has to be
transcribed into RNA and translated into protein before displaying
enzymatic activity, which is the missing link in our study since enzyme
activities were not measured. However, a large number of papers have
justified the strong, positive relationship between enzyme activities
and abundances of their encoding genes (Blackwood, Waldrop, Zak, \&
Sinsabaugh, 2007; Fan, Li, Wakelin, Yu, \& Liang, 2012; Trivedi et al.,
2016). Notably, the positive correlations between \emph{R}
\textsubscript{eco} and fungal functional genes should be interpreted
with caution since GPP was also increased by warming (Table S1).
Therefore, the increase of \emph{R} \textsubscript{eco} could be
explained in terms of the increased photosynthesis and corresponding
increased respiration of new photosynthate. Additionally, higher
bacterial C degradation capacities under warming conditions at this site
have been previously documented (Xue et al., 2016), which could also
account for the \emph{R} \textsubscript{eco} increase. Therefore, future
studies that distinguish plant respiration and heterotrophic soil
respiration are needed.

\textbf{Concluding Remarks}

This study shows that winter warming rapidly alters fungal functional
gene composition and potential fungal interactions in Alaskan tundra
soil, which are closely linked to the changes of soil hydrological
conditions and plant productivity. Higher fungal C degradation
capacities under warming suggest positive feedback of fungal communities
to a warmer climate because an increase in \emph{R}
\textsubscript{eco}is observed, signifying the potential importance of
fungi in de-stabilizing tundra soil C storage. Although it is uncertain
whether these conclusions are ubiquitous in other arctic tundra
ecosystems due to the limited sampling size and small-scale spatial
variability of this study, it represents an important step forward to
understand the fungal functional potentials changed by warming. Overall,
this study demonstrates fungal community vulnerability to climate
warming and provides new insights into the fungi-mediated C storage in
tundra soil. Therefore, incorporating fungi into current Earth system
models may help predict the future fate of global C storage.

\textbf{References}

Banerjee, S., Kirkby, C. A., Schmutter, D., Bissett, A., Kirkegaard, J.
A., \& Richardson, A. E. (2016). Network analysis reveals functional
redundancy and keystone taxa amongst bacterial and fungal communities
during organic matter decomposition in an arable soil. \emph{Soil
Biology \& Biochemistry, 97} , 188-198.

Banerjee, S., Schlaeppi, K., \& van der Heijden, M. G. A. (2018).
Keystone taxa as drivers of microbiome structure and
functioning.\emph{Nature Reviews Microbiology, 16} (9), 567-576.

Blackwood, C. B., Waldrop, M. P., Zak, D. R., \& Sinsabaugh, R. L.
(2007). Molecular analysis of fungal communities and laccase genes in
decomposing litter reveals differences among forest types but no impact
of nitrogen deposition. \emph{Environmental Microbiology, 9} (5),
1306-1316.

Brown, S. P., Rigdon-Huss, A. R., \& Jumpponen, A. (2014). Analyses of
ITS and LSU gene regions provide congruent results on fungal community
responses. \emph{Fungal Ecology, 9} , 65-68.

Buckeridge, K. M., Banerjee, S., Siciliano, S. D., \& Grogan, P. (2013).
The seasonal pattern of soil microbial community structure in mesic low
arctic tundra. \emph{Soil Biology \& Biochemistry, 65} (65), 338-347.

Cheng, L., Zhang, N. F., Yuan, M. T., Xiao, J., Qin, Y. J., Deng, Y., .
. . Zhou, J. Z. (2017). Warming enhances old organic carbon
decomposition through altering functional microbial
communities.\emph{The Isme Journal, 11} (8), 1825-1835.

Christiansen, C. T., Haugwitz, M. S., Prieme, A., Nielsen, C. S.,
Elberling, B., Michelsen, A., . . . Blok, D. (2016). Enhanced summer
warming reduces fungal decomposer diversity and litter mass loss more
strongly in dry than in wet tundra. \emph{Global Change Biology, 23}
(1), 406.

Clemmensen, K. E., Michelsen, A., Jonasson, S., \& Shaver, G. R. (2006).
Increased ectomycorrhizal fungal abundance after long-term fertilization
and warming of two arctic tundra ecosystems. \emph{New Phytologist, 171}
(2), 391-404.

Culleton, H., Mckie, V., \& de Vries, R. P. (2013). Physiological and
molecular aspects of degradation of plant polysaccharides by fungi: What
have we learned from Aspergillus? \emph{Biotechnology Journal, 8} (8),
884-U822.

de Graaff, M. A., Classen, A. T., Castro, H. F., \& Schadt, C. W.
(2010). Labile soil carbon inputs mediate the soil microbial community
composition and plant residue decomposition rates. \emph{New
Phytologist, 188} (4), 1055-1064.

de Vries, F. T., Griffiths, R. I., Bailey, M., Craig, H., Girlanda, M.,
Gweon, H. S., . . . Bardgett, R. D. (2018). Soil bacterial networks are
less stable under drought than fungal networks. \emph{Nature
Communications, 9} , 3033.

De Vries, F. T., Liiri, M. E., Bjornlund, L., Bowker, M. A.,
Christensen, S., Setala, H. M., \& Bardgett, R. D. (2012). Land use
alters the resistance and resilience of soil food webs to
drought.\emph{Nature Climate Change, 2} (4), 276-280.

Deng, Y., Jiang, Y. H., Yang, Y. F., He, Z. L., Luo, F., \& Zhou, J. Z.
(2012). Molecular ecological network analyses. \emph{Bmc Bioinformatics,
13} (1), 113.

Deslippe, J. R., Hartmann, M., Simard, S. W., \& Mohn, W. W. (2012).
Long\selectlanguage{english}-term warming alters the composition of Arctic soil microbial
communities. \emph{Fems Microbiology Ecology, 82} (2), 303-315.

Fan, F. L., Li, Z. J., Wakelin, S. A., Yu, W. T., \& Liang, Y. C.
(2012). Mineral fertilizer alters cellulolytic community structure and
suppresses soil cellobiohydrolase activity in a long-term fertilization
experiment. \emph{Soil Biology \& Biochemistry, 55} , 70-77.

Faust, K., \& Raes, J. (2012). Microbial interactions: from networks to
models. \emph{Nature Reviews Microbiology, 10} (8), 538-550.

Fraser, R. H., Lantz, T. C., Olthof, I., Kokelj, S. V., \& Sims, R. A.
(2014). Warming-induced shrub expansion and lichen decline in the
western Canadian Arctic. \emph{Ecosystems, 17} (7), 1151-1168.

Fuhrman, J. A. (2009). Microbial community structure and its functional
implications. \emph{Nature, 459} (7244), 193-199.

Guimera, R., \& Amaral, L. A. N. (2005). Functional cartography of
complex metabolic networks. \emph{Nature, 433} (7028), 895-900.

He, Z., Gentry, T. J., Schadt, C. W., Wu, L., Liebich, J., Chong, S. C.,
. . . Zhou, J. (2007). GeoChip: a comprehensive microarray for
investigating biogeochemical, ecological and environmental
processes.\emph{The Isme Journal, 1} , 67.

Hinkel, K. M., \& Hurd, J. K. (2006). Permafrost destabilization and
thermokarst following snow fence installation, Barrow, Alaska,
USA.\emph{Arctic Antarctic and Alpine Research, 38} (4), 530-539.

Hultman, J., Waldrop, M. P., Mackelprang, R., David, M. M., Mcfarland,
J., Blazewicz, S. J., . . . Shah, M. B. (2015). Multi-omics of
permafrost, active layer and thermokarst bog soil
microbiomes.\emph{Nature, 521} (7551), 208.

Levy-Booth, D. J., Prescott, C. E., \& Grayston, S. J. (2014). Microbial
functional genes involved in nitrogen fixation, nitrification and
denitrification in forest ecosystems. \emph{Soil Biology \&
Biochemistry, 75} , 11-25.

Liljedahl, A. K., Boike, J., Daanen, R. P., Fedorov, A. N., Frost, G.
V., Grosse, G., . . . Zona, D. (2016). Pan-Arctic ice-wedge degradation
in warming permafrost and its influence on tundra hydrology.\emph{Nature
Geoscience, 9} (4), 312-318.

Liu, S., Wang, F., Xue, K., Sun, B., Zhang, Y., He, Z., . . . Yang, Y.
(2015). The interactive effects of soil transplant into colder regions
and cropping on soil microbiology and
biogeochemistry.\emph{Environmental Microbiology, 17} (3), 566-576.

Malcolm, G. M., Lopez-Gutierrez, J. C., Koide, R. T., \& Eissenstat, D.
M. (2008). Acclimation to temperature and temperature sensitivity of
metabolism by ectomycorrhizal fungi. \emph{Global Change Biology, 14}
(5), 1169-1180.

Mann, P. J., Eglinton, T. I., Mcintyre, C. P., Zimov, N., Davydova, A.,
Vonk, J. E., . . . Spencer, R. G. M. (2015). Utilization of ancient
permafrost carbon in headwaters of Arctic fluvial networks. \emph{Nature
Communications, 6} , 7856.

Mau, R. L., Dijkstra, P., Schwartz, E., Koch, B. J., \& Hungate, B. A.
(2018). Warming induced changes in soil carbon and nitrogen influence
priming responses in four ecosystems. \emph{Applied Soil Ecology, 124} ,
110-116.

Morgado, L. N., Semenova, T. A., Welker, J. M., Walker, M. D., Smets,
E., \& Geml, J. (2016). Long-term increase in snow depth leads to
compositional changes in arctic ectomycorrhizal fungal
communities.\emph{Global Change Biology, 22} (9), 3080-3096.

Mundra, S., Halvorsen, R., Kauserud, H., Bahram, M., Tedersoo, L.,
Elberling, B., . . . Eidesen, P. B. (2016). Ectomycorrhizal and
saprotrophic fungi respond differently to long-term experimentally
increased snow depth in the High Arctic. \emph{Microbiologyopen, 5} (5),
856-869.

Myers-Smith, I. H., Elmendorf, S. C., Beck, P. S. A., Wilmking, M.,
Hallinger, M., Blok, D., . . . Vellend, M. (2015). Climate sensitivity
of shrub growth across the tundra biome. \emph{Nature Climate Change, 5}
, 887-891.

Natali, S. M., Eag, S., \& Rubin, R. L. (2012). Increased plant
productivity in Alaskan tundra as a result of experimental warming of
soil and permafrost. \emph{Journal of Ecology, 100} (2), 488-498.

Natali, S. M., Schuur, E. A. G., Trucco, C., Pries, C. E. H., Crummer,
K. G., \& Lopez, A. F. B. (2011). Effects of experimental warming of
air, soil and permafrost on carbon balance in Alaskan
tundra.\emph{Global Change Biology, 17} (3), 1394-1407.

Natali, S. M., Schuur, E. A. G., Webb, E. E., Pries, C. E. H., \&
Crummer, K. G. (2014). Permafrost degradation stimulates carbon loss
from experimentally warmed tundra. \emph{Ecology, 95} (3), 602-608.

Nowinski, N. S., Taneva, L., Trumbore, S. E., \& Welker, J. M. (2010).
Decomposition of old organic matter as a result of deeper active layers
in a snow depth manipulation experiment. \emph{Oecologia, 163} (3),
785-792.

Rinnan, R., Michelsen, A., Baath, E., \& Jonasson, S. (2007).
Mineralization and carbon turnover in subarctic heath soil as affected
by warming and additional litter. \emph{Soil Biology \& Biochemistry,
39} (12), 3014-3023.

Rovira, P., \& Vallejo, V. R. (2002). Labile and recalcitrant pools of
carbon and nitrogen in organic matter decomposing at different depths in
soil: an acid hydrolysis approach. \emph{Geoderma, 107} (1-2), 109-141.

Schadt, C. W., Martin, A. P., Lipson, D. A., \& Schmidt, S. K. (2003).
Seasonal dynamics of previously unknown fungal lineages in tundra
soils.\emph{Science, 301} (5638), 1359-1361.

Schuur, E. A. G., Crummer, K. G., Vogel, J. G., \& Mack, M. C. (2007).
Plant species composition and productivity following permafrost thaw and
thermokarst in Alaskan tundra. \emph{Ecosystems, 10} (2), 280-292.

Schuur, E. A. G., Vogel, J. G., Crummer, K. G., Lee, H., Sickman, J. O.,
\& Osterkamp, T. E. (2009). The effect of permafrost thaw on old carbon
release and net carbon exchange from tundra. \emph{Nature, 459} (7246),
556-559.

Semenova, T. A., Morgado, L. N., Welker, J. M., Walker, M. D., Smets,
E., \& Geml, J. (2016). Compositional and functional shifts in arctic
fungal communities in response to experimentally increased snow
depth.\emph{Soil Biology \& Biochemistry, 100} , 201-209.

Shi, S. J., Nuccio, E. E., Shi, Z. J., He, Z. L., Zhou, J. Z., \&
Firestone, M. K. (2016). The interconnected rhizosphere: High network
complexity dominates rhizosphere assemblages. \emph{Ecology Letters, 19}
(8), 926-936.

Sturm, M., Racine, C., \& Tape, K. (2001). Climate change - Increasing
shrub abundance in the Arctic. \emph{Nature, 411} (6837), 546-547.

Tarnocai, C., Canadell, J. G., Schuur, E. A. G., Kuhry, P., Mazhitova,
G., \& Zimov, S. (2009). Soil organic carbon pools in the northern
circumpolar permafrost region. \emph{Global Biogeochemical Cycles, 23} ,
GB2023.

Tedersoo, L., Bahram, M., Polme, S., Koljalg, U., Yorou, N. S.,
Wijesundera, R., . . . Abarenkov, K. (2014). Global diversity and
geography of soil fungi. \emph{Science, 346} (6213), 1256688.

Timling, I., Walker, D. A., Nusbaum, C., Lennon, N. J., \& Taylor, D. L.
(2014). Rich and cold: diversity, distribution and drivers of fungal
communities in patterned-ground ecosystems of the North American
Arctic.\emph{Molecular Ecology, 23} (13), 3258-3272.

Trivedi, P., Delgado-Baquerizo, M., Trivedi, C., Hu, H. W., Anderson, I.
C., Jeffries, T. C., . . . Singh, B. K. (2016). Microbial regulation of
the soil carbon cycle: evidence from gene-enzyme relationships.\emph{The
Isme Journal, 10} (11), 2593-2604.

Tuomi, M., Vanhala, P., Karhu, K., Fritze, H., \& Liski, J. (2008).
Heterotrophic soil respiration---comparison of different models
describing its temperature dependence. \emph{Ecological Modelling, 211}
(1-2), 182-190.

Walker, M. (1996). \emph{Community baseline measurements for ITEX
studies} (2nd edn ed.). Copenhagen: Danish Polar Centre.

Wallenstein, M. D., McMahon, S., \& Schimel, J. (2007). Bacterial and
fungal community structure in Arctic tundra tussock and shrub
soils.\emph{Fems Microbiology Ecology, 59} (2), 428-435.

Wang, H., He, Z., Lu, Z., Zhou, J., Van Nostrand, J. D., Xu, X., \&
Zhang, Z. (2012). Genetic linkage of soil carbon pools and microbial
functions in subtropical freshwater wetlands in response to experimental
warming. \emph{Applied and Environmental Microbiology, 78} (21),
7652-7661.

Weber, C. F., Zak, D. R., Hungate, B. A., Jackson, R. B., Vilgalys, R.,
Evans, R. D., . . . Kuske, C. R. (2011). Responses of soil cellulolytic
fungal communities to elevated atmospheric CO\textsubscript{2} are
complex and variable across five ecosystems. \emph{Environmental
Microbiology, 13} (10), 2778-2793.

Weinstein, R. N., Montiel, P. O., \& Johnstone, K. (2000). Influence of
growth temperature on lipid and soluble carbohydrate synthesis by fungi
isolated from fellfield soil in the maritime Antarctic. \emph{Mycologia,
92} (2), 222-229.

Wu, L., Yang, Y., Chen, S., Zhao, M., Zhu, Z., Yang, S., . . . Zhou, J.
(2016). Long-term successional dynamics of microbial association
networks in anaerobic digestion processes. \emph{Water Research, 104} ,
1-10.

Xue, K., Yuan, M. M., Zhou, J. S., Qin, Y., Deng, Y., Cheng, L., . . .
Bracho, R. (2016). Tundra soil carbon is vulnerable to rapid microbial
decomposition under climate warming. \emph{Nature Climate Change, 6} ,
595--600.

Yergeau, E., Kang, S., He, Z., Zhou, J., \& Kowalchuk, G. A. (2007).
Functional microarray analysis of nitrogen and carbon cycling genes
across an Antarctic latitudinal transect. \emph{The Isme Journal, 1}
(2), 163-179.

Yuste, J. C., Penuelas, J., Estiarte, M., Garciamas, J., Mattana, S.,
Ogaya, R., . . . Sardans, J. (2011). Drought-resistant fungi control
soil organic matter decomposition and its response to
temperature.\emph{Global Change Biology, 17} (3), 1475-1486.

Zak, D. R., \& Kling, G. W. (2006). Microbial community composition and
function across an arctic tundra landscape. \emph{Ecology, 87} (7),
1659-1670.

Zhao, M., Xue, K., Wang, F., Liu, S., Bai, S., Sun, B., . . . Yang, Y.
(2014). Microbial mediation of biogeochemical cycles revealed by
simulation of global changes with soil transplant and cropping.\emph{The
Isme Journal, 8} (10), 2045-2055.

Zhou, J. Z., Xue, K., Xie, J. P., Deng, Y., Wu, L. Y., Cheng, X. H., . .
. Luo, Y. Q. (2012). Microbial mediation of carbon-cycle feedbacks to
climate warming. \emph{Nature Climate Change, 2} (2), 106-110.

\textbf{Acknowledgments}

We thank three anonymous reviewers and the editor for constructive
comments and suggestions to make this paper greatly improved. This study
is supported by the Second Tibetan Plateau Scientific Expedition and
Research (STEP) program (2019QZKK0503), the National Science Foundation
of China (41877048 and 41825016), US Department of Energy, Office of
Science, Genomic Science Program under Award Numbers DE-SC0004601 and
DE-SC0010715, the NSF LTER program, and the Office of the Vice President
for Research at the University of Oklahoma. The authors declare no
conflict of interest.

\textbf{Authors' Contributions}

JZ, ES, JRC, and JT conceived the ideas and designed methodology; MY,
ZS, and LW collected the data; MY did sequencing and GeoChip
hybridization; JC and QG analyzed the data; JC and YY wrote the first
draft of the manuscript. All authors contributed to the manuscript and
gave final approval for publication.

\textbf{Data Accessibility}

All data used in this manuscript have been uploaded. GeoChip raw and
normalized signal intensities can be found in Gene Expression Omnibus
(GEO) database by the accession number GSE127899. 28S rRNA gene
sequencing data can be found in Sequence Read Archive (SRA) database by
the project number PRJNA525959. All the environmental data are available
online (LTER Network Data Portal,
https://portal.lternet.edu/nis/home.jsp).

\textbf{Figures and Tables}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

\textbf{Figure 1.} Winter warming increased the relative abundances of
fungal functional genes for C degradation and P utilization. Response
ratios show significant increases of genes encoding xylose reductase,
vanillin dehydrogenase, invertase, and phytase (indicated by **
for\emph{P} \textless{} 0.01 and * for \emph{P} \textless{} 0.05). C
substrates are arranged in the order from labile C to recalcitrant C.
Error bars indicate 95\% confidence intervals of differences of
functional gene data between warmed and control samples.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

\textbf{Figure 2.} Winter warming changes topological characteristics
and key genes of functional molecule ecological networks
(fMENs):\textbf{(a)} the fMEN of control samples; \textbf{(b)} the fMEN
of warmed samples. Each dot or diamond represents a node (fungal
functional gene), and each line represents an edge (potential
interaction). Networks are separated into modules by using the greedy
modularity optimization method and shown as separate circles. Diamond
nodes represent key genes of fMENs. Diamond nodes inside modules
represent module hubs. Diamond nodes outside modules represent
connectors. Different colors of diamond nodes represent different
functional genes, and the number marked on each diamond node represents
fungal species containing the functional gene. Red lines represent
negative links, and black lines represent positive links.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

\textbf{Figure 3. (a)} Winter soil temperature (winter soil T), thaw
depth, soil moisture (soil M) and gross primary productivity (GPP) are
linked to fungal functional gene composition, as revealed by the
canonical correlation analysis (CCA); \textbf{(b)} The variation
partitioning analysis (VPA) shows explanatory percentages of
environmental factors to fungal functional gene composition.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

\textbf{Figure 4.} Significant linear regressions between ecosystem
respiration (\emph{R} \textsubscript{eco}) and relative abundances of
fungal functional genes encoding xylose reductase, vanillin
dehydrogenase, glyoxal oxidase, and aspartate protease.

\textbf{Table 1.} Results of the multi-response permutation procedure
(MRPP) and analysis of similarity (ANOSIM) based on three distance
algorithms. Bold values represent\emph{P} \textless{} 0.05.\selectlanguage{english}
\begin{longtable}[]{@{}lllllllll@{}}
\toprule
Distance method & MRPP & MRPP & MRPP & MRPP & ANOSIM & ANOSIM & ANOSIM &
ANOSIM\tabularnewline
\midrule
\endhead
& Functional genes & Functional genes & OTUs & OTUs & Functional genes &
Functional genes & OTUs & OTUs\tabularnewline
& \selectlanguage{greek}δ & \emph{P} & \selectlanguage{greek}δ & \emph{P} & R & \emph{P} & R &
\emph{P}\tabularnewline
Bray-Curtis & 0.11 & \textbf{0.03} & 0.75 & 0.64 & 0.32 & \textbf{0.01}
& -0.18 & 0.94\tabularnewline
Euclidean & 17.74 & \textbf{0.02} & 5008 & 0.59 & 0.37 &
\textbf{\textless{} 0.01} & -0.15 & 0.89\tabularnewline
Horn & 0.08 & \textbf{0.02} & 0.72 & 0.64 & 0.35 & \textbf{\textless{}
0.01} & -0.13 & 0.83\tabularnewline
\bottomrule
\end{longtable}

\textbf{Table 2.} The similarity percentage (SIMPER) analysis reveals
genes with significant contributions (\emph{P} \textless{} 0.05) to the
differences of functional gene composition between control and warmed
samples.\selectlanguage{english}
\begin{longtable}[]{@{}llll@{}}
\toprule
Gene product & Functional process of the gene product & Contribution
(\%) & \emph{P}\tabularnewline
\midrule
\endhead
Cellobiase & C degradation (Cellulose) & 4.97 & 0.04\tabularnewline
Sulfate transporter & S transportation & 2.93 & 0.03\tabularnewline
Xylanase & C degradation (Hemicellulose) & 2.51 & 0.03\tabularnewline
Glyoxal oxidase & C degradation (Lignin) & 2.42 & 0.03\tabularnewline
Metalloprotease & C degradation (Protein) & 2.01 & 0.03\tabularnewline
Phospholipase A2 & C degradation (Phospholipids) & 1.75 &
0.02\tabularnewline
Phytase & P utilization & 1.50 & 0.01\tabularnewline
PAPS reductase & S assimilation & 0.77 & 0.04\tabularnewline
Superoxide dismutase & Stress tolerance & 0.59 & 0.04\tabularnewline
Vanillin dehydrogenase & C degradation (Aromatics) & 0.53 & \textless{}
0.01\tabularnewline
Xylose reductase & C degradation (Hemicellulose) & 0.50 &
0.01\tabularnewline
Invertase & C degradation (Sucrose) & 0.38 & 0.03\tabularnewline
\bottomrule
\end{longtable}

\textbf{Table 3.} Topological characteristics of functional molecule
ecological networks (fMENs) of control and warmed samples and their
corresponding random networks. Bold values represent \emph{P}
\textless{} 0.001\selectlanguage{english}
\begin{longtable}[]{@{}lll@{}}
\toprule
Network topological characteristics & Control & Warmed\tabularnewline
\midrule
\endhead
fMENs of control and warmed samples & &\tabularnewline
Threshold & 0.98 & 0.98\tabularnewline
Total nodes & 1211 & 1122\tabularnewline
Total edges & 3009 & 2862\tabularnewline
R\textsuperscript{2} of power-law & 0.855 & 0.877\tabularnewline
Average connectivity & 4.969 & 5.102\tabularnewline
Average clustering coefficient & \textbf{0.329} &
\textbf{0.298}\tabularnewline
Average path distance & \textbf{10.325} & \textbf{18.718}\tabularnewline
Modularity & \textbf{0.778} & \textbf{0.703}\tabularnewline
Negative interactions & 36.32\% & 19.78\%\tabularnewline
Random networks & &\tabularnewline
Average clustering coefficient & 0.009\selectlanguage{ngerman}±0.002 &
0.022\selectlanguage{ngerman}±0.003\tabularnewline
Average path distance & 4.235\selectlanguage{ngerman}±0.022 & 3.948\selectlanguage{ngerman}±0.030\tabularnewline
Modularity & 0.437\selectlanguage{ngerman}±0.004 & 0.417\selectlanguage{ngerman}±0.003\tabularnewline
\bottomrule
\end{longtable}\selectlanguage{ngerman}

\selectlanguage{english}
\FloatBarrier
\end{document}

