Structural and functional characteristics of soil microbial community in a Pinus massoniana forest at different elevations

Site description

The Leigong Mountain Nature Reserve (26°15′–26°32′ N; 108°5′–108°24′ E) was first established as a Provincial Nature Reserve in 1982 and then promoted to a National Nature Reserve in 2001. Leigong Mountain with 83% forest coverage belong to the typical mid-subtropical mountain forest ecosystem (Zhang et al., 2020). Leigong Mountain belongs to the humid climatic region of the mid-subtropical monsoon mountainous region. The annual rainfall is 1,300∼1,600 mm; the annual average temperature is about 9.2 °C and 16.3 °C at the top and at the foot of the mountain, respectively, with an annual average temperature drop rate of 0.5 °C per 100 m. The soils are mainly acidic mountain yellow soils and mountain yellow brown soils, with a deep soil profile. The zonal vegetation belongs to the humid evergreen broad-leaved forest in the eastern mid-subtropical zone, where it is distributed below 1,350 m a.s.l., while the mountain evergreen deciduous broad-leaved mixed forest is distributed between 1,350 and 2,100 m a.s.l., above which the alpine shrub is observed. Leishan county also is one of the main tea producing areas in Guizhou province (the organic tea area is 653.2 hm², and the annual tea products is 3,950 t), where tea garden are typical artificial shrub vegetation in mountain ecosystem.

Study site and soil sampling

In June 2019, four sites were selected along the altitudinal gradient established on the P. massoniana forest in the Leigong Mountain at 1,000, 1,200, 1,400, and 1,600 m a.s.l. and named Mp1000, Mp1200, Mp1400, and Mp1600, respectively. An artificial tea garden (Tea-garden) was selected as the control plot at 1,000 m a.s.l. (Fig. 1). One vegetation survey plot with 20 m × 20 m in size was also set at each altitude. Information on site conditions (Table 1), plant communities, and litter was obtained for each plot. Soil samples were excavated from the 0–20 cm topsoil layer at six points located along an “S” path drawn within each plot, and thoroughly mixed into one composite sample. Five soil samples were collected (one from each survey plot plus one from the control plot). For each sample, 500 g of soil was sieved though a 2-mm mesh and air-dried for physicochemical analyses; additionally, 50 g of soil was encapsulated in a plastic tube and kept at −80 °C for microbial DNA extraction.

Soil analysis and microbial DNA extraction

Soil pH was determined at a ratio of 1:2.5 (soil to water, w/w). Soil organic carbon (SOC) was determined using the K₂Cr₂O₇ titration method (Walkley–Black method). Total nitrogen (TN) was measured using the semi-micro Kjeldahl method. Total phosphorus (TP) was determined colorimetrically after wet digestion with H₂SO₄ + HClO₄ (Xu et al., 2014).

Microbial DNA extraction and metagenome sequencing

Total genomic DNA was isolated from 0.50 g of soil using the MoBio Ultraclean Soil DNA Isolation Kit (MoBio Laboratories, Carlsbad, CA, USA) following the instructions of the manufacturer. Extracted DNA was stored at −80 °C until use. Tenfold diluted DNA samples were then used for metagenome analysis on the Illumina MiSeq platform (Illumina Inc., San Diego, CA, USA) according to the standard protocols of Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China).

The MiSeq platform was used to perform parallel mixing and sequencing of DNA multiple samples. The sequence of each sample was introduced with an index tag sequence indicating sample source information. Software “Seqprep” was used to cut the 3′ end and 5′ adapter sequence at the end. Software “Sickle” was used to remove sequences less than 50 bp; the average quality value was below a certain threshold (default 20) and N-base reads; high-quality pair-end reads and single-end reads were retained (Pabinger et al., 2014; Schirmer et al., 2015). We then used “fastp” (https://github.com/OpenGene/fastp) software to perform quality control processing on the original sequencing data to obtain high-quality control data (clean data). Software “Megahit” (https://github.com/voutcn/megahit) was used to assemble sequences with different sequencing depths (Li et al., 2015) and to preserve sequences with assembly contigs >300 bp. MetaGene (http://metagene.cb.ku-tokyo.ac.jp/) software was then used to predict open reading frames (ORF) of contigs ≥100 bp (gene sequences) in clean sequences (Noguchi, Park & Takagi, 2006). Using CD-HIT (http://www.bioinformatics.org/cd-hit/) software, clustering (identity ≥90% and coverage ≥90%) was performed on the predicted gene sequences of all samples, taking the longest gene in each class as a representative sequence to construct a non-redundant gene set (Fu et al., 2012). SOAPaligner software was used to compare the high-quality reads of each sample with the non-redundant gene set (identity ≥95%) and to obtain gene abundance information in the corresponding samples (Luo et al., 2012). For species identification, BLASTp (http://blast.ncbi.nlm.nih.gov/Blast.cgi) software was first used to blast the obtained gene sequences against the non-redundant (NR) protein sequence database (Altschul et al., 1997) and the evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG) database (http://eggnog.embl.de/) with the expected e-value for the functional annotation set to 1e−5 (Jensen et al., 2008). The Kyoto Encyclopedia of Genes and Genomes (KEGG) Database and an e-value of 1e−5 were used for functional annotation (Du et al., 2014). The Carbohydrate-active enzymes (CAZy) database (http://www.cazy.org/, Version 5.0) and an e-value of 1e−5 were used for active carbohydrate enzymes functional annotation (Lombard et al., 2014). Fungal guilds were assigned using the FUNGuild platform (http://www.funguild.org/) (Nguyen et al., 2016).

The number of total and unique species and the functional genes of Archaea, Bacteria, Fungi, and Virus groups in the different soil samples were determined and, using the VennDiagram R package (https://www.r-project.org) a Venn diagram was produced for the numbers of Archaea, Bacteria, Fungi, and Viruses in the soil microbial community. A horizontal community histogram, also produced in R using cluster analysis, was based on the similarity of soil microbial community between samples. Similarity results of the soil microbial community among different soil samples are presented as a heatmap.

The correlation between the soil microbial community matrix and the function matrix among the different samples was also tested. The soil microbial community species-level matrix based on the Bray-Curtis distance algorithm was used for redundancy analysis (RDA). The “networks” software was used for co-occurrence network analysis (Barberan et al., 2012). Soil environmental factors were analysed using One-way analysis of variance (ANOVA) in SPSS 22.5 (https://www.ibm.com/analytics/spss-statistics-software). Finally, Sigma Plot 12.5 and Adobe Illustrator 2020 (https://www.adobe.com/products/illustrator.html) software were used for figure drawing.

Composition and structure of soil microbial communities

A total of 452,526,750 original sequences were obtained from the metagenomic sequencing of soil microorganisms in the P. massoniana forest along the altitudinal gradient established. After quality control, 447,074,316 high-quality sequences were generated and then divided into different soil microbial groups, including Archaea, Bacteria, Fungi, and Viruses. There were 62,208,980 high-quality sequences and each sample produced 9,960,236~16,161,802 high-quality sequences. After metagenomic sequencing, BALSTp on the NR database allowed identifying 13,307 species of soil microorganisms. These included 912 Archaea, 11,576 Bacteria, 417 Fungi, and 402 Viruses (Table 2).

The changes in richness and in Shannon’s index of soil microbial groups at the different altitudes selected at the study site showed different trends (Table 3). In terms of species richness, the abundances of soil Bacteria, Fungi, and Viruses exhibited a “unimodal” pattern with altitude, while Archaea showed a major peak at Mp1200 (752) and a secondary peak at Mp1600 (716). With increasing altitude, Shannon’s index values increased for Archaea and Fungi but decreased for Bacteria; for Viruses, no obvious pattern was found.

The distributions of exclusive and common species of the different soil microbial groups showed obviously different profiles (Fig. 2). Common species of Archaea, Bacteria, Fungi, and Viruses accounted for 58.6%, 75.1%, 50.4%, and 7.5% of the soil microbial groups, respectively, whereas exclusive species of Archaea, Bacteria, Fungi, and Viruses accounted for 12.3%, 6.6%, 15.6%, and 54.5%, respectively. The spatial distribution patterns of the different soil microbial groups also differed across the altitudinal gradient. The ratio of exclusive and common species for Archaea and Fungi were similar, but with obvious differences of the ratio for Bacteria or Viruses.

Euryarchaeota and Thaumarchaeota were the predominant groups of soil Archaea in the P. massoniana forest surveyed across the selected altitudinal gradient. The order of relative abundance of these two groups at the five sites was: Tea-garden (84.1%) > Mp1000 (80.3%) > Mp1200 (75.8%) > Mp1400 (75.6%) > Mp1600 (75.5%), thus showing a decreasing trend from 1,000 to 1,200 m altitude. Euryarchaeota and Thaumarchaeota were followed in relative abundance by Candidatus Bathyarchaeota, Crenarchaeota, Candidatus Thorarchaeota, Candidatus Lokiarchaeota, Candidatus Micrarchaeota, Candidatus Korarchaeota, Candidatus Parvarchaeota, Candidatus Nanohaloarchaeota, and Nanoarchaeota (Fig. 3A).

Proteobacteria and Acidobacteria were the predominant groups of Bacteria. The order of relative abundance of these two bacterial groups at the five sites was: Mp1000 (73.1%) < Tea-garden (74.0%) < Mp1200 (74.3%) < Mp1600 (76.0%) < Mp1400 (76.2%), thus showing an increasing trend with increasing altitude. The bacterial groups following in relative abundance were Actinobacteria, Verrucomicrobia, Chloroflexi, Planctomycetes, Cyanobacteria, Firmicutes, Bacteroidetes, Gemmatimonadetes, and other bacteria (Fig. 3B).

Ascomycota and Basidiomycota were the predominant groups of Fungi. The order of relative abundance of these two fungal groups at the five sites was: Tea-garden (89.4%) < Mp1000 (92.2%) < Mp1200 (93.3%) < Mp1400 (94.2%) < Mp1600 (94.4%), thus showing an increasing trend with increasing altitude. Moreover, ascomycetes increased with elevation and their relative abundance decreased significantly whereas basidiomycetes showed an obvious increase with elevation. The relative abundances of the two groups showed increasing trends with increasing altitude. The fungal groups following in relative abundance were Chytridiomycota, Glomeromycota, other Fungi, and unclassified Fungi (Fig. 3C).

Myoviridae, Podoviridae, and Siphoviridae were the predominant virus groups. The order of relative abundance of these groups at the five sites was: Mp1600 (42.1%) < Tea-garden (48.7%) < Mp1000 (49.6%) < Mp1200 (50.8%) < Mp1400 (54.5%). The groups following in relative abundance were Phycodnaviridae, Tectiviridae, Mimiviridae, unclassified virus, Caudovirales, and other viruses (Fig. 3D).

Clustering of the different soil microbial groups (Archaea, Bacteria, Fungi, and Viruses) at the genus level showed significant differences with altitude. Overall, Candidatus_Nitrosotalea, unclassified_p_Candidatus_Bathyarchaeota (from phylum to genus), Methanosarcina, Nitrososphaera, and unclassified_d_Archaea (from domain to genus) were the predominant archaeal genera (Fig. 4A), while unclassified_f_Acidobacteriaceae (from family to genus), Candidatus_Solibacter, Candidatus_Koribacter, and Bradyrhizobium were the predominant bacterial genera (Fig. 4B). Moniliophthora, Penicillium, Oidiodendron, Colletotrichum, Aspergillus, and Cladophialophora were the predominant fungal genera (Fig. 4C). Unclassified_d_Viruses (from domain to genus), unclassified_f_Myoviridae (from family to genus), unclassified_f_Siphoviridae (from family to genus), and unclassified_f_Podoviridae (from family to genus) were the predominant virus genera (Fig. 4D).

Functional gene composition and structure of soil microbial communities

Based on the eggNOG database, the soil microbial community functional annotation was determined for the P. massoniana forest at different altitudes (Fig. 5). The main groups of functional genes were related to amino acid transport and metabolism (7.51–7.98%), energy production and transformation (6.62–7.46%), signal transduction mechanisms (6.76–7.24%), and DNA replication, recombination, and repair (6.46–7.16%). These dominant functional gene categories are important for energy production and cell metabolism. Genes with unknown functions accounted for 30% (28.58–31.16%) with the most abundant of all functional genes classified.

A high representation of carbohydrate metabolism (12.22–13.00%), global and overview map (11.51–11.62%), amino acid metabolism (10.23–10.67%), and energy metabolism (7.51–7.80%) was observed when the KEGG database was applied to the functional annotation into putative pathways (Fig. 6A). Glycosyl transferase (GT, 34.59–38.58%) and glycoside hydrolase (GH, 28.98–33.59%) were most represented carbohydrate active enzymes (CAZy), followed by carbohydrate esterase (CE, 15.91–16.28%), coenzyme activity (AA, 10.43–10.92%), and polysaccharide lyase (PL, 3.00–3.64%), combined with the carbohydrate Module (CBM, 1.98–2.14%) (Fig. 6B).

Saprophytic and parasitic trophic mode made up the two largest proportions in soil fungal communities in the P. massoniana forest; while Ectomycorrhizal fungi and Ericoid mycorrhizal fungi were dominant groups in the symbiotic trophic mode. The relative abundances of different functional groups in soil fungal communities in the P. massoniana forest showed varying trends with increasing altitudes (Fig. 7), such as undefined saprotroph, plant pathogen, animal pathogen-undefined saprotroph, and ericoid mycorrhizal fungi showed a decreasing trend with increasing altitudes, while wood saprotroph, plant pathogen-undefined saprotroph, and Ectomycorrhizal fungi showed an increasing trend with increasing altitudes (Fig. 7). These results suggest that specific functional group in soil fungal communities may be better indicators of forests in terms of elevation changes.

Similarity of species and functions among soil microbial communities

The similarity and dissimilarity of species as per NR, NOG, KO, and CAZy functional genes among all sample sites were compared to results obtained from similarity clustering analysis and RDA (Fig. 8). The similarity of species and functions of soil microbial communities along the altitudinal gradient in the P. massoniana forest was as described above, i.e., the RDA ordination of soil microbial communities at different altitudes based on functional gene composition showed a good level of consistency with the ordination based on species identification. The distribution of identified species, NOG, KO, and CAZy functional genes of soil microbial community at MP1000, MP1400, and MP1600 were clearly separated from those observed at MP1200 (Fig. 8B). Furthermore, a significant difference was observed between the P. massoniana forest and the tea garden surveyed.

Co-occurrence patterns of the different soil microbial groups

In order to understand the correlation between species of the different soil microbial communities, we detected significant differences among different soil microbial groups (Archaea, Bacteria, Fungi, and Viruses) and their corresponding random networks (Fig. 9) in terms of average degree, modularity, number of communities, and degree centrality were found (Table 4). Modularity ranked as follows: Archaea (0.866) > Bacteria (0.794) > Fungi (0.698) > Viruses (0.539), while degree centrality showed the following ranking: Bacteria (0.181) > Viruses (0.155) > Fungi (0.087) > Archaea (0.048) (Table 4). Overall, the co-occurrence networks of archaeal and fungal species showed a relatively high similarity, while the co-occurrence network of bacterial and viral species showed a relatively high similarity.