Hiroaki Fujita1*
Masayuki Ushio2
Kenta Suzuki3
Masato S. Abe4
Masato Yamamichi5,6
Yusuke Okazaki7
Alberto Canarini1Ibuki Hayashi1
Keitaro Fukushima8
Shinji Fukuda9,10,11,12
E. Toby Kiers13
Hirokazu Toju1*- 1Center for Ecological Research, Kyoto University, Kyoto, Japan
- 2Department of Ocean Science (OCES), The Hong Kong University of Science and Technology (HKUST), Kowloon, Hong Kong SAR, China
- 3Integrated Bioresource Information Division, BioResource Research Center, RIKEN, Tsukuba, Japan
- 4Faculty of Culture and Information Science, Doshisha University, Kyoto, Japan
- 5School of Biological Sciences, The University of Queensland, Brisbane, QLD, Australia
- 6Department of International Health and Medical Anthropology, Institute of Tropical Medicine, Nagasaki University, Nagasaki, Japan
- 7Institute for Chemical Research, Kyoto University, Kyoto, Japan
- 8Faculty of Food and Agricultural Sciences, Fukushima University, Fukushima, Japan
- 9Institute for Advanced Biosciences, Keio University, Tsuruoka, Japan
- 10Gut Environmental Design Group, Kanagawa Institute of Industrial Science and Technology, Kawasaki, Japan
- 11Transborder Medical Research Center, University of Tsukuba, Tsukuba, Japan
- 12Laboratory for Regenerative Microbiology, Juntendo University Graduate School of Medicine, Tokyo, Japan
- 13Amsterdam Institute for Life and Environment, Vrije Universiteit, Amsterdam, Netherlands
Facilitative interactions between microbial species are ubiquitous in various types of ecosystems on the Earth. Therefore, inferring how entangled webs of interspecific interactions shift through time in microbial ecosystems is an essential step for understanding ecological processes driving microbiome dynamics. By compiling shotgun metagenomic sequencing data of an experimental microbial community, we examined how the architectural features of facilitative interaction networks could change through time. A metabolic modeling approach for estimating dependence between microbial genomes (species) allowed us to infer the network structure of potential facilitative interactions at 13 time points through the 110-day monitoring of experimental microbiomes. We then found that positive feedback loops, which were theoretically predicted to promote cascade breakdown of ecological communities, existed within the inferred networks of metabolic interactions prior to the drastic community-compositional shift observed in the microbiome time-series. We further applied “directed-graph” analyses to pinpoint potential keystone species located at the “upper stream” positions of such feedback loops. These analyses on facilitative interactions will help us understand key mechanisms causing catastrophic shifts in microbial community structure.
Introduction
In nature, species form complex webs of interactions, thereby driving various types of community-and ecosystem-level phenomena (May, 1972; Ives and Carpenter, 2007; Allesina and Tang, 2012). Roles of interspecific interactions in sudden community structural shifts are among the most important targets of ecological research (Scheffer et al., 2001; Scheffer and Carpenter, 2003; Ratzke et al., 2020). Theoretical studies have predicted that the architecture (topology) of interaction networks determines consequences of ecological interactions such as species coexistence or community collapse (Thébault and Fontaine, 2010; Rohr et al., 2014; Levine et al., 2017). Although a number of empirical studies on plants and animals have been conducted to uncover the characteristics of species interaction networks (Olesen et al., 2007; Lever et al., 2014; CaraDonna and Waser, 2020), our knowledge of potential relationship between network architecture and community-level consequences have been limited for microbial ecosystems.
In microbial ecology, estimating the architecture of potential interactions has been increasingly common (Faust and Raes, 2012; Friedman and Alm, 2012; Berry and Widder, 2014; Kurtz et al., 2015). Amplicon sequencing (DNA metabarcoding) of prokaryote 16S rRNA gene, for example, have been frequently used to infer structure of networks depicting co-occurrence patterns of microbial species (Barberan et al., 2012; Faust et al., 2012; Berry and Widder, 2014). Nonetheless, those networks obtained with co-occurrence pattern analyzes include pairs of species that merely share environmental preferences, making it difficult to investigate webs of direct facilitative/competitive interactions between species (Warton et al., 2015; Toju et al., 2017; Kurtz et al., 2019; Blanchet et al., 2020). Moreover, although studies on co-occurrence patterns assume bidirectional associations between species, interspecific interactions in nature are not necessarily bidirectional (Sugihara et al., 2012; Ushio et al., 2018; Delmas et al., 2019). In other words, assuming mutually positive (mutualistic) interactions in all pairs of species within such co-occurrence networks (i.e., neglecting commensalistic or other types of interactions) may lead to misunderstanding of ecological community processes. Consequently, reconstructing networks consisting of not only bidirectional but also unidirectional interactions between species [i.e., “directed graphs” (Newman, 2010)] is an essential step for advancing our understanding of microbial community processes.
A promising approach for investigating complex webs of microbial interactions is to estimate potential flows of metabolites between microbial species based on metagenomic datasets (Stolyar et al., 2007; Klitgord and Segrè, 2010; Zomorrodi and Maranas, 2012; Levy and Borenstein, 2013). Because species’ ability to metabolize given chemical compounds is encoded in their genomes, metabolic modeling has been applied to infer potential metabolic interactions between microbes (Zelezniak et al., 2015; Magnúsdóttir et al., 2017). If genomic information is available for a pair of species, potential dependence of a species on the other species can be evaluated in terms of the list of metabolites presumably emitted by the other species (Zelezniak et al., 2015; Magnúsdóttir et al., 2017). By applying such community-scale metabolic modeling (Frioux et al., 2020), we will be able to gain insights into network architecture of facilitative interactions (Sung et al., 2017; Hassani et al., 2018; Gralka et al., 2020). Analyzes on temporal shifts in such metabolic interaction network architecture, in particular, are expected to enhance our understanding of processes or mechanisms causing community collapse [or dysbiosis (Carding et al., 2015; Kriss et al., 2018)]. Nonetheless, there have been, to the best of our knowledge, no study reporting changes in facilitative interaction network architecture through microbial community dynamics.
In this study, we performed an analysis of metabolic interaction networks using the shotgun metagenomic sequencing dataset of an experimental bacterial community (Fujita et al., 2022). Across the 110-day monitoring of a co-culture system of a freshwater bacterial community (Fujita et al., 2023), the previous study examined temporal shifts in the level of ecological niche overlap between species in order to infer dynamics of competitive interactions (Fujita et al., 2022). In this study, we reconstructed networks of facilitative interactions based on the metabolic modeling of shotgun metagenomic data at 13 time points across the time-series. We then evaluated changes in the architectural features of the directed graphs through the time-series. Specifically, we tested the hypothesis that positive feedback loops, which have been predicted to destabilize biological communities, existed prior to a sudden community-compositional shift observed in the microbiome experiment. In addition, we examined the presence of microbial species that could be located at the source or sink positions within the directed graphs of metabolic flows. Overall, the preliminary application of community-level metabolic modeling provides a platform for understanding relationship between dynamics of interaction network architecture and ecosystem-level consequences.
Materials and methods
Time-series data of the microbial experiment
We used the 110-day time-series dataset of the microbial community experiment described elsewhere (Fujita et al., 2022, 2023). In the experiment, the source microbial community was sampled from a pond (“Shoubuike”) near Center for Ecological Research, Kyoto University, Otsu, Japan (34.974 °N, 135.966 °E). The source community was then introduced into the deep-well culture system of an oatmeal broth medium [0.5% (w/v) milled oatmeal (Nisshoku Oats; Nippon Food Manufacturer)] with eight replicates, kept at 23°C for 5 days (Fujita et al., 2023). After the five-day pre-incubation, 200 μl out of 1,000-μL culture medium was sampled from each well of the deep-well plate after mixing (pipetting) every 24 h for 110 days. In each sampling event, 200 μl of fresh medium was added to each well so that the total culture volume was kept constant. For the samples, amplicon sequencing of 16S rRNA was conducted as reported previously (Fujita et al., 2023). Based on the amplicon sequencing of the community compositional dynamics, we selected the replicate community that showed the largest community compositional changes within the time-series (Fujita et al., 2022): a rapid and substantial community compositional change occurred around Day 18 in the replicate community (Figure 1). The extracted DNA samples of the replicate community was subjected to a shotgun metagenomic sequencing analysis, which targeted 13 time points across the time-series (Day 1, 10, 20, 24, 30, 40, 50, 60, 70, 80, 90, 100, and 110; ca. 10 Gb per sample) (Fujita et al., 2022). The analysis described below was performed by compiling the metagenomic sequencing data (Fujita et al., 2022).
Figure 1. Time-series data of the community structure. (A) Community structure inferred with amplicon sequencing. Through the 110-day experiment, community compositions were monitored based on 16S rRNA sequencing. To quantify the speed and magnitude of community shifts through time, the “abruptness” index was calculated through the time-series (blue line). Specifically, an estimate of the abruptness index for time point t was obtained as the Bray-Curtis β-diversity between average community compositions from time points t – 4 to t and those from t + 1 to t + 5 (i.e., dissimilarity between 5-day time-windows). An abruptness score larger than 0.5 indicates that turnover of more than 50% of community compositions occurred between the time-windows. Reproduced from the amplicon sequencing data of a previous study on the microbiome system (Fujita et al., 2023). (B) Community structure inferred with metagenomic sequencing. The 32 metagenome-assembled genomes (MAGs) detected in the shotgun metagenomic sequencing are shown.
Shotgun metagenomic sequencing
The metagenomic sequencing data were processed as detailed previously (Fujita et al., 2022). Briefly, after adaptor trimming with Cutadapt (Martin, 2011) and quality filtering with Fastp0.21.0 (Chen et al., 2018) [in total, the number of output sequencing reads was 1002.49 M (160.08 Gb)], the sequences of each time-point sample were assembled using metaSPAdes 3.15.2 (Bankevich et al., 2012). Binning and quality assessing, were then performed with MetaWRAP 1.3.2 (Uritskiy et al., 2018) and CheckM 1.1.3 (Parks et al., 2015), respectively. The identity between metagenome-assembled genomes (MAGs) were calculated with FastANI 1.33 (Jain et al., 2018) and MAGs with >99% identity were grouped through the time-series. The read-coverage calculation was then conducted using CoverM 0.6.0 (Woodcroft, 2021). Taxonomic annotation and genome annotation were conducted, respectively, with GTDB-Tk 1.6 (Chaumeil et al., 2020; Parks et al., 2022) and Prokka 1.14.6 (Seemann, 2014) 1.14.6. The MAGs with >80% completeness and < 5% contamination were used in the analyzes below. The orthology numbers of Kyoto Encyclopedia of Genomes (KEGG) were retrieved for respective genes using GhostKOALA 2.2 (Kanehisa et al., 2016) and the completeness of metabolic pathways was estimated for each MAG using KEGG decoder 1.3 (Graham et al., 2018). In total, 32 MAGs belonging to 20 genera (16 families; 12 orders) were detected across the time-series (Supplementary Data 1; Fujita et al., 2022).
Metabolic modeling
To explore potential effects of facilitative interactions between microbes within the microbiome, we performed an analyzes of metabolite-exchange interaction networks based on the MAGs detailed above. For each MAG, we reconstructed a metabolic model based on the top-down carving approach of curated “universal models” (i.e., manually curated and simulation-ready metabolic models) (Machado et al., 2018) using CarveMe 1.5.0 (Machado et al., 2018). Potential metabolic interactions between microbial MAGs were then evaluated based on species coupling scores indicating dependency of target species in the presence of others as implemented in SMETANA 1.0.0 (Zelezniak et al., 2015). In each direction in each pair of species (MAGs), this score can vary between 0 (complete independence) and 1 (essentiality) (Zelezniak et al., 2015). In this approach, all potential exchanges of metabolites between species were mapped with the default parameters as implemented in SMETANA. To evaluate the reproducibility of the metabolic modeling, 10 SMETANA runs were performed with different random seed numbers.
To infer the community-scale magnitude of facilitative and competitive interactions, metabolic interaction potential and metabolic resource overlap scores were estimated, respectively, at each time point using SMETANA. The former was calculated as the difference between the minimal number of metabolites required for the growth of all species (MAGs) in a community without interspecific interactions and that with interactions (Zelezniak et al., 2015). The later was estimated as the mean proportion of shared nutritional requirements (metabolites) in pairs of species (MAGs) within a community (Zelezniak et al., 2015).
Network analysis
The architecture of the inferred network was visualized for each time point (day) based on the “backbone” layout algorithm (Nocaj et al., 2015) using R v.4.1.3 (R Core Team, 2020). In each network, metabolic dependencies between pairs of species (i.e., species coupling score) were indicated as arrows (i.e., directed graphs). The metabolites behind the inferred species interactions were shown for each day in the form of heatmaps.
The inferred metabolic interaction network of each time point was then analyzed based on the treeness, feedforwardness, and orderability (Corominas-Murtra et al., 2013). Treeness is a measure of pyramidal (top-down) network structure, in which small numbers of nodes at upper layers have outward links to many other nodes at lower layers. Feedforwardness is a measure of network-scale bias in the direction of links: a high feedforwardness value represents strong upstream-downstream structure within a network. Meanwhile, orderability represents the degree of the lack of feedback loops within directed graphs (networks). As the orderability index is defined as the proportion of nodes outside feedback network loops, it ranges from 0 (loop structure involving all nodes) to 1 (absence of loops).
To evaluate topological positions of respective microbial MAGs within the networks, influence centrality, which was a measure of the degree to which a focal node had influence on the others within a directed graph (Masuda et al., 2009), was calculated. Likewise, PageRank centrality, which was a measure of the degree to which a focal node had links from other nodes with many inward links (Page and Brin, 1998) was calculated as well using the igraph v.1.3.5 package of R. In R codes used for the calculation of influence, treeness, feedforwardness, and orderability are available from the GitHub repository1.
Results
Dynamics of metabolic interaction networks
Microbial MAGs belonging to different taxa were linked with each other within the network of potential facilitative interactions (Figure 2; Supplementary Data 2–3). In particular, microbes in the class Gammaproteobacteria were inferred to provide metabolites to microbes in other taxonomic groups. Likewise, Terracidiphilus bacteria (Acidobacteriae) had links of potential metabolite supply toward some bacteria belonging to Gammaproteobacteria and Alphaproteobacteria at some time points (Figure 2). The number of detectable nodes suddenly decreased between Days 20 and 30, entailing rapid decline of the inferred metabolic interaction networks (Figure 2). The microbial community then reached a quasi-stable state characterized by several bacteria in the genera Hydrotalea, Terracidiphilus, Mangrovibacter, and Rhizomicrobium (from Day 40 to Day 50; Figure 2). Among them, unidirectional facilitative effects from Mangrovibacter to other bacteria were inferred based on the metabolic modeling analysis (Figure 2). The number of detectable MAGs gradually increased from Day 60, resulting in the restoration of an entangled web of potential metabolic interactions on Day 110 (Figure 2).
Figure 2. Inferred network of metabolic interactions between microbes. Based on the shotgun metagenomic sequencing data, genome-scale metabolic modeling was conducted at each of the target time point. The results were used to infer potential flows metabolites between microbial MAGs. Positive effects inferred by metabolic modeling are shown with arrows connecting donor and recipient microbial MAGs. Darker colors of arrows indicate higher species coupling scores inferred in the metabolic modeling analysis. Edges constituting feedback loops are highlighted in red within the networks. Results based on a random seed number are shown: see Supplementary Data 2–3 for the metabolic modeling results of all the 10 SMETANA runs.
The chemicals (metabolites) potentially flowing through the inferred networks are shown in Figure 3. Iron (Fe2+ and Fe3+) and its siderophore (ferrypyoverdine) were among the chemicals most frequently flowing through the facilitative interaction networks. Copper, acetaldehyde, hydrogen sulfide, and phosphate were frequently reused within the networks (Figure 3). In the period in which the network structure was the simplest (Day 40 to 50; Figure 2), arginine and ethanol were provided from Mangrovibacter to other bacteria (Hydrotalea, Terracidiphilus, and Rhizomicrobium) within the microbiome (Figure 3B).
Figure 3. Metabolites presumably flowing through the facilitative network. (A) Absolute frequencies of inferred metabolite transfer between MAGs. The number of MAG combinations in which each metabolite was inferred to be provided was counted. Mean values across the 10 SMETANA runs with different random seed numbers are shown. The top-25 chemicals (metabolites) with the highest frequencies are presented. (B) Standardized frequencies of inferred transfer between pairs of species. The values in the panel A are standardized so that the total equals 1 at each time point.
The treeness, feedforwardness, and orderability of the network of the potential metabolic interactions varied considerably across the time-series (Figure 4). Until Day 20, the network structure was characterized by low treeness, low feedforwardness, and low to moderate orderability (Figure 4B). The facilitative interaction network then showed drastic architectural shift until Day 40 as characterized by the rapid increase of orderability (Figure 4B). This result suggests that the dynamics of the network architecture are characterized by the presence of positive feedback loops (as represented by low orderability) early in the time-series and that such feedback loops collapsed in the microbial community by Day 40 (Figure 3). Through the gradual restoration of network complexity after Day 60, the presence of feedback loops was inferred again on Day 110 (Figure 3) as indicated by lowered network orderability estimate on the day (Figure 4).
Figure 4. Network topology analysis. (A) Schema of network architectural properties. Treeness and feedforwardness represent pyramidal and upstream-downstream structures of directed graphs, respectively. Orderability represents lack of feedback loops within a network. Along the axis of orderability, the nodes and links included in feedback loops are highlighted in red. (B) Dynamics of network characteristics. Changes in network architectural properties are shown in terms of treeness, feedforwardness, and orderability. Networks with low “orderability,” by definition, contain loops of flow of metabolites, while those with maximum orderability (= 1) lack feedback loops. Mean and standard deviation across the 10 SMETANA runs with different random seed numbers are shown with diamonds and bars, respectively. The grey line represents the magnitude of community compositional changes shown in Figure 1A.
Potential keystone species
Within the metabolic interaction networks (Figure 2), some microbial MAGs belonging to the class Gammaproteobacteria were located at the “upper stream” of the network, showing high influence scores (Figure 5; Supplementary Data 1). In particular, a gammaproteobacterial MAG in the genus Mangrovibacter consistently showed the highest influence among the microbes detected at most time points (Figure 5). Meanwhile, microbes located at the sink positions within the inferred metabolic interaction networks (i.e., MAGs with high PageRank scores) represented diverse taxonomic groups (Figure 5). From Day 40 to 50, through which a small number of bacterial taxa represented the microbiome structure, simple source–sink relationship of potential metabolite flow was observed between Mangrovibacter and others (i.e., Hydrotalea, Terracidiphilus, and Rhizomicrobium; Figures 2, 5).
Figure 5. Potential keystone species/taxa within metabolic interaction networks. Within each directed graph of metabolic dependence network (Figure 4), influence-centrality (a measure of the degree to which a focal node has influence on the others within a directed graph) and PageRank-centrality (a measure of the degree to which a focal node has links from other nodes with many inward links) measures of network centrality was calculated for each microbe. Mangrovibacter tended to show high impacts (influence) on other bacteria within the metabolic interaction networks throughout the time-series. Mean and standard deviation across the 10 SMETANA runs with different random seed numbers are shown with symbols and bars, respectively.
Community-scale facilitation and competition
The community-scale magnitude of facilitative interactions, as inferred with the metabolic interaction potential index, was the highest on Day 10, and it dropped sharply until Day 24 (Figure 6A). It then remained low until Day 60 and increased slightly late in the time-series. The community-scale magnitude of competitive interactions, as inferred with the metabolic resource overlap index, was the highest from Day 1 to Day 20 (Figure 6B). The level of potential competition then decreased until Day 50, but it returned to relatively high levels after Day 90 (Figure 6B).
Figure 6. Community-scale inference of facilitation and competition. (A) Facilitation. Metabolic interaction potential was calculated as the difference between the minimal number of metabolites required for the growth of all MAGs in a community without interspecific interactions and that with interactions. (B) Competition. Metabolic resource overlap was estimated as the mean proportion of shared nutritional requirements (metabolites) in pairs of MAGs within a community (Zelezniak et al., 2015). Mean and standard deviation across the 10 SMETANA runs with different random seed numbers are shown with diamonds and bars, respectively.
Discussion
We here showed temporal shifts in the network architecture of facilitative interactions by compiling a metagenomic sequencing dataset of experimental microbiome dynamics. While ecosystem-level profiles of metabolic functions have been intensively investigated (Raes and Bork, 2008), shotgun metagenomic data also allow us to infer potential ecological interactions between species (Figures 3–6). Basic ecological theory predicts that facilitative interactions can destabilize biological communities (Allesina and Tang, 2012). However, recent theoretical studies suggest that such effects of facilitative interactions depend on network architecture of interactions (Bastolla et al., 2009; Thébault and Fontaine, 2010; Fontaine et al., 2011; Morton et al., 2022). Nestedness, for example, have been intensively investigated as a potential key property of facilitative interaction networks in terms of species coexistence (Bascompte et al., 2003; Thébault and Fontaine, 2010; Rohr et al., 2014). Meanwhile, most studies on facilitative ecological interactions have relied on the assumption that all links within a network are bidirectional (i.e., mutualistic). In this study, we explored ways for uncovering the structure of directed graphs of species interactions (Sugihara et al., 2012; Ushio et al., 2018; Delmas et al., 2019) based on a metagenomic analysis of potential metabolic interactions (Zelezniak et al., 2015). By showing that architecture of directed interaction networks can drastically change through time, we hope to spur discussion on potential roles of interaction network structure in biological community dynamics and stability.
Among the directed-graph indices examined in this study, orderability was of particular interest (Figure 4). It has been theoretically predicted that presence of positive feedback loops in facilitative interaction networks can destabilize ecological communities (Coyte et al., 2015; Levine et al., 2017). Specifically, such feedback structure of dependence may magnify cascades of population collapse once balance of population size among constituent species fluctuates within the feedback loops. In a previous study on the examined experimental microbiome, a high level of niche overlap among bacterial species was inferred to have promoted community compositional shifts (Fujita et al., 2022). In particular, niche overlap within the gammaproteobacterial or alphaproteobacterial sub-community (guild) presumably resulted in competitive exclusion of constituent microbial populations (Fujita et al., 2022). Such competition-driven decline of some gammaproteobacterial or alphaproteobacterial species may have triggered a cascade breakdown of species (Rezende et al., 2007) through the positive feedback loop observed in this study (Figure 2). In other words, once competitive exclusion occurs within an ecological guild, species depending on the metabolites of the declining guilds are expected to be negatively influenced by the reduced flow of metabolites through the facilitative interaction network. The relatively high risk of potential resource competition early in the microbiome dynamics was inferred as well in the metabolic modeling of this study (Figure 6), suggesting the importance of simultaneously considering facilitative and competitive interactions.
Treeness and feedforwardness of network architecture give additional important information about propagation of negative effects within networks. If a facilitative interaction network has hierarchical structure represented by high treeness and feedforwardness (Corominas-Murtra et al., 2013), placement of the ecological guilds from which fluctuations are initiated would influence subsequent ecological processes through the network. Specifically, fluctuations occurred in upstream positions may be propagated more rapidly throughout the network, while those derived from downstream positions would entail minimal impacts on the entire community. Albeit the potential roles of such hierarchical structure, the network architecture observed early in the experimental community (until Day 20) was represented by low treeness and feedforwardness (Figure 4B). Thus, the influence of hierarchical network structure on community collapse remains to be examined in future studies on networks with high treeness and feedforwardness.
In parallel with investigations on the entire network structure, directed-graphs reconstructed with metagenomic data provide us with insights into species occupying upstream/downstream positions within networks. Species located at upstream positions within a “supply chain” of metabolites may impose greater impacts on population dynamics of other species within the network than species at downstream positions. In our data, a bacterium in the genus Mangrovibacter continued to occupy upstream positions throughout the community dynamics as indicated by the analysis of network influence scores (Figure 5). Thus, although the Mangrovibacter bacterium was a minor component of the community (Figure 1), it might have disproportionately large impacts on the dynamics of the entire microbiome. This working hypothesis could only be tested by removing the Mangrovibacter bacterium from the experimental system. Nonetheless, such selective removal of specific bacterial species from microbiomes remains a challenge because the use of antibiotics often causes unexpected side-effects on non-target species (Cho et al., 2012; Francino, 2016; Langdon et al., 2016). Technical advances that allow selective removal of potential “keystone species” (Paine, 1966; Power et al., 1996) within microbiomes as well as in silico simulation of those highlighted species will extend the discussion on species’ roles within networks.
While the conclusions drawn from our analysis remain preliminary, we hope that it can showcase the use of metabolic modeling approaches to understand dynamics and consequences of facilitative interactions in ecological communities for future studies. Context-dependency of network architecture, for example, needs to be examined by comparing network dynamics among different experimental settings (e.g., different culture media or different temperature conditions) (Zelezniak et al., 2015; Magnúsdóttir et al., 2017). It is also important to evaluate to what extent network architectural properties inferred with the metabolic modeling approaches are consistent with those estimated with other informatics approaches. In this respect, comparison with recently developed methods for inferring species interactions based on time-series data is of particular interest (Deyle et al., 2016; Ushio et al., 2018; Suzuki et al., 2022). Furthermore, integrating information of facilitative interactions with that of competitive interactions is an essential step for examining how relative balance of multiple interaction types affect community stability (Bastolla et al., 2009; Fontaine et al., 2011; Mougi and Kondoh, 2012; Goldford et al., 2018). Interdisciplinary studies combining genomics and ecological theory will broaden our views on fundamental mechanism driving microbial community dynamics.
Data availability statement
The datasets and codes used in this study are available from the Github repository (https://github.com/hiroakif93/Facilitative-interaction-networks-in-experimental-microbial-community-dynamics-).
Author contributions
HT designed the work with HF. HF performed experiments. HF analyzed the data with HT. HF and HT wrote the paper with all the authors. All authors contributed to the article and approved the submitted version.
Funding
This work was financially supported by JST PRESTO (JPMJPR16Q6), JSPS Grant-in-Aid for Scientific Research (20 K20586), NEDO Moonshot Research and Development Program (JPNP18016), and JST FOREST (JPMJFR2048) to HT, Human Frontier Science Program (RGP0029/2019) to HT and EK, NWO-VICI (202.012) to EK, JSPS Grant-in-Aid for Scientific Research (20 K06820 and 20H03010) to KS, and JSPS Fellowship to HF and AC.
Acknowledgments
We thank Sayaka S. Suzuki and Keisuke Koba for support in the experiment. We are grateful to the reviewers for their constructive comments that improved the manuscript.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2023.1153952/full#supplementary-material
Footnotes
1. ^https://github.com/hiroakif93/Facilitative-interaction-networks-in-experimental-microbial-community-dynamics-
References
Allesina, S., and Tang, S. (2012). Stability criteria for complex ecosystems. Nature 483, 205–208. doi: 10.1038/nature10832
Bankevich, A., Nurk, S., Antipov, D., Gurevich, A. A., Dvorkin, M., Kulikov, A. S., et al. (2012). SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J. Comput. Biol. 19, 455–477. doi: 10.1089/cmb.2012.0021
Barberan, A., Bates, S. T., Casamayor, E. O., and Fierer, N. (2012). Using network analysis to explore co-occurrence patterns in soil microbial communities. ISME J. 6, 343–351. doi: 10.1038/ismej.2011.119
Bascompte, J., Jordano, P., Melia, C. J., and Olesen, J. M. (2003). The nested assembly of plant-animal mutualistic networks. Proc. Natl. Acad. Sci. U. S. A. 100, 9383–9387. doi: 10.1073/pnas.1633576100
Bastolla, U., Fortuna, M., Pascual-García, A., Ferrera, A., Luque, B., and Bascompte, J. (2009). The architecture of mutualistic networks minimizes competition and increases biodiversity. Nature 458, 1018–1020. doi: 10.1038/nature07950
Berry, D., and Widder, S. (2014). Deciphering microbial interactions and detecting keystone species with co-occurrence networks. Front. Microbiol. 5:219. doi: 10.3389/fmicb.2014.00219
Blanchet, F. G., Cazelles, K., and Gravel, D. (2020). Co-occurrence is not evidence of ecological interactions. Ecol. Lett. 23, 1050–1063. doi: 10.1111/ele.13525
CaraDonna, P. J., and Waser, N. M. (2020). Temporal flexibility in the structure of plant–pollinator interaction networks. Oikos 129, 1369–1380. doi: 10.1111/oik.07526
Carding, S., Verbeke, K., Vipond, D. T., Corfe, B. M., and Owen, L. J. (2015). Dysbiosis of the gut microbiota in disease. Microb. Ecol. Health Dis. 26:26191. doi: 10.3402/mehd.v26.26191
Chaumeil, P. A., Mussig, A. J., Hugenholtz, P., and Parks, D. H. (2020). GTDB-Tk: a toolkit to classify genomes with the genome taxonomy database. Bioinformatics 36, 1925–1927. doi: 10.1093/bioinformatics/btz848
Chen, S., Zhou, Y., Chen, Y., and Gu, J. (2018). Fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34, i884–i890. doi: 10.1093/bioinformatics/bty560
Cho, I., Yamanishi, S., Cox, L., Methé, B. A., Zavadil, J., Li, K., et al. (2012). Antibiotics in early life alter the murine colonic microbiome and adiposity. Nature 488, 621–626. doi: 10.1038/nature11400
Corominas-Murtra, B., Goñi, J., Solé, R., and Rodríguez-Caso, C. (2013). On the origins of hierarchy in complex networks. Proc. Natl. Acad. Sci. U. S. A. 110, 13316–13321. doi: 10.1073/pnas.1300832110
Coyte, K. Z., Schluter, J., and Foster, K. R. (2015). The ecology of the microbiome: networks, competition, and stability. Science 350, 663–666. doi: 10.1126/science.aad2602
Delmas, E., Besson, M., Brice, M. H., Burkle, L. A., Dalla Riva, G., Fortin, M. J., et al. (2019). Analysing ecological networks of species interactions. Biol. Rev. 94, 16–36. doi: 10.1111/brv.12433
Deyle, E. R., May, R. M., Munch, S. B., and Sugihara, G. (2016). Tracking and forecasting ecosystem interactions in real time. Proc. R. Soc. B Biol. Sci. 283:20152258. doi: 10.1098/rspb.2015.2258
Faust, K., and Raes, J. (2012). Microbial interactions: from networks to models. Nat. Rev. Microbiol. 10, 538–550. doi: 10.1038/nrmicro2832
Faust, K., Sathirapongsasuti, J. F., Izard, J., Segata, N., Gevers, D., Raes, J., et al. (2012). Microbial co-occurrence relationships in the human microbiome. PLoS Comput. Biol. 8:e1002606. doi: 10.1371/journal.pcbi.1002606
Fontaine, C., Guimarães, P. R. Jr., Kéfi, S., Loeuille, N., Memmott, J., van der Putten, W. H., et al. (2011). The ecological and evolutionary implications of merging different types of networks. Ecol. Lett. 14, 1170–1181. doi: 10.1111/j.1461-0248.2011.01688.x
Francino, M. P. (2016). Antibiotics and the human gut microbiome: Dysbioses and accumulation of resistances. Front. Microbiol. 6:1543. doi: 10.3389/fmicb.2015.01543
Friedman, J., and Alm, E. J. (2012). Inferring correlation networks from genomic survey data. PLoS Comput. Biol. 8, e1002687–e1002611. doi: 10.1371/journal.pcbi.1002687
Frioux, C., Singh, D., Korcsmaros, T., and Hildebrand, F. (2020). From bag-of-genes to bag-of-genomes: metabolic modelling of communities in the era of metagenome-assembled genomes. Comput. Struct. Biotechnol. J. 18, 1722–1734. doi: 10.1016/j.csbj.2020.06.028
Fujita, H., Ushio, M., Suzuki, K., Abe, M., Yamamichi, M., Iwayama, K., et al. (2023). Alternative stable states, nonlinear behavior, and predictability of microbiome dynamics. Microbiome 2023:505041. doi: 10.1101/2022.08.23.505041
Fujita, H., Ushio, M., Suzuki, K., Abe, M. S., Yamamichi, M., Okazaki, Y., et al. (2022). Metagenomic analysis of ecological niche overlap and community collapse in microbiome dynamics. bioRxiv 2022:524457. doi: 10.1101/2023.01.17.524457
Goldford, J. E., Lu, N., Bajić, D., Estrela, S., Tikhonov, M., Sanchez-Gorostiaga, A., et al. (2018). Emergent simplicity in microbial community assembly. Science 361, 469–474. doi: 10.1126/science.aat1168
Graham, E. D., Heidelberg, J. F., and Tully, B. J. (2018). Potential for primary productivity in a globally-distributed bacterial phototroph. ISME J. 12, 1861–1866. doi: 10.1038/s41396-018-0091-3
Gralka, M., Szabo, R., Stocker, R., and Cordero, O. X. (2020). Trophic interactions and the drivers of microbial community assembly. Curr. Biol. 30, R1176–R1188. doi: 10.1016/j.cub.2020.08.007
Hassani, M. A., Durán, P., and Hacquard, S. (2018). Microbial interactions within the plant holobiont. Microbiome 6:58. doi: 10.1186/s40168-018-0445-0
Ives, A. R., and Carpenter, S. R. (2007). Stability and diversity of ecosystems. Science 317, 58–62. doi: 10.1126/science.1133258
Jain, C., Rodriguez-R, L. M., Phillippy, A. M., Konstantinidis, K. T., and Aluru, S. (2018). High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. Nat. Commun. 9:5114. doi: 10.1038/s41467-018-07641-9
Kanehisa, M., Sato, Y., and Morishima, K. (2016). BlastKOALA and GhostKOALA: KEGG tools for functional characterization of genome and metagenome sequences. J. Mol. Biol. 428, 726–731. doi: 10.1016/j.jmb.2015.11.006
Klitgord, N., and Segrè, D. (2010). Environments that induce synthetic microbial ecosystems. PLoS Comput. Biol. 6:e1001002. doi: 10.1371/journal.pcbi.1001002
Kriss, M., Hazleton, K. Z., Nusbacher, N. M., Martin, C. G., and Lozupone, C. A. (2018). Low diversity gut microbiota dysbiosis: drivers, functional implications and recovery. Curr. Opin. Microbiol. 44, 34–40. doi: 10.1016/j.mib.2018.07.003
Kurtz, Z. D., Bonneau, R., and Müller, C. L. (2019). Disentangling microbial associations from hidden environmental and technical factors via latent graphical models. bioRxiv 2019:885889. doi: 10.1101/2019.12.21.885889
Kurtz, Z. D., Müller, C. L., Miraldi, E. R., Littman, D. R., Blaser, M. J., and Bonneau, R. A. (2015). Sparse and compositionally robust inference of microbial ecological networks. PLoS Comput. Biol. 11:e1004226. doi: 10.1371/journal.pcbi.1004226
Langdon, A., Crook, N., and Dantas, G. (2016). The effects of antibiotics on the microbiome throughout development and alternative approaches for therapeutic modulation. Genome Med. 8:39. doi: 10.1186/s13073-016-0294-z
Lever, J. J., van Nes, E. H., Scheffer, M., and Bascompte, J. (2014). The sudden collapse of pollinator communities. Ecol. Lett. 17, 350–359. doi: 10.1111/ele.12236
Levine, J. M., Bascompte, J., Adler, P. B., and Allesina, S. (2017). Beyond pairwise mechanisms of species coexistence in complex communities. Nature 546, 56–64. doi: 10.1038/nature22898
Levy, R., and Borenstein, E. (2013). Metabolic modeling of species interaction in the human microbiome elucidates community-level assembly rules. Proc. Natl. Acad. Sci. U. S. A. 110, 12804–12809. doi: 10.1073/pnas.1300926110
Machado, D., Andrejev, S., Tramontano, M., and Patil, K. R. (2018). Fast automated reconstruction of genome-scale metabolic models for microbial species and communities. Nucleic Acids Res. 46, 7542–7553. doi: 10.1093/nar/gky537
Magnúsdóttir, S., Heinken, A., Kutt, L., Ravcheev, D. A., Bauer, E., Noronha, A., et al. (2017). Generation of genome-scale metabolic reconstructions for 773 members of the human gut microbiota. Nat. Biotechnol. 35, 81–89. doi: 10.1038/nbt.3703
Martin, M. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J 17:200. doi: 10.14806/ej.17.1.200
Masuda, N., Kawamura, Y., and Kori, H. (2009). Analysis of relative influence of nodes in directed networks. Phys. Rev. E Stat. Nonlinear Soft Matter Phys. 80, 1–10. doi: 10.1103/PhysRevE.80.046114
May, R. M. (1972). Will a large complex system be stable? Nature 238, 413–414. doi: 10.1038/238413a0
Morton, D. N., Keyes, A., Barner, A. K., and Dee, L. E. (2022). Merging theory and experiments to predict and understand coextinctions. Trends Ecol. Evol. 37, 886–898. doi: 10.1016/j.tree.2022.06.004
Mougi, A., and Kondoh, M. (2012). Diversity of interaction types and ecological community stability. Science 337, 349–351. doi: 10.1126/science.1220529
Nocaj, A., Ortmann, M., and Brandes, U. (2015). Untangling the hairballs of multi-centered, small-world online social media networks. J. Graph Algorithms Appl. 19, 595–618. doi: 10.7155/jgaa.00370
Olesen, J. M., Bascompte, J., Dupont, Y. L., and Jordano, P. (2007). The modularity of pollination networks. Proc. Natl. Acad. Sci. U. S. A. 104, 19891–19896. doi: 10.1073/pnas.0706375104
Page, L., and Brin, S. (1998). The anatomy of a large-scale hypertextual web search engine. Comput. Netw. 30, 107–117. doi: 10.1016/s0169-7552(98)00110-x
Paine, R. T. (1966). Food web complexity and species diversity. Am. Nat. 100, 65–75. doi: 10.1086/282400
Parks, D. H., Chuvochina, M., Rinke, C., Mussig, A. J., Chaumeil, P.-A., and Hugenholtz, P. (2022). GTDB: an ongoing census of bacterial and archaeal diversity through a phylogenetically consistent, rank normalized and complete genome-based taxonomy. Nucleic Acids Res. 50, D785–D794. doi: 10.1093/nar/gkab776
Parks, D. H., Imelfort, M., Skennerton, C. T., Hugenholtz, P., and Tyson, G. W. (2015). CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 25, 1043–1055. doi: 10.1101/gr.186072.114
Power, M. E., Tilman, D., Estes, J. A., Menge, B. A., Bond, W. J., Mills, L. S., et al. (1996). Challenges in the quest for keystones. Bioscience 46, 609–620. doi: 10.2307/1312990
R Core Team. (2020). R: A Language and Environment for Statistical Computing. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
Raes, J., and Bork, P. (2008). Molecular eco-systems biology: towards an understanding of community function. Nat. Rev. Microbiol. 6, 693–699. doi: 10.1038/nrmicro1935
Ratzke, C., Barrere, J., and Gore, J. (2020). Strength of species interactions determines biodiversity and stability in microbial communities. Nat. Ecol. Evol. 4, 376–383. doi: 10.1038/s41559-020-1099-4
Rezende, E. L., Lavabre, J. E., Guimarães, P. R., Jordano, P., and Bascompte, J. (2007). Non-random coextinctions in phylogenetically structured mutualistic networks. Nature 448, 925–928. doi: 10.1038/nature05956
Rohr, R. P., Saavedra, S., and Bascompte, J. (2014). On the structural stability of mutualistic systems. Science 345:1253497. doi: 10.1126/science.1253497
Scheffer, M., and Carpenter, S. R. (2003). Catastrophic regime shifts in ecosystems: linking theory to observation. Trends Ecol. Evol. 18, 648–656. doi: 10.1016/j.tree.2003.09.002
Scheffer, M., Carpenter, S., Foley, J., Folke, C., and Walker, B. (2001). Catastrophic shifts in ecosystems. Nature 413, 591–596. doi: 10.1038/35098000
Seemann, T. (2014). Prokka: rapid prokaryotic genome annotation. Bioinformatics 30, 2068–2069. doi: 10.1093/bioinformatics/btu153
Stolyar, S., van Dien, S., Hillesland, K. L., Pinel, N., Lie, T. J., Leigh, J. A., et al. (2007). Metabolic modeling of a mutualistic microbial community. Mol. Syst. Biol. 3:92. doi: 10.1038/msb4100131
Sugihara, G., May, R., Ye, H., Hsieh, C., Deyle, E., Fogarty, M., et al. (2012). Detecting causality in complex ecosystems. Science 338, 496–500. doi: 10.1126/science.1227079
Sung, J., Kim, S., Cabatbat, J. J. T., Jang, S., Jin, Y. S., Jung, G. Y., et al. (2017). Global metabolic interaction network of the human gut microbiota for context-specific community-scale analysis. Nat. Commun. 8, 1–12. doi: 10.1038/ncomms15393
Suzuki, K., Matsuzaki, S., and Masuya, H. (2022). Decomposing predictability to identify dominant causal drivers in complex ecosystems. Proc. Natl. Acad. Sci. U. S. A. 119:e2204405119. doi: 10.1073/pnas.2204405119
Thébault, E., and Fontaine, C. (2010). Stability of ecological communities and the architecture of mutualistic and trophic networks. Science 329, 853–856. doi: 10.1126/science.1188321
Toju, H., Yamamichi, M. Jr., Guimarães, P. R. Jr., Olesen, J. M., Mougi, A., Yoshida, T., et al. (2017). Species-rich networks and eco-evolutionary synthesis at the metacommunity level. Nat. Ecol. Evol. 1:0024. doi: 10.1038/s41559-016-0024
Uritskiy, G., DiRuggiero, J., and Taylor, J. (2018). MetaWRAP–a flexible pipeline for genome-resolved metagenomic data analysis. Microbiome 6:158. doi: 10.1186/s40168-018-0541-1
Ushio, M., Hsieh, C. H., Masuda, R., Deyle, E. R., Ye, H., Chang, C. W., et al. (2018). Fluctuating interaction network and time-varying stability of a natural fish community. Nature 554, 360–363. doi: 10.1038/nature25504
Warton, D. I., Blanchet, F. G., O’Hara, R. B., Ovaskainen, O., Taskinen, S., Walker, S. C., et al. (2015). So many variables: joint modeling in community ecology. Trends Ecol. Evol. 30, 766–779. doi: 10.1016/j.tree.2015.09.007
Woodcroft, B.. (2021). CoverM: Program. Available at: https://github.com/wwood/CoverM
Zelezniak, A., Andrejev, S., Ponomarova, O., Mende, D. R., Bork, P., and Patil, K. R. (2015). Metabolic dependencies drive species co-occurrence in diverse microbial communities. Proc. Natl. Acad. Sci. U. S. A. 112, 6449–6454. doi: 10.1073/pnas.1421834112
Keywords: community stability, dysbiosis, ecosystem functions, microbe-microbe interactions, metabolic modeling, microbial functions, mutualism, species interactions
Citation: Fujita H, Ushio M, Suzuki K, Abe MS, Yamamichi M, Okazaki Y, Canarini A, Hayashi I, Fukushima K, Fukuda S, Kiers ET and Toju H (2023) Facilitative interaction networks in experimental microbial community dynamics. Front. Microbiol. 14:1153952. doi: 10.3389/fmicb.2023.1153952
Edited by:
Babak Momeni, Boston College, United StatesReviewed by:
Pandeng Wang, Sun Yat-sen University, ChinaErida Gjini, University of Lisbon, Portugal
Copyright © 2023 Fujita, Ushio, Suzuki, Abe, Yamamichi, Okazaki, Canarini, Hayashi, Fukushima, Fukuda, Kiers and Toju. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Hiroaki Fujita, fujita.h@ecology.kyoto-u.ac.jp; Hirokazu Toju, toju.hirokazu.4c@kyoto-u.ac.jp