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.

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.

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-).

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.

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.

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.

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.