The metagenomics-derived gene catalogs used for the current analysis was previously reported in Ofek-Lalzar et al. (2014) and is publically available (BioProject accession number PRJNA208116). In brief, the DNA data were extracted from two agricultural crops: wheat and cucumber. For each crop plant, samples were taken from rhizosphere – the area under the direct influence of the root, and the more distant soil not under direct effect, termed here root and soil samples, respectively. The replicated experiment included 10 samples in total: root samples were in triplicates (a total of six root samples) and soil samples in duplicates (a total of four samples). Reproducibility between replicas was tested and reported, clearly demonstrating higher variance between root and soil treatments. The data were sequenced, annotated and mapped to taxonomic bins and KEGG ortholog groups. Taxonomic assignments were done using the lowest common ancestor algorithm, MEGAN (version 4.0) (Huson et al., 2007). The gene catalog represents approximately 71% of cucumber root reads, 50% of wheat root reads and 34% of soil reads. DNA reads were mapped to different taxonomic ranks. Approximately 72, 63, and 47% of the reads in soil cucumber and wheat samples were mapped to bacteria, while 22, 16, and 24% were ‘not assigned,’ respectively. Generally, out of the total number of reads mapped to bacteria, over 97% of the reads mapped to the bacterial level were assigned to order level. Overall, this gene catalog provides a description of genes detected in each of the four treatments (cucumber root and soil; wheat root and soil), their relative abundance, functional annotation (e.g., KEGG assignment), and taxonomic origin. This gene catalog was used as a starting point for network construction and subsequent analyses.

KEGG ortholog groups associated with enzymatic functions were detected across the four environments. Differential abundance (root vs. soil) of enzyme-associated reads was determined independently for wheat and cucumber using the EdgeR R package (Robinson et al., 2010) under the set of conditions previously described by Ofek-Lalzar et al. (2014). Differential abundance between environments was considered significant if the difference was greater than two fold and the FDR-adjusted p-value was < 0.01, requiring consistency between replicas. Overall, the differential abundance analysis produced four sets of differentially abundant enzyme sets describing: cucumber soil, cucumber root, wheat soil and wheat root (Supplementary Tables S1, S2 and Figure S1). In each environment, the set of differentially abundant enzymes was used for construction of an environment-specific network, following the procedure outlined in Kreimer et al. (2012). Similarly, a meta-network was constructed, containing all enzymatic functions annotated across the metagenomic data.

Using the NetSeed algorithm (Carr and Borenstein, 2012), through its implantation in NetCmpt (Kreimer et al., 2012), an approximation of the relevant metabolic environment was retrieved for each of the five networks (meta-network and four environment-specific networks). Based on network topology, the algorithm provided a list of metabolites that were predicted to be externally consumed from the environment, termed here ‘source metabolites.’ Since the four environment-specific networks were constructed from differentially abundant enzymes only, they were highly fragmented, leading to prediction of artificial source-metabolites (Supplementary Figure S2). Source-metabolites, identified for each environment specific network, were hence compared to the source-metabolite set identified for the meta-network. Only source-metabolites present in both sets were further considered. Following this filtration, we complied four environment-specific sets of source metabolites providing an approximation of the metabolic content in the corresponding environments (root and soil of wheat and cucumber, a total of four environments, Supplementary Tables S3, S4 and Figure S1).

To predict metabolic activities in each environment we made use of the Expansion algorithm (Ebenhoh et al., 2004). Briefly, the algorithm allows the predicting of an active metabolic network (expanded) given a pre-defined set of substrates and reactions. The algorithm starts with a set of source-metabolites acting as substrates; it scans the reaction bank for feasible reactions for which all the possible substrates exits; all feasible reactions are added to the network, their products being the substrates for the next set of reactions. The network stops expanding when no feasible reactions are found. Thus, the full expansion of the network reflects both the reaction repertoire and the primary set of compounds (source-metabolites). Here, simulations of environmental activity were carried by expanding the meta-network (the full set of reaction detected across all samples) four times – each time using the environmental specific source metabolite set. We made use of the full set of reactions for all simulations, since despite differences in abundance, almost all enzymes were detected in all samples (Ofek-Lalzar et al., 2014). The environment-specific expanded networks are provided at Supplementary Table S5.

All sequenced reads, collected in the metagenomic dataset, were linked to taxonomic groups, in the order level, using mapping from previous analysis (Ofek-Lalzar et al., 2014). Each read was assigned a Gene Id which was used as key parameter. Enzymes were then linked to the taxonomy mapping through gene IDs. The Simpson index (Heip et al., 1998), typically used to determine species dominance in ecological surveys, was newly applied here to determine the dominance of specific taxonomic groups in regards to a function. To this end, instead of looking at the frequencies of species in a sample (as in ecological surveys), for each enzyme (equivalent to a sample), we looked at the distribution the taxonomic affiliations of its associated reads. Accordingly – a low score denoted a function carried by many taxonomic groups; high score denoted a function carried out by single or few groups. To this end, Simpson Indices were calculated for each enzyme in the dataset for each of the original 10 samples (Supplementary Table S6). Replicates show similarity in the dominance/diversity indices (Supplementary Figure S3). Then, an environmental Simpson index value was determined for each enzyme by calculating mean values across corresponding samples. Finally, within each environment, we described a function to be dominated by a taxonomic group if: (i) the environmental Simpson index value was greater than 0.4 – mean dominance value across all enzymes (Supplementary Table S6) typically also associated with low diversity (Supplementary Figure S3); (ii) the same taxonomic group is dominate in all replicate samples. Hence, associations between taxonomy are representative of a treatment (environment) and consistent between replicas. Most (556) of the dominant enzymes were associated with a single taxonomic group (that is, dominance by a single taxonomic groups in different environments; Supplementary Table S7). For each environment, we ranked taxonomic groups according to the number of enzymes they dominate. Top five groups were termed ’key’ taxonomic groups.

For each environment, network expansions were carried six times using the corresponding set of source metabolites. In the first of the six iterations, the reaction set included the full set of metabolic functions (as described above). In each of the subsequent five iterations – all edges (metabolic functions) specifically dominated by one of the key taxonomic groups were removed from the original enzyme set. The impact of the removal of each key taxonomic group was estimated according to differences in the metabolite content (metabolite number) between the network expanded from the truncated enzyme set, and the original meta-network (first iteration, expanded from the full enzymatic set). The removed metabolite vectors, created for each iteration, were mapped to KEGG pathways. A removal effect score was calculated for each pathway as the fraction of metabolites left after the removal out of the original number of metabolites per pathway (counted in the first iteration, considering the un-truncated network).

In order to discover whether there is a synergistic metabolic complementation between combinations of key taxonomic groups, we applied a reverse approach to the removal procedure described above. Starting from a core set of enzymes representing common functions (i.e., functions that are not dominated by a key taxonomic group) we added combinations of specific and unique taxa-dominated enzyme sets. The possible combinations are described in Supplementary Table S8. An enzyme set was automatically created for each possible combination in each environment, according to the Simpson dominance scores. Thus leading to, for each combination tested, a number of enzyme sets as the number of combination members in addition to the enzyme set describing the amalgamate. Each enzyme set was used to expand a network as previously described. Each combination-type was tested in all four environments. All the networks were scanned for complementary metabolites. To identify complementary metabolites, we compared the metabolite content of multi-members networks to these of its corresponding single-member networks. That is, for a combination of a pair of taxonomic groups A and B, we expanded three networks per environment, giving a total of three networks i.e., (1) a core network + enzymes dominated by group A; (2) a core network + enzymes dominated by group B and (3) a core network + enzymes dominated by groups A and B (Supplementary Figure S4). Metabolites that were produced in the joint network (network 3), but not in the individual networks (networks 1 and 2), were termed complementary metabolites. For each such combination, the process was carried in each of the four environments. All complementary metabolites were than mapped to KEGG pathways (Supplementary Table S9).

Enzymes and metabolites were cataloged and mapped to pathways according to the KEGG database (Kanehisa et al., 2014). Testing for significantly enriched pathways (metabolites/enzymes) was done using the results of the hypergeometric enrichment test as in Yadav et al. (2016). In addition, using the R chisq.test function (R version 3.2.2), the X² goodness of fit test was used to evaluate the compliance of subsets with the relative distribution of samples. A pathway was considered significantly enriched if it passed both the hypergeometric distribution (p < 0.05) and the X² goodness of fit (p < 0.05) tests.

Venn diagrams were made using Venny (Oliveros, 2007) then adapted using Microsoft VISIO 2013. PCA analysis was performed using R prcomp function (R version 3.2.2). Heatmap and PCA plots were made using the ggplot2 R package (version 1.0.1) (Ito and Murphy, 2013). Pathways were plotted into a heatmap using the pheatmap R package (version 1.0.8) (Kolde, 2015). All network visualizations were made using Cytoscape (version 3.3.0) (Shannon et al., 2003).

Gene catalogs used here were constructed based on published DNA metagenomic data that were collected from four environments – roots of wheat and cucumber and the corresponding soils (Ofek-Lalzar et al., 2014). A total of 3436 KEGG orthologs, identified across all data, were mapped to 1574 unique enzymes (denoted by a four digits EC numbers) which represent the overall cross-environment compilation. We first identified differentially abundant enzymes in soil vs. root environments, considering each crop independently. Both plant treatments showed an overall similarity in their root vs. soil divergence pattern with most of the differentially abundant enzymes shared by both crops (Figure 1A). The differentially abundant enzyme sets were mapped to 100 KEGG pathways (Supplementary Table S10), with only nine showing significant enrichment in a specific environment (Figure 1B). Enrichment pattern of functional categories corresponds with previous reports: while some of the enriched root associated pathways were found to be involved in lipopolysaccharide metabolism, in agreement with Owen et al. (2007) and Ofek-Lalzar et al. (2014), most of enriched soil associated pathways were mostly of primary metabolism such as the TCA cycle and carbon metabolism. Next step was going beyond the list of discrete genes and integrating data into networks, according to stages framework outlined above. The sets of environment specific enzymes were used for the construction of four corresponding environment-specific networks. Subsequently, computationally based approximations of each metabolic environment were calculated based on the network topology. Similarly to the differentially abundant enzyme groups, most of the predicted source-metabolites in the root and soil environments are shared by both crops (Figure 1A). The source-metabolites were mapped to 90 KEGG pathways (Supplementary Table S11). When comparing the pathway distribution of source-metabolites across the different environments, seven pathways showing significant environment-specific associations were detected (Figure 1B). Most of the significant pathways were identified for the soil environments, with the exception of pathways from the biosynthesis of secondary metabolites category, that were significantly enriched in the cucumber root environment (Figure 1B), possibly reflecting the effect of root exudates.

The description of the predicted source-metabolites, together with the metabolic potential (the enzymes), allowed simulating metabolic activity in each of the four environments. As can be expected in natural robust systems (Freilich et al., 2010), the large majority of basic metabolism was carried despite environmental variations (Figure 1A bottom, and Supplementary Table S12). To delineate environmental induced metabolic activity, network metabolites were mapped to 122 KEGG pathways. Pathway mapping of the expanded networks showed differences in a wide range of secondary metabolism functions (Figure 1B). This stood in contrast to enriched pathways found in the initial enzymes and source-metabolites sets, where differences were found mainly in primary metabolic functions. Most of the divergent pathways found for the environment-specific networks were associated with root environments (Figure 1B). Some of these root-enriched pathways, belong to secondary metabolism categories, including metabolism of terpenoids, polyketides, and anthocyanins. These, are common plant metabolites that are less likely to be abundant with increasing distance from the root (Owen et al., 2007; Megharaj et al., 2011; Jeon and Madsen, 2013; Jha et al., 2015). These root unique network functions support the ecological relevance of the expanded environment-specific networks and their relevance for delineating robust versus unique metabolic capacities.

Considering the triangular relationship between an organism, a function and an environment, we set out to project taxonomy information over the environment-specific networks. In each environment, enzymes were scored according to the taxonomic diversity of their associated reads. A total of 667 (out of 1574) enzymes were paired with dominant taxa (order-level) in at least a single environment. Most differences in the classification profiles of these sets were associated with relative representation in secondary metabolism categories (Supplementary Figure S5). In each environment, the five taxonomic groups with the highest number of dominated enzymes were defined as key taxonomic groups (Supplementary Table S13 and Figure S6). Overall, these key taxonomic groups were similar across the four environments and included Actinomycetales (in soil environments only), Burkholderiales, Pseudomonadales, Rhizobiales, Sphingomonadales, and Xanthomonadales. Conserved vs. dominated functions in the communal metabolic network of the cucumber root are illustrated in Figure 2A. To directly explore the contribution of taxa-dominated functions to community performances, we simulated network activity while eliminating such enzymes, all at once and group by group (Methods, Supplementary Table S14). This iterative removal-expansion process directly explored the impact of the key taxonomic groups on metabolic processes carried in their specific environment. In general, the effect of function removal was found to depend both on its hierarchical positioning in a pathway (for example – an enzyme converting a source-metabolite into a compound accessible for multiple groups in the community will have a high impact), and on the robustness of the pathway (the prospects of finding alternative routes for the production of the corresponding metabolites). The impact is only relevant to network studies, representing a snapshot of community structure at a given time point. Despite the robustness of the expanded metabolic networks, the removal of all key taxonomic group specific functions led to up to approximately 26% reduction in network size (Methods, Figure 2B). The highest removal impact was observed in the cucumber root environment. In general, the Rhizobiales and the Actinomycetales groups have had the highest removal impact in root and the soil environments, respectively (Figure 2B). These taxonomic groups are generally known to be abundant and highly dominant in soil and root environments (Mendes et al., 2013; Tian and Gao, 2014).

The metabolic processes affected by the removal of key taxonomic group (order level) were delineated by mapping the omitted metabolites into KEGG pathways (Methods, Supplementary Tables S15–S18). As expected in ecological systems, the removal effect was observed to be both group specific, environment dependent and differed between root and soil (Figure 3A). Pathway mapping has shown that the large majority of pathways affected were belongs to secondary-metabolism categories (Figure 3B). In the soil, the removal of the order Actinomycetales invoked the highest predicted impact. One example for an affected function is streptomycin biosynthesis (Figure 3B), an antibiotic produced by bacteria from the order Actinomycetales (Nett et al., 2009; Zhou et al., 2012; Mitter et al., 2013; Panov et al., 2013). Actinomycetales’ significant impact was also simulated in various pathways associated with degradation of fluorobenzoate and compounds from the polychlorinated biphenyls (PCBs) pollutant family. This may stem from one of Actinomycetales high dominance enzyme, benzoate 1,2-dioxygenase (EC 1.14.12.10), which catalyzes a variety of ’first-step’ reactions toward the degradation of an array of benzoate analogs (Solyanikova et al., 2015). In the root environment, pathways that were affected by the removal included lipids, terepnoids, and plant induced secondary metabolites categories (Figure 3B). The taxonomic group with the highest impact in the root, in accordance with its key role in the rhizosphere (Berg and Smalla, 2009), was the Rhizobiales (Figure 3B). Out of 14 reactions involved in arachidonic acid metabolism, only a single reaction (aryl-4-monooxygenase) was simulated to be highly dominant by Rhizobiales sequences. However, its upstream location in the pathway suggests that Rhizobiales may be crucial for its metabolism in the surveyed environment. In addition, the simulated removal of Rhizobiales in the root (but not in the soil) affected the metabolism of linoleic acid and geraniol associated pathways. Both compounds are plant exudates that are used as carbon sources in the rhizosphere (Folman et al., 2001; Owen et al., 2007). Similarly, the simulated removal of Sphingomonadales in the root (but not in soil) affected mostly phenylpropanoid and flavonoid-related pathways (Figure 3B). These root-specific effects correspond with the role of plant exudates such as flavonoids, organic acids, and carbohydrates as determinants of the microorganism community structure in the rhizosphere (Narasimhan et al., 2003; Schulz and Dickschat, 2007; Ofek-Lalzar et al., 2014). Bisphenol degradation in the root, was affected by the removal of Actinomycetales, Pseudomonadales and Burkholderiale, but not by Rhizobiales, in correspondence with recent reports (Matsumura et al., 2015). Other pathway categories uniquely affected included these involved in the metabolism of potential regulators of plant–microbe interactions. For example, vitamin B6 was uniquely affected by the removal of the Pseudomonadales taxonomic group (Figure 3B), in accordance with their role in the production of B- group vitamins in the rhizosphere (Marek-Kozaczuk and Skorupska, 2001).

Bacterial communities, or combinations of taxonomic groups, are suggested to perform tasks that no species could perform on their individually (Freilich et al., 2011; Zomorrodi and Segre, 2016). Here, a simulative system was applied for surveying such potential synergistic interactions between the key taxonomic groups. A synergistic interaction was estimated according to complementary metabolites, defined as metabolites that are produced only in the presence of a combination of species, and not by individual members of the combination. These, environment-specific, combination-specific complementary metabolites were mapped to pathways (Figure 4, methods). Many of these predicted processes aligned well with ecological theories, hence providing a functional rational to observed co-occurrence patterns.

The highest number of complementary metabolites was simulated for a combination of Pseuodomonadales and Sphingomonadales in the wheat root environment (20 metabolites, Supplementary Table S19). Bacteria from these two taxonomic groups were demonstrated as having a co-dependent distribution pattern across environments (Pascual-Garcia et al., 2014). Most of the complementary metabolites predicted for the Pseuodomonadales–Sphingomonadales combination were mapped to the phenylpropanoid biosynthesis pathway (Figure 4). This corresponds with the demonstrated activity of Pseudomonadales in the rhizosphere that via the phenylpropanoid pathway assist the plants’ response to biotic stresses by contributing to a consortium that elicits the accumulation of phenolic compounds (Jain et al., 2012; Singh et al., 2013). Complementary metabolites predicted to be produced by a Actinomycetales–Burkholderiales combination included those in the pathway of inositol-phosphate metabolism (Figure 4). This corresponds with the suggested role of rhizobacteria in increasing phosphorus availability in the rhizosphere, hence contributing to plant growth (Unno et al., 2005; Singh and Satyanarayana, 2011; Wang et al., 2013). Burkholderiales and Xanthomonadales showed a simulated synergistic complementary effect related to tryptophan metabolism in the cucumber root environment. Tryptophan, secreted by the root, is converted by rhizobacteria to auxin, an hormone promoting plant growth (Kamilova et al., 2006). In the cucumber rhizosphere, amounts of secreted tryptophan were reported to be low in comparison to other crops (Kamilova et al., 2006). Hence, this cucumber-root specific synergistic complementary predicted effect may indicate a unique adaptation to a low tryptophan environment. Finally, complementation between Actinomycetales and Sphingomonadales was simulated in several pathways associated with PCB degradation and especially in dioxin degradation. As mentioned above, Actinomycetales dominate the ‘first step’ catalyzing enzyme. Bacteria from Sphingomonadales order dominate downstream enzymes. Thus the presence of both taxonomic groups, in turn, may lead to the use of PCBs as a carbon source.