Metabolic interdependencies in thermophilic communities are revealed using co-occurrence and complementarity networks

Domain: earth_science
Published: 2024-09-17
Editorial outcome: major_revision

Abstract
Microbial communities exhibit intricate interactions underpinned by metabolic dependencies. To elucidate these dependencies, we present a workflow utilizing random matrix theory on metagenome-assembled genomes to construct co-occurrence and metabolic complementarity networks. We apply this approach to a temperature gradient hot spring, unraveling the interplay between thermal stress and metabolic cooperation. Our analysis reveals an increase in the frequency of metabolic interactions with rising temperatures. Amino acids, coenzyme A derivatives, and carbohydrates emerge as key exchange metabolites, forming the foundation for syntrophic dependencies, in which commensalistic interactions take a greater proportion than mutualistic ones. These metabolic exchanges are most prevalent between phylogenetically distant species, especially archaea-bacteria collaborations, as a crucial adaptation to harsh environments. Furthermore, we identify a significant positive correlation between basal metabolite exchange and genome size disparity, potentially signifying a means for streamlined genomes to leverage cooperation with metabolically richer partners. This phenomenon is also confirmed by another composting system which has a similar wide range of temperature fluctuations. Our workflow provides a feasible way to decipher the metabolic complementarity mechanisms underlying microbial interactions, and our findings suggest environmental stress regulates the cooperative strategies of thermophiles, while these dependencies have been potentially hardwired into their genomes during co-evolutions.

Introduction
Microbial ecosystems thrive on intricate webs of interactions, often obscured by their vast diversity and hidden by limitations in conventional culture-based methods1,2. While studying pairwise interactions provides valuable insights3,4, understanding the complete picture within natural complex communities remains challenging5. With the advances in sequencing and computational modeling, microbial co-occurrence networks, constructed from patterns of species co-occurrence across environmental samples, have become a powerful tool to infer interspecies connections6–8. The use of co-occurrence network analysis became popular after several ready-to-use workflows with user-friendly visualizations were published9–14. However, as several previous studies argue, co-occurrence alone does not provide conclusive evidence of ecological interactions15–17. To address this limitation, we combine co-occurrence data with evidence of metabolic complementarity to infer more robust interactions and explore potential underlying mechanisms.

Metabolic dependencies have been proposed as a major driver of species co-occurrence18. However, it has also been suggested that species co-occurrence could in turn drive metabolic dependencies2. This highlights the need for a nuanced understanding of the dynamic interplay between metabolic interactions and microbial community assembly. Long-read sequencing and improved assembly enable high-quality metagenome-assembled genomes (MAGs), opening doors to map metabolic cooperation and competition in complex communities through nutrient exchange19,20. Although genome-scale metabolic models (GSMMs), along with flux balance analysis (FBA), have shown potential in predicting metabolite exchange and growth under curated laboratory conditions21,22, there remain substantial potential for further exploration in this area. Integrating these models with co-occurrence networks to validate and refine predictions of metabolic dependencies within uncontrolled natural habitats is an untapped frontier.

Geothermal springs, characterized by extreme temperatures, fluctuating ionic strengths, and enriched minerals, harbor a unique community of thermophiles with specialized physiological and metabolic adaptations23–25. These adaptations often involve genome streamlining, leading to reduced genomes and temperature-adapted proteins26–31. Yet, thermophiles exhibit remarkable genomic plasticity and metabolic flexibility32. Metabolic reconstructions of some geothermal phylas suggest dependence on interspecies exchange of amino acids, vitamins, and cofactors, implying close syntrophy33,34. This is further supported by observations of highly cooperative communities formed by small-genome auxotrophs, where cross-feeding promotes mutual benefit35–37. This suggests that harsh geothermal conditions while restricting growth, may simultaneously drive synergistic interactions. Nevertheless, this hypothesis requires a systemic examination. Addressing this question may yield insights into the ongoing debate regarding the prevalence of either antagonistic or synergistic interactions within natural microbial communities, as well as how environmental stresses affect their prevalence38–42.

Here, we leverage the power of random matrix theory (RMT) for robust co-occurrence and metabolic network constructions while tapping into the rich information contained within MAGs9,43. This synergy allows us to identify species co-occurrence patterns and predict potential metabolic interdependencies based on metabolic pathways. Adopting this powerful tool, we show that metabolic complementarity intensifies with rising temperatures, as thermophiles under heat stress increasingly rely on interspecies exchange of essential metabolites for survival. Additionally, we demonstrate that phylogenetic distance dictates cooperative strategies, with closely related thermophiles competing for similar resources while distantly related ones engage in mutually beneficial metabolic exchanges. Furthermore, our findings reveal that species with small genomes, potentially lacking vital metabolic pathways, depend heavily on metabolic partnerships for survival and growth. Our research offers a deeper understanding of microbial synergy in challenging environments, providing valuable insights into the interplay between environmental stress, metabolic dependencies, and the evolution of cooperative strategies within microbial communities.

Results
Network construction and thermal preferences of thermophiles
To comprehensively analyze both co-occurrence and metabolic complementarity networks among thermophiles, we developed a unique bioinformatic workflow (Fig. 1) and seamlessly integrated it into the iNAP platform (https://inap.denglab.org.cn)44 (Supplementary Fig. 1), as updated iNAP 2.045. This workflow leverages multiple metabolic complementarity indices, infers the network thresholds of co-occurrence and metabolic complementarity using the random matrix theory (RMT), extracts shared interactions within both network inferences, and identifies potential exchangeable metabolites.

A total of 449 million Illumina short reads (223.35 Gb of raw bases, N=40) and 51 million Nanopore long reads (52.01 Gb of raw bases, N=40) were obtained from the hot spring sediment samples across a temperature gradient (63.5°C~85.8°C), with eight samples each having five replicates. The Nonpareil method confirmed the estimated average coverage of each sample above 80%, indicating sufficient sequencing depth (Supplementary Fig. 2). There is an increasing trend of sequencing diversity (Nd) from the high- to low-temperatures (Fig. 2a), indicating that the environmental conditions of extreme thermal temperatures have reduced the sequence diversity of its microbial communities, and the inferred species diversity might be more specialized. After assembly, binning, and CheckM quality control, 401 medium- and high-quality metagenome-assembled genomes (MAGs) were retained for subsequent network analysis.

The taxonomic assignment revealed that 85.78% of MAGs belonged to 38 bacterial phyla, whereas 57 MAGs were affiliated with seven archaeal phyla (Fig. 2b, Supplementary Data 1). The MAG abundance varied significantly across temperature ranges (PERMANOVA, p<0.001) (Fig. 2c), leading to the classification of sampling sites into three groups: extremely thermal (ET, 78.5–85.8°C, 2 samples×5 replicates), highly thermal (HT, 67.5–73.9°C, 3 samples×5 replicates), and moderately thermal (MT, 63.5–65.8°C, 3 samples×5 replicates). Most MAGs exhibited distinct temperature preferences, with only a few evenly distributed across the gradient. For example, 80.52% of the relative abundance assigned to Thermoproteota MAGs originated from the ET group, highlighting this archaea’s strong preference for extremely high-temperature habitats. In contrast, Cyanobacteria and Bacteroidota preferred the cooler group, with 96.28% and 96.98% of relative abundance originating from HT and MT groups (Fig. 2d, Supplementary Data 2).

Archaea also displayed clear temperature preferences despite accounting for only 7.31% of relative abundance across all samples. While Thermoproteota favored the ET group, Micrarchaeota, the second most abundant, exhibited a unique distribution, with nearly equal contributions from HT and MT groups. These findings underscore the diverse thermal adaptations and preferences of thermophile communities within this hot spring ecosystem.

The co-occurrence patterns of thermophiles within three temperature groups
Three co-occurrence networks were constructed using the relative abundance information, employing the RMT-based method (Fig. 3, Supplementary Data 3). Positive edges constituted the majority in the three networks (ET: 98.82%; HT: 75.54%; MT: 77.99%). Post-extraction of positive edges as subnetworks, the R2 values of the power-law model for the co-occurrence networks of the HT and MT groups were 0.59 and 0.53, respectively, but the ET group could not fit well to the power-law distribution, losing its scale-free property. We found that the ET group subnetwork exhibited a higher network density, shorter average path distance (GD), lower harmonic geodesic distance (HD), higher average clustering coefficient, and reduced modularity compared to the other two subnetworks (Table 1). These characteristics suggested a tighter interaction structure in the ET group subnetwork than in the other two subnetworks. A harsh environment (high temperature) induces tight and synergistic interactions among thermophiles. Furthermore, the global topological properties of the three networks were significantly different from the respective randomized networks (generated 100 times), demonstrating that all observed interactions filtered by RMT cutoff were not randomly connected (Table 1, one sample t-test).

The inferred pairwise metabolic complementarity among thermophiles
To assess metabolic dependencies, we reconstructed genome-scale metabolic models for each MAG. We defined the PhyloMint metabolic complementarity index (MIcomplementarity) to quantitatively assess the degree of metabolic dependencies between each pair of MAGs within three temperature groups (see “Methods” section). Specifically, 29 MAGs were unique to the ET group, while 70 MAGs were detected across all temperatures, likely representing heat-tolerant generalists. Both approaches revealed a surprising rarity of synergistic metabolic interactions (Fig. 4), while the RMT-based threshold yielded an even more stringent classification, with less than 3% of interactions deemed significant (ET: 4.13%; HT: 4.68%; MT: 6.46%) (Fig. 4).

Notably, observed synergistic interactions displayed a marked asymmetry of MIcomplementarity values under all three groups. We categorized these into mutualistic (both MIcomplementarity values of pairwise MAGs exceeding the threshold) and commensalistic (only one MIcomplementarity value of pairwise MAGs exceeding the threshold) metabolic interactions. Regardless of the threshold used, mutualistic interactions were extremely rare, constituting less than 0.3% of all pairs in any group (Fig. 4). This suggests that metabolism-based synergy within the hot spring community is primarily driven by unidirectional feeding, with one thermophile benefiting from the metabolic products of another. In particular, the species pairs involved in commensalistic and mutualistic interactions exhibit dramatically different patterns regarding genome size differences. Across the three temperature groups, the average genome size differences (estimated genome size of the giver minus that of the taker) identified as commensalistic interactions by the RMT threshold were ET: 716 Kbp, HT: 1204 Kbp, and MT: 879 Kbp. However, the average genome size differences for mutualistic interactions were close to 0 for all three temperature groups. In all commensalistic interactions, the proportion of interactions where the estimated genome size of the taker is less than 2 Mbp and that of the giver is more than 2 Mbp accounted for as high as 31.29% (87 pairs) in the ET group, while in the HT and MT groups, these values were 37.01% (624 pairs) and 29.87% (1057 pairs), respectively. These results indicated that most metabolic dependency between one streamlined genome and another more comprehensive genome is a more prevalent pattern in commensalistic interactions.

Following RMT threshold determination (Fig. 4), the constructed metabolic networks for the ET group posed more robust scale-free properties than that of the occurrence networks (Table 2). In addition, the topological properties of metabolic networks varied across temperatures. The hottest ET network, despite its smaller size, exhibited the highest average clustering coefficient, highest density, and shortest communication paths (Table 2), suggesting a highly interconnected community in this harsh environment. However, the patterns in average degree, average path distance, and harmonic geodesic distance were opposite to those of the co-occurrence networks. Such findings illustrated that metabolic complementarity and co-occurrence networks do not exhibit the same characteristics of microbial communities, which were then reflected in the topology of the network. The hub node identification using z–P analysis illuminated the key players in the metabolic networks (Supplementary Data 4). The ET network, with 36 hubs, hosted 6 archaeal hubs (3 connector hubs, 1 module hub, and 2 network hubs), including 4 assigned to Thermoproteota, 1 to Asgardarchaeota, and 1 to Aenigmarchaeota. The larger HT and MT networks revealed diverse hubs (137 and 97), while both contained hubs from Thermoproteota, Micrarchaeota, and Methanobacteriota. These results indicated that despite their low abundance, archaea, particularly Thermoproteota and Methanobacteriota, emerge as essential hubs across all temperature networks, demonstrating their unique metabolic capabilities and indispensable roles in the thermophile communities.

To understand how microbial interactions arise, we investigated the relationships between phylogenetic distance (Faith’s PD) and the co-occurrence patterns/ metabolic complementarity. The co-occurrence strength, measured by Spearman’s correlation, negatively correlated to PD in all temperature groups (Fig. 5a–c). However, a fascinating pattern emerged after RMT filtering (cutoff = 0.830). The density distribution of PD exhibited two distinct peaks: one near 0 and another near 2. This suggested two potential drivers of co-occurrence: the close association of similar species with shared resource requirements (PD at peak 0) and the co-occurrence of phylogenetically distant species with complementary metabolic needs (PD at peak 2). Besides, significant positive correlations between MIcomplementarity and PD were observed in all temperature groups (Fig. 5d–f), indicating that distantly related MAGs were more likely to exhibit strong metabolic dependencies. This trend was particularly pronounced in the hottest temperature. Notably, after RMT filtering, those MAG pairs of PD lower than 1 became rare, with most values clustering around 2. It suggested that the analysis of metabolic complementarity omitted co-occurrences based on phylogenetic proximity, highlighting the importance of metabolic complementarity between distantly related species.

Genomic clues to species metabolic synergies
Genome sizes and contents were shown to be crucial connecting links for metabolic interactions within the hot spring community, which also were intricately linked to nutrient requirements. The estimated genome sizes of MAGs distributed at different temperatures underwent a transition from negative to positive correlation with temperature from the ET group (T85 and T78) to lower temperature (Fig. 6a). In harsh environments with extremely high temperatures, small-genome species tend to occupy higher abundances. However, when the environmental temperature drops to a range where most species can adapt, the small-genome species no longer hold this abundance advantage and tend to occupy only a small abundance in the environment. Genome size was also confirmed to correlate with the proportion of genes in the genome that perform different functions (Fig. 6b). Estimated genome size showed a significant positive correlation with the proportion of genes involved in Secondary metabolites biosysthesis, transport and catabolism (COG-Q, Spearman’s Rho=0.670, p<0.001), Carbohydrate transport and metabolism (COG-G, Spearman’s Rho=0.473, p<0.001), Lipid transport and metabolism (COG-I, Spearman’s Rho=0.234, p<0.001) and Inorganic ion transport and metabolism (COG-P, Spearman’s Rho=0.125, p<0.05). Conversely, smaller genomes prioritized housekeeping functions like Translation, ribosomal structures and biogenesis (COG-J, Spearman’s Rho=−0.932, p<0.001), Nucleotide transport and metabolism (COG-F, Spearman’s Rho=−0.634, p<0.001) and Replication, recombination and repair (COG-L, Spearman’s Rho=−0.462, p<0.001). Furthermore, a linear mixed-effect model was constructed to reveal that when genome size was set to be random effect, differences in biological functions such as energy production and conversion, metabolism of nucleotide, amino acid and lipid had a strong positive effect on MIcomplementarity (Fig. 6c). The MIcomplementarity of pairwise MAGs significantly increased with the difference in genome sizes (Fig. 6d–f). This correlation was strongest in the ET group, suggesting that genome size plays more crucial role in shaping metabolic partnerships in harsh environments. These results reflected a trade-off of smaller size genomes, where the genes associated with genetic information storage and processing were retained to a greater degree while genes related to various metabolisms underwent loss, relying on synergistic interactions to compensate for their limited metabolic repertoire. However, it is essential to consider that the absence of specific genes might also be due to incomplete genome reconstruction, given the accepted cutoff for MAG completeness at 50%.

More genomic clues of species synergies were found in the overlapped partnerships between co-occurrence and metabolic complementarity networks. Only a handful partnerships (7, 49, and 58 in the ET, HT, and MT groups, respectively) exhibited strong metabolic complementarity alongside persistent co-occurrence (Supplementary Data 5). Among these partnerships, 274 metabolites were detected as potentially transferable, while 58 were coenzyme A derivatives. Amino acids (with peptides and analogs) and carbohydrates (and carbohydrate conjugates) followed, with 43 and 36, respectively. There were also 30 molecules classified as nucleosides, nucleotides, and analogs deemed to be potentially transferable. Notably, between those paired species, the types of metabolites potentially transferred from one to the other were quite different (Supplementary Data 5).

In exploring the intricacies of species metabolic synergies, our focus was drawn to 17 archaea-bacteria interactions, which were supported by both metabolic networks and co-occurrence networks, indicating that inter-domain synergies were prevalent in harsh environments such as hot springs. In these interactions, there were some MAGs that exhibited the properties of co-occurrence and metabolic complementarity with multiple other species. For instance, an archaeal MAG (T64.bin.92, Micrarchaeota phylum) had seven distinct bacterial partners (Fig. 7). This centrally positioned MAG, characterized by its smaller estimated genome size of 0.857 Mbp, predominantly assumed the role of a taker in these interactions. This was evidenced by its substantial reliance on obtaining several types of metabolites, such as some coenzyme A derivatives and other crucial metabolic need that it cannot fulfill endogenously, from other species. Despite its taker centric position, T64.bin.92 contributed to these commensalistic relationships by providing potential surplus materials like carbohydrates and nucleoside/nucleotide- related substances, embodying a reciprocal dynamic in these interactions. This necessity for external coenzyme A derivatives, vital for biosynthesis reactions involving acyl transfers, underline the significant dependence of T64.bin.92 and its partners. Furthermore, its requirement for certain inorganic substances, including various metal ion and non-metal compounds from its metabolic counterparts, further illustrated the complexity and significance of its role as a taker. However, it is important to note that these conclusions are based on genetic potential inferred from genomic data. To confirm the actual metabolic exchanges and interactions, additional validation with functional omics and metabolomics data is necessary.

The validation of our study on composting system experiment
In hot spring habitats, we noticed that symmetric interactions (mutualism) were rare. Besides, in commensalistic and mutualistic interactions, the genome size differences between species pair varied. In a similar vein, our findings from a composting system where temperature is a key stress factor mirror those from hot springs. A low occurrence of mutualistic interactions were observed across different sampling time (temperature) groups: D00 (Tavg=26.02°C), D05 (Tavg=64.29°C), and D25 (Tavg=41.35°C), with mutualism rates at 0.23%, 0.087%, and 0.23%, respectively. Asymmetrical interactions, suggesting commensalism, were slightly more common at 4.23%, 2.73%, and 4.18%. The genome size differences in commensalistic interactions at D05 were around 1.36 Mbp, similar to findings in hot springs, whereas mutualistic interactions across all groups showed negligible genome size differences (Supplementary Fig. 3).

The correlation between phylogenetic distance and pairwise MIcomplementarity was also observed to be significant in composting system. Specifically, D05, as the highest temperature sample in the experiment, showed the strongest positive correlation. Furthermore, the pattern that the greater the genome size difference between two species corresponds to higher their complementarity potential was confirmed by testing the significance of the correlation between MIcomplementarity and genome size difference (D00: Spearman’s Rho= 0.122, p<0.001; D05: Spearman’s Rho=0.213, p<0.001; D25: Spearman’s Rho=0.090, p<0.001).

The sequencing data from the composting experiment recovered few overly streamlined genomes, suggesting that such genomes might arise from long-term selection rather than short-term temperature changes. Nevertheless, we observed asymmetrical synergistic interactions, with a significant proportion being commensalistic (Supplementary Fig. 4). These interactions showed a considerable difference in genome size, particularly in the highest temperature group (genome size difference: D00, 815 Kbp; D05, 1368 Kbp; D25, 525 Kbp), similar to the hot spring habitats. This implies that high temperatures are a crucial factor driving synergy between genomes of vastly different sizes.

Methods
Sample collection, DNA extraction, and sequencing
Sediment samples for our study were collected from a hot spring in Tengchong, Yunnan Province, China (N24°56′ ~ 25°27′, E98°26′ ~ 98°27′) in June 2020. The temperature range measured at this hot spring was between 63.5°C and 85.8°C, and the water was slightly alkaline (pH range=8.36~8.70) (Supplementary Fig. 5a). Eight sampling sites displaying a gradual temperature decrease were selected along the spring flow, with five replicate samples collected at each site. According to the measured temperature of each sample, eight temperature sites were labeled as T85, T78, T73, T70, T67, T65, T64, and T63 (Supplementary Fig. 5b). After sample collection, they were promptly placed in a liquid nitrogen container and transported to laboratory within 2 days. Subsequently, they were stored at −80°C.

Before storage, sediment samples for DNA extraction were pre-divided to prevent DNA damage due to repeated freeze-thaw cycles in later experiments. After 48 h of lyophilization, total DNA within the microbial community was extracted from 1.5 g of freeze-dried sediment using the grind plus kit method as previously described64.

The acquired DNA was used for Illumina NovaSeq6000 PE250 metagenomic sequencing (DNA library insertion size: 450bp) and Oxford Nanopore sequencing (the PromethION R9.4 flow cells FLO-PRO002, Oxford Nanopore PromethION sequencer). The Illumina sequencing was conducted on all 40 samples by Magigene Biotechnology Co., Ltd. (Guangzhou, China). Replicate samples of each site were mixed and sequenced using Nanopore sequencing by Benagen Technology Co., Ltd (Wuhan, China).

Metagenome assembly, genome binning, and MAG classification
Raw Illumina metagenomic reads were quality trimmed using Trimmomatic (v0.39, LEADING:3 TRAILING:3 MINLEN:50 SLIDINGWINDOW:4:20)65. The remaining reads of each sample were de novo assembled using IDBA-UD (v1.1.3) with default parameters66. Here, the reads of 5 replicates from one sampling site were co-assembled to improve the robustness of assembly and the possibility of obtaining a higher diversity of genomes within the homogeneous environment. The sequencing coverage and read diversity were estimated using the Nonpareil method67,68. Then, contigs of samples from the same sampling site were pooled together and used for co-assembly with Nanopore long reads using OPERA-MS (v0.9.0) with default parameters69. Longer contigs generated by this hybrid assembly, with a minimum length of 1000 bp (1500 bp for metaBAT2), were used for genome binning using metaWRAP (v1.3.2)70, with metaBAT2 (v2.12.1)71, MaxBin2 (v2.2.6)72, and CONCOCT (v1.0.0)73 as the core binning tools. Draft bins were quality controlled using metaWRAP bin_refinement module with parameters -c 50 -x 10, indicating only MAGs with completeness higher than 50% and contamination lower than 10% were retained for the following analysis, declared to be medium- and high-quality MAGs74. The refined bins were then replicated using dRep (v3.5.0)75 with parameters -pa 0.9 -sa 0.99. The estimated genome size was calculated by dividing the genome size by the sum of completeness and contamination76. The taxonomic classification of the selected MAGs was conducted with GTDB-Tk classify_wf workflow against the GTDB genome database (Release 202)77. The phylogenetic trees for bacterial and archaeal MAGs were constructed using multiple sequence alignment results generated by GTDB-Tk workflow. The unrooted bacterial and archaeal phylogenetic trees were rooted using midpoint rooting method, performed by midpoint.root function in phytools R package78. and were then visualized using the online tree display tool iTOL (v6.6)79. The phylogenetic distance was then calculated using the cophenetic.phylo function in R.

Quantification of MAGs and co-occurrence network construction
The refined bins were quantified by using CoverM (v0.7.0, https://github.com/wwood/CoverM, genome mode and coupled reads as input) with parameter: --method relative_abundance. In the output relative abundance table (referred as MAG table), every column representing a sample contains unmapped reads percentage, making it reasonable to compare relative abundance across samples. A non-metric multidimensional scaling analysis (NMDS) was conducted with an encapsulated function metaMDS in the vegan package (R) to investigate the community-level difference among sampling sites. A dissimilarity test based on the Bray–Curtis distance was then applied to divide samples into three groups. MAGs were assigned to three groups (ET, HT, MT) based on their relative abundance in samples. Each group comprised multiple samples with five replicates. A MAG was considered present in a group if it had non-zero abundance in more than half of the samples (≥6 for ET; ≥8 for HT and MT). If a MAG did not meet this criterion in any group, it was assigned to the group with the highest total relative abundance. Additionally, MAGs were ensured to be in the group corresponding to their sample prefix (e.g., T85 and T78 in ET). With the abundance table, the pairwise Spearman’s rank correlation of each two MAGs was calculated. Since the relative abundance table originated from the distribution of hypothetical species in a natural environment, it can be regarded as the matrix that the random matrix theory (RMT)-based approach requires. Before using RMT cutoff tool on iNAP website44, majority_selection tool was used to filter MAGs that had zero abundance in less than half of samples. This step is to exclude biases from correlation calculations due to the effect of too many zeros. The separated noise formed a network representing MAG pairs with a high abundance correlation, inferring their co-occurrence. The network’s general properties were calculated using the igraph R package and visualized using Cytoscape (v3.10.0)80.

Genome-scale metabolic models and interaction network construction
Genome-scale metabolic models (GSMMs) for all MAGs were constructed using CarveMe (v1.4.1) with default parameters81. The input for CarveMe was coding sequences of each MAG, which were predicted using Prokka82. The interactions between pairwise GSMMs were predicted using PhyloMint47. All modeling software mentioned above required the installation of an optimization solver, which for our study was IBM CPLEX Optimizer (ILOG COS 20.10 Linux x86-64 version). PhyloMint metabolic complementarity index (MIcomplementarity) was used to represent the cooperation potential of pairwise GSMMs, or MAGs. For a given MAG pair, two MIcomplementarity values, ranging from 0 to 1, were calculated due to the asymmetry of the index that a MAG could be both the metabolic giver and taker. Therefore, a higher MIcomplementarity value was regarded as the maximum metabolic cooperation ability between pairwise MAGs and was used for the following analysis. MIcomplementarity values were determined using the RMT cutoff tool on iNAP website, and the corresponding pairwise interactions were selected to form the PhyloMint complementarity network. In the network, synergistic interactions are determined based on whether the original PhyloMInt values for all species pairs represented by the edges exceed the RMT screening threshold. If only the unidirectional PhyloMInt value exceeds the RMT threshold, the edge is defined as representing commensalistic interactions; if the bidirectional PhyloMInt values both exceed the RMT threshold, the edge is defined as representing mutualistic interactions. Similarly, the network attributes were also computed. The predicted coding sequences of each MAG were searched against the database of Clusters of Orthologous Genes (COGs, https://www.ncbi.nlm.nih.gov/research/cog/)83 and were categorized into various biological functions. The effects of biological functions on the complementarity potential between pairwise MAGs were estimated using linear mixed-effects models (LMMs) to eliminate the impact of genome size on functions. This step was performed using the R package lme4. Pairwise MIcomplementarity was regarded as the response variable, and gene proportion differences of each biological function were regarded as fixed effects. The genome size difference was termed a random intercept effect.

Extraction of shared interactions and determination of potential exchangeable metabolites
If a specific pairwise MAG interaction possessed both co-occurring properties, as detected by the co-occurrence network and the metabolic cooperation potential detected by the PhyloMint MIcomplementarity network, they were defined as a “dual interaction”. A dual interaction indicated that two MAGs shared a strong correlation between environmental co-occurrence patterns and metabolite profiling complementarity. Using the PhyloMint algorithm, seed metabolites were defined by the strongly connected component and represented substrates acquired exogenously. When considering two MAGs, A and B, a metabolite found in A’s seed set but not in B’s seed set indicated a potential transferability from B to A. This suggests that B could synthesize the metabolite while A could utilize it. This definition aligns with the computation method of the PhyloMint complementarity index. In line with the previously described definitions, the potential of one MAG to utilize the metabolites produced by another was determined by the overlap between one MAG’s seed metabolite set and the other MAG’s non-seed metabolite set. The metabolite profiles were taken directly from the corresponding genome-scale metabolic models generated by CarveMe with the BiGG Models database (Version 1.6)84. The metabolites were categorized using the Human Metabolome Database (HMDB Version 5.0)85. Due to the miscellaneous nomenclature of diverse metabolites, we manually browsed all related metabolites by their name in the BiGG Models database to recheck the classification in HMDB and made some changes to the metabolite classification rules of HMDB as follows (at “super class” level): (i) fatty acyl CoAs were separated from lipids and lipid-like molecules, and were then grouped with other CoA derivatives to form the class “Coenzyme A derivatives”; (ii) compounds identified as amino acids, peptides, and analogs were separated from organic acids and derivatives into a separate class; (iii) Carbohydrates and carbohydrate conjugates were separated from organic oxygen compounds into a separate class; (iv) all inorganic compounds were combined into one class, “Inorganics”. This class included homogeneous metal and non-metal compounds (inorganic compounds that contain solely metal or non-metal elements, respectively); (v) Organic nitrogen compounds, organic oxygen compounds, organoheterocyclic compounds, benzenoids, and organosulfur compounds were combined into one class “Organic compounds with specific atoms or structures”; (vi) Oxidized ferrodoxin, reduced ferrodoxin, protoheme, staphyloferrin B, and pyoverdine P. putida specific were categorized into “Other metabolites”. With these modified rules, all metabolites involved in our study were classified into nine categories (Supplementary Fig. 6, Supplementary Data 5).

Analysis of the composting system
Zhao et al. conducted a 30-day composting experiment in 10 composting piles at a food waste composting facility. As the temperature changed during the experiment, they selected six sampling time points, with three replicates from each composting pile for each sampling, obtaining a total of 180 samples for metagenomic sequencing. After co-assembly and binning, they obtained 159 high-quality MAGs (completeness >90%, contamination <5%) for further study. Three sampling time points were selected, Day 0, 5, and 25 (tagged as D00, D05, and D25), with average temperatures of 26.02°C, 64.29°C, and 41.35°C. These three time points were regarded as three states of temperature change in the composting experiment (Day 0: the beginning of the experiment, the coolest temperature; Day 5: the highest temperature; Day 25: the temperature dropped in the late stage of the experiment). The relative abundance table of MAGs at these three sampling time points was summarized as a proxy for the state of these MAGs at that time point. The metabolic models and metabolic complementarity indices of MAGs were constructed and calculated using the methods described above. In the three groups, MAGs with relative abundance greater than 0 were assigned to the corresponding group for the construction of metabolic complementarity networks. The threshold for identifying metabolic complementarity within networks was established at an RMT cutoff of 0.280.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Review of \"Deciphering Metabolic Interdependencies in Thermophilic Communities Using Co-Occurrence and Complementarity Networks\" by Xi Peng, Shang Wang, Miaoxiao Wang, Kai Feng, Qing He, Xingsheng Yang, Weiguo Hou, Fangru Li, Ye Deng.\n**Manuscript ID:** NCOMMS-24-14471\n---\nIn the manuscript titled \"Deciphering Metabolic Interdependencies in Thermophilic Communities Using Co-Occurrence and Complementarity Networks,\" Xi Peng and colleagues attempt to unravel the complex web of metabolic interdependencies in microbial communities residing in thermophilic environments. Employing a novel integrative approach combining metagenome-assembled genomes (MAGs) with random matrix theory (RMT) to construct both co-occurrence and metabolic complementarity networks, the study highlights how microbial interactions intensify under extreme temperature conditions. Key findings suggest that metabolic exchanges, particularly involving amino acids and coenzyme A derivatives, are crucial for the survival and adaptation of microbial communities in high-temperature ecosystems. This research contributes to our understanding of microbial ecology by linking environmental stressors with metabolic cooperation strategies, potentially informing future biotechnological applications in bioenergy production and environmental remediation.\nThe major issue in the manuscript revolves around the interpretation of microbial interactions, which predominantly suggests a framework of cooperation and synergy that may not accurately reflect the nature of the interactions observed, particularly the prevalence of commensalism. While the innovative methodologies applied significantly advance our understanding of microbial ecosystems in extreme environments, it is critical to refine the conceptual framework presented. The study's tendency to categorize commensal interactions, where benefits accrue to one party without reciprocal interactions, as cooperative or synergistic, is misleading. This distinction is crucial not only for accurate ecological characterization but also for understanding the evolutionary pressures that shape these interactions. Re-framing the discussion in the manuscript to accurately reflect these unidirectional dependencies would provide a clearer and more precise interpretation of the data, enhancing the manuscript's contribution to microbial ecology.\n---\n****************\nMinor Points: *\n****************\n1. **Spelling and Typographical Errors:**\n- **Abstract, Line 40**: \"complementary mechanisms\" should be \"complementarity mechanisms\" to maintain consistency with the terminology used throughout the document.\n- **Introduction, Line 54**: \"cooccurrence\" should be corrected to \"co-occurrence\" to ensure consistent hyphenation as used elsewhere in the text.\n- **Results, Line 110**: \"metabolic complementary indices\" should be \"metabolic complementarity indices\" to match the term used in the Methods section and maintain term consistency across the paper.\n- **Discussion, Line 354**: \"Notably, these synergistic relationships can be asymmetrical\" - The term \"asymmetrical\" should consistently match the terminology used in the rest of the document. If 'asymmetric' is used elsewhere, it should be used here too.\n- **Methods, Line 477**: \"Genome-scale metabolic models (GSMMs) for all MAGs were constructed using CarveMe...\" - Ensure that all acronyms like 'GSMMs' are consistently defined at their first occurrence. If 'GSMMs' is first used here, a prior definition is necessary.\n2. **Punctuation and Formatting:**\n- **Results, Line 195**: \"In a similar vein, our findings from a composting system...\" - Check for consistent use of serial commas in lists throughout the document.\n- Throughout the document, consistency in comma usage before conjunctions in compound sentences should be checked and standardized according to the chosen style guide (e.g., APA, MLA).\n3. **Clarification and Precision:**\n- **Results, Line 210**: The phrase \"This finding aligns with our co-occurrence network analysis...\" could be made clearer. Specify what aspect of the findings aligns to avoid vague references.\n- **Materials and Methods, Line 455**: The term \"back-mapped\" might be too jargon-heavy or unclear for readers unfamiliar with the technique. A brief definition or alternative description could enhance clarity.\n4. **Consistency in Terminology:**\n- **Throughout the Document**: Ensure consistent terminology when referring to 'metagenome-assembled genomes.' If abbreviated as MAGs, introduce this at the first occurrence and use consistently.\n5. **Figure and Table Presentation:**\n- Provide more detailed explanations or clearer legends for complex figures such as Figure 4, where multiple layers of interaction are presented. This could include simplifying the visuals or offering a supplementary file that walks through the interpretation of these data layers step-by-step.\n6. **References:**\n- **Results, Line 109**: The reference \"Feng, K., Peng, X., Zhang, Z., Gu, S., He, Q., Shen, W., ... & Deng, Y. (2022). iNAP: An Integrated Network Analysis Pipeline for Microbiome Studies.” iMeta 1: e13.\" is missing. This reference should be included where the iNAP tools are discussed or results derived from these tools are mentioned to provide proper citation and acknowledge the source.", "reviewer_2": "In this manuscript Peng et al combine the power of co-occurrence network and metabolic networks to look for metabolic interdependencies. The concept is very interesting, and if just focusing on the described results I would consider them novel. However, I worry that the details as written for the methods make it unclear to me if the results are even valid.\nThe authors collected in total 40 samples from 8 different temperature hot spring sediments. The extracted DNA followed by short read and long read sequencing. From there, the authors coassembly samples of the same temperature which they claim increases robustness. A cleaner approach would have used both single sample assembly and co-assembly approach https://link.springer.com/article/10.1186/s40168-022-01259-2 and follow by bin refinement. Independent of that, they selected 405 MAGs, and instead of clearly stating which MAGs was obtained from which temperature sample, they went ahead and mapped and then looked at percentage of reads. Which introduces yet more biases, because at this point, it is not clear to me if the reads were even across all samples. The mapping is also the source of the calculation for the co-occurrence network. SO overall, without clearing the methods section, I am afraid it is hard to validate the results.\nHere I also have more detailed comments on the rest of the manuscript:\nLine 52 – Suggestion: a recent article that came out a few days ago fits really well to this statement DOI: https://doi.org/10.1128/jb.00073-24\nLine 59 - It is good that Blanchet et al is cited in here. However, my concern is that as Blanchet et al states, co-occurrence is not evidence of ecological interactions\nLine 60 – Yes, the title of Zelezniak et al states that metabolic dependencies drive species co-occurrence. And while I strongly believe that metabolic dependencies and co-occurrence patterns are interlinked, the question is, what came first? The metabolic dependencies or the co-occurrence? Would be in that case hard to argue if dependencies drive co-occurrence patterns, or co-occurrence patters drive the evolution of dependencies (See figure 3 in your other citation in this statement Kost et al). Please consider this argument and if you agree I please modify the statement in starting in line 60.\nLine 114 - Please consider including the number of samples\nLine 118 – Figure S2 is currently not indicating the estimated coverage of each sample\nLine 119 – you statement “increasing trend of sequencing diversity (Nd) from the high- to low-temperatures” points to figure 2a, but looking at figure 2a, I see in the poorly labled x axis what I interpret as temperature and in the y axis the sequencing effort in basepair. I do not see how figure 2a displays sequencing diversity with the two parameters displayed.\nLine 122 – You mention the 405 MAGs are high quality, however in your methods you describe completeness >70% and contamination <10%. According to MiMAG quality standards https://www.nature.com/articles/nbt.3893/tables/1 these correspond to medium quality\nQuestion to paragraph line 114 10 1124: Moreover, even after checking coverage of sequencing, once the quality MAGs are selected, how well do they represent the samples?\nLine 127 – Figure 2b. 396 MAGs used to reconstruct the tree. Why did you discard 9 MAGs to reconstruct the tree?\nLine 147 – How many samples per network?\nLine 247 – Figure 6b displays MT first with red color and it is confusing as ET was shown with red color in figure 2 and displayed first.\nLine 272 – only a handful of partnerships exhibited strong metabolic complementarity\nLine 311 – co-occurrence networks cannot infer relationships – remember Blanchet et al which you cite\nLine 372 – It seems unsusual to bring up comparative results that the same authors generated only in the discussion. The results of the composting system maybe need to be mentioned in results as well.\nLine 423 – please briefly describe the method\nLine 455 – was the mapping competitive?\nLine 458 – to my knowledge, CheckM only calculates completeness and contamination of MAGs and not coverage\nLine 459 – If you only mapped MAGs back to the reads of the samples they were binned from, how did you compare the composition of all different samples using NMDS? However if you mapped the MAGs to all samples, then I would expect that you first perform dereplication of all 405 MAGs to eliminate redundancy between samples.\nLine 465 – This rule does not make sense to me. If one sample has a lot more reads than a sample from another temperature, then the results are already biased by this rule.", "reviewer_3": "In this manuscript, the authors introduce a workflow incorporating co-occurrence and metabolic modelling to describe interactions within microbial community inhabiting hot spring sediments across a temperature gradient. Potential synergistic relations are interpreted among different community members, and the influence of temperature range, genome size differences and phylogenetic distances are explored. While this study is overall well-presented and engaging, and highlights aspects of cooperative interactions within complex microbiomes relevant to various ecological habitats, many of the findings and conclusions are already well established. Moreover, the reported results are heavily based on co-occurrence analysis, and though valuable, this method is widely used to interpret microbial community dependencies in existing literature. Although the integration of methods in a single workflow may be new, the individual components are well established, and thus appears as incremental. Moreover, without support of e.g. functional omics, metabolomics, spatial organization analysis etc, many of the statements made in this manuscript remains descriptive and speculative. Although the authors attempt to validate their analysis by exploring a composting system, the experimental details and results are insufficiently described. Overall, while the presented study and workflow is not uninteresting, the novelty of the research is not adequately conveyed.\nL122: Here, the authors states that they have 405 high-quality MAGs, yet the Metods section describes a quality cut-off on 70% completeness and 10% contamination, hence not high quality. Please refer to the MIMAG standards for genome quality.\nL140: This is a repetition of L135.\nL158: This is a repetition on L149.\nL228-229: This is highly speculative without validation of cell aggregation by microscopically observations.\nL286-293: The analysis is based only on the genetic potential, and statements like this would need backup of functional omics or metabolomics for validation. The text should reflect this limitation.\nL448-449: This way of estimating genome size of MAGs is problematic in several ways. For example, completeness and contamination estimates are based on presence/absence of single-copy marker genes and does not indicate the size of the missing fragment. Is the method correctly described in the manuscript, and if so; what is the rationale behind this approach? This should include references to others who have done this calculation previously.\nL420: This sentence reads “All samples were immediately stored at -80C AFTER being transported to the laboratory” but it does not mention how is stored immediately after sampling and during transport, or how long the time between sampling and arriving the laboratory was. Even exposure to ambient temperatures at short periods of time can change induce changes in the microbial community.\nL465-467: While this may indicate the temperature preference of a particular population, it does not provide proves the “source”.\nL513- 514: As far as I can see, this is highly hypothetical and speculative without any functional data. Which should be reflected in the\nL539 (whole chapter): This experiment is not sufficiently explained."}, "author_response": "Responses to the reviewers’ comments\n1. Response to Reviewer #1\n1.1 Comments:\nIn the manuscript titled \"Deciphering Metabolic Interdependencies in Thermophilic Communities Using Co-Occurrence and Complementarity Networks,\" Xi Peng and colleagues attempt to unravel the complex web of metabolic interdependencies in microbial communities residing in thermophilic environments. Employing a novel integrative approach combining metagenome-assembled genomes (MAGs) with random matrix theory (RMT) to construct both co-occurrence and metabolic complementarity networks, the study highlights how microbial interactions intensify under extreme temperature conditions. Key findings suggest that metabolic exchanges, particularly involving amino acids and coenzyme A derivatives, are crucial for the survival and adaptation of microbial communities in high-temperature ecosystems. This research contributes to our understanding of microbial ecology by linking environmental stressors with metabolic cooperation strategies, potentially informing future biotechnological applications in bioenergy production and environmental remediation.\nThe major issue in the manuscript revolves around the interpretation of microbial interactions, which predominantly suggests a framework of cooperation and synergy that may not accurately reflect the nature of the interactions observed, particularly the prevalence of commensalism.\nResponse: Thank you for all your insightful comments and valuable feedback! We appreciate your thorough review and acknowledge the importance of accurately interpreting microbial interactions in our study.\nOur framework was designed to combine the analysis of co-occurrence networks and metabolic complementarity networks. We believe that two network approaches can be a reasonable refinement (or a filter) to each other. For example, those highly related species’ pairs based on co-occurrence only because they have similar responses to some specific environmental variable will not be presented in metabolic complementarity. On the other hand, other species pairs that possess great metabolic cross-feeding potential may not live in the same habitats (just like across the temperature gradients in our study). In this context, we want to determine whether certain co-occurring species have metabolic complementarity or to assess whether species with a strong tendency for metabolic complementarity (such as cross-feeding or species with incomplete metabolic pathways) are present in similar habitats, and that is where our integrative method comes in.\nIn our study, we first noticed that positive interactions were mostly observed in co-occurrence networks. However, statistical correlation cannot be used to decide the direction of interactions. Meanwhile, metabolic trade can have direction attributes. So, we chose positive interactions and metabolic complementarity as our main focus.\n1.2 While the innovative methodologies applied significantly advance our understanding of microbial ecosystems in extreme environments, it is critical to refine the conceptual framework presented. The study's tendency to categorize commensal interactions, where benefits accrue to one party without reciprocal interactions, as cooperative or synergistic, is misleading. This distinction is crucial not only for accurate ecological characterization but also for understanding the evolutionary pressures that shape these interactions.\nResponse: We appreciate your recognition of our methodologies. We fully understand the importance of accurately representing the nature of metabolic complementarity and its direction. In our revised manuscript, we restated our classification and definition of synergistic interactions for microbial species (Fig. R1)1, in the Material and Methods section. In a well-cited publication, Kost and his colleagues broadly divided ecological interactions into 2 types: synergistic interactions and antagonistic interactions. In synergistic interactions, there was commensalism (one partner gains fitness advantages from the interactions, +/0) and mutualism (both partners gain fitness advantages from the interactions, +/+). Here, we used the PhyloMInt workflow to determine the metabolic complementarity potential in two directions. When a partner showed the great ability to give out some metabolites while the other one did not, the relationship was defined as commensalism, and when both partners could significantly potentially exchange different metabolites, the relationship was then defined as mutualism.\n1. Kost, C., Patil, K. R., Friedman, J., Garcia, S. L. & Ralser, M. Metabolic exchanges are ubiquitous in natural microbial communities. Nat Microbiol 8, 2244-2252 (2023).\nFig. R1 Obligate metabolic interactions are frequently synergistic (Kost, Patil, et al, 2023)\n1.3 Re-framing the discussion in the manuscript to accurately reflect these unidirectional dependencies would provide a clearer and more precise interpretation of the data, enhancing the manuscript's contribution to microbial ecology.\nResponse: We appreciate your observation regarding the need to accurately distinguish between different types of metabolic synergies, particularly commensalism and mutualism. To address this, we have revised the paragraph related to this topic in the discussion as follows:\n“Notably, these synergistic relationships can be asymmetrical56. The distinct patterns observed in genome size differences between commensalistic and mutualistic interactions underscore this imbalance. Our results illuminated that a considerable proportion of synergistic relationships are commensalistic (Fig. 4). In these commensalistic interactions, there is a significant tendency for one participant (the taker) to have a streamlined genome, while the other (the giver) possesses a comparatively larger genome. This is particularly evident in the extremely thermal group, where over half of the commensalistic interactions involve a taker with a genome size streamlined to lower than 2Mbp (Fig. 6D). Similar trends, though less pronounced, are observed in the HT and MT groups. These findings suggest an ecological strategy among geothermal microbial communities. The commensalistic relationships often formed between microorganisms where one has undergone genome streamlining, forming loss of function (LOF), suggesting a specialized, especially energy-efficient or metabolism-conserving role, while their partners maintain a larger, potentially more versatile genomic repertoire. This dynamic could reflect an evolutionary optimization, where streamlined genomes reduce metabolic redundancy, relying instead on the metabolic versatility of their partners with more complete genomes and various capabilities. In contrast, mutualistic interactions, where both participants benefit, showed different genome patterns. These mutualistic relationships often involve species with more balanced genomic capabilities, suggesting that both partners contribute to and benefit from their combined metabolic activities. By accurately distinguishing between commensalism and mutualism, our findings provide a clearer understanding of the ecological and evolutionary dynamics within geothermal microbial communities.”\n(Lines 402-425)\nWe explicitly distinguished between commensalistic and mutualistic interactions, describing the genomic patterns and ecological implications of each type. By differentiating between these interaction types, we provided a clear and more precise interpretation of our data, as suggested.\nReviewer #1’s Minor Points\n1.4 **Spelling and Typographical Errors:**\n- **Abstract, Line 40**: \"complementary mechanisms\" should be \"complementarity mechanisms\" to maintain consistency with the terminology used throughout the document.\n- **Introduction, Line 54**: \"cooccurrence\" should be corrected to \"co-occurrence\" to ensure consistent hyphenation as used elsewhere in the text.\n- **Results, Line 110**: \"metabolic complementary indices\" should be \"metabolic complementarity indices\" to match the term used in the Methods section and maintain term consistency across the paper.\n- **Discussion, Line 354**: \"Notably, these synergistic relationships can be asymmetrical\" - The term \"asymmetrical\" should consistently match the terminology used in the rest of the document. If 'asymmetric' is used elsewhere, it should be used here too.\n- **Methods, Line 477**: \"Genome-scale metabolic models (GSMMs) for all MAGs were constructed using CarveMe...\" - Ensure that all acronyms like 'GSMMs' are consistently defined at their first occurrence. If 'GSMMs' is first used here, a prior definition is necessary.\nResponse: Thank you for your correction! We have made changes to correct all spelling errors as follows:\n- We have changed “complementary” to “complementarity” in the related text to maintain consistency with the terminology.\n- We have checked and ensured that hyphens were properly used in terms like “co-occurrence”.\n- We have replaced all “asymmetric” with “asymmetrical” to make the terminology consistent.\n- GSMM, as the abbreviation of genome-scale metabolic models, firstly appeared in the introduction section, so its first occurrence has been moved there (Line 81).\n1.5 **Punctuation and Formatting:**\n- **Results, Line 195**: \"In a similar vein, our findings from a composting system...\"\n- Check for consistent use of serial commas in lists throughout the document.\n- Throughout the document, consistency in comma usage before conjunctions in compound sentences should be checked and standardized according to the chosen style guide (e.g., APA, MLA).\nResponse: Thank you for the heads up on comma usage. We have rewritten the sentence you mentioned (“In a similar vein, our findings from a composting system where temperature is a key stress factor mirror those from hot springs.” Line 319) to avoid the use of consecutive commas.\nWe have also gone through the APA style guidelines and corrected the comma usage before conjunctions.\n1.6 **Clarification and Precision:**\n- **Results, Line 210**: The phrase \"This finding aligns with our co-occurrence network analysis...\" could be made clearer. Specify what aspect of the findings aligns to avoid vague references.\nResponse: Thank you for your suggestions! We re-described our new results, here for the comparison of the topologies of the two networks, revised as follows:\n“In addition, the topological properties of metabolic networks varied across temperatures. The hottest ET network, despite its smaller size, exhibited the highest average clustering coefficient, highest density, and shortest communication paths (Table 2), suggesting a highly interconnected community in this harsh environment. However, the patterns in average degree, average path distance, and harmonic geodesic distance were opposite to those of the co-occurrence networks. Such findings illustrated that metabolic complementarity and co-occurrence networks do not exhibit the same characteristics of microbial communities, which were then reflected in the topology of the network.” (Lines 217-225)\n1.7 **Materials and Methods, Line 455**: The term \"back-mapped\" might be too jargon-heavy or unclear for readers unfamiliar with the technique. A brief definition or alternative description could enhance clarity.\nResponse: Thank you for your suggestion! Due to the change in our method used for the quantification of MAGs, we re-described the method here as follows:\n“The refined bins were quantified by using CoverM (v0.7.0, genome mode and coupled reads as input) with parameter: --method relative_abundance. In the output relative abundance table (referred to MAG table), every column representing a sample contains unmapped reads percentage, making it reasonable to compare relative abundance across samples.” (Lines 514-518)\n1.8 **Consistency in Terminology:**\n- **Throughout the Document**: Ensure consistent terminology when referring to 'metagenome-assembled genomes.' If abbreviated as MAGs, introduce this at the first occurrence and use consistently.\nResponse: Thank you for your correction! The terminology of metagenome-assembled genomes, MAGs, first appeared in the introduction section.\n1.9 **Figure and Table Presentation:**\n- Provide more detailed explanations or clearer legends for complex figures such as Figure 4, where multiple layers of interaction are presented. This could include simplifying the visuals or offering a supplementary file that walks through the interpretation of these data layers step-by-step.\nResponse: Thank you for your suggestion! We also agree that the readability of Figure 4 is not satisfying. In this case, we plotted the pie charts to simply and explicitly describe the classification of all interspecies relationships. We also represented synergistic interactions as green connecting lines in the corresponding network plots.\nFig. R2 (Fig. 4 in the revised manuscript) The distributions of three categories of interactions based on the pairwise MI values and their corresponding complementarity\nMI networks.\ncomplementarity\n1.10 **References:**\n- **Results, Line 109**: The reference \"Feng, K., Peng, X., Zhang, Z., Gu, S., He, Q., Shen, W., ... & Deng, Y. (2022). iNAP: An Integrated Network Analysis Pipeline for Microbiome Studies.” iMeta 1: e13.\" is missing. This reference should be included where the iNAP tools are discussed or results derived from these tools are mentioned to provide proper citation and acknowledge the source.\nResponse: Thanks for reminding us about this citation! We have added the citation in the proper places.\n2. Response to Reviewer #2\n2.1 In this manuscript Peng et al combine the power of co-occurrence network and metabolic networks to look for metabolic interdependencies. The concept is very interesting, and if just focusing on the described results I would consider them novel. However, I worry that the details as written for the methods make it unclear to me if the results are even valid.\nThe authors collected in total 40 samples from 8 different temperature hot spring sediments. The extracted DNA followed by short read and long read sequencing. From there, the authors coassembly samples of the same temperature which they claim increases robustness. A cleaner approach would have used both single sample assembly and co-assembly approach https://link.springer.com/article/10.1186/s40168-022-01259-2 and follow by bin refinement.\nResponse: Thank you for your valuable feedback on our manuscript! We appreciate your insights and have addressed your concerns as follows:\nWe acknowledge your suggestion regarding the use of both single-sample assembly and co-assembly. Our decision to co-assemble samples of the same temperature gradient was based on several considerations:\n1. Improving assembly of low-abundance fragments: Co-assembling similar samples, such as replicate samples in our study, enhances the reconstruction of low-abundance genes and contiguous fragments.\n2. Avoiding fragmentation of genes due to strain variation: By not mixing samples across different conditions, we avoid issues arising from highly similar but distinct strains, which could otherwise lead to fragmented gene assemblies.\n3. Computational constraints: The computational burden, particularly in terms of memory usage, for co-assembling all samples was prohibitively high. Therefore, we opted for a balanced approach by co-assembling replicate samples to achieve the best possible assembly quality within our resource limits.\nTaking these three considerations into account, we believe that the goal of our replicate sample co-assembly strategy is consistent with the mix-assembly in the literature you mentioned. We have clarified these points in the revised Material and Methods section to provide a better understanding of our approach and its rationale.\n2.2 Independent of that, they selected 405 MAGs, and instead of clearly stating which MAGs was obtained from which temperature sample, they went ahead and mapped and then looked at percentage of reads.Which introduces yet more biases, because at this point, it is not clear to me if the reads were even across all samples. The mapping is also the source of the cçalculation for the co-occurrence network. SO overall, without clearing the methods section, I am afraid it is hard to validate the results.\nResponse: Thank you for this concern! We apologize for the lack of clarity regarding the origin of MAGs. To address your concern, we have made the following changes to our current analyzing methods:\n1. We adjusted the completeness and contamination thresholds to 50% and 10%, respectively, for filtering medium- and high-quality MAGs to match the criteria required in MIMAG1. Then, we dereplicated the MAGs using a secondary clustering threshold of 0.99% and retrieved 401 medium- and high-quality MAGs. About this part, we rewrote the corresponding part in Material and Methods as below:\n“Draft bins were quality controlled using metaWRAP bin_refinement module with parameters -c 50 -x 10, indicating only MAGs with completeness higher than 50% and contamination lower than 10% were retained for the following analysis, declared to be medium- and high-quality MAGs72. The refined bins were then replicated using dRep (v3.5.0)73 with parameters -pa 0.9 -sa 0.99.” (Lines 498-502)\n2. Instead of using direct mapped reads to represent the abundance of MAGs, we turned to use the relative abundance method provided in CoverM, to calculate the relative abundance of each MAG. The method has been previously used and cited by several studies2,3. This method calculates the relative abundance of a genome by dividing the number of read bases that align with the genome by its total length2. In the result table, in addition to the relative abundance of MAG in each column, there is a column representing the percentage of unmapped reads for each sample. About this part, we rewrote the corresponding part in Material and Methods as below:\n“The refined bins were quantified by using CoverM (v0.7.0, https://github.com/wwood/CoverM, genome mode and coupled reads as input) with parameter: --method relative_abundance. In the output relative abundance table (referred to MAG table), every column representing a sample contains unmapped reads percentage, making it reasonable to compare relative abundance across samples.” (Lines 514-518)\n3. With the new MAG relative abundance table obtained, we can better perform MAG assignment to each temperature group. Due to the existence of zero in the relative abundance table, we used majority selection to assign MAGs. In detail, a MAG must have a relative abundance that is not zero in more than half of the samples under a given grouping (given that there are 2-3 samples under each grouping and 5 replicates in each) in order for it to be considered to be present under that grouping. Besides, we kept the original binning source in the MAG name such as T85.bin.4, indicating the MAG was binned from T85. In this case, no matter whether it was assigned to ET (extremely thermal, T85, and T78), it will be assigned to ET anyway. With regard to this part, we have provided the corresponding additional explanations in Material and Methods as follows:\n“MAGs were assigned to three groups (ET, HT, MT) based on their relative abundance in samples. Each group comprised multiple samples with five replicates. A MAG was considered present in a group if it had non-zero abundance in more than half of the samples (≥6 for ET; ≥8 for HT and MT). If a MAG did not meet this criterion in any group, it was assigned to the group with the highest total relative abundance. Additionally, MAGs were ensured to be in the group corresponding to their sample prefix (e.g., T85 and T78 in ET).” (Lines 522-528)\n4. When constructing co-occurrence networks, we first performed majority selection of the MAG relative abundance table before calculating the Spearman’s correlation matrix, i.e., we selected as candidates only those MAGs whose relative abundance was not zero in more than 50% of the total number of samples. This was done to exclude vectors with too many zeros that might lead to inaccurate correlation calculations. About this part, we rewrote the corresponding part in Material and Methods as below:\n“With the abundance table, the pairwise Spearman’s rank correlation of each two MAGs was calculated. Since the relative abundance table originated from the distribution of hypothetical species in a natural environment, it can be regarded as the matrix that the random matrix theory (RMT)-based approach requires. Before using the RMT cutoff tool on iNAP website43, majority_selection tool was used to filter MAGs that had zero abundance in less than half of the samples. This step is to exclude biases from correlation calculations due to the effect of too many zeros. In this study, the test method used for the RMT threshold setting was the Kolmogorov-Smirnov test. The separated noise formed a network representing MAG pairs with a high abundance correlation, inferring their co-occurrence.” (Lines 528-537)\n2.3 Here I also have more detailed comments on the rest of the manuscript:\nLine 52 – Suggestion: a recent article that came out a few days ago fits really well to this statement DOI: https://doi.org/10.1128/jb.00073-24\nResponse: It is a good suggestion! We have read this literature carefully, and we agree with the importance of model microbial communities for microbial ecology studies and believe that the knowledge gaps in understanding microbial interactions mentioned in the paper do exist, and the perspective is now cited in our manuscript (Line 67).\n2.4 Line 59 - It is good that Blanchet et al is cited in here. However, my concern is that as Blanchet et al states, co-occurrence is not evidence of ecological interactions\nResponse: Thank you! Your concern here is very necessary and we fully understand it regarding the use of co-occurrence networks as evidence of ecological interactions, as highlighted by Blanchet et al. (2020). We have clarified this point in our revised manuscript to better reflect the limitations of co-occurrence networks and our approach to addressing them. Specifically, we have revised the relevant sentences to acknowledge that co-occurrence alone does not provide conclusive evidence of ecological interactions. Instead, we emphasize our strategy of combining co-occurrence data with evidence of metabolic complementarity to infer more robust interactions and explore the potential underlying mechanisms. We rewrote this sentence as follows:\n“However, as several previous studies argue, co-occurrence alone does not provide conclusive evidence of ecological interactions15-17. To address this limitation, we combine co-occurrence data with evidence of metabolic complementarity to infer more robust interactions and explore potential underlying mechanisms.” (Lines 71-75)\n2.5 Line 60 – Yes, the title of Zelezniak et al states that metabolic dependencies drive species co-occurrence. And while I strongly believe that metabolic dependencies and co-occurrence patterns are interlinked, the question is, what came first? The metabolic dependencies or the co-occurrence? Would be in that case hard to argue if dependencies drive co-occurrence patterns, or co-occurrence patters drive the evolution of dependencies (See figure 3 in your other citation in this statement Kost et al). Please consider this argument and if you agree I please modify the statement in starting in line 60.\nResponse: Thank you! We agree that the direction of causality is not always clear-cut and that both factors likely influence each other. We have revised the statement you mentioned to acknowledge this complexity and highlight the need for further investigation into the dynamic interplay between metabolic interactions and microbial community assembly. The relevant sentences have been rewritten as follows:\n“Metabolic dependencies have been proposed as a major driver of species co-occurrence18. However, it has also been suggested that this driving relationship may also be reversed2. This highlights the need for a nuanced understanding of the dynamic interplay between metabolic interactions and microbial community assembly.” (Lines 76-79)\n2.6 Line 114 - Please consider including the number of samples\nResponse: Thank you for your suggestion! We have specified the number of samples and replicates in this sentence. The revised sentence is:\n“A total of 449 million Illumina short reads (223.35Gb of raw bases, N=40) and 51 million Nanopore long reads (52.01Gb of raw bases, N=40) were obtained for the hot spring sediment samples across a temperature gradient (63.5°C~85.8°C), with eight samples each having five replicates.” (Lines 131-134)\n2.7 Line 118 – Figure S2 is currently not indicating the estimated coverage of each sample\nResponse: Thank you for your kind reminder! We omitted this figure from the supplementary material. We have updated the figure.\n0 100%Coverage\n.\n1\n95%Coverage\ne\n8\na\ng\n0\n. T85(85.8℃)\nr\ne T78(78.5℃)\nv\no T73(73.9℃)\nc 6\n0 . T70(70.8℃)\ng\na T67(67.5℃)\nr\ne\nv T65(65.8℃)\na 4\n0 . T64(64.3℃)\nd\ne\nt\nT63(63.5℃)\na\nm\n2\ni\nt .\ns 0 Actual sequencing effort\nE\n↑ Exponential value of N\nd\n103 105 107 109 1011 1013\nSequence effort (bp)\nFig. R3 (Fig. S2 in the revised manuscript) The Non-pareil curves of 40 geothermal spring metagenome samples.\nThe horizontal dashed lines indicate 100% (light red) and 95% (red) coverage. The empty circles represent the sequencing efforts, and the arrows indicate the Nonpareil index of sequencing diversity.\n2.8 Line 119 – you statement “increasing trend of sequencing diversity (Nd) from the high- to low-temperatures” points to figure 2a, but looking at figure 2a, I see in the poorly labled x axis what I interpret as temperature and in the y axis the sequencing effort in basepair. I do not see how figure 2a displays sequencing diversity with the two parameters displayed.\nResponse: Thank you for pointing out this problem! The axis title of this chart and its description in the manuscript are indeed unclear and misleading. We have modified this figure, and its description as follows:\n- We have changed the metric represented by the vertical axis N (the Nonpareil index d of sequence diversity) values to make this figure more directly convey the trend of increasing N values with decreasing temperature in the samples. N was proposed d d to be a coverage-independent measurement of the diversity of the sampled community4.\n- We have added the title “Sample” to the horizontal axis to avoid confusion.\nAfter revision, the figure is displayed as below.\nFig. R4 (Fig. 2A in the revised manuscript)\nSequence diversity (N ) for each temperature group, quantified using Illumina d sequencing. Each lollipop plot represents a specific temperature group, with the endpoints indicating the maximum and minimum N values observed across five d replicates. The central point on each lollipop plot denotes the average N value for the d respective temperature group, calculated from the same five replicates.\n2.9 Line 122 – You mention the 405 MAGs are high quality, however in your methods you describe completeness >70% and contamination <10%. According to MiMAG quality standards https://www.nature.com/articles/nbt.3893/tables/1 these correspond to medium quality\nResponse: Thank you for your kind reminder! We apologize for arbitrarily selecting the criterion of high-quality MAGs. Our original intention in making such a choice was to find a balance between the number of MAGs (the more MAGs, the more conducive to building microbial interaction networks) and the quality of MAGs (the higher the completeness of MAGs, the more conducive to building better genome-scale metabolic models). However, to ensure the rigor of the study, we re-selected MAGs according to the standards of medium- and high-quality MAGs proposed by MIMAG1 and reproduced all analyses.\n2.10 Question to paragraph line 114 10 1124: Moreover, even after checking coverage of sequencing, once the quality MAGs are selected, how well do they represent the samples?\nResponse: Thank you for this insightful comment! We acknowledge that the representativeness of MAGs is crucial for accurately depicting the microbial communities within environmental samples. The study by Meziti et al. highlights several important considerations, such as completeness/contamination assessment and metagenomic read recruitment5. According to this paper, although the representativeness of MAGs for environmental samples still needs to be cautious. In our study, we believe that its representativeness is acceptable after ensuring the quality of metagenomic sequencing and quality of retrieved MAGs.\n2.11 Line 127 – Figure 2b. 396 MAGs used to reconstruct the tree. Why did you discard 9 MAGs to reconstruct the tree?\nResponse: It is a good question! The tree was inferred using PhyloPhlAn 3.0, a tool that is useful for microbial genome characterization and phylogenetic analysis. Specifically, the tool provides two pre-build universal marker gene databases for prokaryotes, PhyloPhlAn (400 universal marker genes) and AMPHORA2 (136 universal marker genes). When building a phylogeny, PhyloPhlAn stipulates that input genomes or proteomes that map to less than 25% (default value, which can be changed by set the parameter --min_num_markers) of markers will be discarded. That was why some MAGs went missing in our tree.\nIn our revised manuscript, we adjusted the method of constructing MAG phylogeny. We used bacterial and archaeal trees generated by GTDB-Tk (inferring MAG taxonomic information), which were inferred from 120 and 122 single-copy marker genes, respectively. Since the bacterial and archaeal trees were both unrooted, we used the midpoint rooting method to infer the midpoint and put the root there, performed by midpoint.root(tree) function in phytools R package. Two trees were then merged by connecting the midpoint roots. In fact, in the document of GTDB-Tk, researchers have claimed that their opinions on the merging of archaeal and bacterial trees remain open because two trees were inferred from different marker genes.\nThe new tree is shown below:\nFig. R5 (Fig. 2D in the revised manuscript) The phylogenetic tree of 401 medium- and high-quality MAGs.\nThe tree's background color represented the domain of MAGs (light blue for bacteria, light red for archaea). The branches represented the phylum of MAGs. The bar plot in the outer circle represents the relative abundance of MAGs.\n2.12 Line 147 – How many samples per network?\nResponse: Thank you for asking about this detail! Each of the three networks represented three thermal groups (ET, HT, and MT). They had 10, 15, and 15 samples, respectively. We have clarified the numbers of samples of each group in revised manuscript as below:\n“The MAG abundance varied significantly across temperature ranges (PERMANOVA, p<0.001) (Fig. 2C), leading to the classification of sampling sites into three groups: extremely thermal (ET, 78.5-85.8°C, 2 samples × 5 replicates), highly thermal (HT, 67.5-73.9°C, 3 samples × 5 replicates), and moderately thermal (MT, 63.5-65.8°C, 3 samples × 5 replicates).” (Lines 144-148)\n2.13 Line 247 – Figure 6b displays MT first with red color and it is confusing as ET was shown with red color in figure 2 and displayed first.\nResponse: Thank you for this reminder! We did mis-assign the colors to the groups here. In the reworked results, we corrected this. The reworked figures are shown below:\n(A) (D)\n0.20 T85r=−0.498,p<0.0001 ETGroup\nT T 7 7 8 3 r r = = − − 0 0 . . 4 0 9 8 0 9 , , p p < = 0 0 . . 0 0 0 7 0 5 1 6 Meanestimatedgenomesize\nmalized\nAverage\nAbundance\n0\n0\n0\n.\n.\n.\n0\n1\n1\n5\n0\nT\nT T\nT\n6\n6\n6 6\n7\n7\n5\n4 3\n0\nr\nr\nr r\nr\n=\n=\n=\n=\n=\n−\n0\n0\n0 0\n0\n.\n.\n.\n.\n4\n3\n3 3\n.0\n2\n9\n0 0\n9\n2\n1\n4 9\n1\n,\n,\n,\n,\n,\np\np\np p\np\n<\n<\n<\n<\n=\n0\n0\n0 0\n0\n.\n.\n.\n.\n0\n0\n0 0\n.0\n0\n0\n0 0\n6\n0\n0\n0 0\n7\n1\n1\n1 1\n6 1990927.20bp\nNor 0.00\nSpearman’sRho=0.158\np-value<0.0001\n-0.05\n2 4 6\nEstimated Genome Size(Millionbasepairs)\n(B) (E)\nC C F Estimatedgenomesize Meanestimatedg H e T no G me ro si u ze p\nF E 2808581.29bp\nE I\nI J\nJ P\nP H\nH G\nG A\nA K\nK L pvalue Spearman’sRho=0.156\nL Q <0.05\nQ <0.01 p-value<0.0001\n<0.001\n-1 -0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 0.8 1\n(C) Spearman’scorrelation (F)\nEnergy production and conversion(C) *** MTGroup\nNucleotide transport and metabolism(F) *** Meanestimatedgenomesize\nAmino Acid transport and metabolism(E) *** 2908186.45bp\nLipidtransportand metabolism(I) ***\nTranslation,rib a o n s d om bio a g ls e t n r e u s c i t s ur ( e J s ) ***\nns (P) Inorganiciontransportandmetabolism\n** (H) Coenzymetransportandmetabolism\n*** (G) Carbohydratetransportandmetabolism\n*** (A) RNAprocessingandmodification\n*** (K) Transcription Spearman’sRho=0.121\n*** (L) Replication,recombinationandrepair p-value<0.0001\n*** ( tr Q a ) n s S p e o c r o t n a d n a d ry ca m ta e b ta o b li o s l m itesbiosynthesis,\n-0.2 -0.1 0.0 0.1 0.2\nEffectSize\nFig. R6 (Fig. 6 in the revised manuscript) Relationships between MI and\ncomplementarity\nestimated genome size and biological function proportion.\n(A) Correlation between the normalized average abundance of each temperature group and the estimated genome size. (B) Correlation matrix on proportions of twelve COG biological functions in MAGs. The lines between estimated genome size and each COG biological function indicated Spearman’s correlation. The color shade of the lines represented the strength of the correlation, and the thickness represented the confidence level. (C) Effect of COG biological function proportion on MI by linear\ncomplementarity\nmixed-effects models (LMMs). Bar length and error bar indicated the mean values and standard errors of the estimated effect sizes. Significant effects are represented by asterisks: **p<0.001, **p<0.01, ‘ns’ stands for not significant. (D~F) Correlation\nbetween MI and genome size difference. The fitted curve was modeled using\ncomplementarity\nthe generalized linear model. The correlation was tested using Spearman’s rank\ncorrelation test.\n2.14 Line 272 – only a handful of partnerships exhibited strong metabolic complementarity\nResponse: Thank you for your feedback! This comment is a partial sentence without further details or specific concerns. Could you please provide additional information or clarify your concerns regarding this part of our manuscript? We want to ensure that we address your comments thoroughly and make the necessary improvements to our work.\n2.15 Line 311 – co-occurrence networks cannot infer relationships – remember Blanchet et al which you cite\nResponse: Thank you for your reminder! We understand your concern regarding the interpretation of co-occurrence network approaches and their limitations in providing evidence for ecological interactions. We have revised the relevant sentence in our manuscript to better reflect this point:\n“Conversely, co-occurrence network approaches infer species relationships in a series of real natural communities, but they only reflect superficial co-occurrence patterns that might not indicate true ecological interactions16,17,47. Therefore, we aim to combine the strengths of both approaches to address their limitations (Fig. 1).” (Lines 358-32)\nThe revision aims to clarify that while co-occurrence networks provide valuable insights into species relationships, they may not necessarily indicate true ecological interactions. By combining co-occurrence data with evidence of metabolic complementarity, we strive to develop a more comprehensive understanding of microbial interactions.\n2.16 Line 372 – It seems unsusual to bring up comparative results that the same authors generated only in the discussion. The results of the composting system maybe need to be mentioned in results as well.\nResponse: Thank you for this suggestion! Our original intention of citing Zhao’s composting experiment and conducting analysis was to verify the patterns we observed using our newly proposed method in hot springs, an environment where temperature is the main stress factor.\nIn our revised manuscript, we added a section in the Results section to describe in more detail the acquisition of relevant MAG and the related results. All descriptions in the Material and Methods section are also more detailed.\nRevised paragraph in Material and Methods section:\n“Zhao et al. conducted a 30-day composting experiment in 10 composting piles at a food waste composting facility. As the temperature changed during the experiment, they selected six sampling time points, with three replicates from each composting pile for each sampling, obtaining a total of 180 samples for metagenomic sequencing. After co-assembly and binning, they obtained 159 high-quality MAGs (completeness > 90%, contamination < 5%) for further study. Three sampling time points were selected, Day 0, 5, and 25 (tagged as D00, D05, and D25), with average temperatures of 26.02°C, 64.29°C, and 41.35°C. These three-time points were regarded as three states of temperature change in the composting experiment (Day 0: the beginning of the experiment, the coolest temperature; Day 5: the highest temperature; Day 25: the temperature dropped in the late stage of the experiment). The relative abundance table of MAGs at these three sampling time points was summarized as a proxy for the state of these MAGs at that time point. The metabolic models and metabolic complementarity indices of MAGs were constructed and calculated by the above methods. In the three groups, MAGs with relative abundance greater than 0 were assigned to the corresponding group for the construction of metabolic complementarity networks. The threshold for identifying metabolic complementarity within networks was established at an RMT cutoff of 0.280.” (Lines 605-623)\nThe revised paragraph in the Results section:\n“The validation of our study on composting system experiment\nIn hot spring habitats, we noticed that symmetric interactions (mutualism) were rare. Besides, in commensalistic and mutualistic interactions, the genome size differences between species pair varied. In a similar vein, our findings from a composting system where temperature is a key stress factor mirror those from hot springs. A low occurrence of mutualistic interactions across different sampling time (temperature) groups: D00(T =26.02°C), D05(T =64.29°C), and D25(T =41.35°C), with mutualism rates\navg avg avg\nat 0.23%, 0.087%, and 0.23%, respectively. Asymmetrical interactions, suggesting commensalism, were slightly more common at 4.23%, 2.73%, and 4.18%. The genome size differences in commensalistic interactions at D05 were around 1.36Mbp, similar to findings in hot springs, whereas mutualistic interactions across all groups showed negligible genome size differences (Supplementary Table S4). The correlation between phylogenetic distance and pairwise MI was also observed to be significant\ncomplementarity\nin composting system. Specifically, D05, as the highest temperature sample in the experiment, showed the strongest positive correlation. Furthermore, the pattern that the greater the genome size difference between two species, the higher their complementarity potential was confirmed by testing the significance of the correlation between MI and genome size difference.\ncomplementarity\nIn the sequencing data from the composting experiment, little did it recover overly streamlined genomes, suggesting that such streamlined genomes might arise from long-term selection rather than short-term temperature changes. Nevertheless, we observed asymmetrical synergistic interactions, with a significant proportion being commensalistic. These interactions showed a considerable difference in genome size, particularly in the highest temperature group (D05), similar to the hot spring habitats. This implies that high temperatures are a crucial factor driving synergy between genomes of vastly different sizes.” (Lines 319-346)\n2.17 Line 423 – please briefly describe the method\nResponse: Thank you for paying attention to this method! The DNA extraction method employed in our study is based on a protocol previously developed6. It integrates a classic cell disruption technique involving freeze grinding and SDS lysis in extraction buffer, ensuring good DNA recovery and intactness from various microorganisms.\nThe extraction process commenced with the sample being weighed out and ground with sterile quartz sand in a mortar under liquid nitrogen. This freeze-grinding step was repeated to ensure thorough cell disruption. The sample was then transferred to a fresh tube, and an extraction buffer containing CTAB was added. Following incubation with proteinase K and SDS, the mixture was subjected to heating and centrifugation to separate the supernatant.\nThe supernatant was carefully transferred to avoid the white surface layer, and an extraction buffer with CTAB was added to the remaining sand pellet. After another round of incubation and centrifugation, the supernatants were combined. The DNA was precipitated using 2-isopropanol and incubated at sub-zero temperatures to facilitate precipitation.\nPost incubation, the sample was warmed and centrifuged to pellet the DNA. The pellet was then washed with cold 70% ethanol and centrifuged again to remove the ethanol. For purification, we followed the manufacturer's protocol of the MO BIO PowerSoil DNA isolation kit.\nThe purity of the extracted DNA was assessed using a Nanodrop spectrophotometer, targeting a 260/280 ratio of approximately 1.8 and a 260/230 ratio of at least 1.7. The integrity of the DNA was checked via agarose gel electrophoresis.\n2.18 Line 455 – was the mapping competitive?\nResponse: Thank you for noticing this detail! In the manuscript before the revision, we did not use competitive mapping. This is indeed an important factor that can lead to biases in mapping quantification. For this part, CoverM has been used to calculate relative abundance in this study. In future studies, we will consider competitive mapping when quantifying reads.\n2.19 Line 458 – to my knowledge, CheckM only calculates completeness and contamination of MAGs and not coverage\nResponse: It is a good question! CheckM is well-known for assessing the quality of recovered microbial genomes. It also provides a command checkm coverage that can produce coverage profiles for all sequences within a set of genome bins. However, this part is also removed because we turned to use CoverM.\n2.20 Line 459 – If you only mapped MAGs back to the reads of the samples they were binned from, how did you compare the composition of all different samples using NMDS? However if you mapped the MAGs to all samples, then I would expect that you first perform dereplication of all 405 MAGs to eliminate redundancy between samples.\nResponse: Thank you for this suggestion! We apologize for not clearly describing the steps of MAG abundance calculation in our previous manuscript. Previously, we used the number of reads mapped to MAGs in each sample as the abundance of MAGs in the corresponding sample. We could then obtain a MAG table (like an ASV/OTU table in amplicon diversity analysis) to perform NMDS analysis and subsequent co-occurrence network-related analysis.\nIn our revised manuscript, we added dereplication in our study. The relevant descriptions are shown below:\n“Draft bins were quality controlled using metaWRAP bin_refinement module with parameters -c 50 -x 10, indicating only MAGs with completeness higher than 50% and contamination lower than 10% were retained for the following analysis, declared to be medium- and high-quality MAGs72. The refined bins were then replicated using dRep (v3.5.0)73 with parameters -pa 0.9 -sa 0.99.” (Lines 501-505)\n“The refined bins were quantified by using CoverM (v0.7.0, https://github.com/wwood/CoverM, genome mode and coupled reads as input) with parameter: --method relative_abundance. In the output relative abundance table (referred to MAG table), every column representing a sample contains unmapped reads percentage, making it reasonable to compare relative abundance across samples.” (Lines 514-518)\n2.21 Line 465 – This rule does not make sense to me. If one sample has a lot more reads than a sample from another temperature, then the results are already biased by this rule.\nResponse: It is a good reminder! In the MAG relative abundance table generated by CoverM, the unmapped relative abundance is retained in each sample to represent those reads that cannot be mapped to any MAG, and the total relative abundance of each sample is normalized to 100%.\n1. Bowers, R., Kyrpides, N., Stepanauskas, R. et al. Minimum information about a single amplified genome (MISAG) and a metagenome-assembled genome (MIMAG) of bacteria and archaea. Nat Biotechnol 35, 725–731 (2017).\n2. Hessler, T., Huddy, R.J., Sachdeva, R. et al. Vitamin interdependencies predicted by metagenomics-informed network analyses and validated in microbial community microcosms. Nat Commun 14, 4768 (2023).\n3. Mara, P., Geller-McGrath, D., Edgcomb, V. et al. Metagenomic profiles of archaea and bacteria within thermal and geochemical gradients of the Guaymas Basin deep subsurface. Nat Commun 14, 7768 (2023).\n4. Rodriguez-R L.M., Gunturu S., Tiedje J.M. et al. Nonpareil 3: Fast Estimation of Metagenomic Coverage and Sequence Diversity. mSystems 3, e00039-18 (2018).\n5. Meziti A., Rodriguez-R L.M., Hatt J.K. et al. The Reliability of Metagenome-Assembled Genomes (MAGs) in Representing Natural Populations: Insights from Comparing MAGs against Isolate Genomes Derived from the Same Fecal Sample. Appl Environ Microbiol, 87, e02593-20 (2021)\n6. Zhou J., Bruns M.A., et al. DNA recovery from soils of diverse composition. Appl Environ Microbiol 62, 316-322 (1996).\n3. Response to Reviewer #3\n3.1 In this manuscript, the authors introduce a workflow incorporating co-occurrence and metabolic modelling to describe interactions within microbial community inhabiting hot spring sediments across a temperature gradient. Potential synergistic relations are interpreted among different community members, and the influence of temperature range, genome size differences and phylogenetic distances are explored. While this study is overall well-presented and engaging, and highlights aspects of cooperative interactions within complex microbiomes relevant to various ecological habitats, many of the findings and conclusions are already well established. Moreover, the reported results are heavily based on co-occurrence analysis, and though valuable, this method is widely used to interpret microbial community dependencies in existing literature. Although the integration of methods in a single workflow may be new, the individual components are well established, and thus appears as incremental. Moreover, without support of e.g. functional omics, metabolomics, spatial organization analysis etc, many of the statements made in this manuscript remains descriptive and speculative. Although the authors attempt to validate their analysis by exploring a composting system, the experimental details and results are insufficiently described. Overall, while the presented study and workflow is not uninteresting, the novelty of the research is not adequately conveyed.\nResponse: Thank you for your thoughtful and detailed review of our manuscript. We appreciate your positive feedback on the presentation and engagement of our study and your recognition of its relevance to understanding synergistic interactions within complex microbiomes. We also value your constructive criticisms and have addressed each of your concerns as follows:\n1. You mentioned that many of the findings and conclusions are already well-established. We acknowledge that the individual components of our workflow, such as co-occurrence analysis and metabolic modeling, are well-established methods. However, the novelty of our study lies in integrating these methods to create a comprehensive workflow that can infer microbial interactions across a temperature gradient. Besides, the application of random matrix theory noise filtering on metabolic complementarity is also a brand new asset. The integrated approach allows for a more holistic understanding of microbial community dynamics.\n2. You commented that the reported results are heavily based on co-occurrence analysis. We must clarify that in our study, we moved beyond simple co-occurrence because we have noticed that co-occurrence cannot fully be evidence of disentangling complex microbial interactions. Our focus for this study is using metabolic complementarity networks, which provide additional insights into the potential functional interactions between species. Combining these two approaches improves the robustness of synergistic interactions, and digging out the potential proxy for the interactions (e.g., transferable metabolites) is the ultimate goal.\n3. We understand the importance of supporting functional data for comprehending the metabolic and molecular mechanisms for microbial synergies. In our study, we first used two methods to filter microbial synergy to refine the interactions between microbes and explore their similarities and differences in predicting synergy. We believe that in analyzing the refined relationships, tools such as functional omics and metabolomics will be important for further analysis of these relationships.\n4. You mentioned our insufficient description of composting system validation. We appreciate this feedback, and we have provided a more detailed description of this composting experiment and how we used their work to validate the patterns we discovered in our study. We have strengthened the specificity of the sampling setup, data collection, and the results obtained, demonstrating how they support our main findings. The revised part is shown below:\nRevised paragraph in Material and Methods section:\n“Zhao et al. conducted a 30-day composting experiment in 10 composting piles at a food waste composting facility. As the temperature changed during the experiment, they selected six sampling time points, with three replicates from each composting pile for each sampling, obtaining a total of 180 samples for metagenomic sequencing. After co-assembly and binning, they obtained 159 high-quality MAGs (completeness > 90%, contamination < 5%) for further study. Three sampling time points were selected, Day 0, 5, and 25 (tagged as D00, D05, and D25), with average temperatures of 26.02°C, 64.29°C, and 41.35°C. These three time points were regarded as three states of temperature change in the composting experiment (Day 0: the beginning of the experiment, the coolest temperature; Day 5: the highest temperature; Day 25: the temperature dropped in the late stage of the experiment). The relative abundance table of MAGs at these three sampling time points was summarized as a proxy for the state of these MAGs at that time point. The metabolic models and metabolic complementarity indices of MAGs were constructed and calculated by the above methods. In the three groups, MAGs with relative abundance greater than 0 were assigned to the corresponding group for the construction of metabolic complementarity networks. The threshold for identifying metabolic complementarity within networks was established at an RMT cutoff of 0.280.” (Lines 605-623)\nThe revised paragraph in the Results section:\n“The validation of our study on composting system experiment\nIn hot spring habitats, we noticed that symmetric interactions (mutualism) were rare. Besides, in commensalistic and mutualistic interactions, the genome size differences between species pair varied. In a similar vein, our findings from a composting system where temperature is a key stress factor mirror those from hot springs. A low occurrence of mutualistic interactions across different sampling time (temperature) groups: D00(T =26.02°C), D05(T =64.29°C), and D25(T =41.35°C), with mutualism rates\navg avg avg\nat 0.23%, 0.087%, and 0.23%, respectively. Asymmetrical interactions, suggesting commensalism, were slightly more common at 4.23%, 2.73%, and 4.18%. The genome size differences in commensalistic interactions at D05 were around 1.36Mbp, similar to findings in hot springs, whereas mutualistic interactions across all groups showed negligible genome size differences (Supplementary Table S4). The correlation between phylogenetic distance and pairwise MI was also observed to be significant\ncomplementarity\nin composting system. Specifically, D05, as the highest temperature sample in the experiment, showed the strongest positive correlation. Furthermore, the pattern that the greater the genome size difference between two species, the higher their complementarity potential was confirmed by testing the significance of the correlation between MI and genome size difference.\ncomplementarity\nIn the sequencing data from the composting experiment, little did it recover overly streamlined genomes, suggesting that such streamlined genomes might arise from long-term selection rather than short-term temperature changes. Nevertheless, we observed asymmetrical synergistic interactions, with a significant proportion being commensalistic. These interactions showed a considerable difference in genome size, particularly in the highest temperature group (D05), similar to the hot spring habitats. This implies that high temperatures are a crucial factor driving synergy between genomes of vastly different sizes.” (Lines 319-346)\n3.2 L122: Here, the authors states that they have 405 high-quality MAGs, yet the Metods section describes a quality cut-off on 70% completeness and 10% contamination, hence not high quality. Please refer to the MIMAG standards for genome quality.\nResponse: Thank you for your kind suggestion! We adjusted the completeness and contamination thresholds to 50% and 10%, respectively, for filtering medium- and high-quality MAGs to match the criteria required in MIMAG1. Then, we dereplicated the MAGs using a secondary clustering threshold of 0.99% and retrieved 401 medium- and high-quality MAGs. About this part, we rewrote the corresponding part in Material and Methods as below:\n“Draft bins were quality controlled using metaWRAP bin_refinement module with parameters -c 50 -x 10, indicating only MAGs with completeness higher than 50% and contamination lower than 10% were retained for the following analysis, declared to be medium- and high-quality MAGs72. The refined bins were then replicated using dRep (v3.5.0)73 with parameters -pa 0.9 -sa 0.99.” (Lines 501-505)\n3.3 L140: This is a repetition of L135.\nResponse: Thank you for your reminder! Due to the rework, we have re-described the results here, omitting the superfluous description of Thermoproteota to avoid repetition, as follows:\n“Archaea also displayed clear temperature preferences despite accounting for only 7.31% of relative abundance across all samples. Thermoproteota favored the ET group, while Micrarchaeota, the second most abundant, exhibited a unique distribution, with nearly equal contributions from HT and MT groups. These findings underscore the diverse thermal adaptations and preferences of thermophile communities within this hot spring ecosystem.” (Lines 159-164)\n3.4 L158: This is a repetition on L149.\nResponse: Thank you for your reminder! We have deleted the sentence (formerly Line 158) to avoid repetition.\n3.5 L228-229: This is highly speculative without validation of cell aggregation by microscopically observations.\nResponse: Thank you for raising our attention to these descriptions! Our original intention of using the word “aggregation” here was to convey that species with close phylogenetic distance (PD at peak 0) and species with distance phylogenetic relationship but potentially having complementarity metabolic needs (PD at peak 2) were likely to show strong positive correlation in relative abundance (co-occurring). We understand that the term “aggregation” could be misinterpreted as referring to physical cell clustering rather than co-occurrence patterns in relative abundance, so we revised this sentence for clarification as follows:\n“However, a fascinating pattern emerged after RMT filtering (cutoff = 0.830). The density distribution of PD exhibited two distinct peaks: one near 0 and another near 2. This suggested two potential drivers of co-occurrence: the close association of similar species with shared resource requirements (PD at peak 0) and the co-occurrence of phylogenetically distant species with complementary metabolic needs (PD at peak 2).” (Lines 238-243)\n3.6 L286-293: The analysis is based only on the genetic potential, and statements like this would need backup of functional omics or metabolomics for validation. The text should reflect this limitation.\nResponse: Thank you for your suggestion! We appreciate your point regarding the need to acknowledge the limitations of our findings. We have revised this section of the manuscript to clarify that our findings were based on genetic potentials and highlighted the necessity for further validation with tools such as functional omics or metabolomics. Here is the revised paragraph for this part:\n“In exploring the intricacies of species metabolic synergies, our focus was drawn to 17 archaea-bacteria interactions, which were supported by both metabolic networks and co-occurrence networks, indicating that inter-domain synergies were prevalent in harsh environments such as hot springs. In these interactions, there were some MAGs that exhibited the properties of co-occurrence and metabolic complementarity with multiple other species. For instance, an archaeal MAG (T64.bin.92, Micrarchaeota phylum) had seven distinct bacterial partners. This centrally positioned MAG, characterized by its smaller estimated genome size of 0.857 Mbp, predominantly assumed the role of a taker in these interactions. This was evidenced by its substantial reliance on obtaining several types of metabolites, such as some coenzyme A derivatives and other crucial metabolic needs that it cannot fulfill endogenously from other species. Despite its taker-centric position, T64.bin.92 contributed to these commensalistic relationships by providing potential surplus materials like carbohydrates and nucleoside/nucleotide-related substances, embodying a reciprocal dynamic in these interactions. This necessity for external coenzyme A derivatives, vital for biosynthesis reactions involving acyl transfers, underlines the significant dependence of T64.bin.92 and its partners. Furthermore, its requirement for certain inorganic substances, including various metal ion and non-metal compounds from its metabolic counterparts, further illustrated the complexity and significance of its role as a taker. However, it is important to note that these conclusions are based on genetic potential inferred from genomic data. To confirm the actual metabolic exchanges and interactions, additional validation with functional omics and metabolomics data is necessary.” (Lines 296-317)\n3.7 L448-449: This way of estimating genome size of MAGs is problematic in several ways. For example, completeness and contamination estimates are based on presence/absence of single-copy marker genes and does not indicate the size of the missing fragment. Is the method correctly described in the manuscript, and if so; what is the rationale behind this approach? This should include references to others who have done this calculation previously.\nResponse: Thank you for your reminder! We did miss the references on adjusting genome size for estimated completeness and contamination. Our initial idea was that completeness must affect the actual genome size, so we searched for relevant studies and found that in the article by Parks D. H., et al2. The original description is: “The genome size of the UBA genomes was adjusted to account for its estimated completeness and contamination: adjusted genome size = (genome size)/(completeness + contamination).” We have updated the references and citations in our manuscript.\n3.8 L420: This sentence reads “All samples were immediately stored at -80C AFTER being transported to the laboratory” but it does not mention how is stored immediately after sampling and during transport, or how long the time between sampling and arriving the laboratory was. Even exposure to ambient temperatures at short periods of time can change induce changes in the microbial community.\nResponse: It is a great question! After swift sample collection, we promptly placed them in a liquid nitrogen container and then transported them to the laboratory within 2 days. Subsequently, they were stored at -80°C. Before the storage, sediment samples for DNA extraction were pre-divided to prevent DNA damage due to repeated freeze-thaw cycles in later experiments. We believe that through the above operations, the loss of DNA and related errors have been reduced to an optimal level. The relevant section in our manuscript has been revised as follows:\n“After sample collection, they were promptly placed in a liquid nitrogen container and transported to the laboratory within 2 days. Subsequently, they were stored at -80°C. Before storage, sediment samples for DNA extraction were pre-divided to prevent DNA damage due to repeated freeze-thaw cycles in later experiments. After 48 hours of lyophilization, total DNA within the microbial community was extracted from 1.5g of freeze-dried sediment using the grind plus kit method as previously described63.” (Lines 474-479)\n3.9 L465-467: While this may indicate the temperature preference of a particular population, it does not provide proves the “source”.\nResponse: It is a good reminder! The previous definition of “potential source” did have some controversies. In the revised manuscript, we turned to using CoverM to calculate the relative abundance of MAGs. We have described this MAG assignment rule to other reviewers (With the new MAG relative abundance table obtained, we can better perform MAG assignment to each temperature group. Due to the existence of zero in the relative abundance table, we used majority selection to assign MAGs. In detail, a MAG must have a relative abundance that is not zero in more than half of the samples under a given grouping (given that there are 2-3 samples under each grouping and 5 replicates in each) in order for it to be considered to be present under that grouping. Besides, we kept the original binning source in the MAG name such as T85.bin.4, indicating the MAG was binned from T85. In this case, no matter whether it was assigned to ET (extremely thermal, T85 and T78), it will be assigned to ET anyway). The relevant paragraph in our manuscript has been revised as follows:\n“MAGs were assigned to three groups (ET, HT, MT) based on their relative abundance in samples. Each group comprised multiple samples with five replicates. A MAG was considered present in a group if it had non-zero abundance in more than half of the samples (≥6 for ET; ≥8 for HT and MT). If a MAG did not meet this criterion in any group, it was assigned to the group with the highest total relative abundance. Additionally, MAGs were ensured to be in the group corresponding to their sample prefix (e.g., T85 and T78 in ET).” (Lines 522-528)\n3.10 L513- 514: As far as I can see, this is highly hypothetical and speculative without any functional data. Which should be reflected in the\nResponse: It is a good question! However, this comment is a partial sentence without further details or specific concerns. We guess that we should emphasize the hypothetical property when using this method. Here, we followed the principle of the PhyloMint method (the same method we used to calculate the complementarity index) and used the functions in the Python code they provide to mine metabolites in metabolic complementarity models (metabolic networks) of pairwise MAGs. The method itself is based on multiple techniques, such as genome-scale metabolic modeling and flux balance analysis, rather than solely relying on genetic potential. However, we still acknowledge that further omics analysis can provide more evidence for our findings, and this is also an important direction for our future research.\n3.11 L539 (whole chapter): This experiment is not sufficiently explained.\nResponse: Thank you for pointing out this flaw. The primary purpose of citing the composting experiment and performing our analysis was to validate the patterns we observed with our newly developed method in hot springs, where temperature acts as the main stress factor. In the revised manuscript, we have added a more detailed section in the Results to elaborate on the acquisition of relevant MAGs and the associated findings. Additionally, we have enhanced the descriptions in the Materials and Methods section to provide greater detail. (See Reponses for 3.1)\n1. Bowers, R., Kyrpides, N., Stepanauskas, R. et al. Minimum information about a single amplified genome (MISAG) and a metagenome-assembled genome (MIMAG) of bacteria and archaea. Nat Biotechnol 35, 725–731 (2017).\n2. Parks, D.H., Rinke, C., Chuvochina, M. et al. Recovery of nearly 8,000 metagenome-assembled genomes substantially expands the tree of life. Nat Microbiol 2, 1533–1542 (2017).\nREVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nThe authors have satisfactorily addressed the points I previously raised. The article is now ready for publication.\nReviewer #2 (Remarks to the Author):\nThe authors successfully addressed all my comments.\nI have a few comments about the figures;\nFigure 2: The phylogenetic tree should get more space to be able to appreciate it.\nFigure 2, 5 and 6: the text in the axis is too small to read and in some panels is missing and needs to be added.\nReviewer #3 (Remarks to the Author):\nI appreciate that the Authors have taken most of the concerns and suggestions raised into consideration in the revised version. I also acknowledge the Authors’ effort to reason for the novelty of the manuscript. Some good points are raised, and I recognize that I may have been too strict on this matter. This version of the manuscript shows improvement; however, some issues still need to be resolved before publication. The line numbers below correspond to the lines in the revised version (merged pdf).\nL149-L153: Using accumulative percentage (>100%) here does not make any sense (to me).\nL284: Genetic loss, or missing from the genome reconstruction? With the accepted cut-off for MAG completeness at 50%, it is important to also consider that absence of a gene in theory could be due to missing part of the genome, rather than an actual gene loss. This should be considered (and mentioned) when talking about specific genes.\nL513-515: The details around the merging of the bacterial and the archaeal trees (which I honestly am a bit skeptical to) should be described here.\nFinally, some of the sentences appear difficult to grasp and/or wording resemble spoken language, and I suggest some minor modifications to improve the clarity:\nL77-78: It is very good that this important note has been included in the revised manuscript, however the sentence could benefit from some rephrasing to clarify the point.\nL85: Consider rephrasing “there is still blank room for further exploration in this area” to e.g. “there remains substantial potential for further...”\nL99: Consider using “new” or “novel” instead of “fresh”\nL132: Change “for” to “from”\nL324-327: Is a word (e.g. “was observed”) missing?\nL335-339: A long and heavy sentence, please consider restructuring to ease the reading.\nL340-342: This sentence could also benefit from rephrasing to improve the readability.\nL372: “cope with the scorching heat” could be implied as mechanisms that are directly involved in e.g. heath resistance. Consider rephrasing.\nL558-560: Consider to re-structure the sentence to improve the readability, e.g. “In the network, synergistic interactions are determined based on whether the original PhyloMInt …” etc\nL623: Consider changing “above” to “abovementioned” or “methods described above”.\nFigure 2: The bars around the phylogenetic tree are very difficult to see. I would suggest making panel B (the tree) larger.\nFigure 3: The figure legend is very minimalistic. Please add more information to ease the interpretation for the reader. Minor note: lines in the figure key are not aligned for B and C.\nFigure 4: Same as for figure 4; The figure legend is not sufficient.\nFigure 5: Legend does not follow the same style as the other (A,B,C vs a-c,d-f)\nResponses to referees\n1. Response to Reviewer #1:\n1.1 The authors have satisfactorily addressed the points I previously raised. The article is now ready for publication.\nResponse: Thank you very much for your comments and recognition. Your comments have greatly improved our manuscript!\n2. Response to Reviewer #2:\n2.1 The authors successfully addressed all my comments.\nResponse: Thank you very much for your comments and recognition.\n2.2 I have a few comments about the figures;\nFigure 2: The phylogenetic tree should get more space to be able to appreciate it.\nResponse: Thanks for the good advice. We adjusted the layout of this graph to make the tree easier to observe.\n2.3 Figure 2, 5 and 6: the text in the axis is too small to read and in some panels is missing and needs to be added.\nResponse: Thank you for reminding us! We have added the missing coordinate axes in Figure 6d-f. We also have tried to optimize the fonts and sizes of texts for these figures.\n3. Response to Reviewer #3:\nI appreciate that the Authors have taken most of the concerns and suggestions raised into consideration in the revised version. I also acknowledge the Authors’ effort to reason for the novelty of the manuscript. Some good points are raised, and I recognize that I may have been too strict on this matter. This version of the manuscript shows improvement; however, some issues still need to be resolved before publication. The line numbers below correspond to the lines in the revised version (merged pdf).\nResponse: Thank you for your recognition of our revision! We should keep improving according to your advice.\n3.1 L149-L153: Using accumulative percentage (>100%) here does not make any sense (to me).\nResponse: It is a good point! We have reconsidered this result and decided to remove it because the sum of relative abundance of one MAG or a phylum across samples are meaningless (because they are normalized by samples, let alone the unmapped read proportion). The removal of this result would not affect any other results and discussion.\n3.2 L284: Genetic loss, or missing from the genome reconstruction? With the accepted cut-off for MAG completeness at 50%, it is important to also consider that absence of a gene in theory could be due to missing part of the genome, rather than an actual gene loss. This should be considered (and mentioned) when talking about specific genes.\nResponse: Thanks for reminding us of this. As you commented, after lowering the completeness criteria to 50%, more MAGs with more possible incomplete gene profiles came in. So, we pointed it out in our latest manuscript to address this inevitable concern.\nHere is the revision:\nThis correlation was strongest in the ET group, suggesting that genome size plays more crucial role in shaping metabolic partnerships in harsh environments. These results reflected a trade-off of smaller size genomes, where the genes associated with genetic information storage and processing were retained to a greater degree while genes related to various metabolisms underwent loss, relying on synergistic interactions to compensate for their limited metabolic repertoire. However, it is essential to consider that the absence of specific genes might also be due to incomplete genome reconstruction, given the accepted cutoff for MAG completeness at 50%.\n3.3 L513-515: The details around the merging of the bacterial and the archaeal trees (which I honestly am a bit skeptical to) should be described here.\nResponse: Thanks for your reminder! We carelessly forgot to mention this in Methods. We added in our latest manuscript as below:\nThe phylogenetic trees for bacterial and archaeal MAGs were constructed using multiple sequence alignment results generated by GTDB-Tk workflow. The unrooted bacterial and archaeal phylogenetic trees were rooted using midpoint rooting method, performed by midpoint.root function in phytools R package77. and were then visualized using the online tree display tool iTOL (v6.6)78. The phylogenetic distance was then calculated using the cophenetic.phylo function in R.\n3.4 L77-78: It is very good that this important note has been included in the revised manuscript, however the sentence could benefit from some rephrasing to clarify the point.\nResponse: Thanks for the advice. We rephrased the sentences as below:\nMetabolic dependencies have been proposed as a major driver of species co-occurrence18. However, it has also been suggested that species co-occurrence could in turn drive metabolic dependencies2.\n3.5 L85: Consider rephrasing “there is still blank room for further exploration in this area” to e.g. “there remains substantial potential for further...”\nResponse: Thanks for the suggestion. We rephrased the sentence as your example.\n3.6 L99: Consider using “new” or “novel” instead of “fresh”\nResponse: Thanks for the suggestion! However, the editor suggested we should refrain from using words such as new/novel/first when referring to results or scientific findings. So, we deleted related words as the editor mentioned.\n3.7 L132: Change “for” to “from”\nResponse: Thanks for the notice. The error has been corrected.\n3.8 L324-327: Is a word (e.g. “was observed”) missing?\nResponse: Yes. We made a mistake here and we corrected it.\n3.9 L335-339: A long and heavy sentence, please consider restructuring to ease the reading.\nResponse: Thanks for your suggestion! We rephrased this sentence as below:\nFurthermore, the pattern that the greater the genome size difference between two species corresponds to higher their complementarity potential was confirmed by testing the significance of the correlation between MI and genome size difference\ncomplementarity\n(D00: Spearman’s Rho = 0.122, p<0.001; D05: Spearman’s Rho = 0.213, p<0.001; D25: Spearman’s Rho = 0.090, p<0.001).\n3.10 L340-342: This sentence could also benefit from rephrasing to improve the readability.\nResponse: Thanks for your suggestion! We also rephrased this sentence as below:\nThe sequencing data from the composting experiment recovered few overly streamlined genomes, suggesting that such genomes might arise from long-term selection rather than short-term temperature changes.\n3.11 L372: “cope with the scorching heat” could be implied as mechanisms that are directly involved in e.g. heath resistance. Consider rephrasing.\nResponse: Thank you for this suggestion! We used “adapt to” instead of “cope with” to avoid misunderstanding.\n3.12 L558-560: Consider to re-structure the sentence to improve the readability, e.g. “In the network, synergistic interactions are determined based on whether the original PhyloMInt …” etc\nResponse: Thanks for the suggestion! We rephrased this sentence as your example.\n3.13 L623: Consider changing “above” to “abovementioned” or “methods described above”.\nResponse: Thanks for this reminder. We have corrected this mistake as you suggested.\n3.14 Figure 2: The bars around the phylogenetic tree are very difficult to see. I would suggest making panel B (the tree) larger.\nResponse: Thanks for pointing out this problem. We adjusted the layout of this graph to make the tree easier to observe.\n3.15 Figure 3: The figure legend is very minimalistic. Please add more information to ease the interpretation for the reader. Minor note: lines in the figure key are not aligned for B and C.\nResponse: Thanks for your carefulness. We have adjusted the layout of the figure. Besides, we also added information for every figure’s legend/caption according to the editor’s request.\n3.16 Figure 4: Same as for figure 4; The figure legend is not sufficient.\nResponse: Thanks for this comment! We think the added information in legend/caption should solve the problem.\n3.17 Figure 5: Legend does not follow the same style as the other (A,B,C vs a-c,d-f)\nResponse: Thanks for your carefulness again! According to the figure requirement, we have changed all legends from uppercase letters to lowercase letters."}]

Decision letter