Tong Bao · Shilin Deng · Kaiyue Yu · Weiyi Li ·Airong Dong

Abstract The effect of seasons on the soil microbiome in a Larix gmelinii forest of Mohe, China, where winter temperatures are generally below ? 40 °C, was evaluated with metagenomics analysis. Taxonomic prof iling using sequencing information revealed that Proteobacteria, Actinobacteria,Acidobacteria and Verrucomicrobia were the dominant phyla in spring, summer, and fall, as were Bradyrhizobium, Chthoniobacter, Streptomyces, Acid Candidatus Koribacter at the genus level. Some species that were abundant in spring and fall greatly diminished in abundance in summer. Clusters of orthologous groups (COG) of proteins, carbohydrate-active enzymes(CAZy), Kyoto Encyclopedia of Genes and Genomes (KEGG)and NCBI databases were used to elucidate the function of diverse proteins and metabolites of the microbial community of L. gmelinii forest. COG analysis showed that fewer genes were detected in spring than in fall and summer, indicating that many soil microbes in the L. gmelinii forest were not tolerant to cold. Based on KEGG analysis, some pathways in the soil microbes were activated in spring and autumn and deactivated in summer. CAZy analysis revealed that most CAZy were more active in summer than in spring or autumn and were severely inhibited in the spring. Many functional pathways, proteins,and CAZy involved in the community changes were concerned with cold or heat resistance. Therefore, the soil in the L. gmelinii forest can be a valuable resource for further research on heat and cold tolerance of soil microbes.

Keywords Soil microbial community · Microbial function · Metagenomic · Seasonal variation

Introduction

The Great Xing?an Mountains is covered with the largest modern, state-owned forest in northern China. Mohe County,in the northern area of these mountains in Heilongjiang Province, is in a cold temperate zone with a continental climate and is at the highest latitude with the coldest winter in China. In winter, the temperature is generally below ? 40 °C.The extreme cold directly affects native tree species and the microbiome in their rhizosphere. Larix gmelinii is a native species and forms a representative forest in Mohe.Extremophiles grow best in extreme environments (Rothschild and Mancinelli 2001), which can have a variety of conditions, and they have specific metabolic abilities and physical structures to survive (Rampelotto 2016). Liquid water is the essential substance for life, so water at boiling and freezing points can limit life. At subzero temperatures,there can still be liquid water present in soils in the form of a thin liquid film on the surface of soil particles (Rivkina et al.2000; Steven et al. 2006). Temperature and water interact and affect the diversity of soil microbes, but numerous extremophiles have adapted to temperatures below 0 °C (Chen et al. 1993; Xu et al. 1985). Therefore, it is of great economic and scientific value to study the seasonal variation of soil microbial diversity in L. gmelinii forest of Mohe County.

Microorganisms that live in extreme environments such as high temperature, cold, acid and alkali are more diffi-cult to isolate, cultivate and thus study than those living in ordinary environments, metagenomic technology provides a great advantage in the study of these extremophiles (Wu et al. 2017). Therefore, we used a metagenomics analysis to study seasonal changes in the soil microbiome of L. gmelinii forest of Mohe, China.

Materials and methods

Study area, field sampling and metagenomic sequencing

The study area in Mohe County, Heilongjiang Province,China has a mean annual temperature of ? 4 °C and mean annual precipitation of about 400 mm, with most precipitation from June to August (Wang et al. 2015). The forest coverage is over 90% (Gao and Zhang 2018). Soil sampling sites were selected randomly in the forest, which is dominated by Larix gmelinii forest. Soil samples were collected in three seasons (May, July and September, 2018)from the same locations described in Table 1. Samples taken on May 2, 2018 were designated as D1 (Yang and Li 1985). May was characterized by a freezing-thawing period when ice and snow at the sampling sites started to melt (Kikani et al. 2017). Thus, May, with cold features of winter and the revival signs of spring, was the most suitable choice for the cold period. Samples taken on July 11,2018, the hottest time of the year in Mohe County, were designated as D2. Samples taken on September 23, 2018,the most suitable choice for the mild environment and reference, microbiome were designated as D3. Briefly, leaflitter and any animal remains were removed first, then soil 0-10 cm below the soil surface was collected with a sterile trowel, sifter through a 2-mm-mesh sieve and kept in a sterile self-sealing bag. Three types of sample sites were selected each season, and each sample site was sampled three times in an S-shape. Three soil samples were taken from each repeat. Each soil sample contained a mixture of soil from five to eight sampling sites and weighed at least 50 g (Salam et al. 2017). All samples were placed in a freezer right after collection and kept frozen until used.

Metagenomic sequencing, which has been widely used,was adopted to analyze the species and functional diversity of microorganisms in the soil samples (Zhang 2014).

Data assessment, quality control and assembly

The raw image data file generated by Illumina Hiseq(ST1000-3) was turned into raw sequenced reads by CASAVA base calling analysis, which we called raw data or raw reads and contains sequence information of the sequenced reads and its corresponding sequencing quality information. The raw sequenced reads were evaluated by FastQC (version 0.11.2) and filtered by Trimmomatic (version 0.36) to obtain relatively accurate and effective data (Bolger et al. 2014;Remove sequences with N bases; Remove the joint in the reads sequences (foward: AGA TCG GAA GAG CAC ACG TCT GAA C; reverse: AGA TCG GAA GAG CGT CGT GTA GGG A); low quality bases were removed from reads 3′-5′(Q < 20); low quality bases were removed from reads 5′-3′(Q{N0} < 20); the bases with the reads tail mass value below 20 (the window size was 5 bp) were removed by sliding window method; reads less than 35 nt in length and their paired reads were removed.). The high quality reads were assembled by IDBA_UD (version 1.1.2) based on the principle of De Bruijn graph. According to overlapping between reads, contigs were obtained (Peng et al. 2012). Then several assembly results of K-mer were comprehensively evaluated, and best K-mer assembly results were selected.

Gene prediction and species taxonomic analysis

Prodigal (version 2.60) was used to predict the open reading frames (ORFs) of the splicing results. Genes ≥ 100 bp wereselected, and the amino acid sequences were predicted (Liu et al. 2013). For the gene prediction results of each sample,CD-HIT (version 4.6) (Li and Godzik 2006) software was used to remove redundancy and obtain a non-redundant gene set. Bowtie 2 (version 2.1.0; Langmead and Salzberg 2012)was used to compare clean reads with the non-redundant gene set. The number of reads was obtained using SAMtools(version 0.1.18; Li et al. 2009), and gene abundances in each sample were calculated by combining the gene length.

Table 1 Trees species, location, and weather characteristics where soil samples were collected in Mohe County, Heilongjiang Province, China

where Sobsis the total number of genes, niis the number of sequences contained in theith gene, N is the total number of sequences) and the number of genes were used to characterize microbial genetic diversity (Shannon 1948).

DIAMOND (version 0.8.20) was used for a blastp (NCBI Blast+ version 2.28) (Buchf ink et al. 2015) homology comparison between the predicted protein sequences and the Nrdatabase to obtain functional annotation and homologous species information (Altschul 1997). Screening standard was set to E-value < 10?5, Score > 60. The NCBI non-redundant protein sequences is a non-redundant protein sequences database collected by NCBI, which includes all non-redundant GenBank CDS translation sequences, Protein DataBank(PDB), Swiss-Prot and other databases such as Protein Information Resource (PIR) and Protein Research Foundation(PDF). Taxonomy, gene annotation and their relative abundances were obtained based on NCBI database.

Functional analysis

The nucleotide sequences of the genes were compared to the information NR, KEGG, eggNOG, and CAZy databases to acquire taxonomic information and functional annotation of the genes. Function and species abundances were obtained based on the abundances of the genes. Then the species and functional analysis, and other multi-directional statistical analysis and exploration were conducted.

The database of clusters of orthologous groups (COG) is an attempt on a phylogenetic classification of the proteins encoded in 21 complete genomes of bacteria, archaea and eukaryotes (Tatusov et al. 2000). The database of EggNOG(evolutionary genealogy of genes: non-supervised orthologous groups) annotated the constructed clusters of orthologous groups of proteins by Smith-Waterman comparison algorithm. DIAMOND (DIAMOND Version 0.8.20) (Buchf ink et al. 2015) was used to compare the gene set protein sequence with eggNOG (Huerta-Cepas et al. 2016) database to obtain the COG corresponding to the gene. Screening conditions: E-value < 1e?5, Score > 60. On this basis, the functional annotation and classification of the gene set were carried out, and the abundances of each COG functional level in each sample were calculated.

KEGG is a relatively complete database of biological systems, integrating chemical, genomic and system functional information (Kanehisa and Goto 2000). GhostKOALA(GhostKOALA Version 1.0) (Kanehisa et al. 2016) was used to compare the gene set protein sequence with the KEGG database (http://www.kegg.jp), and the KO number corresponding to the sequence was obtained. The pathway and module annotation information of the sequence was obtained based on the link among KO, pathway and module. The abundances of KEGG functional levels in each sample were calculated.

The carbohydrate-active enzymes database (CAZy; Lombard et al. 2014) includes a family of enzymes that catalyze carbohydrate degradation, modification, and biosynthesis(Lombard et al. 2014). HMMER3 (version 3.1b1) (Eddy 2009) was used to compare the gene sequences with those in the CAZy database to generate the corresponding carbohydrate active enzyme annotation. The screening condition was e-value < 1e?5, and the abundance of the carbohydrate active enzyme at each functional level was calculated.

Results

Overview of assembly

From the metagenomic sequencing, 28.2 GB of raw reads were obtained. DNA sequencing from the samples from the three soil depths (D1, D2 and D3) by Illumina resulted in 188,279,256 raw reads with a total of 28,241,888,400 bp(average: 150 ± 0 bp). After quality control was implemented, the total reads decreased to 180,744,014 sequence reads with 26,185,371,307 bp (average: 144 ± 4 bp). The sequences from soil microorganism from D1, D2 and D3 soils were respectively assembled into 8727, 51,806 and 46,490 contigs (> 500 bp), consisting of 7,905,287,54,363,597, and 45,259,034 bp (average: 905.84, 1049.37 and 973.52 bp), and GC content of 56.98%, 58.49%, and 61.91%. Table 2 shows that the amount of data was enough for this study (Albertsen et al. 2012). Prodigal predicted 12,225 ORFs from the D1 soil, 84,168 from D2 and 73,209 ORFs from D3. Genes larger than or equal to 100 bp were selected for predicting amino acid sequences. D2 had more ORFs than did D1 and D3 soils, which equaled the contig distribution.

Table 2 Overview of the metagenomic sequencing

COG annotation and analysis

In the analysis of the gene set among all samples using egg NOG (http://eggno gdb.embl.de/), except for functions unknown, replication, recombination and repair, energy production and conversion, amino acid transport and metabolism and signal transduction mechanisms were the main functions represented from 25 categories (Fig. 1). Of the known functions, replication, recombination and repair was predominant in three samples, and DNA recombination is reported to be an important part of DNA replication and cell survival (Klein and Kreuzer 2002). However, the gene number of replication, recombination and repair in D1 was notably lower than in D2 and D3. The gene number of the next three categories (amino acid transport and metabolism,energy production and conversion and signal transduction mechanisms) above in D2 was higher than that of D1 and D3. In the same way, the gene number of the three categories above in D1 was notably lower. The fewest genes were found for RNA processing and modification, cytoskeleton,and extracellular structures (each ≤ 0.01%) among the soil samples from the three seasons. Genes for nuclear structure and for general function were not found.

Fig. 1 COG functional annotation by percentage for the microbiome in the three soil depths (D1-D3). The vertical axis shows the full names of the 25 functional categories, and the horizontal axis shows the percentage of genes in the corresponding functional categories

These sequencing data indicated that replication, recombination and repair and signal transduction mechanisms were the dominant functions of soil microorganisms in the L. gmelinii forest, especially for soils D2 and D3. Recombination allows replication to succeed proceed through the blocks presented by damaged or folded replication forks,and signals and controls switches between replication and recombination (Klein and Kreuzer 2002). Cold temperatures were accompanied by a significant reduction in the number of genes for replication, recombination and repair and signal transduction Mechanisms. These results indicate that many molecular mechanisms adapted to low temperature have been developed by soil microorganism for cold survival (Fan and Sun 2008). Weather warming stimulated microbial metabolism, shown by the rise in the number of metabolic genes in D2, which is expected since amino acid transport and metabolism and energy production and conversion are closely connected with the growth and development of microorganisms (Feng et al. 2018).

Diversity analysis of the microbial communities

Taxonomic annotation and abundances of the species were obtained using Nr(NCBI non-redundant protein sequences;(http://ncbi.nlm.nih.gov/), and clustered among three samples. “Norank” designates the lack of a clear classification or name at a certain taxonomic level, and “noname” designates the species could not be determined. “Other” was used for species with less than 0.10% abundance (Figs. 2, 3).

Alpha diversity is an indicator for evaluating the diversity within communities. Shannon-Wiener index belongs to Alpha diversity, and the diversity of gene number also reflects Alpha diversity (Xing et al. 2007; Wang et al. 2002).In Table 2, the number of genes and Shannon-Wiener index in summer (D2) were higher than in the spring (D1) and autumn (D3). Further evidence that genetic diversity and species diversity of soil microorganisms in L. gmelinii forest is higher in the spring than in the spring and autumn.

For diversity of species composition among the three seasons, 96 phyla (excluding norank) were identified in the three samples. Proteobacteria, Actinobacteria , Acidobacteria, and Verrucomicrobiaas were dominant phyla in all three samples,accounting for more than 70% of the total number of phyla(Fig. 2). In the three samples, Proteobacteria was the most abundant, with D1 accounting for 49.93%, D2 for 44.33% and D3 for 45.43% of the total abundance by phyla. These bacteria are the most abundant of all microorganisms and are crucial in many biochemical processes in soil (Miah et al. 2010). Alphaproteobacteria, Betaproteobacteria, Gammaproteobacteria,Deltaproteobacteria and Epsilonproteobacteria in the three samples were belonged to Proteobacteria. And the abundance of Alphaproteobacteria was the highest. Actinobacteria was the second most abundant phylum in D1 (14.29%) and D3 (26.76%),but lower in D2 (4.52%). Therefore, taxonomic prof iling showed Proteobacteria, Actinobacteria, Acidobacteria, and Verrucomicrobia played important roles in the soil microbial communities.

Fig. 2 Relative abundances and cluster tree for phyla represented in the three seasonal soil samples (D1-D3). On the left is the similarity tree. The smaller the difference between samples, the samples will be in the same branch with similar color groups information differentiation. The histogram on the right shows the distribution of species in the sample. Different colors represent different phyla

Fig. 3 Relative abundances and cluster tree for major genera represented in the three seasonal soil samples (D1-D3). On the left is the similarity tree. The smaller the difference between samples, the sam-ples will be in the same branch with similar color groups information differentiation. The histogram on the right shows the distribution of species in the samples. Different colors represent different species

At the genus level, 1830 genera (except norank) were identified from the three samples. Among three samples,1644 genera were in D1, 1645 genera in D2, and 1692 in D3.Figure 3 shows that in one of the three samples, each major genus accounted for more than 1.00% of the total genera.Bradyrhizobium was the dominant genus among the entire soil community in the L. gmelinii forest, at 21.03% of all genera present in D1, 12.96% in D2 and 20.98% in D3. Apart from the “noname” group, Chthoniobacter (4.26%) and Streptomyces (2.25%) were the two most abundant genera in D1. Chthoniobacter (4.93%) and Acid Candidatus Koribacter (4.52%) were the two most abundant genera in D2,but the relative abundance of Streptomyces was only 0.75%.In D3 samples, Streptomyces (4.32%) and Chthoniobacter(3.75%) were the two most abundant genera (Fig. 3).

According to the cluster tree for phylum and genera in soil samples, species differences between D1 and D3 were small. As seen in Figs. 2 and 3, D1 and D3 were in the same branch close to each other because the seasonal conditions such as low temperature for D1 and D3 were more similar compared with D2.

KEGG function annotation and analysis

To analyze the influence of seasonal changes on microbial functional diversity of L. gmelinii forest soil, GhostKOALA (version 1.0) (Kanehisa et al. 2016) was used to compare the gene set protein sequence with the KEGG database (http://www.kegg.jp), and we obtained the relevant annotations. According to VennDiagram (version 1.6.16)(Chen and Boutros 2011), there were 71 KOs exclusive to the spring soil sample (D1), 98 KOs exclusive to the summer soil sample (D2), 63 KOs exclusive to the autumn soil sample (D3), and 2142 KOs common to all three soil samples (Fig. 4). Moreover, 2451, 2376, and 2437 KOs were identified in D1, D2, and D3, respectively (Fig. 4).These data show that the diversity of pathways, modules,and orthologies in summer was lower than in spring and autumn despite the large number of exclusive KOs and genes (Table 3). Therefore, some KEGG pathways of soil microorganisms might be repressed in summer. Next, the Table 3 further indicates that in the soil microorganisms in the spring (D1), the diversity of pathways, modules, and orthologies were likely to have a corresponding functional increase to adapt to the cold environment, despite a small quantity of exclusive KOs and genes.

Fig. 4 Venn diagrams of KEGG Orthology groups (KOs) in three seasonal soil samples. Different samples are shown in different colors, and the Numbers in the figure represent the number of specific or common components. The similarity and overlap of sample components can be shown intuitively

Table 3 KEGG functional annotation statistics

We def ined the dominant pathways as those with the relative abundances ranking in the top 30 of the total observation pathways (Fig. 5). The pathway with the highest abundances in the three samples was carbon metabolism(PATH: ko01200) (Fig. 5, Table 4), essential for life (Wang and Bai 2006). Citrate cycle (map00020) and carbon fixation pathways in prokaryotes (map00720) were the major metabolic pathways associated with soil microorganisms in the L. gmelinii forest (https ://www.kegg.jp/kegg-bin/show_pathw ay?ko012 00). The total abundances of the top 30 pathways accounted for about 54.57% of the 343 pathways in D1, 53.46% of the 325 pathways in D2, and 54.28%of the 347 pathways in D3. Although the relative abundances of each pathway of the three seasonal soil samples did not differ significantly overall, the relative abundances of ABC transporters (PATH: ko02010) and quorum sensing (PATH:ko02024) in D2 were significantly lower than in D1 and D3. This result indicated that some soil microbial KEGG pathways were activated in spring and autumn, but some were weaker in summer.

To analyze information for all 354 pathways, the 354 pathways were split into 507 various functional modules.Among them, 30 major high-ranking modules were selected according to relative abundance (Fig. 6). M00178 (ribosome, bacteria) and M00179 (ribosome, archaea), the two modules with the highest abundance in the sample were part of the map03010 pathway (Table 5). This indicated that the dominant soil microorganisms in the Larix gmelinii forest were mainly bacteria and archaea. In addition to M00178 and M00179, the relative abundance of M00009(citrate cycle) was also quite high. In the three samples, the sum of the relative abundances of the 30 dominant modules accounted for more than 40% of the all modules. M00178,M00179, M00009, M00237, M00173, M00183, M00011,M00374, M00239, M00036, M00570, M00335, M00019,and M00017 were the dominant modules in the soil microbial metabolism of spring. The total abundances of the 14 major modules accounted for 47% of D1, higher than those of D2 and D3, indicating that the 14 main modules had latent roles in microbial cold resistance. The relative abundances of the other 16 main modules in D2 and D3 were higher than D1, indicating that some major modules were restrained in the cold of spring. The results of KEGG module showed that the relative abundances of some main modules obviously rose during the spring.

CAZy function annotation and analysis

Besides phylogeny and KEGG functional analysis, the degradation of plant organic residues (especially cellulose) has a very important effect on the circulation of soil nutrients.Therefore, we pay special attention to the variation of carbohydrate active enzyme gene abundances. In the CAZy analysis, the most abundant enzymes in order of most to least abundant were glycoside hydrolases (GHs), polysaccharide lyases (PLs), glycosyl transferases (GTs), carbohydrate esterases (CEs), auxiliary activities (AAs), and carbohydrate-binding modules (CBMs) (Fig. 7). The gene numbers for these five groups were significantly higher for D2 than for D1 and D3 and significantly lower for D1 than for D2 and D3. These results indicated that the levels of these enzymes were extremely elevated during the summer and severely inhibited in the spring. Therefore, we speculate that these five groups of carbohydrate-active enzymes could resist high temperature (daily mean maximum temperature > 26 °C) and help soil microbes in Larix gmelinii forest resist high temperature.

Fig. 5 Relative abundances of the top 30 KEGG pathways represented in the three seasonal soil samples (D1-D3). On the left is the similarity tree. The smaller the difference between samples, the samples will be in the same branch with similar color groups information differentiation. The histogram on the right shows the distribution of KEGG pathways in the samples. Different colors represent different KEGG pathways

At the level of the CAZy family, 198 carbohydrate-active enzymes were identified from the metagenomic. Using Circos version 0.69, we chose the top 10 most abundant carbohydrate-active enzymes (Fig. 8). In the three seasonal samples, the collective abundances of the top 10 carbohydrate-active enzymes accounted for more than 50% of the total. GT41 was the most abundant carbohydrate-active enzyme category in these three samples, which is made up of UDP-GlcNAc: peptide β- N-acetylglucosaminyltransferase(EC 2.4.1.255); UDP-Glc: peptide N- β-glucosyltransferase(EC 2.4.1.-). AA3 was the second most abundant, which is made up of cellobiose dehydrogenase (EC 1.1.99.18);glucose 1-oxidase (EC 1.1.3.4); aryl alcohol oxidase (EC 1.1.3.7); alcohol oxidase (EC 1.1.3.13); pyranose oxidase(EC 1.1.3.10) (Zhao et al. 2015). The relative abundances of GT41 accounted for 15.47% in D1, 18.95% in D2, and 13.64% in D3. Apart from AA3 and GH13 in these 10 enzymes, the relative abundances in D2 were greater than that in D1 and D3. At the same time, we found the relative abundances of these eight enzymes all accounted for more than 40% in the three seasonal samples. In addition, most carbohydrate-active enzymes in the soil microorganisms of Larix gmelinii forest were more active in summer than in spring or autumn, and a few enzymes such as AA3 and GH13 were inhibited in summer by high temperature.

Discussion

We expected our metagenomics sequencing to show seasonal variation in the community composition and potential function of soil microorganisms in the Larix gmelinii forest. Previous molecular ecological methods studies reported that in soil microflora libraries on average, Proteobacteria accounted for 39% (from 10 to 77%), Actinobacteria accounted for 13% (from 0 to 34%), Acidobacteria 20%(from 5 to 46%), and Verrucomicrobia 7% (range from 0 to 21%) (Janssen PH 2006). In our study, the relative abundances of Proteobacteria, Actinobacteria, and Verrucomicrobia were in the same range. The gene abundance of these genera obtained by metagenomic sequencing suggested the microbial community composition changed over the threeseasons. For instance, the relative abundance of Streptomyces in D1 was 2.25% and 4.32% in D3, but its relative abundance in D2 was 0.75%; the relative abundance of Chthoniobacter in D1 was 4.26% and 3.75% in D3, but in D2 was 4.93% (Fig. 3). These results revealed that high temperatures in summer raised the relative abundances of Chthoniobacter,and Chthoniobacter might be resistant to heat. Currently,there are few the functional analyses for Chthoniobacter (Shi et al. 2004; Sangwan et al. 2004), so our findings need to be further verified. At the same time, the relative abundances of other soil microorganisms that were not resistant to heat,such as Streptomyces also dropped (Shi et al. 2004).

Table 4 Top 30 KEGG pathways annotations in the three soil samples

The core microorganisms were mainly Alphaproteobacteria, usually living in symbiosis with plants or as plant pathogens (Batut et al. 2004), and among the various metabolic pathways are photosynthesis, nitrogen fixation, ammonia oxidation, and methylotrophy (Williams et al. 2007). Rhodospirillum are photosynthetic bacteria (Gi 2004), Bradyrhizobium are nitrogen fixing bacteria (Chun et al. 1994), Nitrosomonas are ammonia-oxidizing bacteria(Dong et al. 2008), and Methylobacterium can use plantderived methanol as an energy substrate (Abanda-Nkpwatt et al. 2006). All these functions can contribute to growth of the L. gmelinii forest and response to seasonal changes.

The COG annotations were primarily concerned with microbial metabolism (Fig. 1) (Feng et al. 2018); the relative abundances of replication, recombination and repair(L) were the highest among the 25 categories, but the gene numbers of most proteins in D1 were notably lower than in D2 and D3. A decrease in the gene number of most proteins in D1 suggested that the corresponding soil microorganisms in L. gmelinii forest did not possess the proteins of potential tolerance to cold and slowed down or stopped typical metabolism in the cold environment. Among all the 354 pathways of the soil microorganisms, the ko01200 pathway,assigned to carbon metabolism had the highest abundance,not surprisingly, since carbon metabolism is fundamental to life, and soil microorganisms participate in the soil carbon cycle and regulate soil carbon interception capacity, carbon mineralization and ecosystem productivity. The citrate cycle and carbon fixation pathways in prokaryotes are major metabolic pathways associated with soil microorganisms in the L. gmelinii forest. The current research on microbial carbon metabolism showed that most plantderived carbon (C) was lost in the form of CO 2 after passing through microbial metabolic pathways, or finally in the form of stable soil organic matter, and some soil microbes use C substrates for energy production and biosynthesis (Dijkstra et al. 2011). Plant-soil-microorganism interactions largely determine net C sequestration in forest ecosystems (Lange et al. 2015). The main source of forest C emission is from decomposition of soil organic carbon by microorganisms(?if?áková et al. 2016). Soil organic C can be divided into labile C (LC) and recalcitrant C (RC) (Rovira and Vallejo 2002). The richness of carbon metabolis genes in the soil microbial community determines the ability of the community to decompose C substrates (Bardgett and Wh 2014).Soil organic matter (SOM) and C content have the greatest impact on the soil C cycle, which not only directly affects the functional activity of soil microbes, but also indirectly affects their by regulating soil pH and soil enzymatic activities (urease and invertase). And with an increase in SOM,the abundances of genes that decompose RC increase more than those that decompose LC (Yu et al. 2015). We found the relative abundance of the carbon metabolism pathway(PATH: ko01200) in summer was lower than in spring and autumn. In spring, the soil temperature is lower with fewer plant species actively growing; the quantity and activity of soil microbes are also low; thus, the gene numbers of the five soil microbial carbohydrate-active enzymes are lower, but the relative abundances of the dominant modules are mostly higher. The soil C pool is likely to be dominated by LC, so microorganisms mainly decompose the LC at the surface,and the functional abundance of microbial C metabolism rises. Microbial C metabolism functions are also dominant in autumn. In summer, the weather is hot, the forest is dense,and the surface is humid; gene numbers for the five soil microbial carbohydrate-active enzymes are higher, the SOM and organic C content are lower than in other seasons, and RC decomposition gradually becomes dominant, resulting in reduced soil microbial activity and C metabolism. This result is similar to that of Ma et al. (2016). In the spring, the abundance of the most important soil microbial pathways rises: ABC transporters (Smart and Fleming 1996), quorum sensing (Bhardwaj et al. 2013; Ziuzina et al. 2015), and 14 dominant modules such as M00335 (Soderberg et al. 2004)and M00017 (Lu et al. 2017). These metabolic pathways are concerned with the cold resistance of microorganisms(Fig. 5, Table 4). Some KEGG pathways in the soil microbial community are stimulated in the spring and autumn(Smart and Fleming 1996; Lu et al. 2017), but other KEGG pathways such as quorum sensing are inhibited in the summer (Lee et al. 2018). The KEGG analysis showed that some soil microorganisms had potential cold resistance, which provides a new direction for identifying microorganisms with bioremediation potential.

Fig. 6 Relative abundances of the top 30 KEGG modules for the soil samples from D1 to D3 soil depths. See Table 5 for the module annotations. On the left is the similarity tree. The smaller the difference between samples, the samples will be in the same branch with similar color groups information differentiation. The histogram on the right shows the distribution of KEGG modules in the samples. Different colors represent different KEGG modules

Table 5 KEGG modules annotation (top 10) in the three seasonal soil samples

Fig. 7 Functional classification of carbohydrate-active enzymes based on CAZy, glycoside hydrolases (GHs), polysaccharide lyases (PLs), glycosyl transferases (GTs), carbohydrate esterases (CEs), auxiliary activities (AAs), and carbohydratebinding modules (CBMs) The vertical axis shows the full names of the 25 functional categories, and the horizontal axis shows the percentage of genes in the corresponding functional categories

Fig. 8 Circos determination of the levels of CAZy family enzymes in the three seasonal soil samples (D1-D3). Different samples are shown in different colors, and the Numbers in the figure represent the number of specific or common components. The similarity and overlap of sample components can be shown intuitively

At the level of the CAZy family, 198 different carbohydrate-active enzymes were identified from the metagenomic sequencing. The CAZy analysis revealed that most carbohydrate-active enzymes in the soil microorganisms of L. gmelinii forest were more active in summer than in spring or autumn, but severely inhibited in the spring. Carbon metabolism in summer is significantly higher than in other seasons, and significantly lower in winter than in other seasons (Yan et al. 2014). At low temperatures, only a finite amount of the microbiome continues to utilize glucose via both glycolysis and pentose phosphate pathways of carbon metabolism (Apostel et al. 2015). Therefore, we speculate that these five carbohydrate-active enzymes can resist high temperature (daily mean maximum temperature > 26 °C) and help soil microbes tolerate high temperature.

In the metagenomics analysis of the soil samples in different seasons, 96 phyla and 1830 genera were identified among the soil microbes in the L. gmelinii forest, many of which had numerous of potential functions. (Janssen 2006) also reported Proteobacteria, Actinobacteria, Acidobacteria, and Verrucomicrobia were the most abundant phyla in soils (Janssen 2006), and they made up over 70% of the total abundances of every sample in the present study. Taxonomic prof iling showed that the community composition of L. gmelinii forest soil microorganisms apparently change with the season. In addition, the main functional pathways of soil microbes varied as the community composition changed, and many of these functions were concerned with microbial reactions to seasonal changes. These results suggest that seasonal soil should be used for research on cold- or heat-resistant microorganisms.The major group found in the present study was Alphaproteobacteria, belonging to Proteobacteria, with the majority concerned with carbon metabolism and seasonal variation. We also found that in summer, the number of dominant species in the soil declined. In a study based on traditional cultural techniques, when the temperature in summer is high and vegetation grows vigorously, with strong root cell activity, microbial reproduction is also fast and the populations are relatively large (Feng et al. 2009). Such differences in results likely are due to the difference in methodology. Soil microorganisms and some of their functions are critical to the ability of the soil microbial community to adapt to seasonal variations.

Conclusions

The research showed that the variation in seasons caused changes in the soil microbial community structure in the Larix gmelinii forest. The soil microbial main pathways,proteins and CAZy also varied with the seasons, and many were related to cold or heat resistance. The cold and heat resistance of the soil microbes warrants further study in different seasons.