Soil metabolomics and bacterial functional traits revealed the responses of rhizosphere soil bacterial community to long-term continuous cropping of Tibetan barley

Site description and rhizosphere soil sampling

The study was carried out at an alpine Tibetan barley experiment site located in the Qinghai–Tibetan Plateau (37°21′N, 101°43′E; altitude, 3,700 m), Beishan Township, Menyuan County, Qinghai Province, China. This region has a plateau continental climate with a mean annual temperature of 1.3 °C and mean annual precipitation of 530–560 mm. According to the Food and Agriculture Organization (FAO) classification, the soil type is classified as Kastanozems (Schad & Spaargaren, 2006). Kunlun 14 barley (Hordeum vulgare L.) was used as the experiment material. The study site consisted of independent plots over many years of the successive monoculture of Tibetan barley. Tibetan barley seeds were cultivated in April each year after the annual application of blended fertilizer and harvested in August following previously described methods (Yao et al., 2020). Briefly, in July 2020, we selected three 5 m × 5 m quadrats including 2, 5, and 10 years of successive monocropping of Tibetan barley, which we named CCY02, CCY05, and CCY10, respectively. Each quadrat was originally divided into five rectangular 5 m × 1 m sub-quadrats. In each sub-quadrat, we collected a total of 15 rhizosphere soil samples (closely adhered to the roots) in a “Z” pattern and all rhizosphere soil samples within a quadrat were mixed together as a single biological sample. The same steps were repeated four times to make four biological rhizosphere soil duplicates per continuous monocropping year. Finally, we obtained a total 12 rhizosphere soil samples for the subsequent experiments. Soil samples were then immediately transferred on ice to the laboratory and stored at −20 °C for subsequent analysis.

High-throughput sequencing of bacterial 16S rRNA gene

The genomic DNA was extracted from each rhizosphere soil sample using an E.Z.N.A^®. Soil DNA kit (Omega Bio-Tek, Norcross, GA, USA) following the manufacturer’s instructions. The bacterial V3-4 region was then amplified with the primer set 338F (5′-ACT CCT ACG GGA GGC AGC AG-3′) and 806R (5′-GGA CTA CHV GGG TWT CTA AT-3′) (Zeng & An, 2021). PCR amplification and purification were performed following previously described methods (Yang, Liu & Zhang, 2019). Briefly, each sample consisted of a 30 µL mixture in triplicate, comprised of 10 µL TaKaTa EX Tag PCR premix (TaKaRa Bio Inc., Kusatsu, Japan), 0.5 µL of each primer (10 µM), 1 µL DNA template (20 ng), and 18 µL PCR-grade water, which was amplified according to the stetted procedure (95 °C for 5 min, followed by 25 cycles at 95 °C for 30 s, 55 °C for 30 s, 72 °C for 30 s, and a final extension at 72 °C for 10 min). The replicate PCR reactions for each triplicate were pooled and purified using a QIAquick Gel Extraction Kit (Qiagen, Chatsworth, CA, USA). A single composite sample for sequencing was prepared by combining approximately equimolar amounts of PCR products from each sample. The prepared PCR products were then paired-end sequenced (2 × 250 bp) on an Illumina MiSeq platform (Majorbio Biotech CO., Ltd., Shanghai, China) following the manufacturer’s instructions.

Soil non-targeted metabolomic detection and analysis

Twelve rhizosphere soil samples were sent to Majorbio (Shanghai, China) on dry ice for metabolite extraction, detection, and analysis. In detail, 1,000-mg soil aliquot of each sample was homogenized with 1,000 µL methanol/water (4:1, v/v) solution, including 0.02 mg/mL internal standard (L-2-chlorophenylamine acid) for 6 min at −10 °C and 50 kHz using a Geno-grinder 2,000 (SPEX, Metuchen, NJ, USA) and spun down for 30 min at 5 °C and 40 kHz. After resting for 30 min at −20 °C, each material was centrifuged at 13,000 g (relative centrifugal force) for 15 min at 4 °C, then the supernatant was transferred and concentrated by a Termovap Sample Concentrator (DC-24, Anpel Laboratory Technologies, Shanghai, China). The dry residue was derivatized by adding 50 µL acetonitrile/water (1:1), homogenized for 30 s at 5 °C and 40 kHz, and centrifuged at 13,000 g (relative centrifugal force) for 10 min at 4 °C. Finally, the supernatant was subjected to LC-tandem MS (MS/MS) analysis; 20 µL supernatant of each sample was mixed for the quality control sample.

The UHPLC-Triple TOF system (AB Sciex, Foster City, CA, USA) equipped with an ACQUITY UPLC HSS T3 (100 mm × 2.1 mm i.d., 1.8 µm; Waters, Milford, CT, USA) was applied to chromatographic separating the metabolites with two mobile phases, phase A (95% water and 5% acetonitrile (with 0.1% formic acid)) and phase B (5% water (with 0.1% formic acid), 47.5% acetonitrile, and 47.5% isopropanol). The whole system was integrated to a quadrupole time-of-flight mass spectrometer (Triple TOFTM 5600th; AB Sciex, Foster City, CA, USA) equipped with an electrospray ionization source and operated in positive and negative modes.

The online analysis platform (Majorbio, Shanghai, China) was applied to the raw data following the manufacturer’s instructions. Briefly, the ProgenesisQI (Waters Corp., Milford, CT, USA) was used for baseline filtering, peak recognition, integration, retention time correction, and peak alignment, which produced the retention time, mass charge ratio, peak intensity, and data matrix. Then, the characteristic peak was detected and the information from the MS and MS/MS analyses was mapped with a metabolic specific database; the mass error was set to <10 PPM. Finally, the metabolites were identified according to the matching scores of the secondary MS. The data were normalized with Pareto scaling and log-transformed before further analysis.

QMEC

Twelve qualified soil DNA samples were sent to Guangdong Magigene Biotechnology Co., Ltd. (Guangzhou, China) for QMEC analysis, wherein a high-throughput quantitative PCR (qPCR)-based chip was applied to assess the microbial functional potential, including 71 functional genes related to C-, N-, P-, and S-cycling (Zheng et al., 2018). QMEC manipulation was conducted following previously described methods (Chen et al., 2020a). Briefly, DNA templates and qPCR reagents were added to the sample source-plate; each primer and qPCR reagent was added to a separate source-plate. In each run, using a non-template reaction as the negative control, 100-nL reactions were mixed in parallel three times using an automated high-throughput sample preparation device, and then added to the nanopore of the qPCR chip on a Wafergen Smart-Chip Real-time PCR platform (Wafergen, Fremont, CA, USA). Specifically, a threshold cycle (CT) of 31 was used as the detection limit and multiple melting peaks; amplification efficiencies outside the range of 0.9–1.1 were discarded.

Data and statistical analyses

The 16S rRNA gene raw reads were processed by the QIIME2 v2020.08 pipeline (Bolyen et al., 2019). Briefly, paired-end reads were joined by FLASH v1.2.11 (Magoč & Salzberg, 2011), then sequences were demultiplexed by the q2-demux plugin. Afterwards, plugin q2-dada2 was used to conduct quality control, chimeric sequence removal, and sequence clustering. Taxonomic analyses were conducted using the plugins q2-feature-classifier and SILVA v132 database (Quast et al., 2012). The α diversity of the bacterial communities was represented by the Chao1, Shannon, ACE, and Faith’s phylogenetic diversity (PD) indices, which were calculated from the amplicon sequence variants (ASVs) table rarefied to 20,349 sequences in QIIME2 for each sample.

Beta diversity metrics differ widely in the types of differences they detect, therefore, four distance-based principal coordinate analysis (PCoA), including the Bray–Curtis, Jaccard, Weighted Unifrac, and Unweighted Unifrac were applied to visualize the differentiation among samples. Non-parametric multivariate statistical tests, Kruskal-Wallis or Adonis, were implemented to test for significant differences between the variances of bacterial communities with a P < 0.05 significance threshold. The similarity of percentages analysis (SIMPER) was applied to identify the most affected genus between different continuous cropping years using the R package, Vegan (Oksanen et al., 2013). Procrustes included in the Vegan R package was used to compare the congruence of the microbiome and metabolism or functional profiles based on the P-values and goodness of fit (m₂). The co-occurrence network construction and parameter statistics were performed using Cytoscape v3.9.0 software (Shannon et al., 2003). The network was visualized using the interactive platform, Gephi (Cherven, 2013). In each network, the node size (i.e., degrees) was proportional to its number of connections; the thickness of each connection between two nodes (i.e., edges) was proportional to the Spearman’s correlation coefficient (|r| = 0.6–1).

Sequence accessions

The bacterial 16S rRNA gene sequencing data are publicly available in the NCBI Short Read Archive (SRA) under Bioproject accession No.: PRJNA759342.

Overall measures of diversity

We obtained 503,330 raw sequences and an average of 41,944 sequences per sample. After quality control and sample normalization, a total of 244,188 sequences were filtered and clustered into 6,745 unique ASVs (Edgar, 2018). The coverage of samples was >93% and the rarefaction curve of each sample approached a saturation plateau (Fig. S1), indicating that the current sequencing depth reflected the microbial composition. The Chao1, ACE, and PD indices were all significantly higher in CCY02 and CCY10 compared with CCY05 (all P < 0.05) (Table 1). The Shannon index was significantly higher in CCY10 when compared to CCY05, but no significant differences existed between CCY10 and CCY02, or CCY02 and CCY05, indicating that the bacterial diversity measured by the species richness with evenness had been reduced, but was not significant, during short continuous cropping years; after longer monocropping years, the whole bacterial diversity seemed to recover.

According to four different distance-based PCoA, we found that samples from the same cropping year clustered tightly, while samples from different years clearly separated based on both Jaccard (Fig. S2A) and Unweighted UniFrac PCoA (Fig. S2B), indicating that there clear differences existed in the absence or presence of ASVs among groups. Meanwhile, we found by Bray–Curtis PCoA (Fig. S2C) that, when considering the presence, absence, and abundance of ASVs, CCY05 obviously separated from CCY02 and CCY10. Additionally, there were no obvious differences detected among communities per the Weighted-Unifrac distances (Fig. S2D), which considered the presence, absence, abundance, and up-weighting differences in distant evolutionarily-related ASVs. Furthermore, the soil microbial compositions under different treatments were analyzed by Adonis statistical tests, which revealed significant dissimilarities among different cropping years (Table S1).

Microbial community composition

Overall, Actinobacteria (average ± standard deviation (SD), 28.84 ± 1.50%) and Proteobacteria (26.88 ± 1.50%) were the dominant phyla across all rhizosphere soil samples. The genera with relative abundances >1% included Blastococcus (2.57 ± 0.30%), Nocardioides (2.55 ± 0.39%), Sphingomonas (2.44 ± 0.28%), Bacillus (1.09 ± 0.30%), and Solirubrobacter (1.02 ± 0.04%), but these genera did not significantly differ among continuous cropping years. After comparing the relative abundances of each known taxa among the three successive monoculture years, we found that six phyla, five classes, 11 orders, 23 families, and 47 genera exhibited a straight increasing trend over time, while five phyla, eight classes, 14 orders, 19 families, and 32 genera exhibited the opposite trend. A total of 20 taxa exhibited significant continuously strengthened or weakened tolerance to the long-term continuous cropping of Tibetan barley based on the ANOVA results (Fig. 1), but most of these taxa were low in abundance with a mean relative abundance <1%. For example, Fibrobacteres is an important phylum of cellulose-degrading bacteria (Ransom-Jones et al., 2012); the genus Devosia possesses bioremediation potential (Talwar et al., 2020); and the denitrifying bacteria Noviherbaspirillum (Wu et al., 2021) were significantly upregulated over time. The following significantly decreased over time: Firmicutes, which are plant growth-promoting bacilli (Kumar, Khare & Dubey, 2012); Bacillales, which have effective biological control and biodegradation potential (Barathi et al., 2020); Nitrosomonadales, which are related to sulfate and iron reduction; and Desulfuromonadales, which are an order capable of iron and sulfate reduction (Wunder et al., 2021).

A SIMPER analysis was used to detect the main genera that drove the compositional shifts of bacterial communities over time. The results showed that 22, 24, and 25 genera were responsible for >50% cumulative dissimilarity of the microbial community shifts between CCY02 and CCY05, CCY02 and CCY10, and CCY05 and CCY10, respectively (Fig. 2). Fourteen of these overlapping genera were mainly responsible for the differences between each pair of cropping years (Table S2), and the genera Nocardioides had the largest dissimilarity contribution.

Core taxa among different cropping years

The consistency of crop type, planting strategies, and research site may lead to a cohort of rhizosphere soil-shared microbial communities from Tibetan barley, while continuous cropping over time could lead to year-unique taxa in the rhizosphere soil. In this study, the taxa that were simultaneously present in all samples were defined as core or shared taxa, and when their mean relative abundances (at the genus level) were >or <1%, they were defined as dominant or rare core taxa and subsequently classified at the phylum, class, order, family, and genus levels (Gobet, Quince & Ramette, 2010). At the genus level, the core shared microbiome consisted of nearly 50.87% of all taxa (taxonomic richness) and this percentage reached 60% at the phylum level, however, a proportion of unshared taxa seemed to be sample-specific (Fig. 3A). Surprisingly, the total relative abundances of the core taxa at the genus level occupied >98% of all sequences on average (Fig. 3A), while unshared taxa, or sample-unique taxa, only occupied 1.66 ± 0.31% of all sequences. A total of 164 core taxa were successfully annotated at the genus level and the relative abundances ranged from 0.02 ± 0.01% to 2.57 ± 0.61% (Table S3). The dominant shared genera among samples were Blastococcus, Nocardioides, Sphingomonas, Bacillus, and Solirubrobacter (Fig. 3B).

Bacterial functional genes

QMEC was implemented to determine the abundances of a wide spectrum of functional genes related to C-, N-, P-, and S-cycling (71 genes in total) in a high-throughput manner. A total of 60 genes were successfully detected by QMEC across 12 samples, including 17 genes related to C-degradation, 13 related to C-fixation, 19 related to N-cycling, 6 related to P-cycling, and 5 related to S-cycling (Table S4). Samples from CCY02 and CCY10 displayed a clear clustered pattern based on the principal components analysis (PCA), while samples from CCY05 were scattered (Fig. S3). The abundances of each gene in 12 samples are presented in Fig. S4. Briefly, we found that the genes involved in hemicellulose degradation (abfA and xylA) were most prevalent in C-degradation. The gene involved in the Wood–Ljungdahl pathway (acsA) and the gene involved in ammonization (UreC) were most abundant in C-fixation and N-cycling, respectively. In P-cycling, the genes involved in organic P mineralization (phoD and phnK) were the dominant genes across samples, while the gene related to S reduction (apsA) and the genes (YedZ and SoxY) related to S oxidation were the main genes in S-cycling.

We employed Z-values to compare differences in cycling genes among different cropping years (Table S5). The long-term continuous cropping of Tibetan barley significantly increased the abundances of genes related to C-degradation (F = 9.25, P = 0.01) and P-cycling (F = 5.35, P = 0.03) (Figs. S5A and S5B). The abundances of N-cycling genes increased from CCY02 to CCY05 and decreased by CCY10, but no significant differences were detected between groups (Fig. S5C). We also investigated the relationship between the abundances of functional genes and bacterial diversity (Shannon) (Table S6). Results revealed that only N-cycling significantly correlated with bacterial diversity (r = −0.71, P = 0.01).

Relationships between taxa and functional genes

The Procrustes analysis showed that there were significant correlations between bacterial composition and functional gene profiles after long-term continuous cropping (P = 0.012) (Fig. 4A), but not after short-term continuous cropping (P = 0.111) (Fig. 4B). A co-occurrence network was constructed to visualize the relationships between the functional cycling genes and bacterial genera (Fig. 5). Only the interactions between the genera and functional genes were kept, of which 173 nodes and 481 edges (196 positive and 285 negative associations) were obtained. The network had good modularity (0.41) and contained 6 modules; each included tightly connected nodes. Module I, which was the largest, contained 65 nodes, including 43 functional genes and 22 genera; most of the C-, N-, P-, and S-cycling genes were tightly connected with a small group of genera. Nine key genera with a degree value >10 were identified, including Saccharothrix, Planomicrobium, Edaphobaculum, Nitrosospira, Angustibacter, Sanguibacter, Brevibacillus, Caulobacter, and Pantoea. Interestingly, Saccharothrix, Edaphobaculum, Sanguibacter, and Pantoea were positively connected with the functional genes, while Planomicrobium, Nitrosospira, Angustibacter, Brevibacillus, and Caulobacter were mainly negatively connected (Table S7).

We also examined the correlations between the cycling genes and 14 top contributors with the dissimilarities among different cropping years. Results revealed that Ohtaekwangia (r = 0.64, P = 0.02), Pedobacter (r = 0.63, P = 0.03), and Gaiella (r = 0.63, P = 0.03) significantly positively correlated with gmGDH, which is related to C degradation (Fig. 6). With regard to C fixation, Devosia significantly positively correlated with cdaR (r = 0.66, P = 0.02), while Gaiella significantly positively correlated with frdA (r = 0.58, P = 0.04) and rbcL (r = 0.68, P = 0.01). For N-cycling, Devosia and Sphingomonas significantly negatively correlated with amoA2 and nxrA (P < 0.05), while Saccharibacillus significantly positively correlated with these genes (P < 0.05). For P-cycling, Blastococcus significantly negatively correlated with dsrB (r = −0.58, P = 0.04), while Marmoricola significantly positively correlated with phnK (r = 0.60, P = 0.04). No significant correlations were detected between the five major categories of functional genes and 14 top genera (Fig. S6).

Soil metabolic profiles

By LC-MS/MS-based non-targeted metabolomics, a total of 14,380 metabolites in the soils were obtained, in which, 822 metabolites were successfully annotated in the HMDB and KEGG databases. Based on the metabolic profiles, we implemented a PCA to uncover the different compositions of metabolites in the soil samples from different continuous cropping years. Results revealed that the first two principal components (PC1 and PC2) explained 62.93% of the total variance and samples from different groups obviously differed from each other (Fig. S7). The varying metabolic compositions among different groups were further confirmed by the heatmap analysis of the top 100 metabolites (Fig. S8).

VIP values in the PLS-DA model were calculated to examine the changes in the soil metabolites in greater detail. A total of 440 metabolites with VIP >1.0 and P < 0.05 (ANOVA) were considered significantly affected by the continuous cropping of Tibetan barley. Further linear regression analyses between the acquired metabolites and enhanced continuous cropping years were conducted. A total of 126 metabolites significantly responded to continuous cropping with an increase in 83 metabolites and decrease in 43 metabolites (Table S8). The changes in metabolites mainly occurred in lipids and lipid-like molecules, organic acids and derivatives, and organoheterocyclic compounds. Lipids and lipid-like molecules had the largest number of metabolites affected by continuous cropping. Moreover, nucleosides, nucleotides, and analogues were mainly upregulated over time, while alkaloids and derivatives were mainly downregulated.

Relationships between taxa and metabolites

First, we applied Procrustes tests to depict the correlations between soil metabolic profiles and the bacterial community structure. Significant correlations were detected between specific metabolites and the bacterial community structure in rhizosphere soils during the long- (P = 0.001) and short-term continuous cropping of Tibetan barley (P = 0.027) (Fig. 7). To elucidate which microbial taxa were responsible for the changes in soil microbial metabolism, an interactive network linking the microbes with significant differences between groups (a total of 100 microbial taxa) and differential metabolites (a of total 126 metabolites) was constructed. Only the interactions between metabolites and microbial taxa were kept. The orders Desulfuromonadales (degree 57) and Nitrosomonadales (degree 19), families Archangiaceae (degree 34), Nocardiaceae (degree 23), and Sanguibacteraceae (degree 11), and genera Terrimicrobium (degree 35), Caldinitratiruptor (degree 34), Dehalogenimonas (degree 21), Rubellimicrobium (degree 15), Aquimonas (degree 13), and Sanguibacter (degree 11) co-occurred most frequently with several differential metabolites, including 4-hydroxyvalsartan, methyl-4-pentanoate, pentigetide, and glaucarubolone-15-O-beta-D-glucopyranoside (Fig. 8). Of these metabolites, the frequently co-occurring ones were downregulated over time, except 4-hydroxyvalsartan.

Effects of continuous cropping on microbial community structure and diversity

Long-term continuous cropping may significantly alter bacterial community structure and diversity. For instance, Yan et al. (2021) reported that soil bacterial diversity significantly decreased (PD and Shannon) after the continuous cropping of Nicotiana tabacum. Chen et al. (2020b) found that the long-term consecutive monocropping of peanuts led to a general, albeit insignificant, loss in bacterial diversity. Liu et al. (2020) found that the bacterial diversity (PD and number of operational taxonomic units (OTUs)) was significantly higher after the long-term (13 year) continuous monoculture of soybeans, while no significant differences were detected after the short-term monoculture periods (3 and 5 years). Zhao et al. (2020) showed that bacterial richness (ACE) significantly decreased after 5–10 years of continuous cropping and recovered after 30 years, while bacterial diversity (Shannon) significantly increased after continuous cropping (5, 10, and 30 years). In the present study, bacterial diversity (Shannon) did not differ during short continuous cropping years, but bacterial richness (Chao1 and ACE) and PD significantly declined after short-term continuous cropping of Tibetan barley, which is consistent with the findings of Yao et al. (2020), but surprisingly recovered after long-term continuous cropping (Table 1). This phenomenon may be due to the fact that all environmental filters require specific adaptation strategies for survival; stronger selection acting on bacteria often leads to a new cohort of microbiota that adapt to the environment (Fierer, Bradford & Jackson, 2007; Meola, Lazzaro & Zeyer, 2014). Even if the microbial composition is sensitive to environmental disturbances, the community may be resilient and stabilize over time (Allison & Martiny, 2008). However, we believe that the application of more powerful methods (e.g., metagenome and metagenome binning) to deeply investigate the reasons for this phenomenon in future studies.

After 10 years of continuous cropping of Tibetan barley, Actinobacteria and Proteobacteria were the two main dominant phyla that were enriched in many rhizosphere soils after continuous cropping, such as in sugarcane and cotton (Pang et al., 2021; Xi et al., 2019). The widely distributed bacterial phyla, Actinobacteria, possess high proportions of CAZymes and exhibited a key important eco-physiological role in plant residue decomposition in a previous study (Bao et al., 2021). Another abundant phylum, Proteobacteria, plays a central role in the cycling of several key elements, including N- (Conthe et al., 2018), C- (Chan et al., 2013), and S-cycling (Zhou et al., 2020). In this study, a total of 20 low abundant taxa (<1%) significantly increased after continuous cropping, during which planting-promoting bacteria or bacteria with bioremediation abilities significantly decreased.

In the current study, we analyzed the top genera that contributed to the dissimilarities among microbial communities from different continuous cropping years. Of these, five were identified as the dominant shared genera across all soil samples, including Blastococcus, Nocardioides, Sphingomonas, Bacillus, and Solirubrobacter (Table S2). These genera have been frequently detected in the rhizosphere soils after the continuous cropping of different crops, including maize (Zhao et al., 2021), cucumber (Li et al., 2021), cotton (Han et al., 2017), and peanuts (Chen et al., 2018). Solirubrobacter and Blastococcus exhibit heavy metal tolerance and advantages in stress resistance, indicating the promising potential for alleviating polluted soil ecosystems (Hou et al., 2021; Wang et al., 2021). A recent study suggested that Sphingomonas possesses multifaceted functions ranging from the remediation of environmental contamination to producing highly beneficial phytohormones that promote plant growth (Asaf et al., 2020). Bacillus has multiple medical, environmental, and industrial applications (Khurana et al., 2020). Hence, it is proposed that they may play important roles in the successive monoculture of Tibetan barley.

Functional genes and their interactions with microbes

For the first time, QMEC was applied to determine the relative abundances of key element cycling genes after long-term continuous cropping. We found that overall bacterial composition significantly correlated with functional gene profiles (Fig. 4A), indicating that after long-term continuous cropping, and alterations in community functions were mainly affected by alterations in bacterial composition. However, after short-term continuous cropping (CCY02 and CCY05), the functional gene profiles and bacterial compositions had poor consistency (Fig. 4B). Thus, we proposed that the functional gene composition was not mainly mediated by changes in microbial composition after short-term continuous cropping. This may be a result of functional redundancy (Rosenfeld, 2002), which contributed to bacterial community tolerance and overall community functioning after the short-term continuous monoculture of Tibetan barley. However, significant fluctuations in the taxonomic composition after the long-term continuous monoculture may potentially reduce the stability of the community, which would thereby become more vulnerable to continuous disturbances (Sheng et al., 2015), ultimately affecting bacterial functional compositions.

Genes related to C-degradation and P-cycling increased over time (Table S5). Hemicellulose degradation (abfA and xylA) was most prevalent in C degradation, indicating the potential loss of soil C storage (Chen et al., 2020a). Organic P mineralization (phoD and phnK) increased after continuous cropping, playing an important role in soil P bioavailability. Meyer et al. (2018) indicated that microbes compensated for single nutrient deficiencies by accelerating P- or N-cycling and may have increased SOC turnover in co-limited subsoils with acquirable P reserves. Collectively, these results indicated that, along with the continuous cropping of Tibetan barley, the activity of hemicellulose degradation increased and accelerated C degradation, while microbes may have acquired greater C or N supplies by promoting organic P mineralization.

Furthermore, a set of genera were identified that showed tight correlations with most of the functional genes, indicating their functional potential in the cycling of key elements underlying the continuous monocropping of Tibetan barley (Fig. 5). In a previous study, Nitrosospira was the primary ammonia-oxidizing bacteria and was mainly responsible for nitrous oxide production (Lourenço et al., 2018). Pantoea fixes N or induces N uptake, thereby promoting N availability in plants (Loiret et al., 2004; Singh et al., 2020). Additionally, Brevibacillus was reported to have desulfurization activity (Nassar et al., 2013).

A relevant study pointed out that the continuous cropping of Tibetan barley significantly increased the bacteria associated with chemo-heterotrophy, aromatic compound degradation, and nitrate reduction (Yao et al., 2020). This study showed that the dominant shared genera, Sphingomonas, exerted significant negative effects on nitrification gene abundances (amoA and nxrA), leading to the imbalance of nitrification and denitrification. Additionally, a previous study showed that Sphingomonas in vegetable systems utilized contaminants originating from pesticide residues as a growth and energy source (Busse et al., 2003).

Soil metabolites and their interactions with microbes

Chemical dialog undergoing rhizosphere active zones promotes the interactions between plant roots and microbial communities (Van Dam & Bouwmeester, 2016). Shen et al. (2020) observed distinct differences in the metabolite composition over different continuous cropping seasons in tobacco. Our study explored the alterations in soil metabolic profiles, which were significantly mediated by soil microbial communities, indicating that certain metabolic pathways were altered for soil microbial communities to adapt to environmental stress (Fig. 7). Obvious differences in the patterns of metabolite compositions over time were observed (Figs. S7 and S8), in which the primary metabolites related to lipid and nucleoside metabolism were significantly upregulated, while secondary metabolites, such as alkaloids and derivatives, were significantly downregulated. A previous study pointed out that secondary metabolites greatly affect the microbial community, where some metabolites were antibiotic and pharmaceutically relevant, while others were involved in disease interactions (Fox & Howlett, 2008).

We found that the hub taxa Desulfuromonadales, which was reported to utilize different electron exchange pathways (Rotaru et al., 2012), mainly co-occurred with metabolites in either negative or positive correlations (Fig. 8), indicating potential important interactions between them. Meanwhile, a previous study suggested that positive co-occurrences indicated that the metabolites may have be mainly originated from microbes, while negative co-occurrences may represent the specific microbial consumption or degradation (Devi et al., 2017). These findings indicated that soil metabolomics can be used to assess the adaptations of soil microbial communities to continuous cropping strategies at the molecular level.