Different mechanisms driving increasing abundance of microbial phosphorus cycling gene groups along an elevational gradient.

Microbes play an integral role in forest soil phosphorus (P) cycling. However, the variation of microbial P-cycling functional genes and their controlling factors in forest soils is unclearly. We used metagenomics to investigate changes in the abundance of genes involved in P-starvation response regulation, P-uptake and transport, and P-solubilization and mineralization along the five elevational gradients. Our results showed the abundance of three P cycling gene groups increasing along the elevational gradient. Acidobacteria and Proteobacteria were the dominant microbial phyla determining the turnover of soil P-solubilization and immobilization. Along the elevational gradient, soil substrates are the major factor explaining variation in P-starvation response regulation genes. Soil environment is the main driver of P-uptake and transport and P-solubilization and mineralization genes. This study provided insights into the regulation of P-cycling from a microbial functional profile perspective, highlighting the importance of substrate and environmental factors for P-cycling genes in forest soils.

Introduction

Forests, accounting for about 30% of the earth’s land area, play a key role in phosphorus (P) cycling. P plays a more important role in forest ecosystem processes than in other ecosystems because a forest has higher productivity, litter fall, and decomposition. (Takyu et al., 2003). Microorganisms drive P cycle by mineralizing organic P, immobilizing inorganic P, synthesizing new organic P, and affecting the solubility of P minerals in terrestrial ecosystems (Willey et al., 2014; Kehler et al., 2021). Such P-cycling processes are mainly mediated by three microbial gene groups, i.e., genes involved in inorganic P-solubilization and organic P-mineralization, P-uptake and transport, and P-starvation response regulation (Dai et al., 2020). Microbes carrying genes involved in inorganic P-solubilization and organic P-mineralization can release organic anions to solubilize inorganic P or expressed enzymes to mineralize organic P (Elias et al., 2001). Microbes carrying genes involved in P-uptake and transport could efficiently utilize P and immobilize P into microbial biomass (Richardson and Simpson, 2011). The genes involved in P-starvation response regulation enable directly connect with P-uptake and transport system, and also indirectly associated with P-solubilization and mineralization (Eder et al., 1996; Helfenstein et al., 2018). However, previous studies were most focus on P-cycle in grassland and agricultural land because P was more directly related with grass and crop growth. In forest ecosystem, soil P-limitation and the P-demand of plant and microbes were totally different from those in other ecosystems. Thus, the microbial gene potential of three P-cycling processes might also be different. That was urgent to clarify the mechanisms of how three microbial gene groups drive P cycle in forest soil.

Environmental factors could alter the soil microbial potential participating in soil P cycle through regulating the expression of different P-cycling functional gene (Vershinina and Znamenskaya, 2002; Penuelas et al., 2020). For example, some research showed that soil N:P stoichiometry could regulate microbial P-solubilization capacity in forest soil (Allison and Vitousek, 2005; Heuck et al., 2015). In addition, recent research found that soil pH significantly regulated P cycle in agro-ecosystems (Dai et al., 2020; Wan et al., 2021a). Although influential factors of microbial P cycle have been detected in various ecosystems, the effect of environmental factor for genes involved in P-cycling, especially in forest, is not clear. Elevational gradients provide various environmental gradients, such as climate, soil properties (e.g., pH), and soil substrate (e.g., carbon and nutrients), which can be particularly useful for quantifying microbial responses to changes in climate and soil characteristics (Feng et al., 2021). For example, recent evidence indicates that phoD genes encoding enzymes involved in organic P-mineralization was higher in low altitudes than high altitudes, and phoD-harboring bacteria showed border environmental thresholds (Wan et al., 2021b). However, we are still unclear how the three groups of P-cycling functional genes respond to elevation-induced environmental changes and mechanisms behind it. Understanding the drivers for three groups of P-cycling functional genes along the elevational gradient is important for understanding and revealing changes in microbial-driven P cycle under climate change.

To fill this knowledge gap, we used metagenomics to explore the changes of P-cycling functional genes and their influence on the P cycle at five sites along an elevational gradient in forests. We hypothesize that soil substrate mostly influences genes involved in P-starvation response, because the P-starvation response is sensitive to P limitation. Likewise, the genes involved in P-uptake and transport affected by microbes, because microorganisms were required to act as carrier. As for P-solubilization and mineralization, soil environment may play a more important role because solubilization and mineralization process are more closely associated with soil physiochemical properties (e.g., pH and BD). The objectives of this study were to: (1) Identify patterns of P-related functional genes in forest soils along the elevational gradient; (2) reveal the primary factors driving the P-cycle related functional genes along the elevational gradient, thereby affecting soil P cycle in forest soils.

Results

Changes of P-cycle genes along the elevational gradient

In this study, 1,323,164,024 million Read sequences from the 15 metagenomes were obtained and identified from 82,193,308 to 100,991,328 sequences per sample, with an average of 88,210,935 in each sample (Table S2). Overall, P-cycling related genes increased along the elevational gradient to different extents. The abundance of genes involved in P-starvation response regulation increased with elevation (Figure1), but only reached the significant level at mid-high elevation (589.47) and high elevation (667.81) sites (p < 0.05) (Figure 1; Table S4). Particularly, phoB, belonging to P-starvation response regulation genes, at high elevation site was 31.51%, 21.40%, 15.05%, and 13.14% higher than at low elevation, low-mid elevation, mid elevation, and mid-high elevation sites, respectively (Figure 2; Table S3). The abundance of genes involved in P-uptake and transport exhibited significant increasing trends with elevation (p < 0.05; Figure 1), with the lowest value (1280.64) at low elevation site and the highest value (1890.94) at high elevation site. In addition, the abundance of genes involved in P-solubilization and mineralization was significantly higher at low-mid elevation (1795.76), mid elevation (1734.94), mid-high elevation (1714.57), and high elevation (1702.58) sites (p < 0.05) than at low elevation (1449.04) site (Figure 1; Table S4), and showed no significant differences among low-mid elevation, mid elevation, mid-high elevation, and high sites (p > 0.05).

Figure 1

Changes in the abundance of phosphorus cycle functional genes along elevational gradient

The 40 phosphorus-cycle functional genes were grouped into 3 main groups (P-starvation response and regulation, P-uptake and transport, P-solubilization and mineralization) and the 3 groups of phosphorus genes varied in elevation. Error bar means ± SD.

Figure 2

Key genes for soil available phosphorus (AP) screened from functional phosphorus cycle genes by random forest analysis

(∗∗∗ when p<0.001, ∗∗ when p<0.05 and ∗ when p<0.01).

Changes of P cycling taxa along the elevational gradient

In this study, 83 phyla microorganisms were annotated from P-functional genes (Table S5). The phyla network showed that Proteobacteria, Acidobacteria, and Deinococcus-Thermus betweenness centrality was 4, 3, and 3, respectively (Figure S3; Table S6). We finally chose Proteobacteria and Acidobacteria as dominant phyla participating in P-cycle, the relative abundance of which was the largest, accounting for 38.17 and 30.89%, on average, of the total phyla, respectively (Table S5). In those two phyla, the abundance of Proteobacteria ranged from 18.96 to 25.60% along the elevational gradient, with the highest abundance detected at low elevation site (25.60%), and the lowest abundance detected at mid elevation site (15.36%). Although the abundance of Acidobacteria ranged from 15.43 to 29.66% along the elevational gradient, with the highest abundance detected at low elevation site (29.66%), and the lowest abundance also detected at mid elevation site (14.74%) (Figure 3). The lowest Acidobacteria: Proteobacteria was observed at low-mid elevation site (0.66), whereas the highest Acidobacteria: Proteobacteria was found at high elevation site (0.91) (Figure 3).

Figure 3

Variation in Proteobacteria and Acidobacteria relative abundance along elevational gradient

Effects of soil environment, soil substrates, and microbial community properties on P-cycle functional genes

Our study showed that genes involved in P-starvation response regulation was most affected by soil substrates (SOC, N: P) (p < 0.001), followed by microbial community (α-diversity) (p < 0.01). Both genes involved in P-uptake and transport and genes involved in P-solubilization and mineralization were significantly influenced by soil environment (BD, SM, pH) (p < 0.01). In total, a large proportion of variation in P-cycle functional genes (86.76% for genes involved in P-starvation response regulation, 92.94% for genes involved in P-uptake and transport, and 63.45% for genes involved in P-solubilization and mineralization) was explained by the variance partitioning analysis (VPA) (Figure5). Soil substrates were the most influential drivers for genes involved in P-starvation response regulation, which explained 24.63% of the variation in genes involved in P-starvation response regulation, followed by microbial community (11.42%), and then soil environment (3.47%) (Figure 5A). Genes involved in P-uptake and transport and P-solubilization and mineralization were predominated by soil environment, with 10.02 and 39.56% of their variation being explained, respectively, followed by soil substrates (0.44 and 4.58%, respectively), and then microbial community (0.61 and 2.47%, respectively) (Figures 5B and 5C).

Figure 5

Proportion of the effect of soil environment (BD, pH, SM), soil substrates (SOC, N: P), and microbial community (α-diversity) on the interpretation of the influence of phosphorus cycling functional genes

VPA analysis was performed on the factors influencing phosphorus cycling functional genes, the numbers between the circles show the intersections of the two types of variables on either side and number in the center of the triangle represents intersections of all three types of variables.

Discussion

Changes of P-related genes abundance along the elevational gradient

Our results showed that the abundance of genes involved in the regulation of P-starvation response, P-uptake and transport, and P-solubilization and mineralization increased with elevation (Figure 1). Our result was inconsistent with previous study indicating the decreasing abundance of P-solubilization and mineralization (only phoD) gene from low to high elevation (Wan et al., 2021b). These discrepancies can be explained by the fact that phoD was significantly positively correlated with soil phosphorus component (Wan et al., 2021b). Our result also confirmed that phoD was positively related to soil available phosphorus (p < 0.05; Figure S3), indicating that the abundance of phoD would increase with the increasing available phosphorus content from low to high elevation. Meanwhile, because the increasing soil available phosphorus has positive correlation with the affinity of transporter proteins (Dai et al., 2020), P-uptake and transport genes (high-affinity pstSCAB and low-affinity pit) also increased with elevation in our study (Figure 1; Figure S4). In addition, phoB (P-starvation response) increased with soil N:P with elevation (Table S1), which was consistent with Dai et al. (2020) suggesting that phoR (also along to P-starvation response) has positive relationship with N:P. Furthermore, previous studies usually focus on single genes involved in P-cycling (Wan et al., 2021b; Dai et al., 2020), whereas P-cycling functional gene groups were considered in this study.

We found that Acidobacteria and Proteobacteria are the key bacterial phyla for P cycle along the elevational gradient in forest soils for their important status and high relative abundance in microbial community (Figure 3). This is consistent with previous findings showing that Actinobacteria, Proteobacteria, and Acidobacteria play a dominant role in the P cycle, with 37.50% genera from Proteobacteria being dominant in oxidative phosphorylation (Ma et al., 2021). Specifically, we found that Bradyrhizobium is one of the major genus from Proteobacteria in this study (relative abundance range of 22.69–11.36% along low-elevation to high-elevation gradient). Previous study showed that Bradyrhizobium harbors phoD gene in steppe soil, which attributed to microbial immobilization by promoting the conversion of the non-labile organic P pool into the labile organic P pool (Zhu et al., 2021). In addition, Acidobacteria were favored in acidic soils (Rousk et al., 2010). In our study, soil pH ranged from 6.06 to 5.92. low pH might provide adequate cations (e.g., Fe3+) for Acidobacteria (e.g., Acidiphilium) making chemotrophic heterotrophic to metabolize (Willey et al., 2014).

In addition, both the relative abundance of Acidobacteria and Proteobacteria were detected highest in low elevation site and lowest in mid elevation (Figure 3). That would be explained by two expects. First, pH value in our study showed highest in low elevation and lowest in mid elevation. As described previously, high pH might preserve less Fe3+ and Al3+, thus maintain more PO43− for microbial use (Willey et al., 2014). Second, the mid-domain effect might induce the lowest abundance of Acidobacteria and Proteobacteria because the microbial diversity would appear peak value in mid elevation (Colwell and Lees, 2000), thereby dampening the dominance of individual species. Moreover, the lowest Acidobacteria: Proteobacteria was observed at low-mid elevation site (0.66), whereas the highest was found at low elevation site (0.91) (Figure 3). That was consistent with that Acidobacteria: Proteobacteria have been suggested as proxies of K-strategists in soil microbial communities (deVries and Shade, 2013; Sun et al., 2021). That was consistent with that the available phosphorus concentration was in relatively low value (mean 0.55 mg/kg<

Factors affecting P-related functional genes in forest soils

Our results indicate that P-solubilization and mineralization genes were strongly affected by soil environment (e.g., pH, SM, and BD) (Figure 4). Consist with our result, Wan et al. (2021a) reported that pH is the key factor in determining organic P-mineralizing-related gene abundance and Ragot et al. (2017) found that soil pH was one of the main factors influencing the composition of the alkaline phosphatase synthesis community in forests. This was likely because soil pH can directly influence P-mineralizing-related microbial growth, thus affecting organic P mineralization (Wan et al., 2021a). This result also confirmed by the VPA analysis (Figure 5C). For example, phoD gene belonging to P-solubilization and mineralization was significantly and negatively correlated with pH (r = −0.769, p < 0.001, Figure S5). Three putative genes of phoA, phoD, and phoX (belonging to P-solubilization and mineralization) in bacterial cells can encode alkaline phosphatase (Maruyama et al., 2016). However, such results contested traditional viewpoints that the regulation of P-solubilization and mineralizing capacity are primarily affected by ‘N:P stoichiometry’ (Heucket al., 2015) The discrepancy may attribute to, first, the abundance of P-solubilization gene responded stronger to changes in pH than N:P stoichiometry. For example, the gcd gene mediates inorganic P solubilization and is responsible for regulating the solubilization of unavailable mineral phosphorus (Neal et al., 2017), which is the most important gene in P-solubilization and mineralization but not determined by nutrient availability (Dai et al., 2020). In this study, the gcd gene was also significantly and negatively correlated with pH (r = −0.817, p < 0.001, Figure S5). Second, soil N:P stoichiometry was one of influencing factors driving the change of soil pH (Rousk et al., 2010), which indicates that pH would be a direct factor, but N:P stoichiometry would be an indirect factor in affecting P-solubilization and mineralization. Third, it might be because of the difference in site conditions among studies. Forest soils were not highly managed in agricultural soils. However, considering soil acidification because of microbial activity and P mineralization coupled with organic carbon use, it is difficult to exclude the possibility that P mineralization leads to lower pH. Thus, this association should be focused in future work.

Figure 4

Direct or indirect effects of elevation gradient (MAT), vegetation characteristics (Tree Shannon), soil environment (BD, pH, SM), soil substrates (SOC, N:P), and microbial community (α-diversity) on functional genes of the phosphorus cycle

PLS-PM analysis was performed on the factors influencing phosphorus cycling functional genes, with the width of the arrows proportional to the strength of the path coefficients. Black arrows indicate positive correlations and red arrows indicate negative correlations. Low transparency indicates non-significance. R2 indicates the proportion of variance explained. Goodness of fit indicates a good model (GOF > 0.7). Significance is indicated by ∗∗∗ when p<0.001, ∗∗ when p<0.05 and ∗ when p<0.01.

Our results showed that soil substrates were the major driver for P-starvation response genes in forest soils (Figure 5A), which was consistent with previous studies (Bergkemper et al., 2016). When P is limiting and organic matter (SOC) substrates are excess, microorganisms will undergo P starvation (Drebot et al., 1990; Zheng et al., 2019). phoR were more abundant in the P-rich soil whereas phoB and phoU have relatively higher abundance in P-depleted soils (Bergkemperet al., 2016). Meanwhile, previous studies have shown that the relative abundance of 35 P-cycling functional genes is significantly or highly correlated with substrates such as SOC, NH4+-N, NO3−-N, and SON (Tang et al., 2020). Our results also showed significant positive correlation between SOC and phoB (r = 0.842, p < 0.001, Figure S5) and between N: P and phoB (r = −0.736, p < 0.01, Figure S5). Moreover, we found that microbial community characteristics also significantly affected P-starvation response genes (Figures 4, Figure 5A). In line with our result, Dai et al. (2020) reported that microbial community has an impact on the P functional genes, which might be because of microbial demands of energy to obtain P by regulating the genes involved in P-starvation response regulation (Hsieh and Wanner, 2010). The influence of microbial diversity and composition in P-cycling genes had been already discovered in P-solubilization and mineralization (Ragotet al., 2017). In addition, our study showed that P-uptake and transport genes was influenced by soil environment (e.g., pH) (Figures 4, Figure 5B), which was also suggested in Dai et al. (2020). The result confirmed that pstA and pH were significantly and negatively correlated (r = −0.871, p < 0.001, Figure S4), whereas pstA and BD were positively correlated (r = −0.661, p < 0.01, Figure S4). This may be because environment variables such as pH, SOC, and enzyme activity regulate encoding transporter proteins and binding components of membranes (Liu et al., 2018), which are mostly responsible for P-uptake and transport genes.

This study used metagenomics to investigate the microbial functional gene regulation of P cycle at five forest sites along an elevational gradient. We found that the abundance of genes from the regulation of P-starvation response, P-solubilization and mineralization, and P-uptake and transport increased with elevation. Further analysis showed that Acidobacteria and Proteobacteria were the dominant microbial groups determining P cycle. Moreover, the genes involved in P-starvation response regulation were significantly affected by soil substrates whereas the genes involved in P-uptake and transport and P-solubilization and mineralization were influenced by soil environment. The present study provides insights into the regulation of P cycle along the elevational gradient from a microbial functional perspective, highlighting different mechanisms for microbial P-cycling functional genes in response to substrate and environmental factors in forest soils.

Limitations of the study

It is essential to highlight that the detailed microbial communities should be more focused on in each elevational gradient so that we can further detect the variation of P-cycling related microbial communities along the elevational gradient. However, because of the difficulty of sampling in forest soil, the lack of numbers of sample lead that we cannot provide a network of microbial communities in each elevation. Thus, further evaluations on single sites should be undertaken in the future, which can provide detailed information about microbial communities along the elevational gradient.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Fazhu Zhao (zhaofazhu@nwu.edu.cn).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

Our study does not use experimental models typical in the life sciences and also has not strains used for in vitro experiments.

Method details

Study area and field sampling

The study area is located in Taibai Mountain (33°45'- 34°10′N and 107°19′-107°58′E, 3767 m above sea level), which is an ideal location to study the variation characteristics of P-cycling functional genes and their influencing factors along an elevation gradient because of Mount Taibai’s steep topography (Figure S1). In this area, the annual mean temperature is 11.4°C and precipitation is 910.6 mm (Zhang et al., 2019). The advantaged tree species were Quercus aliena var. acutiserrata (low elevation) (1503 m), Quercus wutaishanica (low-mid elevation) (1915 m), Betula albosinensis (mid elevation) (2451 m), Abies fargesii Franch (mid-high elevation) (2753 m), Larix chinensis Beissn (high elevation) (3182 m). We randomly selected three sampling points with similar geographic features 20 ×20 m at each elevation gradient, and the method of investigation of vegetation features was described in our previous studies (Zhao et al., 2015, 2020), in which top soil (0–10 cm) was collected and transported by frozen carbon dioxide. Fresh soil was sieved (2 mm) and homogenized to remove fine heel and other plant debris. Approximately 60 g dry weight of soil (on a fresh soil basis) was used to analyze the basic characteristics. The remaining soil (20–30 g) from each composite sample was stored at 4°C for subsequent incubation experiments. Another sample was stored at −80°C and used for DNA analysis. Vegetation characteristics such as tree Shannon and diameter at breast height (DBH) were measured when we set up the experiment. Soil moisture (SM) was measured by drying to constant weight mass at 105°C with the oven (Zhao et al., 2018), and soil bulk weight (BD) was determined by the cutting ring method (Xu et al., 2016). Soil pH was measured in a deionized water suspension (water-soil ratio of 2.5:1) using a pH meter (Mettler Toledo FE28) (Zhang et al., 2016). Soil organic carbon (SOC) was determined by the K2Cr2O7 oxidation method (Vance et al., 1987), total nitrogen (TN) was determined by the Kjeldahl method (Zhao et al., 2019), and total phosphorus (TP) was determined by the colorimetric method (UV spectrophotometer; SHIMADZU UV-2600i/2700i) after wet digestion with HClO4-H2SO4 (Ren et al., 2016), Available phosphorus (AP) was extracted by the NaHCO3 and then measured by UV Spectrometer Subsystem (SHIMADZU UV-2600i/2700i) (Zhao et al., 2020).

DNA extraction, DNA sequencing, and data processing

Soil DNA was extracted according to our previously described method (Renet al., 2016) and quality-checked was determined with a NanoDrop 2000spectrophotomete (Shanghai Personal Biotechnology Co. Ltd: http://www.personalbio.cn/). To obtain sufficient DNA for the shotgun metagenomic sequencing and to guarantee the representation of forest soil, 6 replicates were conducted from each soil sample. The metagenome libraries were sequenced on an Illumina HiSeq 2000 to generate 150bq paired-end reads at greater sequencing depth. The reads aligned to the human genome were removed, and the lengths were trimmed with Sickle. All DNA sequencing can be found on the National Center for Biotechnology Information (NCBI) Website, with the accession number PRJNA780252 (https://www.ncbi.nlm.nih.gov/sra/PRJNA780252).

Metagenomics analysis

The sequencing reads were filtered following Zhang et al. (2017) to enhance the reliability and quality of subsequent analysis. The adapter sequences, trimming the reads, and discarding the quality-trimmer reads which below than 50 bq or containing N (ambiguous bases) were removed. Meanwhile, the megahit software (https://hku-bal.github.io/megabox/) was used to assemble of high-quality reads (Li et al., 2015). The MetaGeneMark (http://exon.gatech.edu/GeneMark/metagenome) was used to predict the genes in the contigs (longer than 200bq), and per-base coverage depth across all contigs was calculated by mapping raw reads from each sample (Allison et al., 2008).

According to the results of KEGG database, the functional annotation and taxonomic assignment from each sample were obtained for further analysis. Based on the definition of P-cycling functional genes in KEGG and previous studies (Dai et al., 2020; Liang et al., 2020), we classified P-cycling functional genes into three categories: P-starvation response regulation, P-uptake and transport, P-solubilization and mineralization. The conceptual diagram showed the function of three P-cycling genes groups (Figure S2). In terms of gene abundance, we determined the trans per million values [TPM: (Reads Number/Gene Length)_Relative] × 100,0000) for each sample. Finally, data and assembly status were showed in Tables S3 and S4. To identify the dominant species participating in P-cycle, we used microbial network to build the structure of microbial community. According betweenness centrality and the relative abundance determine the dominant phyla.

Quantification and statistical analysis

One-way ANOVA was performed using "SPSS" software to assess the effect of elevation gradient on P-cycling functional genes, abundance of key species groups, and soil physiochemical properties. A linear fit analysis was performed using "Origin" software for P-cycling functional genes and soil physiochemical properties such as AP, pH, and SOC at a significance level of 0.05. α-diversity was calculated by the “Shannon” methods in the R package vegan; β-diversity was calculated by Dissimilarity based on Bray-Curtis Distance (Oksanen et al., 2011). P-cycling functional genes were screened by random forest, and random forest analysis was performed with the ‘Random Forest’ package (Liaw and Wiener, 2002) to confirm the key genes for available phosphorus in version R4.0.2. The network of species carried P-cycling genes was performed with ‘Hmisc’ and ‘spearmanCI’ packages. The Chordal graph plotted by the ‘circlize’ package (Gu, 2013). To identify the key environmental factors of P-cycle functional genes, we used Partial least squares pathway model (PLS-PM) to quantify the contributions of soil environment, soil substrates, and microbial community to the abundance of predominant P-cycling functional genes (P-starvation response regulation, P-uptake and transport and P-solubilization and mineralization), which was based on the ‘plspm’ packages. We performed VPA analysis using the ‘ggplot2’ and ‘vegan’ packages to quantify the driving effect of soil physiochemical properties on P-cycling functions.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Biological samples
Forest soil	From Taibai mountain	N/A
Critical commercial assays
OMEGA Soil DNA Kit	OMEGA	Omega Bio-Tek, Norcross, GA, USA; RRID: M5635-02
Deposited data
Metagenomedata of soil microbial community	This paper	PRJNA780252
Software and algorithms
fastp version 0.20.0	Chen et al., 2018	https://github.com/OpenGene/fastp
FLASH version 1.2.7	Magoc and Salzberg, 2011	http://ccb.jhu.edu/software/FLASH/
UPARSE version 7.1	Edgar, 2013	http://drive5.com/uparse
RDP Classifier version 2.2	Wang et al., 2007	http://rdp.cme.msu.edu/
the library sickle	Joshi and Fass, 2011	https://github.com/najoshi/sickle
MEGAHIT version 1.1.2	Li et al., 2015	https://github.com/voutcn/megahit
MetaGene	Noguchi et al., 2006	http://metagene.cb.k.u-tokyo.ac.jp/
CD-HIT	Fu et al., 2012	http://www.bioinformatics.org/cd-hit/
SOAPaligner	Li et al., 2008	http://soap.genomics.org.cn/
BLASTP version 2.2.28+	Altschul et al., 1997	http://blast.ncbi.nlm.nih.gov/Blast.cgi
Mash	Ondov et al., 2016	https://github.com/marbl/mash
Vegan (R package)	Oksanen et al., 2011	=https://CRAN.R-project.org/package=vegan
Random Forest (R package)	Liaw and Wiener, 2002	http://cran.r-project.org/package=Random Forest
circlize (R package)	Gu, 2013	https://CRAN.R-project.org/package=circlize
car (R package)	Fox and Weisberg, 2019	https://CRAN.R-project.org/package=car
stats (R package)	R Core Team, 2020	https://www.R-project.org/
Other
Illumina sequencing	Illumina	MiSeq; HiseqXten