A framework for interpreting and characterizing the network of metabolic interactions of root associated bacteria.

Going from obscure rhizosphere (upper circle) to the elucidated rhizosphere (lower circle): Microbial genomic DNA extracted from rhizosphere soil serves the construction of a rhizosphere-derived metagenome and the recovery of rhizosphere community Metagenome Assembled Genomes (MAGs). For each MAG, a Genome Scale Metabolic Model (GSMM) is built, representing a specific species in the microbial community. Apple-root exudate profiles based on metabolomics data are used for the construction of a simulation environment. Via Constraint Based Modelling (CBM), applied on GSMMs, interactions are characterized in the simulated root environment, yielding predictive information regarding potential trophic exchanges within the native rhizosphere microbial community.

Apple root exudates adopted from the works of Leisso et al. 48,49, included in the rhizosphere environment.

Phylogenomic surveys of the genomic collection assembled from sequence data derived from apple rhizosphere samples.

A. A phylogenomic tree constructed from 395 de-replicated MAGs extracted from the metagenomes. The tree is based on concatenated marker proteins according to GTDB-Tk44. B. Relative abundance of bacterial taxa at the phylum level as inferred from different phylogenomic classification approaches applied for the same samples: Binned-SG (MAGs derived from shotgun metagenomics); Non-binned SG (contigs derived from shotgun metagenomics); 16S rRNA amplicon derived from the same samples. Taxonomic classifications of rRNA amplicon data and non-binned-SG sequences are taken from Somera et al.41 and Berihu et al.42, respectively. Shotgun metagenomics data from Berihu et. al was used here for construction of MAGs catalogue. “Treatments 1 – 6” relate to six specific combinations of rootstock and soil amendment as detailed in Supp. Table 1. Relative abundance was calculated as the average value in five replicates conducted for each treatment.

Microbial Community Succession Module (MCSM), and characterization of trophic dynamics in the community along iterations.

A. An illustration of the iterative microbial community growth module, representing the growth of community members along iterations starting in the “Rhizosphere environment” and updated with microbial secretion outputs. B. Growth rate characterization of GSMMs in the community along iterations. Each column in the heat-map represents a different GSMM (only models which have grown in the rhizosphere environment are presented); growth rate is indicated by the color bar. Iterations are represented by rows. Blank spaces indicate models not growing at that specific iteration. N is the total number of species that grew after each iteration. C. Distribution of GSMMs growing along iterations at the phylum level. D. Distribution of organic compounds secreted along iterations, classified into biochemical groups. Root exudates bar (far left) represents the classification of organic compounds in the initial “Rhizosphere environment”. Numbers on top of bars in both C, and D (designated by N) denotes the number of new entries in a specific iteration, with respect to the previous iteration.

Trophic interactions based on community exchange fluxes predicted along iterations in the simulated rhizosphere environment.

A. Network representation of potential metabolite exchanges between rhizosphere community members. Edges in the network are directional; arrowhead from a grey node (metabolite) pointing towards a colorful node (GSMM) indicates uptake; arrowhead from a colorful node (GSMM) pointing towards a grey node (metabolite) indicates secretion. Node colors correspond to differential abundance classification of GSMMs in the different plots; H, S, NA are Healthy, Sick, Not-Associated, respectively. Metabolites found in the center of the network are of a higher connectivity degree (i.e., are involved in more exchanges). Only organic compounds are included in the network. Grey rectangles illustrate a zoom-in to specific 5 and 3 partite sub-networks (C, D). B. Pie chart distribution of GSMMs classified according to differential abundance of reads mapped to the respective MAGs; N is the total number of GSMMs in the network. C, D. Examples of specific sub-network motifs derived from the community network; PM (a three component sub-network motif, upper); PMM (a five component sub-network motif, lower), respectively. The full list of PM and PMM sub-network motifs is found in Supplementary data 9.

Characterization of PM and PMM sub-network motifs’ features associated with differentially abundant (DA) GSMMs.

A. Count distribution of PMM sub-networks determined by GSMMs in the first microbial position classified as either healthy (H), sick (S), or non-associated (NA). Boxes extend from first quantile to the third quantile, middle line represents the median; dots outside whiskers indicate outliers. Distinction of groups was determined using ANOVA followed by the Tukey post-hoc test. Asterisks indicate significance of test (*** <= 0.005 ** <= 0.01 * <= 0.05) B. Bar plot indicating the number of exudates significantly associated with H or S-classified PM sub-networks (Hypergeometric test; FDR <= 0.05; green: healthy-H, red: sick-S). C. Bar plots indicate the number of secreted compounds in PM sub-networks, which are significantly associated with H-classified (upper, colored green), or S-classified (lower, colored red) (Hypergeometric test; FDR <= 0.05). Compound name identifiers were taken from the BIGG database (X-axes). Venn diagram on the right represents the intersection of secreted compounds derived from both S and H classified PM sub-networks.

Metabolic pathway completion analysis of community de-replicated MAGs.

The presented cluster map links the MAGs with different metabolic pathways based on completion values of a specific metabolic pathway. Values are concluded based on the ratio between the count of EC numbers present in a MAG, which correspond to a specific pathway (derived from annotated genes), and the total EC numbers in that specific KEGG pathway. Left axis represent the KEGG KO of a set of specific functions, which are indicated on the right axis. Top and bottom color bars on the upper axis represent the taxonomic annotation of the MAGs (phylum and class, respectively), which are indicated on the bottom axis.

Distribution of GSMM attributes at the phylum level.

Subplots of reactions, metabolites and exchanges presenting the distribution of the respective attributes for each group of bacteria, affiliated with a specific phylum.

Bacterial growth rates in different environments.

Each row represents the growth rate (biomass increase hour-1) of species in a specialized medium; from top to bottom: poor medium, poor medium supplemented with root exudates, root environment medium and rich medium. Each column represents a GSMM. Models are sorted by their growth score in the root environment medium. Poor, Rich, and Poor + exudates media are model specific, thus sustain growth of all models. The horizontal color bar on top of the plot represents the phyla of the corresponding GSMM. Actual growth values are provided in Supp. Data 5.

FVA performances of GSMMs in different environments (Supp.

Fig. 3, Supp. Data 5). A. Distribution of potentially active exchange reactions (non-zero minimum FVA flux) in the different environments. Solid line inside each violin indicates the interquartile range (IQR). White point in IQR indicates the median value. Whiskers extending from the IQR indicate the range within 1.5 times the IQR from the quartiles. Violin width at a given value represents the density of data points at that value. B. Loopless FVA scores compared to regular FVA for models in the 3 different environments. Bars indicate the count of active fluxes (non-zero minimum FVA flux). Only a subset of models was used for this analysis.

Phylogenetic distribution of models growing along iterations at the order level.

Bars represent the taxonomic distribution of the community (i.e., overall models that grew along iteration, with respect to the corresponding iteration). Colors are indicative of different orders. Orders with insufficient nomenclature were classified as “Others”.

Uptake and secretion degrees of biochemically-classified metabolite nodes, as inferred from the communal interaction network.

Values indicate mean amount of links between community GSMMs and nodes classified to the corresponding category. Biochemical classification was conducted both based on BRITE database and manually.

Application of MCSM over the process of cellulose decomposition as described by Kato et al1. 5-partite network exhibiting the uptake of cellulose oligomers (4 and 6 units of connected D-glucose) by primary decomposers, through secretion of intermediate compounds and their metabolization by secondary decomposers to CO2. Distribution of phyla of primary and secondary decomposers is denoted by pie charts. Though MAGs were notconstructed for the original species as in Kato et al., among the primary consumers, species corresponding to the Acidobacteria (Acidobacteriales)2, Actinobacteria3, Bacteriodetes4, Proteobacteria (Xanthomonadales)5 and Verrucobacteria6 groups are found to be capable of degrading cellulose compounds via enzymatic mechanisms.

Community GSMM-centric characterization of plant-microbe PM) and plant-microbe-microbe (PMM) motifs. A, B, C. Frequency of community GSMMs in) PM, PMM (first/second position) in subnetwork motifs, respectively. D. correlation of GSMMs position frequency in PMM subnetworks; Pearson, P value=0.009.

Summary table of the mapping of reads to MAGs across treatments. Treatment types are fully descriebd in Berihu et al7.

Table of compounds used throughout iterations in the MCSM

Initial medium compounds used in MCSM Cellulose degradation process