Skip to main content

Advertisement

We’d like to understand how you use our websites in order to improve them. Register your interest.

The microbiome as a biosensor: functional profiles elucidate hidden stress in hosts

Abstract

Background

Microbial communities are highly responsive to environmental cues, and both their structure and activity can be altered in response to changing conditions. We hypothesized that host-associated microbial communities, particularly those colonizing host surfaces, can serve as in situ sensors to reveal environmental conditions experienced by both microorganisms and the host. For a proof-of-concept, we studied a model plant-soil system and employed a non-deterministic gene-centric approach. A holistic analysis was performed using plants of two species and irrigation with water of low quality to induce host stress. Our analyses examined the genetic potential (DNA) and gene expression patterns (RNA) of plant-associated microbial communities, as well as transcriptional profiling of host plants.

Results

Transcriptional analysis of plants irrigated with treated wastewater revealed significant enrichment of general stress-associated root transcripts relative to plants irrigated with fresh water. Metagenomic analysis of root-associated microbial communities in treated wastewater-irrigated plants, however, revealed enrichment of more specific stress-associated genes relating to high levels of salt, high pH and lower levels of oxygen. Meta-analysis of these differentially abundant genes obtained from other metagenome studies, provided evidence of the link between environmental factors such as pH and oxygen and these genes. Analysis of microbial transcriptional response demonstrated that enriched gene content was actively expressed, which implies contemporary response to elevated levels of pH and salt.

Conclusions

We demonstrate here that microbial profiling can elucidate stress signals that cannot be observed even through interrogation of host transcriptome, leading to an alternate mechanism for evaluating in situ conditions experienced by host organisms.

This study is a proof-of-concept for the use of microbial communities as microsensors, with great potential for interrogation of a wide range of host systems.

Video Abstract

Background

Advances in sequencing have propelled the field of microbiology and shifted focus from analysis of microbial isolates or low diversity ecosystems to analysis of environments with highly diverse microbial communities. Global surveys of microbial community structure have been conducted in a wide range of natural environments [1,2,3], also reviewed in [4], and many of these studies have focused on host-associated microbiomes. Such host-associated microbial environments include plant-associated communities [5,6,7,8,9] though the greatest effort has been placed on the human-associated microbial communities [10,11,12,13]. Studies examining plant host-associated microbial communities frequently focused on soil microorganisms that are enriched in the rhizosphere—the soil surrounding and affected by the roots [8]. In the rhizosphere, soil type [6, 8, 9] and plant host type [14, 15] have been identified as the main forces determining rhizosphere and root microbiomes. The selection of rhizosphere-competent organisms from soil has been well established, with specific plants and different growth stages of plants each selecting for different microbial communities from among the high diversity of microorganisms in soil [16, 17]. Rhizosphere microorganisms are further enriched to form sub-populations colonizing root surface [9], as plants shape the soil-plant continuum in a gradient-depended manner [6, 8, 9] mainly through carbon flux to the root environment [18]. Functional profiling of microbial communities associated with different plants has demonstrated that these microbiomes differ in their metabolic activities and has suggested the presence of niche conditions associated with a wide range of factors, of which oxygen concentration is one [19].

More broadly, factors influencing plant root microbiome include geographic location [5, 9], plant developmental stage [15, 20, 21], nutrient (e.g., N or P) availability [7, 22] and redox status [23]. Numerous agricultural practices, which modify many of the above, have been shown to have an impact on root microbiome. These include fertilization [24], compost amendment [25] and irrigation with water of lower quality [26]. Each of these practices alters a wide range of environmental variables, thus confounding the ability to identify the most consequential abiotic factor influencing the plant system and modifying the root microbiome.

Microorganisms sense minor changes in environmental conditions and respond rapidly through transcriptional changes, as well as through microbial amplification—the dynamic modification of the abundance of microbial taxa; these changes occur on a time-scale that is much shorter than for the host [27]. Thus, interrogation of the microbial community may be used as a means to understand environmental conditions on short to long time scales, as well as small to large physical scale. In this study, the root surface is used as a model to test the hypothesis that microorganisms can be sensitive in situ detectors of environmental conditions. The root zone has a number of favorable features for such interrogation, including (a) the presence of a high percentage of microorganisms that are transcriptionally active; (b) high microbial competition for access to root exudates, and therefore likely rapid turn-over if environmental conditions change; and (c) access to high microbial diversity in the soil. Thus, both the composition of the root-associated microbial community and the transcriptional activity of the microbial community can be informative regarding root environmental conditions.

In this study, we followed the root surface microbiome functional response as a micro-sensor to identify stresses imposed by irrigation with water of lower quality, such as treated wastewater (TWW). Our study employs the basic assumptions that the microbiome inhabiting the root surface is exposed to the same environmental conditions as its host, and that the response of the microbiome (i.e., alteration of community structure, associated gene abundance and transcriptional profiles) to stress can identify the specific stress or stresses in the root environment. To examine these assumptions and our general hypothesis, we performed deep DNA and RNA sequencing of plant roots grown in soil irrigated with fresh water (FW) or TWW. Plant transcriptional profiling was examined together with microbial community taxonomic and functional gene content characterization (shotgun metagenome sequencing) and microbial transcriptional profiling (shotgun metatranscriptome sequencing). We observed that the host responded to TWW irrigation in a highly general manner, whereas the microbial response was specific to stresses present in TWW, including elevated salinity and elevated pH.

Results

Here, we assess the use of root-associated microorganisms as an indicator tool to reveal environmental conditions and stresses affecting plant hosts. Our model system was a long-term anthropogenic disturbance caused by soil irrigation with water of lower quality (i.e., treated wastewater, TWW) as compared to irrigation with fresh water (FW). We characterized plant-host response and root microbiome composition and response using deep sequencing of RNA and DNA extracted from roots. Shotgun metagenomic (DNA-based) and metatranscriptomic (RNA-based) analyses were performed on root systems from two host plant types (tomato and lettuce) with two different water treatments (FW or TWW), to describe taxonomic shifts and functional responses associated with long-term root irrigation with water of differing quality. The two plant host tested represent common agricultural crops, with reasonable (tomato) to limited (lettuce) deciphered genetic information, facing an actual environmental stress.

Metagenomic predicted genes were mapped to the SEED [28] and KEGG (Kyoto Encyclopedia of Genes and Genomes) [29] databases for functional predictions, and approximately 22% of non-redundant gene list were annotated by each database. Metatranscriptome reads were mapped to the metagenomic-predicted genes, or to the available plant genome, producing comparable sequencing depth between host and associated microbiome (Table S1).

Host functions actively associated with TWW irrigation

Plant host physiological response to irrigation with water of lower quality has been previously reported, with significantly reduced yield of both plants under TWW irrigation (e.g., [26] Table S2). In this study, we performed deep sequencing of plant transcripts to identify stresses that TWW irrigation imposes on plant roots. A total of 45 tomato genes and 645 predicted lettuce genes were significantly differentially expressed between irrigation treatments. Of the 645 lettuce transcripts, only 141 could be annotated by comparison with known tomato genes (Table S3).

To identify the most robust effects of TWW treatment, tomato and lettuce differentially abundant transcripts were analyzed together. A network analysis of enriched transcripts was performed to predict interactions and highlight clusters of associated genes (Fig. 1). The FW-enriched gene network consisted of 97 nodes, indicating the number of enriched genes that were identified by the STRING protein-protein interaction network database. Similarly, the TWW-enriched gene network consisted of 86 nodes. The FW-enriched gene network, however, was linked by only seven edges representing predicted protein interactions (direct physical interactions, as well as predicted functional association). In contrast, the TWW-enriched gene network was linked by 69 edges, with a significantly higher number of interactions then expected (p value < 0.0001, by Random Graph with Given Degree Sequence (RGGDS)) (Fig 1a, b). The TWW gene network of both plants was enriched (aggregate fold change, permutation-based test) primarily with various heat-shock transcripts, including Hsp20, Hsp70, and DnaJ. Heat shock proteins are prevalent in plants and are active during normal growth [30]. Such genes also show a stress response, and can be activated in response to many stress cues, including heat, cold, water stress, salinity, osmotic stress, and oxidative stress [30, 31]. In addition, tubulin and ‘FKBP-type peptidyl-prolyl cis-trans isomerase’ genes were also significantly enriched under TWW exposure (Fig. 1b, c). Tubulin reorganization has been shown under salt stress, cold shock, aluminum exposure, interaction with pathogens and more [32,33,34,35]. Overall, the plant response to irrigation with TWW, as detected by transcriptome analysis, was largely restricted to highly general stress response genes that are expressed under a wide range of environmental conditions.

Fig. 1
figure1

Plant host functional response to TWW irrigation revealed by tomato and lettuce transcriptome analysis. Differentially expressed plant genes between FW or TWW irrigation were determined using the software package EdgeR, with significance set at FDR p < 0.05. Lettuce transcripts were annotated by comparison to the tomato genome. STRING protein interaction networks are presented for a FW and b TWW irrigation enriched gene list. Only nodes linked by edges are shown. The size of each node is proportional to the log10 (p value) of the enriched gene. Colored are PFAM protein domains significantly enriched within the network landscape. The log10 (p value) of the enriched PFAM domains are presented in c FW and d TWW enrichment analyses

Shifts in microbiome associated with TWW irrigation (DNA-based metagenomics)

Functional profiling demonstrates the extent to which root microbiomes respond to environmental factors

The taxonomic affiliation of root-associated microbial communities was determined by analysis of annotated genes from the metagenomes (Fig. 2f; Table S4, S5). The vast majority of annotated genes were derived from bacteria (96.7% of all mapped reads), while the percentage of reads derived from Fungi (1.2%), Archaea (0.5%), and viruses (0.13%) was much lower. Despite a prior mapping step to remove host reads, 1.2% of annotated gene counts could still be mapped to plant genomes. The root microbiome was primarily composed of bacteria from the phyla Proteobacteria (44% of all mapped reads) and Actinobacteria (33%).

Fig. 2
figure2

Effect of plant host type and irrigation treatment on root-associated metagenome. a VENN diagram of differentially-abundant genes between FW-irrigated tomato and lettuce (DESeq2 FDR p < 0.05). VENN diagrams of differentially-abundant genes by irrigation type in b tomato or c lettuce, and d in both hosts combined. e Canonical correspondence analysis (CCA) of all SEED-annotated genes with DOC, EC, and pH as constrained variables. f Microbial composition predicted by a least common ancestor (LCA) pipeline (MEGAN6) of the predicted gene catalogue. Taxonomic groups are displayed in the inner ring, and differentially abundant taxonomic groups between the two tested plant types are highlighted in the middle ring. The outer ring highlights the taxonomic groups that are significantly differentially abundant between irrigation treatments

The relative abundance of taxa was compared across experimental conditions of plant host type (tomato vs. lettuce) and irrigation water quality (FW vs. TWW) using DESeq2 method for comparing differential abundant count data [36] (Fig. 2f). Broadly, 34% of all taxonomic groups (with highest available taxonomic resolution, based on MEGAN6 least common ancestor) were significantly (Wald test, FDR corrected p value < 0.05) more abundant in tomato roots, as compared to 31%, significantly associated with lettuce roots. Many taxa from the phlya Actinobacteria and Bacteroidetes/Chloroflexi were significantly more abundant in tomato roots relative to lettuce, while Betaproteobacteria and Planctomycetes were strongly and significantly associated with lettuce roots. Irrigation water quality mostly affected Proteobacterial taxa (10% of microbial taxa were significantly more abundant in TWW-irrigated roots, as compared to 11% of microbial taxa enriched in FW-irrigated roots. 60% of all significantly abundant taxonomic group were identified as Proteobacteria). Acidobacteria and Betaproteobacteria were significantly more abundant in FW-irrigated roots and Gammaproteobacteria significantly more abundant in TWW-irrigated roots (all data available at supplementary Table S4).

Root microbial metagenomes from lettuce and tomato were annotated and mapped to the SEED database to identify functional genes significantly associated with plant host type (tomato and lettuce) and irrigation water quality (FW and TWW). A comparison of differentially abundant functional genes between tomato and lettuce root microbiomes demonstrated strong host specificity in microbiome gene content (Fig. 2a), consistent with our prior analyses of root microbiomes of different plant species grown in identical soils [19]. In this study, greater than 50% of SEED annotated genes (from a total of 2625 “functional role” [28]) were significantly (based on Wald test, adjusted p value < 0.05) more abundant in either tomato or lettuce roots (26% in tomato relative to lettuce and 27% in lettuce relative to tomato; Fig. 2a).

Irrigation water quality also affected root-associated microbiome functional profile (Fig. 2b–d). Initially, the effect of irrigation type was examined in tomato and lettuce systems independently by comparing SEED-annotated gene abundance using DESeq2 method. Irrigation type determined 36% of the tomato root metagenome (15% of all annotated gene list in tomatoes were significantly more abundant in FW-irrigated plants compared to 21% more abundant in TWW-irrigated plants, Wald test, Fig. 2b) similarly to 34% of the lettuce root metagenome (FW 16%, TWW 18%, Fig. 2c). To identify commonalities in response to irrigation, irrigation effects were examined in a dataset of both plant hosts combined. This combined analysis identified microbiome genes that were positively associated with FW irrigation (11% of all SEED annotated genes) or with TWW irrigation (14% of all SEED annotated genes) (Fig. 2d). Overall, the combined root microbiome functional profiles were significantly associated with plant host (tomato vs. lettuce, F = 10.1, p value = 0.001) and irrigation treatment (FW vs. TWW, F = 6.6, p value = 0.004) as determined by a PERMANOVA test based on the Bray-Curtis dissimilarity index of the SEED annotated gene counts (Figure S1). For further analyses, unless otherwise indicated, ecosystem comparisons of plants grown in FW and TWW were performed on data combined from both plant hosts.

We previously measured significant increases in pH, dissolved organic matter (DOC) and electrical conductivity (EC) in TWW-irrigated soils relative to FW-irrigated soils [26](also at Table S6). Similar patterns were observed in this study through canonical correspondence analysis (CCA) (Fig. 2e). The CCA presents the relationship of the measured soil parameters and root microbiome functional gene profile (all SEED-annotated gene counts). An ANOVA permutation test was significant (F = 3.2, p value = 0.001 for the full model), and the constrained variables (i.e., pH, DOC, EC) accounted for 54.8% of the variance. DOC (F = 2.4, p value = 0.038, loading primarily on the CCA2, deciphering the irrigation treatments) and pH (F = 5.7, p value = 0.002, loading primarily on the CCA1 axis separating plant hosts) were found to significantly explain portions of the variance associated with the observed microbiome functional profile, while EC was not significant (F = 1.57, p value = 0.163).

SEED and KEGG functional categories enriched or depleted in metagenomes of TWW-irrigated roots relative to FW-irrigated roots

SEED-annotated genes highly significantly (p < 0.01) associated with irrigation water quality across both plants were examined, and in total 438 genes were identified (Fig. 3). Of these genes, 286 were enriched in TWW-irrigated roots and 152 enriched in FW-irrigated roots. These genes were clustered into general categories (SEED level 1, based on the hierarchical clustering available on MEGAN6): e.g., carbohydrates, amino acid derivatives, membrane transports, respiration and regulation of cell signaling (Fig. 3a). Rare categories (n < 5 genes) were removed from the analysis. To compare category enrichment, we examined the proportion of genes enriched for each category (Fig. 3b). For some gene categories (e.g., cell division, cell cycle or carbohydrates), a similar number of genes were enriched in both TWW- or FW-irrigated roots (i.e., no specific effect of irrigation treatment), while others were more strongly skewed to either FW or TWW. For example, membrane transport and transposable element genes were substantially enriched in TWW-irrigated roots. Conversely, the gene categories of nucleosides and nucleotides and sulfur metabolism, were substantially enriched in FW-irrigated roots.

Fig. 3
figure3

Genetic profile of the 438 significantly differentially abundant microbial genes between irrigation treatments. a Heatmap of genes enriched or depleted (FDR p < 0.01) in the metagenomes of TWW-irrigated roots (displaying trimmed mean of M values TMM). Gene abundance was normalized by scaling each row separately. The gene list was clustered to high hierarchy SEED categories. b The proportion of genes enriched or depleted in TWW-irrigated root metagenomes compared to the total abundance of that category in all metagenome analyses. Enriched category (TWW, magenta), deprived (FW, light blue), or the proportion within the full gene catalogue (marked as “all”, colored by orange), are highlighted. The proportion was calculated based on the number of genes assigned to the different categories with-in each data set

An enrichment analysis was also conducted for gene subsystem enrichment and depletion by irrigation method (level 2, based on the SEED hierarchical clustering, tested using Wallenius non-central hypergeometric distribution) (Fig. 4a, Table S7). The most strongly enriched category (log2FC = 1.4; relative abundance = 0.07%, p value < 0.0001) was the Na(+)-translocating NADH-quinone reductase (NQR), a membrane complex that utilizes the respiratory chain to generate a sodium gradient in place of a proton gradient in high pH and sodium conditions [37]. In addition, enrichment of multiple membrane-associated subsystems in TWW-irrigated roots was observed, including (i) sodium-hydrogen antiporter, a common membrane transporter that supports sodium balance in exchange for proton motive force (log2FC = 1.67; p value = 0.002) , (ii) pH adaptation potassium efflux system (log2FC = 1.87; p value = 0.007), (iii) mannose-sensitive hemagglutinin, type 4 pilus (MSHA4) (log2FC = 1.03; p value = 0.0009), (iv) alginate metabolism membrane complex (log2FC = 0.23; p value = 0.01). In addition to membrane-associated subsystems, other subsystems were enriched in TWW-irrigated roots, including arginine degradation (log2FC = 0.56; p value = 0.04). Five genes enriched within this subsystem catalyze the complete arginine to glutamate pathway (Table S7). Soil Na+ concentration, K+ concentration and pH were correlated with observed gene abundance patterns (Figure S2, Table S8). The relative abundance of TRAP transporters (Pearson's RDOM = 0.81, p valueDOM = 0.001; RNa = 0.78, PNa = 0.003; RK+ = 0.81, PK+ = 0.001) and sodium hydrogen antiporters (RDOM = 0.86, PDOM = 0.0003; RNa = 0.78, PNa = 0.003; RK+ = 0.8, PK+ = 0.001) correlated to organic matter and Na+/ K+ concentrations. The relative abundance of potassium antiporter genes was correlated with Na+ and K+ concentrations (RNa = 0.9, PNa < 0.0001; RK+ = 0.86, PK+ = 0.0003), and MSHA4 gene relative abundance was correlated with pH (RpH = 0.88, PpH = 0.0001). NQR gene abundance was significantly correlated with salt concentration (RNa = 0.8, PNa = 0.002; RK+ = 0.81, PK+ = 0.001) and pH (RpH = 0.79, PpH = 0.002).

Fig. 4
figure4

Analyses of SEED subsystems and KEGG pathways significantly enriched or depleted in metagenome of TWW-irrigated roots. a Dot plot of log2(fold-change) relative abundance of enriched or depleted SEED subsystem. Significantly (FDR p < 0.05, represented by more than two gene families) enriched or depleted gene abundance was computed using the goseq software package, with correction for read counts. Symbols are proportional to the sub-system relative abundance and colored based on the enrichment/depletion log10(p value). Circles indicate TWW-enrichment while triangles indicated TWW-depleted categories. b Heatmap of TWW-enriched or depleted (p value < 0.05) KEGG pathways (characterized by keggProfiler). Genes of interest, significantly enriched or depleted in TWW-irrigated root metagenomes, are highlighted and colored in pink c.

Analysis of enriched KEGG pathways (Fig. 4b, c, Table S9) and modules (Figure S3) revealed additional biological processes enriched in TWW-irrigated root microbiomes relative to FW-irrigated root microbiomes. Two-component systems were significantly enriched (40 KEGG genes were enriched, p value < 0.0001) in TWW-irrigated root microbiomes relative to FW-irrigated root microbiomes, including pathways involved in misfolded proteins, flagellar assembly, iron acquisition, and Mg2+ starvation (Fig. 4c). In addition, the denitrification gene module was significantly (4 genes, p value = 0.006) enriched in TWW-irrigated roots (Figure S3). Conversely, ABC transporter gene pathways associated with sugar (maltose, galactose, and oligogalacturonide), peptide (dipeptide and glutamate/aspartate), and nutrient (cobalt or nickel) transport were enriched (43 genes, p value < 0.0001) in FW-irrigated root microbiomes relative to TWW-irrigated roots. A type 3 secretion system (T3SS) gene module was also enriched (12 genes, p value < 0.0001) in FW-irrigated root communities.

Meta-analysis of selected gene counts relative to environmental variables from publicly available metagenomes

A meta-analysis was conducted to establish a global link between metagenome functional gene content and measured environmental variables. We focused on a subset of prominent genes from this study that were strongly positively correlated with pH (NQR, Na+–H+ antiporter) or with dissolved organic matter- DOC (periplasmic nitrate reductase, napAB, nitric oxide reductase- norBC, nitrous oxide reductase Z, nosZ). Metagenomes available at the Joint Genome Institute’s (JGI) Integrated Microbial Genomes and Microbiomes repository (n = 14,596) were screened. Environmental pH measurements were available for a subset of these metagenomes (n = 1588), and of this subset, 160 metagenomes had a total number of predicted genes greater than 100,000 (Table S10). Within these 160 metagenomes, the relative abundance of genes annotated as NQR (pairwise Wilcoxon rank test, Bonfferoni correction p valuepH<7: pH7–8 = 0.02; p valuepH7–8:pH > 8 < 0.0001; p valuepH < 7:pH > 8 = 0.002) and Na+–H+ antiporter (p valuepH < 7: pH7–8 < 0.0001; p valuepH < 7:pH > 8 < 0.0001) were strongly associated with higher measured pH values (Fig. 5a, b). Previously, it was suggested that TWW deliver higher levels of organic matter and this may lead to localized oxygen depletion [38]. Using the same filtering criteria, 257 metagenomes were identified with oxygen measurements and greater than 100,000 predicted genes. In these metagenomes, the abundance of napAB, norBC, and nosZ, in the denitrification pathway, were enriched in samples with lower measured oxygen (Fig. 5d–g). No such trend was observed for housekeeping genes such as gyrase B (gyrB, Fig. 5h). Salinity or Na+ concentrations were measured only in small subset of available metagenomes and were not analyzed further.

Fig. 5
figure5

Meta- analysis of selected pH and oxygen responsive genes. Publicly available (by JGI) subset of 160 environmental metagenomes with pH measurements, were screened. Box plot of gene counts in acidic pH (< 7), neutral (7–8), or alkaline (> 8) pH for a nqr operon (6 subunits), b Na-H antiporter operon (7 subunits), and cgyrB as control. Detailed pattern for each subunit is available in Figure S6, S7. Outliers are not displayed (1.5 × 0.25–0.75 quantiles). Significant differences in gene counts, by Wilcoxon rank sum test (Bonfferoni correction, P < 0.05), are marked in letter report (ac) Oxygen measurements were available for a subset of 257 environmental metagenomes. In these metagenomes, the abundance of napA (d), napB (e) norBC (f), and nosZ (g) was compared to oxygen levels. hgyrB was used as control. All gene counts are in proportion to rpoB gene

Microbial gene expression patterns associated with TWW irrigation (RNA-based metatranscriptomics)

An enrichment analysis of the root-associated microbial metatranscriptomes was performed to identify SEED-annotated genes and subsystems that were significantly differentially transcribed between plants irrigated with TWW relative to those irrigated with FW (Fig. 6, Table S11). In total, 10.1% of SEED-annotated genes were significantly differentially expressed in roots of TWW-irrigated plants relative to FW-irrigated plants (Wald test, FDR corrected p value < 0.05). Specifically, 7.2% of such genes had higher expression and 2.9% lower expression in TWW-irrigated roots relative to FW-irrigated roots. SEED-annotated genes were clustered into 761 SEED subsystems (level 2, based on SEED hierarchical clustering), and of these, 8 were over-represented in TWW-irrigated root microbial communities while only a single subsystem was significantly over-represented in the FW-irrigated root transcriptomes. The most highly and significantly overexpressed gene sub-systems in TWW-irrigated roots were NQR, TRAP transporters, sodium-hydrogen antiporters, alginate metabolism genes, and MSHA4 (Fig. 6a). All of these genes were also significantly enriched in metagenomic analysis of TWW-irrigated roots relative to FW-irrigated roots. Genes involved in alginate metabolism were only slightly enriched in metagenomes of TWW-irrigated roots (log2FC = 0.23, p value = 0.01) but were strongly over-expressed in TWW-irrigated metatranscriptomes relative to FW-irrigated roots (log2FC = 1.7, p value = 0.0004). Conversely, the overall expression level of hydrogenase subsystem genes was significantly higher in FW-irrigated roots relative to TWW-irrigated roots (log2FC = 1.6, p value > 0.0001), though their relative abundance at the DNA level was not substantially affected by irrigation treatment (Fig. 6b).

Fig. 6
figure6

Expressed functions associated with irrigation treatment. a Dot plot of Log2(fold-change) SEED subsystems enriched or depleted in TWW irrigated root metatranscriptomes. Significantly (FDR p < 0.05, represented by more than two gene families) enriched or depleted transcript abundance was computed using the goseq software package, with corrections for read abundance. Symbols are proportional to the sub-system relative abundance and colored based on the enrichment or depletion log10 (p value). Circles indicate TWW-enriched categories and triangles indicate TWW-depleted categories. b Differential metagenome enrichment (TWW/FW fold change) compared to differential metatranscriptome expression level. Highlighted (colored) categories significantly enriched or depleted in the metatranscriptome analysis. “ns” = “not significant.” Symbols are proportional to the log10(p value) enrichment in the metatranscriptome analysis. Numbers label the enriched category, as marked in a. c KEGG pathway-level enrichment or depletion in TWW-irrigation root metatranscriptomes (p value < 0.05), based on keggProfiler enrichment analysis. Significantly enriched or depleted gene clusters in TWW-irrigated roots are highlighted and colored in pink

Enrichment analysis of KEGG pathways (Fig. 6c) and modules (Figure S4) was performed (Table S12). The transcriptome of the microbial community of TWW-irrigated roots was significantly enriched in 2.75% of all KEGG-annotated genes. In addition, 8 KEGG pathways and 2 KEGG modules were significantly enriched in TWW-irrigated root metatranscriptomes relative to FW-irrigated root metatranscriptomes. In contrast, the transcriptome of the microbial community of FW-irrigated roots was significantly enriched in 1.8% of all KEGG-annotated genes. In addition, 3 KEGG pathways and 3 KEGG modules were significantly enriched in FW-irrigated root metatranscriptomes relative to TWW-irrigated root metatranscriptomes. This analysis revealed higher microbial relative expression of two-component systems, including the C4 dicarboxylate gene cluster, in TWW-irrigated roots relative to FW-irrigated roots (Fig. 6c). Moreover, in TWW-irrigated roots, higher relative expression of arginine and proline metabolism genes, particularly those in the arginine-to-spermidine pathway, was observed (Figure S5). The relative expression level of ABC transporter genes, including glutamate and galactose transporters, were higher in the metatranscriptomes of FW-irrigated roots relative to TWW-irrigated roots (10 enriched genes, p value = 0.0004). The level of expression of type 6 secretion system (T6SS) genes was most highly expressed under FW-irrigation conditions (7 enriched genes, p value < 0.0001).

Discussion

We previously studied the effect of TWW irrigation on soil and root microbial community structure and composition [26]. In that study, irrigation water quality and soil type were major explanatory variables for the observed soil microbial community structure and were of a similar magnitude. Similarly, the effect of irrigation water quality on root microbial community structure was of a similar magnitude to the plant host effect [26], demonstrating the responsiveness of the microbial community to both host and environmental factors. In the current study, we have attempted to harness the rhizoplane microbiome—existing at the interface between the plant and the surrounding soil—as a sensor for detecting in situ environmental conditions at the plant-soil interface, including factors leading to host stress. The main incentive in using the host-associated microbiome as a biosensor lies in the fact that high resolution is desired for accurate definition of the factors contributing to host physiological status [39]. Comparing the differences in the relative abundance of microbial genetic features (i.e., metagenome analysis) or expression of microbiome genes (i.e., metatranscriptome analysis) can aid in the identification of long-term stressors imposed on the host under these conditions as well as short-term stressors revealed by expression of genes processing environmental cues at time of sampling. Analyses can be performed at different levels of hierarchical gene annotation and can be performed using gene level annotation (e.g., SEED database) and enriched pathways or modules (e.g., KEGG annotation).

Most commonly used methods for studying root- soil interface employs microelectrodes [40], or specific dyed root imaging in “rhizoboxes” [41, 42]. Both methods measure only pre-defined variables, eliminating the possibility of discovering novel or unsuspected stressors. Moreover, experimental design forces manipulating natural environment by growing plants in designed cells or by removing plants from soil for further experimental procedure. Furthermore, in studies where transcriptional response is examined, plant host response is often tested under severe stress in unnatural short term experimental design [43]. We sought to be able to assess environmental factors leading to crop plants physiological status under more natural conditions. Here, we used non-model plant cultivars in semi-controlled environment, located under field conditions and exposed to moderate integrated stress, e.g., water quality, to demonstrate a possible application of microbiome as a bio-sensor revealing hidden environmental stresses experienced by the host.

A secondary motivating factor for the use of the microbiome as a biosensor lies in the observed low-resolution response of the host organism. In this study, upregulation of stress response genes was identified in the transcriptome of host roots irrigated with TWW relative to those roots irrigated with FW. However, the specific nature of the stresses remained unresolved, with transcriptome analysis revealing only the differential expression of genes involved in a highly general stress response associated with heat shock proteins [44, 45]. In fields, plants are exposed to myriad fluctuating biotic and abiotic environmental conditions, which force plants to tailor their gene transcriptional profiles. Therefore, individual abiotic stress responses extrapolated to plant experiencing multiple stress conditions should be treated with caution. The exact nature of the stress can rarely be predicted based on experimental profiling of individual stress response under regulated conditions [46, 47]. Moreover, other types of cellular regulation can also mediate plant stress responses [48, 49]. Such regulation is not as easily measured as gene expression. Higher resolution of sequences participating in plant transcriptional response (by increasing sequencing depth) can improve characterizing plant stress response, but here we demonstrated within a single sample both plant transcriptional response and equivalent microbiome response, compared at a similar sequencing resolution.

In contrast to the host, the genetic diversity of the host-associated microbiome is much greater [27], and the extraordinarily high microbial diversity in soils provides the plant with a wide selection of organisms competing for access to root exudates [18]. While the plant host can alter gene expression profiles in response to changing environmental conditions, both the membership of the root-associated microbial community and the expression patterns of the root-associated microbial community can be altered. Thus, the microbiome provides us with a highly dynamic and sensitive target with the potential for both short-term responsiveness (i.e., metatranscriptome) and long-term responsiveness (i.e., metagenome). In this study, we observed that in response to long-term irrigation with TWW, both the metagenome and metatranscriptome were significantly altered. Statistical analysis of microbial features lead to the identification of significantly differently abundant genes, gene transcripts, and pathways. Some genes of interest, being most significantly enriched, or with known and informative supported data, are presented in Fig. 7. Critically, the identification of the differentially abundant or expressed microbial features was consistent with the known key stresses imposed by TWW irrigation on the microbial community and the host.

Fig. 7
figure7

Conceptual model of hypothetical bacteria harboring physiological features (i.e., genes, pathways and modules) associated with water quality. Features enriched in a TWW- or b FW-irrigated root microbiomes. White symbols indicate features that are significantly enriched at the DNA level (metagenomes), grey features are highly expressed (metatranscriptomes), and green features are significantly abundant and expressed in one treatment relative to the other

Level of pH and salinity

Cytoplasmic pH in microorganisms must be at a range suitable for maintaining protein integrity. Most non- extremophiles bacterial cytoplasmic pH lies at a pH range of 7.4–7.8 [50]. In alkaline environments, organisms deploy various mechanisms to maintain intracellular pH and preserve electrochemical gradient in the presence of low proton concentration. To prevent proton loss in alkaline environments, an increase in cytoplasmic pH is achieved by reducing the activity of the proton pumping machinery of the cell respiratory chain. Under such conditions, some bacteria form a transmembrane sodium gradient, alternatively or concomitantly with a proton gradient [51]. In our plant system, the abundance and expression of Na+-transporting NADH:ubiquinone oxidoreductase (nqr) genes was significantly enriched under TWW irrigation relative to FW irrigation. In many alkaline environments, NQR constructs the primary sodium efflux system through oxidation of NADH and reduction of quinone. This process creates an electrochemical gradient of net negative charge in the cytosol [52], and the gradient is used by cation/proton antiporters (e.g., Na+/H+, and K+/H+) to exchange non-balanced movement of positive charge (H+) to the cell (more protons enter the cytoplasm as compared to the efflux of sodium or potassium ions [50, 52,53,54]. In a global mapping of soil bacterial communities, cation/proton antiporters were observed to be key genes overrepresented in dryland soil, presumably due to the high levels of salt and pH in arid soils [55]. Additionally, we measured an increase in the relative abundance and expression of tripartite ATP-independent periplasmic (TRAP) transporters in TWW roots, recently demonstrated to use a membrane-associated sodium gradient to facilitate transport of ligands [56,57,58]. Similarly, the increase in abundance in flagella assembly genes under TWW irrigation conditions may also suggest the use of sodium motive force for flagella performance [59, 60]. The signs for high pH stress obtained by the metagenome and metatranscriptome analysis are consistent with soil chemical analysis, as elevated pH conditions can result from long-term TWW irrigation [61].

Bacterial adaptation to alkaline conditions is frequently dependent on salt concentration, and elevated salt levels are also found in TWW and TWW-irrigated soils [62]. Elevated soil salinity can develop through long-term TWW irrigation, and can adversely affect protein and cell membrane stability. Commonly, microorganisms stabilize salt concentration in the cell by regulating cation proton antiporter activity [63]. Efflux of sodium ions by NQR activity and the activation of cation/proton transporters demonstrate that the TWW-irrigated root microbiome and the plant roots are indeed exposed to elevated salinity as compared to the FW-irrigated roots. This finding is consistent with the measurement in this study of higher levels of Na+ in the leaves of TWW-irrigated tomato and lettuce plants relative to FW-irrigated plants (Table S2) and in other plants [64].

The association between our significantly differentiated genes to processing specific environmental conditions is established in vitro in numerous studies [52, 65]. Here, we hypothesized that differentially abundant genes could be used as predictive markers of environmental cue, and this hypothesis is supported by meta-analyses demonstrating a link between single isolate studies and microbial communities in vivo.

Oxygen levels

Microbial gene content and expression patterns have great potential for identifying oxygen conditions in situ, as previously demonstrated [19]. In this study, we observed an enrichment in denitrification genes in TWW-irrigated root metagenomes relative to FW-irrigated metagenomes, possibly suggesting a lower oxygen concentration under TWW-irrigation [63]. However, the expression level of denitrification genes was not significantly higher in TWW-irrigated roots relative to FW-irrigated roots. This difference in enrichment between metagenome and metatranscriptome could be due to enrichment in the TWW-irrigated rhizoplane of facultative denitrifying microorganisms, with rhizosphere selection based on other physiological capabilities (e.g., motility). Conversely, the lack of enrichment in denitrification genes in TWW-irrigated metatranscriptomes could be a result of time of sampling. The root microbiome functional profile is expected to fluctuate by diurnal or hydration-dehydration cycles [66,67,68]. Therefore, gene abundance is indicative of the chronic, long term exposure to stress imposed by TWW, while expression levels may represent transient conditions. Plants in this study were subject to twice-daily irrigation and irrigation conditions in TWW deliver higher levels of organic matter and this may lead to localized oxygen depletion [38]. However, at the time of sampling, oxygen levels may have increased. Further short-term longitudinal analysis will be required to demonstrate diel- and irrigation-derived shifts in denitrification gene expression patterns.

Bacterial lifestyle

We observed a significant enrichment of genes associated with surface attachment, colonization and biofilm formation in TWW-irrigated root microbiome. These enriched microbial features included genes encoding for flagella and MSHA type 4 pili; both features have been previously demonstrated to facilitate near-surface motility and bacterial attachment [69]. Furthermore, an increase in the relative abundance and expression of alginate producing genes, which catalyze the formation of extracellular polysaccharide matrix in biofilms of many bacterial clades [70], was observed. These results further point to the critical importance of both motility and attachment for rhizoplane microorganisms, as has been previously indicated [71]. The reason for significant enrichment of these genes in TWW-irrigated roots relative to FW-irrigated roots is not entirely clear but may be due to overall elevated organic matter in TWW [72]. While the focus in this study has been largely on genes enriched in TWW-irrigated systems, several key ABC transporters were depleted in the metagenomes of roots of both plant types irrigated with TWW relative to those irrigated with FW. These genes include transporters for oligogalacturonide, maltose, general sugar, and glutamate. The change in relative abundance of these transporters was also observed in other plants [19], and may indicate differential pattern of root deposits. Plant glutamate secretion patterns have been shown to be mediated by external cues such as salinity, oxidative stress, and availability of nutrients [73]. Glutamate ABC transporters were depleted, while glutamine was enriched in TWW-irrigated roots.

Conclusions

We hypothesized that the microbial functional gene profiles and expression patterns can serve as in vivo sensors of environmental factors affecting hosts and host-associated microbial communities. As the host and its microbiome are similarly exposed to environmental conditions, genetic profiling, and expression analysis of the microbiome may be used as a predictive tool to identify stresses affecting hosts. It is beyond the scope of the current study to determine which of the abiotic factor is directly influencing the host (and subsequently reducing plant yield), but the suggested methodology can better elucidate the environmental conditions encountered by the host. Environmental surveys or host-associated microbiome analyses frequently yield contradictory or context-dependent results, making the predictive power of such observations inconclusive. Studying the microbiome as a functional unit reacting to a specific environment, however, constitutes a non-deterministic approach, thereby eliminating the need for marker features (e.g., genes, pathways or specific taxa) associated with specific conditions. In this study, we employ a well-defined plant host-microbiome system under experimental treatment, but this approach may be used to define microscale conditions in other host systems, potentially reveal other host physiological stressors.

Materials and methods

Experimental design: mesocosm scale experiment

Lysimeters (0.5 m3), located at Kiryat-Gat-Lachish agricultural research station, northern Negev, Israel (31.605760, 34.791179), were filled with loamy sand soil collected from agricultural fields in the western Negev, Israel (31.351722, 34.403471). Prior to this experiment, the lysimeters were repeatedly irrigated for 8 summers with fresh water (FW) or tertiary treated wastewater (TWW), derived from the Kiryat Gat wastewater treatment plant (WWTP), and rain-fed during winters. Tomato (Solanum lycopersicum-Heinz 4107) and lettuce (Lactuca sativa-Romaine-Assaph) were grown in the lysimeters for 98 and 42 days over two consecutive summers (2014, 2015), respectively, and irrigated with either fresh water or TWW. Each replicate comprised composite sample, consisting of 2–4 plants collected from a single lysimeter, with three replicates (individual lysimeters) per treatment. An exception was FW-irrigated lettuce plants which were collected from only two lysimeters, yet the three replicates were still composites of two plants each. At harvest, roots were collected, vigorously washed, dried, and frozen on site for further procedures. Detailed procedures, including soil and plant measurements, were described previously [26].

DNA and RNA extraction

Standard phenol-chloroform nucleic acid extraction protocol was employed for DNA and RNA isolation [26, 74]. In brief, 0.2 gr of homogenized roots were moderately bead beaten for 45 s at low speed (4.5 m/s) by Fast Prep FP120 (Savant Instruments Inc., Holbrook, NY, USA) with phenol, phosphate buffer pH 8 (with additional 10 μl ml−1 β-mercaptoethanol -Sigma-Aldrich, St Louis, MO, USA) and 1.25% CTAB (hexadecyltrimethylammonium bromide, Sigma Aldrich). Following phenol-chloroform wash, nucleic acids were precipitated with polyethylene glycol (PEG) and ethanol. Nucleic acids were split for DNA and RNA isolation. RNA samples were treated with RQ1-DNase (Promega, Madison, WI, USA) and complete DNA removal from RNA samples was validated by real-time reverse transcription PCR. RNA integrity was evaluated with Agilent TapeStation (Santa Clara, CA, USA). Ribosomal RNAs were removed using the Ribo-Zero rRNA Removal Kit (Illumina, San Diego, CA, USA), combining bacteria and plant probes. Double- strand complementary DNA (cDNA) synthesis was conducted by Maxima H Minus Reverse Transcriptase (Thermo Fisher Scientific, Waltham, MA, USA).

Library prep and sequencing

Shotgun metagenome libraries were generated using a Nextera XT library preparation kit according to the manufacturer’s instructions (Illumina). Complementary DNA for transcriptome analysis was sheared using a Covaris S2 acoustic device, and libraries were generated using a Accel-NGS 1S Plus DNA Library Kit (Swift Biosciences, Ann Arbor, MI) according to the manufacturer’s instructions. Libraries were pooled sequenced using high-output flow cells with paired-end 2 × 150 base reads on an Illumina NextSeq500 sequencer. Library preparation and sequencing was performed at the University of Illinois at Chicago Sequencing Core (UICSQC).

Bioinformatic analysis

Quality control of raw double-strand FASTQ sequences was evaluated by FASTQC software [75], and adjusted by Trimmomatic [76] with customized parameters set to: SLIDINGWINDOW:4:15 MINLEN:100 CROP:145 HEADCROP:15.

Metagenome analysis

For shotgun metagenome analyses, 25–49 million DNA sequences (paired-end) were generated per sample (supplementary Table S1). In tomato roots, 13 to 28% of the reads mapped to the host genome, with the remaining reads were attributed to the microbiome. In lettuce, 55–69% mapped to the host with the remaining reads attributed to the microbiome. Host sequences were removed by comparing quality checked reads to host genome (tomato or lettuce) with bowtie2 [77], and subsequently removing the reads with SAMtools [78]. Metagenomics reads from all three replicates were de novo assembled together with metaSPAdes [79]. This assembly yielded 1,760,490 contigs larger than 500bp, and this assembly had an N50 value greater than 1600 bases. Gene prediction was performed on scaffolds using the software package Prodigal [80], yielding 6,422,376 predicted genes. The predicted genes from all samples were combined, and a non-redundant gene catalog was established with 5,359,885 genes (based on 95% similarity), using CD-HIT [81].The gene catalog was aligned to the NCBI non-redundant protein database using the software package DIAMOND in sensitive mode [82]. Sequence annotation (theSEED [28] and Kyoto Encyclopedia of Genes and Genomes-KEGG [29]) and predicted taxonomy were achieved with MEGAN V6 [83]. To attain count data (number of mapped reads for each gene), quality checked reads (after host read removal) were aligned to the annotated gene catalog by bowtie2, while analogous read annotated terms were summed using a custom python script

Metatranscriptome analysis

Shotgun metatranscriptome analyses generated 27–33 million sequences (paired-end) per sample (Table S1). Quality checked RNA reads were aligned to the gene catalog established from the metagenomics analysis, in a similar fashion to metagenomics count data (microbial transcriptome, with 11–51% of the sequences mapped to the microbiome, and an average of 27%).

Host RNAseq

Metatranscriptome reads were mapped to the plant genome (host transcriptome analysis, with 7–41% of the sequences mapped, and an average of 32%). An estimation of transcript abundance for tomato root samples was obtained by aligning quality checked sequences (prior to host reads removal) to the predicted Solanum lycopersicum transcripts with Trinity RSEM transcript quantification method [84]. As the lettuce genome is not fully annotated, sequence data generated from lettuce plants were screened for orthologs of known tomato genes. Lettuce transcripts first predicted by Tophat and cufflinks for transcript prediction [85], than infered to ortholog tomato genes by OrthoFinder. Transcript quantification was done following similar analysis as for tomato samples.

Statistical analysis

Metagenome and metatranscriptome statistical enriched gene list (SEED or KEGG annotated) or taxonomic groups were obtained by DESeq2 [36] and compared using the VennDiagram R package [86]. DESeq2 assume most genes are shared between the compared treatments. Therefore, the plant host effects might be an overestimation. Nevertheless, the main focus remained the irrigation derived functional shift, with values at the range suitable for this analysis. Taxonomic trees were visualized using the interactive tree of life [87] and applying the least common ancestor MEGAN algorithm. SEED subsystems enrichment analysis was conducted with the ‘R’ goseq package [88], normalizing to SEED counts. KEGG pathway and module enrichment were analyzed by clusterProfiler package in R [89]. Statistical test (MANOVA, ANOSIM) were conducted in R package ‘vegan’ [90], and figures were plotted with R ‘ggplot2’ [91] or ‘pheatmap’ [92]. Count read data was normalized by trimmed mean of M values (TMM) using the EdgeR ‘R’ package [93]. Differentially expressed host transcripts were obtained using the EdgeR ‘R’ package [93], followed by annotation and visualization using the STRING network [94, 95] Cytoscape integrated application [96] for both plant hosts combined. The minimum required interaction score was customized to medium confidence (0.4), and PFAM protein domain enrichment was set to a false discovery rate p value of 0.05.

Availability of data and materials

The sequences generated and analyzed during the current study were uploaded to the NCBI Sequence Read Archive (SRA) data depository, with the project number PRJNA602301. Additional data can be shared upon request.

References

  1. 1.

    Sunagawa S, et al. Structure and function of the global ocean microbiome. Science. 2015;348(6237):1261359.

  2. 2.

    Fierer N, Jackson RB. The diversity and biogeography of soil bacterial communities. Proc Natl Acad Sci U S A. 2006;103(3):626–31.

  3. 3.

    Mendes LW, et al. Soil-borne microbiome: linking diversity to function. Microb Ecol. 2015;70(1):255–65.

  4. 4.

    Fierer N. Embracing the unknown: disentangling the complexities of the soil microbiome. Nat Rev Microbiol. 2017;15(10):579–90.

  5. 5.

    Peiffer JA, et al. Diversity and heritability of the maize rhizosphere microbiome under field conditions. Proc Natl Acad Sci U S A. 2013;110(16):6548–53.

  6. 6.

    Lundberg DS, et al. Defining the core Arabidopsis thaliana root microbiome. Nature. 2012;488(7409):86–90.

  7. 7.

    Mendes LW, Kuramae EE, Navarrete AA, van Veen JA, Tsai SM. Taxonomical and functional microbial community selection in soybean rhizosphere. ISME J. 2014;8(8):1577–87.

  8. 8.

    Bulgarelli D, et al. Revealing structure and assembly cues for Arabidopsis root-inhabiting bacterial microbiota. Nature. 2012;488(7409):91–5.

  9. 9.

    Edwards J, et al. Structure, variation, and assembly of the root-associated microbiomes of rice. Proc Natl Acad Sci U S A. 2015;112(8):E911–20.

  10. 10.

    Turnbaugh PJ, et al. The human microbiome project. Nature. 2007;449(7164):804–10.

  11. 11.

    Dewhirst FE, et al. The human oral microbiome. J Bacteriol. 2010;192(19):5002–17.

  12. 12.

    Goodrich JK, et al. Human genetics shape the gut microbiome. Cell. 2014;159(4):789–99.

  13. 13.

    Grice EA, Segre JA. The skin microbiome. Nat Rev Microbiol. 2011;9(4):244–53.

  14. 14.

    Ofek M, Voronov-Goldman M, Hadar Y, Minz D. Host signature effect on plant root-associated microbiomes revealed through analyses of resident vs. active communities. Environ Microbiol. 2014;16(7):2157–67.

  15. 15.

    Walters WA, et al. Large-scale replicated field study of maize rhizosphere identifies heritable microbes. Proc Natl Acad Sci U S A. 2018;115(28):7368–73.

  16. 16.

    Torsvik V, Øvreås L, Thingstad TF. Prokaryotic diversity--magnitude, dynamics, and controlling factors. Science. 2002;296(5570):1064–6.

  17. 17.

    Tringe SG, et al. Comparative metagenomics of microbial communities. Science. 2005;308(5721):554–7.

  18. 18.

    Sasse J, Martinoia E, Northen T. Feed your friends: do plant exudates shape the root microbiome? Trends Plant Sci. 2018;23(1). https://doi.org/10.1016/j.tplants.2017.09.003.

  19. 19.

    Ofek-Lalzar M, et al. Niche and host-associated functional signatures of the root surface microbiome. Nat Commun. 2014;5(1):4950.

  20. 20.

    Edwards JA, et al. Compositional shifts in root-associated bacterial and archaeal microbiota track the plant life cycle in field-grown rice. PLoS Biol. 2018;16(2):e2003862.

  21. 21.

    Chaparro JM, Badri DV, Vivanco JM. Rhizosphere microbiome assemblage is affected by plant development. ISME J. 2014;8(4):790–803.

  22. 22.

    Castrillo G, et al. Root microbiota drive direct integration of phosphate stress and immunity. Nature. 2017;543(7646):513–8.

  23. 23.

    Pett-Ridge J, Firestone MK. Redox fluctuation structures microbial communities in a wet tropical soil. Appl Environ Microbiol. 2005;71(11):6998–7007.

  24. 24.

    Staley C, et al. Urea amendment decreases microbial diversity and selects for specific nitrifying strains in eight contrasting agricultural soils. Front Microbiol. 2018;9:634.

  25. 25.

    Green SJ, Inbar E, Michel FC, Hadar Y, Minz D. Succession of bacterial communities during early plant development: transition from seed to root and effect of compost amendment. Appl Environ Microbiol. 2006;72(6):3975–83.

  26. 26.

    Zolti A, Green SJ, Ben Mordechay E, Hadar Y, Minz D. Root microbiome response to treated wastewater irrigation. Sci Total Environ. 2019;655:899–907.

  27. 27.

    Zilber-Rosenberg I, Rosenberg E. Role of microorganisms in the evolution of animals and plants: the hologenome theory of evolution. FEMS Microbiol Rev. 2008;32(5):723–35.

  28. 28.

    Overbeek R, et al. The subsystems approach to genome annotation and its use in the project to annotate 1000 genomes. Nucleic Acids Res. 2005;33(17):5691–702.

  29. 29.

    Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

  30. 30.

    Wang W, Vinocur B, Shoseyov O, Altman A. Role of plant heat-shock proteins and molecular chaperones in the abiotic stress response. Trends Plant Sci. 2004;9(5):244–52.

  31. 31.

    Swindell WR, Huebner M, Weber AP. Transcriptional profiling of Arabidopsis heat shock proteins and transcription factors reveals extensive overlap between heat and non-heat stress response pathways. BMC Genomics. 2007;8(1):125.

  32. 32.

    Wang C, Li J, Yuan M. Salt tolerance requires cortical microtubule reorganization in Arabidopsis. Plant Cell Physiol. 2007;48(11):1534–47.

  33. 33.

    Abdrakhamanova A, Wang QY, Khokhlova L, Nick P. Is microtubule disassembly a trigger for cold acclimation? Plant Cell Physiol. 2003;44(7):676–86.

  34. 34.

    Sivaguru M, Pike S, Gassmann W, Baskin TI. Aluminum rapidly depolymerizes cortical microtubules and depolarizes the plasma membrane: evidence that these responses are mediated by a glutamate receptor. Plant Cell Physiol. 2003;44(7):667–75.

  35. 35.

    Takemoto D, Hardham AR. The cytoskeleton as a regulator and target of biotic interactions in plants. Plant Physiol. 2004;136(4):3864–76.

  36. 36.

    Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

  37. 37.

    Juárez O, Barquera B. Insights into the mechanism of electron transfer and sodium translocation of the Na+-pumping NADH:quinone oxidoreductase. Biochim Biophys Acta Bioenerg. 2012;1817(10):1823–32.

  38. 38.

    Henderson SL, et al. Changes in denitrifier abundance, denitrification gene mRNA levels, nitrous oxide emissions, and denitrification in anoxic soil microcosms amended with glucose and plant residues. Appl Environ Microbiol. 2010;76(7):2155–64.

  39. 39.

    Jones DL, Nguyen C, Finlay RD. Carbon flow in the rhizosphere: carbon trading at the soil–root interface. Plant Soil. 2009;321(1–2):5–33.

  40. 40.

    Pang JY, Newman I, Mendham N, Zhou M, Shabala S. Microelectrode ion and O2 fluxes measurements reveal differential sensitivity of barley root tissues to hypoxia. Plant Cell Environ. 2006;29(6):1107–21.

  41. 41.

    Blossfeld S, Schreiber CM, Liebsch G, Kuhn AJ, Hinsinger P. Quantitative imaging of rhizosphere pH and CO2 dynamics with planar optodes. Ann Bot. 2013;112(2):267–76.

  42. 42.

    Blossfeld S, Gansert D, Thiele B, Kuhn AJ, Lösch R. The dynamics of oxygen concentration, pH value, and organic acids in the rhizosphere of Juncus spp. Soil Biol Biochem. 2011;43(6):1186–97.

  43. 43.

    Munns R. Genes and salt tolerance: bringing them together. New Phytol. 2005;167(3):645–63.

  44. 44.

    Nishizawa A, et al. Arabidopsis heat shock transcription factor A2 as a key regulator in response to several types of environmental stress. Plant J. 2006;48(4):535–47.

  45. 45.

    Bhardwaj AR, et al. Global insights into high temperature and drought stress regulated genes by RNA-Seq in economically important oilseed crop Brassica juncea. BMC Plant Biol. 2015;15(1):9.

  46. 46.

    Atkinson NJ, Urwin PE. The interaction of plant biotic and abiotic stresses: from genes to the field. J Exp Bot. 2012;63(10):3523–43.

  47. 47.

    Mittler R. Abiotic stress, the field environment and stress combination. Trends Plant Sci. 2006;11(1):15–9.

  48. 48.

    Sunkar R, Zhu J-K. Novel and stress-regulated microRNAs and other small RNAs from Arabidopsis. Plant Cell. 2004;16(8):2001–19.

  49. 49.

    Vlad F, et al. Protein phosphatases 2C regulate the activation of the Snf1-related kinase OST1 by abscisic acid in Arabidopsis. Plant Cell. 2009;21(10):3170–84.

  50. 50.

    Padan E, Bibi E, Ito M, Krulwich TA. Alkaline pH homeostasis in bacteria: new insights. Biochim Biophys Acta Biomembr. 2005;1717(2):67–88.

  51. 51.

    Krulwich TA, Sachs G, Padan E. Molecular aspects of bacterial pH sensing and homeostasis. Nat Rev Microbiol. 2011;9(5):330–43.

  52. 52.

    Vorburger T, et al. Role of the Na+-translocating NADH:quinone oxidoreductase in voltage generation and Na+ extrusion in vibrio cholerae. Biochim Biophys Acta Bioenerg. 2016;1857(4):473–82.

  53. 53.

    Ito M, Guffanti AA, Oudega B, Krulwich TA. mrp, a multigene, multifunctional locus in Bacillus subtilis with roles in resistance to cholate and to Na+ and in pH homeostasis. J Bacteriol. 1999;181(8):2394–402.

  54. 54.

    Yamaguchi T, Tsutsumi F, Putnoky P, Fukuhara M, Nakamura T. pH-dependent regulation of the multi-subunit cation/proton antiporter Pha1 system from Sinorhizobium meliloti. Microbiology. 2009;155(8):2750–6.

  55. 55.

    Delgado-Baquerizo M, et al. A global atlas of the dominant bacteria found in soil. Science. 2018;359(6373):320–5.

  56. 56.

    Mulligan C, Kelly DJ, Thomas GH. Tripartite ATP-independent periplasmic transporters: application of a relational database for genome-wide analysis of transporter gene frequency and organization. J Mol Microbiol Biotechnol. 2007;12(3–4):218–26.

  57. 57.

    Mulligan C, et al. The substrate-binding protein imposes directionality on an electrochemical sodium gradient-driven TRAP transporter. Proc Natl Acad Sci U S A. 2009;106(6):1778–83.

  58. 58.

    Mulligan C, Fischer M, Thomas GH. Tripartite ATP-independent periplasmic (TRAP) transporters in bacteria and archaea. FEMS Microbiol Rev. 2011;35(1):68–86.

  59. 59.

    Atsumi T, McCartert L, Imae Y. Polar and lateral flagellar motors of marine vibrio are driven by different ion-motive forces. Nature. 1992;355(6356):182–4.

  60. 60.

    Yorimitsu T, Homma M. Na+-driven flagellar motor of vibrio. Biochim Biophys Acta Bioenerg. 2001;1505(1):82–93.

  61. 61.

    Tarchouna LG, Merdy P, Raynaud M, Pfeifer HR, Lucas Y. Effects of long-term irrigation with treated wastewater. Part I: evolution of soil physico-chemical properties. Appl Geochem. 2010;25(11):1703–10.

  62. 62.

    Assouline S, Narkis K (2013) Effect of long-term irrigation with treated wastewater on the root zone environment. Vadose Zo J 12(2):0.

  63. 63.

    Oren A. Life at high salt concentrations. The Prokaryotes (springer New York, New York, NY), pp 263–282. 2006.

  64. 64.

    Paudel I, et al. Treated wastewater irrigation: soil variables and grapefruit tree performance. Agric Water Manag. 2018;204:126–37.

  65. 65.

    Dalsgaard T, et al. Oxygen at nanomolar levels reversibly suppresses process rates and gene expression in anammox and denitrification in the oxygen minimum zone off northern Chile. MBio. 2014;5(6):e01966.

  66. 66.

    Staley C, et al. Diurnal cycling of rhizosphere bacterial communities is associated with shifts in carbon metabolism. Microbiome. 2017;5(1):65.

  67. 67.

    Naylor D, DeGraaf S, Purdom E, Coleman-Derr D. Drought and host selection influence bacterial community dynamics in the grass root microbiome. ISME J. 2017;11(12):2691–704.

  68. 68.

    Šťovíček A, Kim M, Or D, Gillor O. Microbial community response to hydration-desiccation cycles in desert soil. Sci Rep. 2017;7(1):45735.

  69. 69.

    Utada AS, et al. Vibrio cholerae use pili and flagella synergistically to effect motility switching and conditional surface attachment. Nat Commun. 2014;5(1):4913.

  70. 70.

    O’Toole G, Kaplan HB, Kolter R. Biofilm formation as microbial development. Annu Rev Microbiol. 2000;54(1):49–79.

  71. 71.

    Rodríguez-Navarro DN, Dardanelli MS, Ruíz-Saínz JE (2007) Attachment of bacteria to the roots of higher plants. FEMS Microbiol Lett 272(2):127–136.

  72. 72.

    Jimenez-Sanchez C, Wick LY, Cantos M, Ortega-Calvo J-J. Impact of dissolved organic matter on bacterial tactic motility, attachment, and transport. Environ Sci Technol. 2015;49(7):4498–505.

  73. 73.

    Forde BG, Lea PJ. Glutamate in plants: metabolism, regulation, and signalling. J Exp Bot. 2007;58(9):2339–58.

  74. 74.

    Angel R, Claus P, Conrad R. Methanogenic archaea are globally ubiquitous in aerated soils and become active under wet anoxic conditions. ISME J. 2012;6(4):847–62.

  75. 75.

    Andrews S (2010) FastQC: a quality control tool for high throughput sequence data. Available at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.

  76. 76.

    Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.

  77. 77.

    Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9.

  78. 78.

    Li H, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.

  79. 79.

    Nurk S, Meleshko D, Korobeynikov A, Pevzner PA. metaSPAdes: a new versatile metagenomic assembler. Genome Res. 2017;27(5):824–34.

  80. 80.

    Hyatt D, et al. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics. 2010;11(1):119.

  81. 81.

    Fu L, Niu B, Zhu Z, Wu S, Li W. CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics. 2012;28(23):3150–2.

  82. 82.

    Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2015;12(1):59–60.

  83. 83.

    Huson DH, Auch AF, Qi J, Schuster SC. MEGAN analysis of metagenomic data. Genome Res. 2007;17(3):377–86.

  84. 84.

    Haas BJ, et al. De novo transcript sequence reconstruction from RNA-seq using the trinity platform for reference generation and analysis. Nat Protoc. 2013;8(8):1494–512.

  85. 85.

    Cantarel BL, et al. MAKER: an easy-to-use annotation pipeline designed for emerging model organism genomes. Genome Res. 2008;18(1):188–96.

  86. 86.

    Chen H, Boutros PC. VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinformatics. 2011;12(1):35.

  87. 87.

    Letunic I, Bork P. Interactive tree of life (iTOL) v3: an online tool for the display and annotation of phylogenetic and other trees. Nucleic Acids Res. 2016;44(W1):W242–5.

  88. 88.

    Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):R14.

  89. 89.

    Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. Omi A J Integr Biol. 2012;16(5):284–7.

  90. 90.

    Dixon P. VEGAN, a package of R functions for community ecology. J Veg Sci. 2003;14(6):927–30.

  91. 91.

    Wickham H (2016) ggplot2: elegant graphics for data analysis (Springer-Verlag New York) Available at: https://ggplot2.tidyverse.org.

  92. 92.

    Kolde R. Pheatmap: pretty heatmaps [software]. 2015.

  93. 93.

    Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

  94. 94.

    Szklarczyk D, et al. The STRING database in 2017: quality-controlled protein–protein association networks, made broadly accessible. Nucleic Acids Res. 2017;45(D1):D362–8.

  95. 95.

    Szklarczyk D, et al. STRING v10: protein–protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015;43(D1):D447–52.

  96. 96.

    Shannon P, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

Download references

Acknowledgements

The authors thank Dr. Minz lab members, mainly Dr. Nesli Tovi and Dr. Sammy Frenk, from the Agricultural Research Organization–Volcani Center. The authors thank Dr. Jonathan Friedman from the Faculty of Agriculture, Hebrew University, for his critical reading and feedback. The authors would also like to thank Dr. Maya Ofek Lalzar from Haifa University, for her useful suggestions.

Funding

This research was supported by research grant no. IS-4662-13 from the Binational Agricultural Research& Development Fund (BARD), research grant no. 821-0142-13 from the Israel Ministry of Agriculture and Rural Development and USAID-MERC research grant no. M34-011.

Author information

Affiliations

Authors

Contributions

DM and YH conceived and designed the work. SJG lead all library preparation and sequencing procedure and contributed to plan, design and choose all the methods used. SJG substantively revised the manuscript. NS provided bioinformatics support.

AZ preformed sampling, lab procedures, bioinformatic analysis, statistical analysis, and data visualization. All co-authors contributed in writing and editing and approved the final manuscript.

Corresponding author

Correspondence to Dror Minz.

Ethics declarations

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Competing interests

The authors declare that they have no competing interests

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Figure S1. Non-metric multidimensional scaling (nMDS) based on the Bray-Curtis dissimilarity index of the SEED annotated gene counts for (a) metagenome, and (b) metatranscriptome. Figure S2. Correlation between metagenomics enriched SEED categories and soil environmental parameters, for (a) subsystems enriched in TWW irrigated roots, or (b) enriched in FW irrigated roots. Figure S3. Metagenomic enriched KEGG modules in (a) FW irrigated roots, or (b) TWW irrigated roots. The bar chart present the log10(P value) of the enrichment analysis calculated by KeggProfiler ‘R’ package. Figure S4. Metatransctiptome enrichement KEGG modules in (a) FW irrigated roots, or (b) TWW irrigated roots. The bar chart present the log10(P value) of the enrichment analysis calculated by KeggProfiler ‘R’ package. Figure S5. Microbial gene expression patterns revealed TWW irrigation enrichment of arginine and proline metabolism genes, highlighting those in the arginine-to-spermidine pathway (colored in pink). Figure S6. Meta- analysis of the full nqr operon in environmental metagenomes. Box plot demonstrates the nqr operon (a-f subunits) gene abundance in acidic pH (<7), neutral (7-8), or alkaline (>8) pH. Figure S7. Meta- analysis of the full mnh operon in environmental metagenomes. Box plot demonstrates the mnh operon (a-g subunits) gene abundance in acidic pH (<7), neutral (7-8), or alkaline (>8) pH.

Additional file 2: Table S1. Sequencing data: number of reads, contigs, N(50), percentage of mapped reads.

Additional file 3: Table S2. Plant variables: yield, leaves cation concentration.

Additional file 4: Table S3. Plant RNAseq analysis with detailed EdgeR table.

Additional file 5: Table S4. Taxonomic affiliation: taxonomy table. Including: count table, TMM normalized counts table, Deseq2.

Additional file 6: Table S5. Tree file, based on MEGAN6 LCA algorithm.

Additional file 7: Table S6. Soil environmental variables: Na, K+, pH, EC, DOC, TN.

Additional file 8: Table S7. TheSEED metagenomics (DNA based) analysis. Including: count table, TMM normalized table, DEseq2 result, goseq enrichment analysis.

Additional file 9: Table S8. Correlation between enriched SEED subsystems and environmental variables (pH, DOC, Na+, K+, DN).

Additional file 10: Table S9. KEGG metagenomic (DNA based) analysis. Including: count table, TMM normalized table, DEseq2 result, pathway and module enrichment analysis.

Additional file 11: Table S10. Meta-analysis: JGI detailed table- including: projects id, habitat, pH concentration, oxygen concentration, gene counts.

Additional file 12: Table S11. TheSEED metatranscriptome (RNA based) analysis. Including: count table, TMM normalized table, DEseq2 result, goseq enrichment analysis.

Additional file 13: Table S12. KEGG metatranscriptome (RNA based) analysis. Including: count table, TMM normalized table, DEseq2 result, pathway and module enrichment analysis.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Zolti, A., Green, S.J., Sela, N. et al. The microbiome as a biosensor: functional profiles elucidate hidden stress in hosts. Microbiome 8, 71 (2020). https://doi.org/10.1186/s40168-020-00850-9

Download citation

Keywords

  • Microbiome
  • Plant microbe interaction
  • Host microbe interaction
  • Biosensor