The exacerbation of soil acidification correlates with structural and functional evolution of the soil microbiome upon long-term agricultural intensification

Background: Agricultural intensification induces prolonged impacts on soil acidification (SA), whereas linking the mechanisms therein to structural and functional evolution of soil microbiota remain unclear, especially over longer time scales. Methods: To investigate the short- and long-term effects of agricultural intensification on soil microbiome and its association with SA severity, we performed metagenomic sequencing of soil samples from typical rice-vegetable rotations for 0 (control, V0), 10 (V10), 20 (V20) years. The objectives were to clarify biomarkers that indicative of SA and its putative functions of acid production. Results: Rice-vegetable rotations for 0, 10, 20 years yielded three well-defined soil clusters with differential acid saturation (V0: 4.5%; V10: 11.8%; V20: 40.7%). Soil pH declined significantly (p < 0.05) from the 6.38 (V0) to 4.82 (V20). Acid cations (H + , Al 3+ ) dominating V20’s exchangeable cation pool, suggesting that soils in V20 were highly acidic and acid-sensitive, which recruited distinct microbiomes. The increased acid cations and NO 3 - -N concentrations are driving forces that lead to a higher abundance of Rhodanobacter , Gemmatirosa , Sphingomonas and Streptomyces in V20 samples, which acting as aciduric biomarkers that significantly and positively correlated with soil acidity. Furthermore, functional modules of microbial genes revealed V20 samples contained more genes associated with“oxidative phosphorylation” and “two-component system”, particularly enriched multiple pathways of cytoplasmic pH homeostasis. More specifically, we found increased genetic abundances of active efflux of protons and cytoplasmic proton consumption in V20 samples over than in V0, which are both involved in the microbial dehydrogenation and acid production, and favoured the growth of Rhodanbacter , Gemmatirosa , Sphingomonas , Streptomyces . Conclusions: Collectively, the identified aciduric biomarkers was closely associated with SA severity in intensive agroecosystems. An overall increase with the ongoing insights the exacerbation of SA.

was the key reaction of acid production process [22][23][24][25]. Recent studies on aciduric ammoniaoxidizing archaea (AOA) have promoted scientific interest and discussion on mechanisms of agricultural soil acidification. Among which, Candidatus Nitrosotalea and Nitrosopumilus have been confirmed as dominant microbiotas that are responsible for microbial dehydrogenation in the wider range of acidified cropland, grassland, and forest soils, leading to further H + loads and aluminum toxicity with the leaching losses of base cations [26][27][28]. In addition, variations in microbial communities largely originate from the inherent differences in microbial ability to cope with extreme environmental stressors [29]. Most archaea have functions maintaining their cytoplasmic pH within a narrower range than the pH outside the cell, termed "pH homeostasis", while pH homeostasis in bacteria may be specific to particular species under acidic stressors, thereby bacteria tend to be more be more responsive than archaea to acidic stressors [30]. Previous studies have suggested that acidification independently caused by several different pathogenic bacteria, including Helicobacter pylori, Escherichia coli, Streptococcus mutans [31,32], while in intensive agricultural systems, the information of exact aciduric biomarkers and its roles in acid production process remain largely unexplored. Consequently, soil microbial adaptation versus selection may result in a highly enrichment of aciduric biomarkers, which promotes microbial dehydrogenation and development of soil acidification, and therefore the soil with abundant aciduric biomarkers would be expected to have an increased risk of developing acidification compared to those within non-aciduric biomarkers.
More recently in metagenomic profiling, the availability of dominant soil microbial complete or nearcomplete genomes and their functions make it possible to disentangle the association of soil microbiome variations with changes in soil acidity, and to identify indicative functions responsible for these variations [17,30]. Pioneer studies revealed the potential roles of pH homeostasis genes in microbial dehydrogenation and acid production [32]. Prolonged periods of low pH stressors favour the growth of aciduric bacteria, leading to an increase of bacterial pH homeostasis genes [22,33]. A major strategy for bacterial pH homeostasis is the use of transporters that catalyze active proton transport, including primary proton pumps like proton-pumping respiratory chain complexes and nicotinamide adenine dinucleotide phosphate (NADH) dehydrogenases, and secondary active transporters like cation-proton antiporters, which energize active proton uptake in exchange for cytoplasmic cations such as Na + or K + [32]. Another pH homeostasis strategy is related to metabolic consumption of cytoplasmic protons, particularly the catabolism of amino acids generates alkaline amines by decarboxylases, such as the lysine and arginine decarboxylases [34,35]. Moreover, Tomb [31] and Chi A [36] et al. suggested proteomes of Acidithiobacillus ferrooxidans and Helicobacter pyloria could provide a passive adjunct to the active mechanisms for pH homeostasis, which formed a transient proton repellent to counteract external acidification. Overall, diverse mechanisms for pH homeostasis driven by functional genes impact microbial adaptation versus selection under acidic stressors. The multiple dehydrogenation pathways that mediating pH homeostasis are usually accompanied by the enrichment of environmental protons, thus the variation of acid production may be linked to diverse pH homeostasis gene-rich soil microbiomes.
Currently, urbanization has been increasing worldwide with rapid transition of populations from rural to urban, which means sharply increasing demands for grain and vegetables [37], and causing large amounts of conventional paddy fields have been converted from low nitrogen input "cereal cultivation" to high nitrogen input "vegetable cultivation" [38,39]. These conversions have inevitably entailed extensive soil acidification and changes in microbial communities. Our previous studies have confirmed that intensive rotations would aggravate nitrogen enrichment and soil acidification [40,41]. However, the prolonged effects of these intensive practices on soil microbial communities, as well as their contributions to soil acid production have not yet received enough attention.
We hypothesize that the exacerbation of soil acidification related to the decades-scale response of soil microbial community composition and its functional attributes following agricultural intensification in soils. Therefore, herein performed a metagenomic study of the soil microbiome in different intensive cropping history to thoroughly investigate the microbial components associated with the soil acidification. The objectives were to (i) clarify the temporal variation in soil microbial composition, abundance and functions upon agricultural intensification, especially its association with soil acidity, (ii) determine the aciduric biomarkers and indicative gene modules that are closely associated with acid production pathways, (iii) explain the interactions between the soil microbiome and the exacerbation of soil acidification.

Results
Altered soil acidity and high-order taxonomic composition upon short-versus long-term agricultural intensification A highly soil acidification pH < 5 was found in long-term rice-vegetable rotations (V20) ( Table 1). As shown in Fig. 1, soils in short-(V10) versus long-term rice-vegetable rotations were at acid-sensitive stage. Acid cations saturation ratios (proportion of H + , Al 3+ ) were increased significantly with the ongoing of intensive cropping history (p < 0.05), which rose from 4.5% (V0) to 11.8% (V10) and 40.7% (V20). Concentrations of H + and Al 3+ in V20 were 3.0 and 7.6 times higher than V0, respectively. Furthermore, soil base cations saturation ratios (proportion of Ca 2+ , Mg 2+ , K + , Na + ) reached the amount to 95.5% in control treatment (V0), among which Ca 2+ and Mg 2+ dominated exchangeable soil cations pool and accounted for approximate 82.1% and 9.4%, respectively. By comparison, base cations saturation ratios significantly decreased to 88.2% in V10 and by 59.3% in V20.
To identify whether soil acidity can explain changes in microbial structure and functions among different intensive cropping history, we performed shotgun metagenomic sequencing on soil samples.
A total of 549 million 150-bp paired-end reads were generated, with an average of 60.98 ± 2.73 (s.e.) million reads for each sample (Additional file 1: Table S1). After quality control, we obtained 543 million high-quality clean reads, with an average of 60.34 ± 2.69 (s.e.) million reads per sample (Additional file 1: Table S2). The clean reads were assembled, predicted and then aligned to the reference genomes from the National Center for Biotechnology Information (NCBI).
We found that soil microbial communities were disturbed after converting from conventional paddy (V0) to rice-vegetable rotations (V10 and V20) (Fig. 2a). Krona charts indicated that these changes in taxonomic composition were evident in high-order levels, particularly at class rank, which community structure shifted from a more balanced cluster comprising Betaproteobacteria, Deltaproteobacteria and Alphaproteobacteria in V0 to a less-balanced community comprising Alphaproteobacteria, Actinobacteria and Grammaproteobacteria dominated community in V10 and V20 (Fig. 2a, b). One-way ANOVA suggested, Betaproteobacteria, Deltaproteobacteria and Nitrospira were more abundant in V0 than in V10 and V20 (p < 0.05), whereas the relative abundances of Grammaproteobacteria, Gemmatimonadetes, Acidobacteriia, Ktedonobacteria and Sphingobacteriia were significantly increased in V20 than V0 (p < 0.05; Fig. 2c).
Shifts in fine-scale taxonomic composition and determine the aciduric biomarkers indicative of soil acidification Furthermore, it is interesting to notice differentially fine-scale taxonomic compositions (Genus level) of soil microbiome in distinct intensive cropping history ( Fig. 3 and Fig. 4a). V10 and V20 samples shared a similar microbial composition but showed a greater distance with V0 samples in both metagenome and 16S rRNA gene (Fig. 3a, c). Unconstrained principal coordinates analysis (PCoA) of Bray-Curtis distance revealed that microbial communities of V10 and V20 samples formed significantly distinct clusters with V0 in both metagenome (ANOSIM: statistic = 0.9835, p < 0.05) and 16S rRNA gene (ANOSIM: statistic = 1, p < 0.05), which separated along the first coordinate axis ( Fig. 3b, d).
To characterize the core microbial composition and abundance, heat map analysis was performed with the top 50 most abundant genera (Fig. 4a). Increased abundance was evident at some dominant genera (e.g. Gemmatimonas, Bradyrhizobimu, Solirubrobacter, Gemmatirosa, Mycobacterium, Rhodanobacter, Arthrobacter and Preudolabrys, etc) of V10 and V20. The largest increase in relative abundance response to distinct intensive cropping history was by the genus Rhodanobacter, with the abundance increasing from an average of 0.06% [± 0.0086, s.e.] in V0 to 4.12% (± 0.6263, s.e.) of those in V20, followed by the genera Gemmatirosa, Gemmatimonas and Sphingomonas. The largest decrease was by the genus Nitrospira, with the abundance decreasing from an average of 1.65% [± 0.1325, s.e; V0] to 0.93 [± 0.0661, s.e; V20], followed by the genus Candidatus_Entotheonella and Geobacter.
Next, we analyzed whether these dominant members can be used as biomarkers to differentiate V0, V10 and V20. LDA scores (log 10) of 3.5 or greater are listed (Fig. 4b) to identify biomarkers for V0 and V20. The results showed that Nitrospira and Candidatus_Entotheonella were significantly enriched in V0, whereas Rhodanobacter, Gemmatirosa, Sphingomonas, Streptomyces, Haliangium, Candidatus_Solibacter and Thermogemmatispora showed higher relative abundance in V20 than V0.
Furthermore, SparCC network analysis was constructed to investigate the possible interactions between the biomarkers. Two clusters were apparent from the network (Fig. 4c)

Identifying differentially abundant KEGG modules in distinct intensive cropping history
To evaluated whether the shifts in soil microbial composition and abundance also had functional consequences concerning genes, we annotated the functions of genes that were specifically enriched in V0, V10 and V20 to Kyoto Encyclopedia of Genes and Genomes Orthology (KEGG) Module database.
LEfSe analysis was performed and LDA scores (log 10) of 2 or greater are listed to identify indicative modules of V0, V10 and V20 respectively, which were related to 94 modules (Additional file 2: Table  S3 ~ S5). Notably, most modules enriched in V0 samples were involved in ABC transporters, carbon metabolism, whereas modules enriched in V20 samples were related to oxidative phosphorylation, two-component system, DNA replication and nitrogen metabolism.
Second, using one-way ANOVA, we validated that modules were differentially enriched among sample groups (Tukey-Kramer test) (Fig. 5b). We found that gene clusters associated with functions supporting ABC transporters, e.g. branched-chain amino acid- ( To investigate the role of gene-annotated modules in the recruitment of soil microbiome, we carried out species and functional contribution analyses targeting the above-mentioned nine acid-sensing biomarkers ( Fig. 4b) and eleven pH homeostasis modules (Fig. 6a). Regression analysis of Bray-Curtis distance in pH homeostasis modules and biomarkers community composition revealed that soil acidsensing microbiomes of V20 samples were separated from V0 and V10 along the first two coordinate axes. Notably, these gene modules provided the largest contribution to Rhodanobacter (total abundance for eleven modules, Additional file 3: Table S6) and contributed (One-way ANOVA, p <

Discussion
Soil acidification is a major environmental issue of agricultural intensification, which is an acid production process driven by proton-generating reactions, such as microbial nitrification, dissociation of carbonic and other organic acids, accumulation of organic matter, and acid deposition, etc [9,42].
When proton loads are elevated, base cations are exchanged out from soil sorption complex to buffer the input of acids [42]. As the leaching loss of base cations, soil pH becomes buffered by aluminum (Al 3+ ) mobilized from soils [9]. Although the practices of agricultural intensification have been proven to reduce the soil pH in our previous findings [41] as well as the general pattern of acidification observed in many other studies [13,43], there has yet been no systematic study of the extent of soil acidification. Herein, we tested the impact of short-versus long-term rice-vegetable rotations on the composition pattern of cation pool in paddy soil. Figure 1 showed that soil exchangeable cation pool was dominated by acid cations (H + , Al 3+ ) with pH dropped to 4.82 (Table 1) upon long-term ricevegetable rotations. More importantly, we further found that long-term rice-vegetable rotations increased the magnitude of pH decline, where soils were at Al 3+ -buffering stage with highly acidified and acid-sensitive. Overall, rice-vegetable rotations treatments, especially the V20 aggravated the soil acidification and acid-sensitivity. Under the gradient of acid saturation, soils in V0 could be characterized as base-enriched taxa (4.5%), and in V10 (11.8%) and V20 (40.7%) could be identified as acid-enriched taxa.
Given that these declines in soil pH may alter microbial biodiversity, and are intimately to soil functioning and in turn agroecosystems integrity [16,17], we argue that changes of microbial biodiversity should be taken into account when assessing the impact of agricultural intensification on soil acidification. We revealed that base-enriched (V0) and acid-enriched (V10 and V20) soils recruited distinct microbiomes (Fig. 2, Fig. 3 and Fig. 4a). Gammaproteobacteria, Gemmatimonadetes and Acidobacteriia dominated the high-order taxonomic composition (class level) of V20 with dramatically declining in abundance of Betaproteobacteria and Deltaproteobacteria than V0, suggesting certain microbiota have the unique selectivity and thriving in acid-enriched soils. Microbiotas in intensive agricultural soils have been reported to be high sensitive to proton loads and subsequent soil acidification, and high acid saturation caused strong selective pressures on microbial community [44].
Selection minimizes microbial diversity and leads to the survival of dominant species, finally assembles the unique community structure, some of which shape biomarkers to adapt specific soil acidity [45,46]. By LEfSe analysis (Fig. 4b) at fine-scale taxonomic composition (genus level), we found nine biomarkers: Nitrospira, Candidatus_Entotheonella, Rhodanobacter, Gemmatirosa, Sphingomonas, Streptomyces, Haliangium, Candidatus_Solibacter and Thermogemmatispora. Further SparCC network and RDA analysis (Fig. 4c, d) showed a strong influence of the acidity parameters on dominant microbial structure and distribution, which may explain the largest source of variation in soil microbiomes was differentially soil acidity along the intensive cropping history. In addition, we revealed potential aciduric characteristics of Rhodanobacter, Gemmatirosa, Sphingomonas, Candidatus_Solibacter, Streptomyces and potential alkaliphilic characteristics of Nitrospira, Candidatus_Entotheonella. Intriguingly, Rhodanobacter, Gemmatirosa, Sphingomonas, Candidatus_Solibacter, Streptomyces, Nitrospira and Candidatus_Entotheonella were abovementioned biomarkers (Fig. 4b, c), thus they have potential for serving as acid-sensing biomarkers.
Actually, previous studies have demonstrated the acid tolerance of Rhodanobacter [47], Sphingomonas [48], Candidatus_Solibacter [49], Streptomyces [50]. A similar enrichment has been observed in acid-enriched soils of our study cohort as well, particularly for the dominant Rhodanobacter (Fig. 4a, b). One should note that Rhodanobacter belongs to the family Xanthomonadaceae, order Xanthomonadales and class Gammaproteobacteria of phylum Proteobacteria [51], and strains of which were isolated from a nuclear legacy waste site where cocontaminated with large amounts of acids, nitrate, metal radionuclides and other heavy metals, and have been found highly abundant and active in acidic, nitrate rich environments and high metal (e.g. uranium) concentrations [52]. Moreover, similar conclusion has been raised by highly abundant of Rhodanobacter in acidified agricultural soils, at pH near 4 [53], which is consistent with our findings that Rhodanobacter showed the survival potential in acid-enriched soils (pH 4.82). In our findings, the high correlation of acid cations and acid-sensing biomarkers further highlighted Rhodanobacter, Gemmatirosa, Sphingomonas, Candidatus_Solibacter and Streptomyces's essential role in serving as aciduric biomarkers, and call for further necessity in exploring whether aciduric biomarkers are the causal conducers to dehydrogenation and acid production.
Considering that environment and community structure drive the variations in ecosystem functions [54], any functional evolution that increases resistance represents an important strategy that may benefit stressful environment for microbiome. Concomitant with the alteration of soil acidity and microbial composition, we observed a functional evolution in microbiome. KEGG Module annotations with LEfSe analysis revealed that 94 KEGG modules were differentially enriched among the three groups (Additional file 2: Table S6), 15 and 19 modules that specifically enriched in V0 and V20 were listed in Fig. 5. Functional genes associated with ABC transporters, such as branched-chain amino acid, iron, tungstate and phosphonate, general L-amino acid, manganese/zinc/iron transport system were dramatically repressed in V20 samples. In contrast, oxidative stress, signal transduction and DNA repair genes were found actively enriched in V20, including oxidative phosphorylation, twocomponent system, DNA replication and nitrogen metabolism, etc. Moreover, spearman correlation analysis revealed the concentrations of H + , Al 3+ and NO 3 − -N were positively correlated with abundance of 11 modules (Fig. 6a). Intriguingly, six out of eleven modules were above-mentioned ( Fig. 5b) dehydrogenation pathways, including proton-pumping: NAD biosynthesis (M00115), cytochrome complex oxidase (M00151, M00417), pentose phosphate pathway (M0004) and cytoplasmic proton consumption: GABA biosynthesis (M00135), assimilatory nitrate reduction (M00531), suggesting that these modules might be involved in pH homeostasis of the soil microbiome in distinct acid saturated soils. Notably, these dehydrogenation pathways acting as a H + consumer and transporter, which have been reported to drive the difference in microbial pH homeostasis of distinct acidic media [30,33]. The alteration in proton regulators suggested that aciduric microbiotas are likely to exhibit active dehydrogenation pathways for pH homeostasis and acid production in acidenriched soils.
Ultimately, we provide evidence that pH homeostasis gene modules had the greatest contribution to Rhodanobacter in V20 than V0 samples (Fig. 6b, 7). Regression analysis in pH homeostasis modules and acid-sensing biomarkers community composition revealed that that pH homeostasis modules might have the role in establishment of the acid-sensing microbiomes (Fig. 6c). Previous studies showed that the typical function of Rhodanobacter is their capability to undergo complete denitrification [53], which probably ascribed to the upregulation of nitrate reduction module in acidenriched soils (Fig. 7d). Denitrification and associated denitrifying microorganisms are known to be strongly inhibited below pH 5, but Rhodanobacter thrives in conditions of low pH (pH 3-4), high nitrate (10-100 mmol) [47]. Indeed, it has been confirmed that dehydrogenation reactions triggered by pH homeostasis genes are the pivotal feature of acidophilic pathogens, such as Helicobacter pylori, Escherichia coli, Streptococcus mutans, and has been identified as an important contributor to acid production [31,32]. Consequently, we described some supplementary clues of microbial acid tolerance, and hinted at a potential role of pH homeostasis gene modules in triggering microbial dehydrogenation and acid production reactions, while are likely responsible for recruiting a large proportion of aciduric biomarkers, especially Rhodanobacter.
Additionally, we found soils in V20 samples noticeably enriched genes associated with vancomycin resistance (M00651) and multidrug resistance (M00648) (Fig. 5b), illustrating the antibiotic resistance was abundant in V20 microbiome. It confirmed the concern that some of these resistance mechanisms may be exported from the "underground" world to the genomes of human pathogens [55]. Actually, alterations in pH homeostasis have been implicated in antibiotic resistance that low pH with upregulated proton pump could accelerate drug efflux [56,57]. In our study, proton pumps modules as well as acid cations saturation in V20 shown highly coincided with those antibiotic resistance modules. We claim that better investigation of possible soil antibiotic resistance, particularly for inhibition of proton pumps, may be a clue used by aciduric microbiota to avoid drug efflux and acid accumulation in intensive agroecosystems. The main sampling sites were from the same intensive farm (Qingjiang Community Farm), which were subdivided into three different intensive cropping history: (i) traditional rice-wheat rotations (nonconversion control treatment, V0); (ii) short-term (V10) versus (iii) long-term (V20) rice-vegetable rotations that were separately converted from rice-wheat rotations 10 and 20 years ago. In choosing sample plots, we aimed for consistency in soil texture (medium loam) and fertility (moderate). Each land-use is replicated (three plots per type) in a random block design, a total of 9 plots with 50 × 50 m (Table 1). Sampling was carried out on 16 April 2018, after growing season. At each plot, soil (0-20 cm) were sampled in 5 points, and the distance between each point was above 10 m [58]. Once sampled, all the samples from each site were manually pooled into one sample mixing [59], and were placed in aseptic-refrigeration boxes, then transported to the laboratory immediately. Each sample was divided into two parts, one was stored at 4 ℃ for physical and chemical analyses, whereas the other was stored at -80 °C until DNA extraction.

Metagenomic Annotation And Abundance Profiling
The taxonomic assignment (predicted genes clustered into metagenomic species) was carried out using BLASTP alignment against the integrated non-redundant (NR) database of the NCBI [66]. The abundance of a taxonomic specie was calculated by summing the abundance of genes annotated to a feature and the abundance of species in each sample was calculated at the taxonomic levels of domain, kingdom, phylum, class and genus, so as to construct the abundance table at the corresponding taxonomic level. High-order taxonomic composition for domain, kingdom, phylum and class level were visualized as Krona charts [67]. Fine-scale taxonomic composition for genus level were visualized as heat map with "vegan" package in R [68].
To assess the functional assignment, The NR gene catalog was aligned against the Kyoto Raw fastq files were quality-filtered by Trimmomatic and merged by FLASH with the following criteria: (i) The reads were truncated at any site receiving an average quality score < 20 over a 50 bp sliding window. (ii) Sequences whose overlap being longer than 10 bp were merged according to their overlap with mismatch no more than 2 bp. (iii)Sequences of each sample were separated according to barcodes (exactly matching) and Primers (allowing 2 nucleotide mismatching), and reads containing ambiguous bases were removed. Operational taxonomic units (OTUs) were clustered with 97% similarity cutoff using UPARSE(version 7.1 http://drive5.com/uparse/) with a novel'greedy'algorithm that performs chimera filtering and OTU clustering simultaneously. The taxonomy of each 16S rRNA gene sequence was analyzed by RDP Classifier algorithm (http://rdp.cme.msu.edu/) against the Silva (SSU128) 16S rRNA database using confidence threshold of 70% [72].

Statistical analysis
Differentially soil physicochemical properties, abundant species and KEGG modules were analyzed by ANOVA and ranked according to Tukey-kramer test (p < 0.05 or 0.01 or 0.001). The samples distances dendrogram and heat map were created with dissimilarity hierarchical clustering of metagenomic and bacterial16S rRNA genes based on beta diversity Bray-Curtis distance matrices that obtained from QIIME [73]; Principal Co-ordinates Analysis (PCoA) and Analysis of Similarity (ANOSIM) with Bray-Curtis distance tested whether sample groups are significantly different; a redundancy analysis (RDA) to test the relationships between acidity parameters, microbial community and samples, and permutest was employed to analyze the statistical significance of acidity parameters, and above-mentioned analysis were visualized in R "vegan" package [68].

Species And Functional Classifiers
To build potential classifiers like acid-sensing biomarkers and indicative KEGG modules for acidified identification, we performed linear discriminant analysis (LDA). Here, we applied on-parametric factorial Kruskal-Wallis sum-rank test with more strictly multi-group comparison strategy (all-against-all) to identify taxonomic biomarkers and indicative KEGG modules that characterize statistical differences between soils with and without acidification; then wilcoxon rank sum test was used to investigate biological consistency of differences in subgroups, and LEfSe (the logarithmic LDA score threshold set at 2.0) was selected to estimate the effect of each component abundance on the differences (http://huttenhower.sph.harvard.edu/galaxy/root?tool_id=lefse_upload) [74] Correlation And Contribution Analysis Sparse correlations for compositional data (SparCC) network [75] was performed to construct the association network of above-mentioned biomarkers. Spearman's rank correlation was used to find correlations between acid-sensing biomarkers and identify potential aciduric biomarkers. Only significant correlations (p < 0.05) were linked in the network, and visualized in Python with "networkx" package.
Using correlation heat map analysis to identify indicative modules most likely correlate to acid cations: spearman's rank correlation was used to find correlations of acidity parameters and indicative modules, and visualized in R "pheatmap" package [76].
Statistical composition based on relative contribution (%) of indicative KEGG modules assigned to aciduric biomarkers were done using the method described by Ofek-Lalzar et al [77].

Additional Files
Additional files 1: Statistics of metagenomic data sequence. Table S1. Characteristics of raw data sequence.

Competing interests
The authors declare that they have no competing interests.

Funding
This work was supported by National Key Research and Development Program of China (2017YFD0301701), Table 1 Due to technical limitations, Table 1