Modeling the Role of the Microbiome in Evolution

There is undeniable evidence showing that bacteria have strongly influenced the evolution and biological functions of multicellular organisms. It has been hypothesized that many host-microbial interactions have emerged so as to increase the adaptive fitness of the holobiont (the host plus its microbiota). Although this association has been corroborated for many specific cases, general mechanisms explaining the role of the microbiota in the evolution of the host are yet to be understood. Here we present an evolutionary model in which a network representing the host adapts in order to perform a predefined function. During its adaptation, the host network (HN) can interact with other networks representing its microbiota. We show that this interaction greatly accelerates and improves the adaptability of the HN without decreasing the adaptation of the microbial networks. Furthermore, the adaptation of the HN to perform several functions is possible only when it interacts with many different bacterial networks in a specialized way (each bacterial network participating in the adaptation of one function). Disrupting these interactions often leads to non-adaptive states, reminiscent of dysbiosis, where none of the networks the holobiont consists of can perform their respective functions. By considering the holobiont as a unit of selection and focusing on the adaptation of the host to predefined but arbitrary functions, our model predicts the need for specialized diversity in the microbiota. This structural and dynamical complexity in the holobiont facilitates its adaptation, whereas a homogeneous (non-specialized) microbiota is inconsequential or even detrimental to the holobiont's evolution. To our knowledge, this is the first model in which symbiotic interactions, diversity, specialization and dysbiosis in an ecosystem emerge as a result of coevolution. It also helps us understand the emergence of complex organisms, as they adapt more easily to perform multiple tasks than non-complex ones.

Introduction

It has been firmly established during the last decade that the microbiota of a multicellular host strongly influences its evolution and adaptation (Ley et al., 2008; Zilber-Rosenberg and Rosenberg, 2008; Rosenberg et al., 2010; Brucker and Bordenstein, 2013; Andrew et al., 2016; Rosenberg and Zilber-Rosenberg, 2016; Sharpton, 2018). In turn, the host's ability to interact with other organisms and modify its environment to its advantage can guide the composition of its microbiota (Mai, 2004; Spor et al., 2011; Yatsunenko et al., 2012; Moeller et al., 2014). For instance, the human microbiota plays an important role in many fundamental physiological functions, such as the development of the immune system (Hooper et al., 2012), degradation of fiber and metabolization of fats and carbohydrates (Krajmalnik-Brown et al., 2012), regulation of bone density (McCabe et al., 2015), metabolization of drugs (Wilson and Nicholson, 2017) and control of infections by pernicious bacteria like Clostridium difficile (Rupnik et al., 2009; Van Nood et al., 2013; Seekatz and Young, 2014). A state of imbalance in the human microbiota, known as dysbiosis, has been correlated with diseases (Cho and Blaser, 2012) such as obesity (Ley et al., 2006b; Ley, 2010), inflammatory bowel disease (Morgan et al., 2012; Halfvarson et al., 2017), cancer (Farrell et al., 2011; Zackular et al., 2013; Francescone et al., 2014; Sears and Garrett, 2014; Contreras et al., 2016; Yang et al., 2017) and even neurological disorders as schizophrenia and autism (Gonzalez et al., 2011; Rogers et al., 2016). The system consisting of the host and its microbiota, known as holobiont, exhibits the unequivocal existence of symbiotic relationships between microbes and multicellular organisms (Theis et al., 2016). Given these complex interdependencies, the holobiont has been proposed to function as a single evolutionary unit (Gilbert et al., 2012; Gordon et al., 2013; Guerrero et al., 2013; Bordenstein and Theis, 2015; Theis et al., 2016; Roughgarden et al., 2018). This is because environmental changes may impose selective pressures on the host which in turn will affect its microbiota (for a recent review see Roughgarden et al., 2018). In the light of these findings, it has been suggested that evolutionary theories have to be either reformulated or expanded in order to account for the adaptability of the holobiont as an evolutionary unit (Laland et al., 2014; Ereshefsky and Pedroso, 2015; Van Opstal and Bordenstein, 2015; Sandoval-Motta et al., 2017). Whether the holobiont is or is not an evolutionary unit is still a matter of debate (Moran and Sloan, 2015; Douglas and Werren, 2016; Doolittle and Booth, 2017; Doolittle and Inkpen, 2018). However, here we show that selective pressures applied to the host and its associated microbes taken as whole, can help us explain how symbiotic relationships in holobionts arose and are currently maintained.

Examples of the influence that microorganisms have had on the adaptation of their hosts range from cases in which microbes help the host to perform specific non-essential functions, to cases in which microbes have completely substituted essential functions of the host (Sagan, 1967; Zilber-Rosenberg and Rosenberg, 2008; Queller and Strassmann, 2016; Roughgarden et al., 2018). Nevertheless, the specific mechanisms by which this influence is carried on are not yet known. Particularly, what are the general benefits that the microbiota provides to the host during its evolution is still an open question. A possible answer to it is that the adaptation time of the host to face new environmental challenges is considerably reduced due to the great diversity and plasticity of its microbiota (Zilber-Rosenberg and Rosenberg, 2008; Rosenberg and Zilber-Rosenberg, 2016). This hypothesis assumes that the emergence of strong symbiotic relationships between the host and its microbiota occurs at the genetic and metabolic levels, for only in this way changes occurring in the microbiota can rapidly propagate to the host's metabolism and affect its adaptability. Indeed, recent evidence shows that the microbiota can regulate metabolic pathways and gene expression patterns of its host, and due to this interaction the host can properly perform cell differentiation, tissue formation, nutrition and other important functions (Hooper et al., 2001; Rawls et al., 2004; Bates et al., 2006; Shin et al., 2011; Nicholson et al., 2012; Camp et al., 2014).

It has been proposed that natural selection operating at the Host-level promotes stable and redundant microbial societies, whereas selection operating at the microbial level promotes functional specialization of their component species (Ley et al., 2006a). Despite all the knowledge we have now on human associated microbial communities, we still do not fully understand the evolutionary forces behind the diversity observed in our microbiota. On the one hand, the most abundant ecological relationship between microbial species is competition (Foster and Bell, 2012; Coyte et al., 2015; Moran and Sloan, 2015; Douglas and Werren, 2016), which often leads to uniform microbial communities where just a few species dominate the whole environment. On the other hand, it has been shown that purely mutualistic interactions lead to unstable communities as their diversity increases. These observations are at odds with the great diversity and stability observed in the microbiota of most plants and animals. Maintaining this diversity is fundamental for the survivability of the host, as it is known that a loss in the microbiota's diversity may produce severe dysbiosis that can result in host diseases or even death (Blaser and Falkow, 2009; Turnbaugh and Stintzi, 2011; Cho and Blaser, 2012; Fernández et al., 2013; Lloyd-Price et al., 2016; Blaser, 2017). An observation that circumvents this caveat is that multicellular organisms have developed different mechanisms to maintain the equilibrium between its diverse microbial communities. These mechanisms tend to compartmentalize the microbes in separate niches while reducing the interactions between microbes in the same niche (Grice and Segre, 2011; Donaldson et al., 2015; Deines et al., 2017; Tropini et al., 2017; Roughgarden et al., 2018). Understanding why microbial diversity is necessary for the evolution and adaptation of the host, and why disease arises when such diversity is lost, is a fundamental question with still no definitive answer.

To address these questions, we adopt the hypothesis that the holobiont constitutes a unit of selection in evolution and explore its consequences. We present an evolutionary population model in which the biological functions of organisms are encoded in the Boolean dynamics of regulatory networks. In our model, a host is represented as a Boolean network that needs to evolve in order to adequately perform a predefined task (or function). This is equivalent to the host acquiring a new phenotype in order to cope with a new environmental challenge. A population of such host networks is evolved in a way that each host network can establish regulatory interactions with a set of microbial species, each one represented also by a network. The main difference between the microbial and host networks is that due to the faster duplication rates of microbes, the generation of mutants is at least one order of magnitude larger in the microbial networks than in the host. Mutants, as explained in detail in the Materials and Methods section (M&M), are simulated by rewiring the connections of their network, or by altering their functionality. As we are dealing with evolutionary dynamics, it is important to mention that we will only consider host-microbe interactions that can be transmitted across generations. This is based on the fact that in many species, parents directly transmit their microbiota to their offspring or they construct environments with a stable microbial composition that bias the microbial composition of their progeny (Rosenberg et al., 2010; Fitzpatrick, 2014). Another important assumption in our model is the persistence across generations of the host-microbe interactions developed throughout the evolution of the holobiont, which is a necessary condition for natural selection to operate (Doolittle and Booth, 2017). Additionally, we implement the “It's the Song not the Singer” approach proposed by Doolittle (Taxis et al., 2015; Doolittle and Booth, 2017; Doolittle and Inkpen, 2018) by preserving throughout the evolution of the holobiont, those regulatory connections that contribute to the host's adaptation to perform a predefined but otherwise arbitrary dynamical function. The conservation of the dynamical function across generations occurs regardless of the specific host-microbe network interactions that are contributing to the adaptation process.

Our evolutionary model is based on the Boolean network model introduced by S. Kauffman (presented in the M&M section) to describe gene regulation and cell differentiation processes (Kauffman, 1969a,b). During the last 20 years, it has been shown that this model adequately captures the main aspects of gene regulation dynamics. For instance, Boolean networks are able to reproduce gene expression patterns and metabolic pathways experimentally observed in organisms such as Arabidopsis thaliana (Espinosa-Soto et al., 2004), Drosophila melanogaster (Albert and Othmer, 2003), yeast (Li et al., 2004; Davidich and Bornholdt, 2008), human epithelial cells (Huang et al., 2005) and murine blood progenitor cells (Hameya et al., 2017) among others. Additionally, Huang et al. experimentally showed that the dynamical attractors of a Boolean network correspond to different cell types or cell fates (Huang et al., 2005). Because of this evidence, we use Boolean networks to represent the gene regulation networks of both the hosts and their microbes. Since we are interested in general principles about the emergence of symbiotic interactions, we use random networks instead of carefully constructed ones corresponding to specific organisms. Although the gene regulatory network of an organism greatly determines its phenotype (Davidson and Levine, 2008; Oliveri et al., 2008), it is known that several functions depend more on the general structure of the network than on the specific genes involved (Wagner, 2007). Therefore, using random Boolean networks in our population model has the advantage of determining the capability of the network to acquire new functions throughout its evolution regardless of its detailed composition (Davidson, 2010). This function-centered approach is consistent with the fact that a core microbiome is more likely to be identified based on functionalities rather than on the particular phylogenetic details of its species (Consortium, 2012; Taxis et al., 2015; Doolittle and Booth, 2017; Doolittle and Inkpen, 2018). We describe in detail the Boolean network model in M&M section.

Simulations of this evolutionary model show that the adaptation of the host network is greatly enhanced when it interacts with the microbial networks, which are the ones that absorb most of the mutations without changing their own adaptation. Additionally, the host network can improve its adaptation to perform multiple functions only if the set of microbial networks is partitioned into specialized subsets (niches), each one participating in the host's adaptation to a small number of functions. This specialization provides the holobiont with a structural and dynamical complexity that facilitates its evolution, whereas non-specialized microbiota is shown to be either inconsequential or detrimental to the holobiont's adaptation. Once the holobiont is adapted, the disconnection of one or more of these specialized niches leads a global incompetence to perform the required set of imposed tasks. This is reminiscent of the dysbiosis observed in real organisms when their microbiota's diversity is reduced. To our knowledge, this is the first model in which symbiotic interactions, diversity, specialization and dysbiosis in an ecosystem emerge as a result of coevolution. It also helps us understand the emergence of complex organisms, as they adapt more easily than unstructured ones.

Model and Results

Task Assignment

Following the work by Stern (1999), in order to define a task for the Boolean network we start by arbitrarily selecting a subset of N_s nodes that we call signal nodes, {σ_s1, σ_s2, …, σ_sNs}, from which we extract the output signal R(t) defined as (see Figure 1A)

$\begin{array}{ll}R\left(t\right)={\displaystyle \sum _{i=1}^{N_s}{\sigma }_{s_i}}\left(t\right).& (1)\end{array}$

Assigning a task to the Boolean network consists in requiring that the output signal R(t) approximates as much as possible a predefined target function (or task) F(t) (see Figure 1B). In our model F(t) is an arbitrary function such that 0 < F(t) < N_s for 1 ≤ t ≤ t_m, where t_m = 15 is the number of time steps of the assigned task. We set t_m = 15 because this is the average number of time steps it takes for the network to stabilize its dynamics (see Figure S1). In biological terms, the task F(t) would represent an expression pattern some genes must acquire in order for the organism to efficiently respond to a particular environmental challenge (like yeast responding to a heat shock). Since the networks are randomly constructed, it is expected that initially none of them have this response (their output signal R(t) and the task F(t) are usually different at the start of the simulation, see Figure 1B). Therefore, it is necessary to evolve the networks so that R(t) approaches F(t) as much as possible, as in Figure 1C. It is only through a series of mutations and adaptations that the phenotype R(t) will approach F(t) in some individuals, and then be transmitted to their offspring.

FIGURE 1

Figure 1. Boolean Network Task. (A) Network with N = 12 nodes and N_s = 4 signal nodes (represented as green squares), which generate the output signal R(t) that has to converge to the task F(t). (B) Initially, the untrained network produces an output signal R(t) (green curve with squares) very different from the target function F(t) (blue curve with circles). (C) After the evolutionary process, the network is well adapted to perform the task. In this example, only one point of the output signal R(t) acquires a wrong value in the time interval 1 ≤ t ≤ 15 in which the task F(t) is evaluated. (D) Average population error ${\overline{\xi }}^H$ as a function of the number of generations g, for a population with P = 100 networks, each having N = 50 nodes and N_s = 12 signal nodes. Note that ${\overline{\xi }}^H$ decreases and crosses the adaptation threshold δ_A = 1 approximately at generation g = 350, after which most of the networks in the population become well adapted to the task F(t).

Throughout this work we use networks with N = 50 nodes, average connectivity K = 2 and N_s = 12 signal nodes (except in some figures where smaller networks are presented for illustrative purposes). The reason for this choice of parameters is the following. It has been observed that genetic networks of several real organisms are structured in functional modules, each one consisting of a few dozen genes or nodes (Resendis-Antonio et al., 2012). For instance, adaptive resistance to antibiotics in Escherichia coli is mediated by the MarA-AcrAB-TolC system which, when activated, produces efflux pumps that pump toxic molecules in the intracellular fluid out of the cell, keeping the internal antibiotic concentration below lethal levels. Activation of this system is controlled by a regulatory network consisting of about 15 nodes (Motta et al., 2015). Analogously, the cAMP-dependent protein kinase regulatory network (PKA-RN), which regulates (among other things) the stress response in Saccharomyces cerevisiae, consists of 15 nodes (Pérez-Landero et al., 2015). There are many more examples showing that specific cellular functions (such a the response to a given environmental challenge) are controlled by network modules composed of a few dozen nodes (Guo et al., 2016; Ma et al., 2017). Since in this work we are not considering any specific organism, we will assume in a generic way that the task F(t) that the network has to acquire is encoded in N_s = 12 nodes, which in turn are embedded in a module of N = 50 nodes.

Finally, we perform all our simulations using Kauffman networks with connectivity K = 2 for two main reasons. First, these networks are trained faster than networks with connectivity smaller or larger than K = 2 (see Figure S2). A second, more fundamental reason is that networks with K = 2 exhibit critical dynamics, which means that their dynamical behavior is at the brink of a phase transition between order and chaos (Derrida and Pomeau, 1986; Aldana, 2003). Dynamical criticality confers the system interesting properties such as evolvability (i.e., the coexistence of robustness and adaptability) (Aldana et al., 2007; Torres-Sosa et al., 2012), faster information storage, processing and transfer (Langton, 1990; Nykter et al., 2008), and collective response to external stimuli without saturation (Kinouchi and Copelli, 2006), (or shorter training times, as in our case, see Figure S2). There is solid evidence indicating that gene regulatory networks of real organisms are dynamically critical or close to criticality (Shmulevich et al., 2005; Serra et al., 2007; Balleza et al., 2008; Daniels et al., 2018). Therefore, by choosing K = 2 we are working with a representative ensemble of networks that have an important dynamical property observed in real organisms.

Host Network Evolution

We consider a population of P = 100 networks, represented as {H₁, H₂, …, H_P}, which have to perform the same task F(t). We will refer to these networks as the host networks (HNs). At the start of the simulation all the HNs are identical replicas of one randomly constructed network. To make the output signal of the HNs approach the task F(t) we implement a traditional evolutionary algorithm in which the networks are mutated, selected and replicated. Variability in the population is implemented by mutating the HNs with a mutation rate μ_H = 0.001 per node per network per generation. Once a node σ_n of a given network H_i has been chosen for mutation, we perform any of the following changes with equal probability: (i) Randomly rewire one of the input or output connections of σ_n. (ii) Add a new input (or output) connection to σ_n from (or to) a randomly chosen node in the network. (iii) Remove one input or output connection of σ_n. (iv) Change one of the entries of the logical function f_n associated to σ_n.

The mutations described above can make each network H_i get closer to the task F(t) or get away from it. To measure the adaptation of the HNs to the task we denote as R_i(t) the output signal of the network H_i and define its adaptation error ${\xi }_i^H$ as

$\begin{array}{ll}{\xi }_i^H={\left(t_m\right)}^{-1}{\displaystyle \sum _{t=1}^{t_m}{\left(R_i(t)-F(t)\right)}^2}.& (2)\end{array}$

Clearly, if ${\xi }_i^H=0$ then the network H_i is perfectly trained (adapted) to perform the task F(t), whereas large values of ${\xi }_i^H$ indicate a poor adaptation. Therefore, when a mutation occurs such that ${\xi }_i^H$ decreases, the adaptation of H_i increases and viceversa. We will say that the network H_i is well adapted to its task when ${\xi }_i^H\le {\delta }_A$, where δ_A is the adaptation threshold. We set δ_A = 1, which means that at most one node out of the N_s signal nodes is allowed to deviate one unit from the correct value at every time step during the interval 1 ≤ t ≤ t_m over which the task F(t) is evaluated (see Figure 1C). The average population error, defined as ${\overline{\xi }}^H=\frac{1}{P}{\displaystyle {\sum }_{i=1}^P{\xi }_i^H}$, measures the adaptation of the entire population to the task F(t).

In each generation we mutate the HNs in the population with the mutation rate μ_H. Then, we choose the 10 best networks (those whose errors ${\xi }_i^H$ have the lowest values) to get through the next generation while the other 90 networks are removed from the simulation. These 10 networks are replicated by making 9 copies of each one in order to restore the population to its original size P = 100. This evolutionary process is repeated until the population crosses the adaptation threshold. A “generation” consists in a full round of mutation, selection and replication processes. Figure 1D shows that the average population error ${\overline{\xi }}^H$ decreases throughout generations. This is expected since at each generation we select the networks that minimize the error. From Figure 1D we see that it takes about 350 generations for the average population error of hosts networks to cross δ_A and become well adapted to the task (see also Movie S1). We have performed simulations with smaller values of the adaptation threshold: δ_A = 0.5 and δ_A = 0.2, and the results are qualitatively the same. The only difference is that the smaller the value of δ_A, the longer the computing time for the average population error ${\overline{\xi }}^H$ to cross this threshold (see Figure S3). The results presented in Figure 1D correspond to a population of HNs evolving by themselves, i.e., without interacting among them or with other networks. We refer to this case as the control case.

Interaction With the Microbiota: Holobiont Evolution

To model the interaction between the host organism and its microbiota we allow the training of each host network H to be assisted by a set of P_M other networks, ${B}=\left\{M_1,M_2,\dots ,M_{P_M}\right\}$, each one representing a microbial network (MN). We will refer to the set ${B}$ as the microbiota, and to the set ${L}=\left\{H,{B}\right\}=\left\{H,M_1,\dots ,M_{P_M}\right\}$, as the holobiont.

Each microbial network $M_j\in {B}$ also has to perform a predefined task $F_j^M(t)$, which is an arbitrary function constructed in the same way as the host-network task F(t). The microbial tasks $F_1^M(t),\dots ,F_{P_M}^M(t)$ are different from each other and from F(t). Before the training of H begins, each $M_j\in {B}$ is previously trained to be well adapted to its own task $F_j^M(t)$. This means that all the microbial errors ${\xi }_j^M$ satisfy, from the very beginning, the well-adapted condition ${\xi }_j^M<{\delta }_A$ (as in Figure 1D; the microbial error ${\xi }_j^M$ is defined similarly as in Equation (2); see the M&M section for the precise definition). Thus, at generation g = 0 the holobiont consists of the untrained host network H and a set of well adapted MNs. The evolution of the holobiont then proceeds with the adaptation of H to its task and allowing it to interact (as described below) with MNs that already have their own interests. The rationale behind this initial setup is twofold. First, allowing the training of H to be assisted by well-adapted MNs captures the fact that at any moment during its evolution, the host organism can recruit from the environment microbial populations already adapted to their environments and able to carry out some functions by their own. Second, we want to determine whether evolutionary conflict emerges between the host and microbial networks when the holobiont evolves as a unit of adaptation, as has been pointed out in Moran and Sloan (2015) and Douglas and Werren (2016). Such a conflict would be apparent in our simulations if a reduction of the host-network error ξ^H occurs with a simultaneous increase in the average microbial error ${\overline{\xi }}^M$, or viceversa.

The interaction between H and its microbiota ${B}$ is implemented as follows (see Figures 2A,B). Consider the case where a given node σ_n of H has been chosen for mutation such that a new input (output) connection is to be added. Then this new connection can be selected with equal probability either within H itself or from any of the microbial networks $M_j\in {B}$. Likewise, when a given node of a microbial network M_j is mutated so as to receive a new connection (either input or output), the new connection can be established within M_j itself, with H or with any other microbial network $M_k\in {B}$. This allows the emergence of regulatory interactions between all the networks that constitute the holobiont.

FIGURE 2

Figure 2. Network Coevolution. (A). Schematic representation of the holobiont, which in this case consists of one host network H (green) and one microbial network M (red), each with N = 20 nodes and N_s = 4 signal nodes (represented by squares). At generation g = 0 the host network H is not adapted to its task (${\xi }^H>{\delta }_A$), whereas the microbial network M is well adapted (${\xi }^M<{\delta }_A$). The mutation rates μ_H and μ_M of the host and microbial networks, respectively, satisfy μ_M = 10μ_H. (B) At generation g = 150 regulatory interactions between H and M have been established (dashed lines). The highlighted nodes in each network have regulators in the other network. H has become well adapted to its task (${\xi }^H<{\delta }_A$) while the microbial error ξ^M has decreased almost to zero. (C) Population average host error ${\overline{\xi }}^H$ as a function of generations for the holobiont case (H and M evolving together, blue dashed curve, P_M = 1) and for the control case (H evolving by itself, green solid curve, P_M = 0). In the holobiont case the adaptability threshold δ_A is reached faster (g ≈ 100) than in the control case (g ≈ 350). Also, at the end of the simulation (g = 500) the error for the holobiont case is about five times smaller than for the control case. (D) Evolution of the average microbial error ${\overline{\xi }}^M$. At generation g = 0, ${\overline{\xi }}^M$ already satisfies the well-adapted condition, ${\overline{\xi }}^M<{\delta }_A$, but it further decreases as the evolution of the holobiont goes on. (E) Probability of adaptation P_A(g) across generations for the holobiont (blue dashed curve) and control (green solid curve) cases. Note that P_A(g) increases and saturates faster in the holobiont case. (F) Average number ${\overline{\Omega }}^H(g)$ of accumulated mutations in the host network H during its adaptation process for the holobiont case (blue dashed curve) and the control case (green solid curve). Interacting with the microbial network halves the number of mutations H has to undergo in order to adapt to its task. The numerical simulations to generate the graphs (C) to (F) were carried out using networks with N = 50, N_s = 12 and populations of P = 100 networks.

For the adaptation of H to its task we consider the evolution of a population of P = 100 holobionts. Throughout the evolutionary process all the networks in each holobiont undergo the same kind of random mutations described in the section Host Network Evolution (with the possibility of interactions across networks, as mentioned in the previous paragraph). However, in our simulations the mutation rate μ_M for the MNs is ten times larger than the mutation rate μ_H for the host network, namely μ_M = 10μ_H. This captures the fact that bacterial colonies, due to their high reproduction rates, develop mutants at least ten times faster than populations of eukaryotic cells in multicellular organisms (Lynch, 2010; Lynch et al., 2016). It is important to emphasize that in our model each network in the holobiont has to be considered not as representing a single cell, but an entire cell population. In each generation, holobionts are ranked according to their error ξ^L, and the ten with the smallest errors are selected for reproduction (see M&M). In all further simulations the unit of selection is the holobiont, as in each generation we select the ten best holobionts (based on the error ξ^L that takes into account the host and microbial errors) and replicate them.

Figure 2C shows the average population error ${\overline{\xi }}^H$ of the host network H across generations for holobionts as well as for the control case (host networks evolving by themselves without interacting with microbial networks). In the simulations reported in Figure 2C each holobiont consists of one host network H and one microbial network M (P_M = 1). It is clear that interacting with only one microbial network M already makes H to adapt much faster to its task than evolving on its own. In the holobiont case, the error ${\overline{\xi }}^H$ crosses the adaptability threshold δ_A in about one fourth of the generations required for the control case to do it (see Movie S2 and compare with Movie S1). Furthermore, the final error after 500 generations is considerably smaller for the holobiont case (${\overline{\xi }}^H\approx 0.2$) than for the control case (${\overline{\xi }}^H\approx 0.95$). Note that an error ${\overline{\xi }}^H\approx 0.2$ means that, on average, at most 3 points of the output signal R(t) deviate one unit from the task F(t) in the whole interval 1 ≤ t ≤ 15, which represents a percent error 100 × 3/(N_s × 15) ≈ 1.6%. This is almost a perfect adaptation hard to achieve in the control case. For the control case, it takes about 3000 generations to reach a similar error of ${\overline{\xi }}^H\approx 0.2$ (See Figure S3). The average microbial error ${\overline{\xi }}^M$ also decreases, as shown in Figure 2D. At generation g = 0, ${\overline{\xi }}^M$ is already below the adaptability threshold δ_A, but it decreases even further as the evolution of holobionts proceeds. Therefore, the adaptation of the holobiont takes place with no conflict of interest between H and its microbial networks.

In Figure 2E we report the probability of adaptation P_A(g), defined as the fraction of holobionts in which the host-network error ξ^H crosses the adaptation threshold δ_A at generation g. It is apparent from Figure 2E that this probability for the holobiont case increases and saturates much faster than for the control case. About 80% of the holobionts are well adapted after only 120 generations, whereas host networks evolving by themselves never reach 75% of adaptation during the whole simulation time. In addition to speeding up and increasing the adaptation of H, the interaction between H and M also considerably reduces the number of mutations H has to accumulate in order to adapt to its task, as Figure 2F reveals. This is not a trivial result, for only the mutations in both H and M that increase the adaptation of the holobiont are selected and fixed in the population. Thus, even though M mutates ten times faster than H, not all of those mutations are beneficial to the adaptation of the holobiont and consequently, not all of them become fixed in the population. Actually, from Figure 2E we observe that the average number ${\overline{\Omega }}^H$ of accumulated mutations in H to reach the adaptation threshold δ_A is not ten, but only two times larger for the control case than for the holobiont case. However, it is true that because μ_M is larger than μ_H the adaptation of H is improved (our simulations show that there is no significant difference between the holobiont and control cases when μ_M = μ_H, see Figure S4).

Symbiosis and Dysbiosis

To show that symbiotic relationships emerge between the host and microbial networks, once the holobiont is well adapted (after 500 generations as in Figure 2C,D), we remove the connections between H and M (the dashed lines in Figure 2B) and compute the errors of each network at performing their respective tasks while disconnected. This can be thought of as an antibiotic administration where several bacterial species are removed from the microbial population, or as trying to cultivate these symbiotic microbes without their respective host. Thus, a set of microbial species, represented by M, are removed from the holobiont and then the fitness of the host is evaluated without them. At the same time, we determine the survivability of these microbial species M in the absence of their host. Since M starts the evolutionary process already well adapted to its task, one can expect that its error does not significantly increase after the connections between H and M are removed. However, Figure 3 shows a typical example in which removing the connections between H and M increases both errors ξ^H and ξ^M to values that correspond to untrained networks. Thus, in the example shown in Figure 3, after the holobiont has been adapted as a whole, none of the networks it consists of can perform their respective tasks when separated (see Figure S5 for population statistical averages).

FIGURE 3

Figure 3. Emergence of symbiosis and dysbiosis. A holobiont consisting of one host network H and one microbial network M₁ evolves for 500 generations. Then H and M₁ are disconnected and the errors at performing their respective tasks evaluated (H and M₁ are represented by hexagons at the bottom of the bar chart). At generation g = 0, when the training of H begins, the host-network error ξ^H (green bar) is large whereas the microbial-network error ξ^M (red bar) is already below the adaptation threshold δ_A. After H and M₁ have coevolved for 500 generations both ξ^H and ξ^M are quite below δ_A, which indicates the adaptation of the entire holobiont. Then, H and M₁ are disconnected and their respective errors evaluated. Note that after disconnection both errors ξ^H and ξ^M in this example increase to levels corresponding to completely untrained networks. The simulations were performed with networks having N = 50 nodes and N_s = 12 signal nodes.

Multitasking and Microbial Diversity

So far we have presented results in which the host network H has to perform only one task. Interaction with one MN significantly improves the adaptation of H (and of the holobiont) and reduces the number of mutations it has to undergo in order to become well adapted. It could be expected that adding more MNs to the microbiota would further enhance the adaptation of H. However, this is not the case. Adding more MNs either has no effect or can even worsen the adaptation of the host (see Movie S3 and Figure S6). This result is in contradiction with the great diversity observed in the microbiota of real organisms and the ability of the holobiont to adequately perform multiple tasks.

For this reason, we now consider the case in which the host network H is trained to perform T multiple tasks F₁(t), F₂(t), …, F_T(t), each being an arbitrary function constructed as described in the Task Assignment section. Since Boolean Network dynamics are deterministic, depending on the initial condition the dynamics of H will be set to follow a specific task F_τ(t). We can measure the adaptation errors ${\xi }_{\tau }^H$ and ${\xi }_{j,\tau }^M$ of H and the microbial network $M_j\in {B}$, respectively, when H is being trained to perform the particular task F_τ(t) (see the M&M section for a precise definition of ${\xi }_{\tau }^H$ and ${\xi }_{j,\tau }^M$). This allows us to compute the adaptation of the holobiont separately for each task. Averaging ${\xi }_{\tau }^H$ and ${\xi }_{j,\tau }^M$ over all the tasks gives us the total adaptation errors ξ^H and ${\xi }_j^M$ for the host and microbial networks, respectively (see the M&M section).

We implement two ways in which the MNs can assist the adaptation of H to perform many different tasks. First, there is the non-specialized interaction case in which all the MNs can interact among each other and with H. Also, all the MNs can participate in the adaptation of H to all of its tasks (see Figure 4A). In each generation the networks are mutated, allowing new interactions to appear between any two networks within the holobiont. This means that new incoming or outgoing connections can be established either between H and any of its MNs, or between any two MNs. We consider again a population of P = 100 holobionts. After the networks in each holobiont have been mutated (with the mutation rates μ_H and μ_M for the host and microbial networks, respectively), the ten best holobionts are selected and replicated (see the M&M section for a definition of the holobiont error ξ^L in the multitasking case). Figure 4B shows the population average ${\overline{\xi }}^H$ of the host-network error for the case in which H has to perform 10 different tasks. It is clear from this figure that adding more than one microbial network to the microbiota has no effect on the adaptation of H to its 10 different tasks. Therefore, in the non-specialized case increasing the diversity of the microbiota does not help the adaptation rate of the host.

FIGURE 4

Figure 4. Multitasking and microbial specialization. The host network H (green circle) has to perform T different tasks F₁, F₂, …, F_T (represented as ellipses) assisted by P_M microbial networks M₁, M₂, …, M_PM (small circles). (A) Schematic representation of a holobiont with non-specialized interaction in the microbial networks. Each microbial network can interact with all the other ones and participate in the adaptation of H to all of its tasks. The dashed lines represent possible interactions. (B) Evolution of the host-network error ${\overline{\xi }}^H$ in the non-specialized case for T = 10. Each curve corresponds to a different value of P_M, starting from P_M = 0 (H evolving by itself) up to P_M = 10. In none of these cases ${\overline{\xi }}^H$ crosses the adaptation threshold δ_A even after 500 generations. Furthermore, adding more than one microbial network to the holobiont has a negligible effect on the adaptation of H. (C) Schematic representation of the holobiont with specialized interaction in the microbial networks. The set of P_M microbial networks is divided into T different niches, G₁, …, G_T, each one participating in the adaptation of H to only one task. The microbial networks can interact among them only if they belong to the same niche. (D) Host-network error ${\overline{\xi }}^H$ across generations for the specialized interaction and T = 10. Again, each curve corresponds to a different number P_M of microbial networks. Note that in this case the more microbial networks in the holobiont the better the adaptation of H. Note also that for P_M = 10 the host-network error ${\overline{\xi }}^H$ crosses the adaptation threshold δ_A in less than 300 generations. (E) Adaptation probability P_A at generation g = 500 as a function of the number P_M of microbial networks participating in the host's adaptation. The blue and green bars correspond to the non-specialized and specialized cases, respectively. Note that in the non-specialized case the adaptation probability remains low regardless of the number of microbial networks, whereas in the specialized case, the more microbial networks the better the adaptation of the host. The simulations were carried out for populations of P = 100 holobionts and networks with N = 50 nodes.

As a second alternative we implement a specialized interaction in the microbial networks. In this case the microbiota ${B}=\left\{M_1,\dots ,M_{P_M}\right\}$ is divided into P_G disjoint non-empty subsets, or “niches”, {G₁, G₂, …, G_PG}. The set of tasks ${F}=\left\{F_1,\dots ,F_T\right\}$ is also partitioned, as evenly as possible, into P_G non-overlapping subsets, $\left\{{{T}}_1,{{T}}_2,\dots ,{{T}}_{P_G}\right\}$. The maximum number of niches is P_G = T, for in this case each subset ${{T}}_{\tau }$ contains only one task. To each niche G_τ we associate a subset ${{T}}_{\tau }$ of tasks (see Figure 4C). The host network H is still trained to perform the T different tasks F₁, …, F_T. However, the training of H to the tasks in the particular set ${{T}}_{\tau }$ is assisted only by the networks in the corresponding niche G_τ. For each niche G_τ we compute an error ${\xi }_{\tau }^G$ that measures the adaptation of the holobiont when H is being trained to perform the tasks in the specific subset ${{T}}_{\tau }$ (see the M&M section for a precise definition of the niche error ${\xi }_{\tau }^G$). During the adaptation of H to the tasks in ${{T}}_{\tau }$, only the mutations in H or in the microbial networks belonging to G_τ that reduce the corresponding error ${\xi }_{\tau }^G$ are selected. The important point to note here is that the adaptation of H to its tasks can be measured separately for each niche. The holobiont error ξ^L is computed as the average of the errors ${\xi }_{\tau }^G$ over all the niches (see the M&M section).

During the training of H, the MNs in one niche can develop interactions between them, but they cannot interact with the MNs in a different niche, as Figure 4C indicates. This is consistent with the observation that multicellular organisms maintain the stability of its microbiota by reducing microbial interactions (Deines et al., 2017). If microbes interacted with no organization, the loss of one microbial species would affect the fitness of all the others, increasing the risk of extinction cascades (Coyte et al., 2015). Therefore, in the specialized interaction case we compartmentalize the MNs allowing interactions among them only if they belong to the same niche (all the MNs in all the niches can, of course, interact with the host network H).

Figure 4D shows the evolution of the average host-error ${\overline{\xi }}^H$ for simple case where each niche has one MN (P_G ≤ T). It is clear from Figure 4D that, contrary to the non-specialized case, adding more MNs to the microbiota in a specialized way considerably improves the adaptation of H to its multiple tasks. Furthermore, in Figure 4E we report the probability of adaptation P_A(500) at generation g = 500 as a function of P_M for both the non-specialized and specialized cases. In the former case P_A remains low and never improves as P_M increases, whereas in the later case P_A monotonously grows with P_M. This clearly shows that both diversity and specialization of the MNs are necessary for H to adapt to multiple tasks.

The results presented for the specialized interaction case also hold when every niche is populated with more than one MN (P_M ≥ 2T). In Figure 5A (see also Movie S4) we report the evolution of a holobiont with T = 10 different tasks and the same number of niches, and P_M = 25 microbial networks (each niche contains either two or three MNs). At generation g = 0, H is poorly adapted to all of its tasks (represented in red), whereas all the MNs are already well adapted (represented in blue). As the evolution proceeds H becomes more adapted to all of its tasks. Furthermore, the microbial networks also become more adapted to their own tasks. Note that the adaptation of H to its different tasks occurs at different rates, as can be seen from the color-code of the tasks through the holobiont evolution. This is consistent with the observation that different symbiotic relationships between the host and its microbial communities emerge at different rates (Doolittle and Booth, 2017). In the example shown in Figure 5A the adaptation of the whole holobiont crosses the adaptation threshold δ_A in <300 generations. Interestingly, the same results are obtained in the specialized case when many more microbial networks are introduced into the holobiont (see Figure S7).

FIGURE 5

Figure 5. Holobiont's complexity and loss of microbial diversity. (A) Snapshots of the evolution of a holobiont in the specialized case. The holobiont consists of T = 10 tasks (ellipses) and P_M = 25 microbial networks (small circles). All errors larger than 3 are colored in red. At generation g = 0 the host network H has a large error in all of its tasks, while the bacterial networks are well adapted. As the evolution proceeds, the error of the entire holobiont decreases, but the adaptation of H to its different tasks takes place at different rates. The number at the center is the value of the host error ${\overline{\xi }}^H$ averaged over all of its tasks. (B) After the holobiont is well adapted at generation g = 500, Δn microbial networks are disconnected and the resulting host-network error ${\overline{\xi }}^H$ is computed. This simulates the effect on the holobiont's adaptation of a reduced microbial diversity. Note that ${\overline{\xi }}^H$ gradually increases as more microbial networks are disconnected, indicating that the host becomes less adapted as the microbial diversity decreases.

The specialized interaction scheme allows us to compute the robustness of the holobiont under loss of microbial diversity. For this, once the holobiont is well adapted, we disconnect Δn microbial networks from it and compute the resulting host-network error ξ^H averaged over all the host's tasks. Figure 5B shows that ξ^H gradually increases as more microbial networks are disconnected from the holobiont. Therefore, a loss in microbial diversity clearly reduces the adaptation of the host.

Discussion

Multicellular organisms and microbes have coevolved in many different ways, not only as holobionts being units of adaptation (Theis et al., 2016). However, the persistence across generations of regulatory interactions between the host and its microbes is a necessary condition for natural selection to operate at the holobiont level (Doolittle and Booth, 2017). These regulatory interactions have to preserve the holobiont's functionality regardless of the specific microbial species that generate them. This is the “It's the song not the singer” (ITSNTS) approach to evolution proposed by Doolittle (Doolittle and Booth, 2017; Doolittle and Inkpen, 2018) and exemplified by Taxis et al. in ruminal ecosystems (Taxis et al., 2015). In this work we have incorporated into a single evolutionary model both the concept that the holobiont is a unit of selection and Doolittle's ITSNTS approach. We have done so by requiring in our simulations that only the best adapted holobionts at each generation are the ones able to go throughout the selective filter, pass to the next generation and replicate all of its constituent networks. However, in our model selection acts on a dynamical property of the holobiont, which is the host's output signal, in order to bring it close to the functions the host needs to perform. In this scheme, it does not matter what nodes or microbial networks participate in the regulation of the dynamical functions. What is important is the preservation across generations of the dynamical functions themselves, and this must hold for both the host and the microbial networks (which also have to perform and preserve their own functions). Thus, although the holobiont might not be the only unit of selection (Theis et al., 2016), we have concentrated on those important host-microbe co-interactions that transmit functionality across multiple generations (Roughgarden et al., 2018). Our model does not assume that the microbial networks reside inside the host, but only that they interact with it and that the host-microbe interactions are transmitted across generations. This propagation can happen in various ways other than vertical transmission from parents to offspring as, for instance, when the host constructs its environment with a stable microbial composition (Fitzpatrick, 2014).

We have shown that the host network can actually be trained to perform one task without the help of any microbial networks, as Figure 1D and Figure S3 illustrate. However, allowing interactions between the host and microbial networks greatly speeds up and improves the adaptation of the entire holobiont. This is because the host network does not only adapt to its tasks faster, better and with less mutations when it is allowed to interact with microbial networks, but the microbial networks themselves considerably improve their own adaptation to their respective tasks. Furthermore, adaptation of the host network to perform multiple tasks is improved only when it is allowed to interact with a diverse and specialized microbiota, as Figure 4D shows.

In light of these results, we observe that the microbiota does not only help the host to adapt to its tasks. There is mutual benefit in which both the host and its microbial communities contribute to each other's adaptation. It is in this sense that the holobiont can be considered as an evolutionary unit.

It is important to mention that in our model the holobiont cannot just be considered as one “big network” evolving to perform a set of tasks. There are two essential aspects that have to be emphasized. First, the rate at which mutants are generated μ_M in the microbial networks is considerably larger than that μ_H of the host networks. Second, the set of microbial networks must be partitioned into disjoint (i.e., non-interacting) niches for the host network to efficiently adapt to multiple tasks, where each niche specializes in the adaptation of the host to one specific subset of tasks. These two aspects provide the holobiont with a complex internal dynamical structure that prevents us from viewing it as just one big network (see Figure 5A and Figure S7). A holobiont for which μ_M = μ_H and the microbial networks are not partitioned into specialized niches, could be considered as a single large homogeneous network. But in such a homogeneous case the holobiont evolution benefits neither from the host-microbe interactions nor from the microbiota's diversity (see Figures S4, S8). Rather, the structural and dynamical internal complexity of the holobiont, embodied in the functional modularity and specialization of the microbial niches as well as in the difference between the host and microbial mutation rates, is required to facilitate and improve the holobiont's adaptation to perform multiple tasks. Hence, our results show that complexity, modularity and functional specialization are necessary properties that naturally facilitate the evolution and adaptation of the holobiont as well as the diversification of the microbiota (Sachs et al., 2014), whereas structural and functional homogeneity is either inconsequential or even detrimental to the holobiont's evolution.

One may wonder whether the relationship μ_M = 10μ_H between the microbial and host network's mutant-generating rates accurately reflects reality. We have explored a wide range of values of the ratio γ = μ_M/μ_H, ranging from γ = 1 to γ = ∞. The latter case corresponds to μ_H = 0, which means that the adaptation of the host to its task does not occur across generations, but within the host's lifespan. In this extreme case, the adaptation of the host to its task occurs due to mutations in the microbiota but not in the host itself. Our simulations show that the adaptation of the host network is almost equally accelerated and improved for γ = 10 than for γ = ∞ (see Figure S9).

In our model the host network interacts with microbial networks which, from the very beginning, are already well adapted to perform their own functions. The reason for this is to determine whether or not the well-adapted condition imposed on the microbial networks represents a restriction that could generate evolutionary conflict within the holobiont. It has been pointed out that the emergence of symbiotic relationships between organisms requires the symbionts to be highly cooperative and show very little conflict (Morris et al., 2012; Sachs and Hollowell, 2012; Sachs et al., 2014; Queller and Strassmann, 2016). Our simulations show that the evolution of the holobiont can very well take place with no evolutionary conflict between its constituent networks, as long as the microbiota is partitioned into specialized niches. This modularization and division of labor are essential to prevent microbial competition and within-group conflict in the holobiont (West et al., 2015) (see Figure S6). Additionally, modularization of the microbiota allows the holobiont to acquire new functions without affecting the ones already present.

Interestingly, similar results regarding the adaptation of the host network to its tasks are obtained when the well-adapted condition is not imposed on the microbial networks. Our simulations show that the adaptation of the host network is equally accelerated and improved when it interacts with microbial networks that do not have to perform any task (see Movie S5). However, even when the microbial networks are free of any selective pressure, their dynamics are stabilized when they coevolve with the host network (see Movie S5). This is important because it can be interpreted as the holobiont acquiring, at any moment, microbes from the environment and coevolving with them, generating intergenomic epistasis that reduces within-group conflict and promotes the adaptation of the entire holobiont (Bordenstein and Theis, 2015).

Finally, we would like to mention that, although there exist many qualitative taxonomical studies showing the existence of a great variety of host-microbe symbiotic interactions, there are very few mathematical and computational models aiming to explain the general mechanisms responsible for the emergence of such interactions and the need for diversification and specialization of the microbiota (Manor et al., 2014). We have not explicitly considered competition or parasitism in our model. However, by integrating the ITSNTS approach with the hologenome hypothesis (the holobiont as a unit of selection Rosenberg and Zilber-Rosenberg, 2016; Roughgarden et al., 2018), we were able to reproduce many of the observed behaviors in the evolution of holobionts, such as reduction of evolutionary conflict, division of labor, emergence of symbiotic interactions and dysbiosis when the microbiota diversity is reduced. Our model may thus lay the foundations for a comprehensive understanding of the long-lasting coevolution of multicellular organisms and microbes.

Materials and Methods

Boolean Network Model

The Boolean network consists of a set of N nodes {σ₁, σ₂, …, σ_N}, each acquiring the values 0 or 1 that represent two possible states of activity: “active” or “inactive.” The value of each node σ_n is determined through a logical function f_n that depends on a set of k_n other nodes in the network denoted as $I_n=\left\{{\sigma }_1^n,{\sigma }_2^n,\cdots ,{\sigma }_{k_n}^n\right\}$. The nodes in the set I_n are known as the inputs or regulators of σ_n. In the context of genetic networks these regulators together with the logical function f_n mimic the effect of k_n transcription factors (synthesized by the regulators) acting on the expression of σ_n. For networks of real organisms both the logical function f_n and the set of regulators I_n associated to each gene are carefully constructed according to the activating and inhibitory nature of the regulatory interactions between the genes. Here, the k_n regulators of each node σ_n are randomly chosen from anywhere in the network. The logical functions f_n are also randomly chosen from the ensemble of all possible logical functions with k_n variables. In this work we start the simulations with networks for which k_n = K = 2, which means that every node in the initial networks has exactly K = 2 regulators (chosen randomly). Note that this is a directed network because if node σ_m is a regulator of σ_n, it does not necessarily happen that σ_n is also a regulator of σ_m (although it may happen). Note also that throughout the evolution of the networks some input and output connections are added to, or removed from, different nodes. Therefore, the final networks do not have a constant connectivity K = 2 for every node. The final networks will have a connectivity distribution similar to the one observed in the Erdös-Rényi topology with an average around K = 2 (see Movies S1, S2). Once each node in the network has been provided with a set of inputs and a logical function, the network dynamics is given by

$\begin{array}{ll}{\sigma }_n\left(t+1\right)=f_n\left({\sigma }_1^n\left(t\right),{\sigma }_2^n\left(t\right),\cdots ,{\sigma }_{k_n}^n\left(t\right)\right).& (3)\end{array}$

Starting the dynamics from an initial condition Σ = {σ₁(0), σ₂(0), …, σ_N(0)}, the network transits throughout a series of states until a periodic pattern is reached, which is known as a dynamical attractor. There is a great body of work showing that the dynamical attractors of the network correspond to different cell types or cell fates (or more generally, to different functional states of the cell). Here we are not interested in the dynamical attractors, but in training the network to perform a specific task.

Microbial and Holobiont Errors: One Task

Let us consider a holobiont ${L}=\left\{H,M_1,\dots ,M_{P_M}\right\}$. When the host network H has to perform only one task F(t), its error ξ^H is given in Equation (2). We similarly define the microbial error ${\xi }_j^M$ corresponding to the j^th microbial network M_j as ${\xi }_j^M=\frac{1}{t_m}{\displaystyle {\sum }_{t=1}^{t_m}{\left(R_j^M(t)-F_j^M(t)\right)}^2}$, where $R_j^M(t)$ and $F_j^M(t)$ are the output signal and target function of M_j, respectively. Different tasks are assigned to the different microbial networks. The error ξ^L of the entire holobiont is then computed as

$\begin{array}{ll}{\xi }^L=\frac{1}{1+P_M}\left({\xi }^H+{\displaystyle \sum _{j=1}^{P_M}{\xi }_j^M}\right).& (4)\end{array}$

In our simulations the population contains P = 100 holobionts, ${{L}}_1,{{L}}_2,\dots ,{{L}}_P$. For each holobiont ${{L}}_i$ we compute its error ${\xi }_i^L$ as in Equation (4), which is then used at each generation to select the best holobionts in the population.

We also preformed numerical simulations using the following definition for the holobiont error:

$\begin{array}{ll}{\xi }^L=\frac{1}{2}\left({\xi }^H+\frac{1}{P_M}{\displaystyle \sum _{j=1}^{P_M}{\xi }_j^M}\right).& (5)\end{array}$

The difference between Equations (4 and 5) is the contribution of the microbiota to the holobiont's error. In Equation (4) the contribution of the total microbial error could be very large as compared to the contribution of the host-network error, especially if there are many microbial networks in the holobiont. Contrary to this, in Equation (5) the microbiota and the host have the same contribution regardless of the number of microbial networks in the microbiota. Our simulations show that both definitions produce qualitatively the same results (see Figure S6B). This is because at each generation in the evolutionary process we are selecting the best holobionts in the population, and selecting the best holobionts eventually leads to the same type of individuals regardless of the way in which the contribution of each particular network is weighted. In this work we present results using the definition given in Equation (5).

Host and Microbial Errors: Multitasking Non-specialized Case

Let us consider a holobiont ${L}=\left\{H,M_1,\dots ,M_{P_M}\right\}$ where now the host network H has to perform T different tasks F₁(t), F₂(t), …, F_T(t). For each task F_τ(t) the network starts its dynamics from a predefined initial condition ${\Sigma }_{\tau }=\left\{{\sigma }_1^{\tau }(0),{\sigma }_2^{\tau }(0),\dots ,{\sigma }_N^{\tau }(0)\right\}$. Let us denote as R_τ(t) the output signal of H when it starts its dynamics from the initial condition Σ_τ. The error ${\xi }_{\tau }^H$ corresponding to the task F_τ(t) is computed as

$\begin{array}{ll}{\xi }_{\tau }^H=\frac{1}{t_m}{\displaystyle \sum _{t=1}^{t_m}{\left(R_{\tau }(t)-F_{\tau }(t)\right)}^2}.& (6)\end{array}$

This allows us to measure the adaptation of the host network to each of its tasks separately. The total adaptation error ξ^H of H is computed by averaging the corresponding errors over all the T tasks that H has to perform: ${\xi }^H=T^{-1}{\displaystyle {\sum }_{\tau =1}^T{\xi }_{\tau }^H}$.

To define the error ${\xi }_{j,\tau }^M$ of the microbial network M_j when the host network is being trained to perform the task F_τ(t) we have to consider the fact that H may be regulating some of the nodes of M_j (throughout the evolution of the holobiont such regulations may appear). Therefore, the output signal of M_j depends on the initial condition Σ_τ used to start the dynamics of H. Let us denote as $R_{j,\tau }^M$ the output signal of M_j when H started its dynamics from the initial condition Σ_τ. The corresponding microbial error ${\xi }_{j,\tau }^M$ is then defined as

$\begin{array}{ll}{\xi }_{j,\tau }^M=\frac{1}{t_m}{\displaystyle \sum _{t=1}^{t_m}{\left(R_{j,\tau }^M(t)-F_j^M(t)\right)}^2},& (7)\end{array}$

where $F_j^M$ is the task assigned to the microbial network M_j (this network was already well adapted to its tasks and has to remain so during the evolutionary process). The total microbial error ${\xi }_j^M$ corresponding to M_j is then computed by averaging ${\xi }_{j,\tau }^M$ over all the tasks: ${\xi }_j^M=T^{-1}{\displaystyle {\sum }_{\tau =1}^T{\xi }_{j,\tau }^M}$. The holobiont error ξ^L is computed using Equation (4), where now ξ^H and ${\xi }_j^M$ are the host and microbial errors averaged over all the tasks as described above.

Niche Error: Multitasking Specialized Case

Let us consider the niche G_τ = {M_τ1, M_τ2, …, M_τpτ}, containing p_τ microbial networks. This niche is helping H to adapt to the tasks in the set ${{T}}_{\tau }=\{F_{{\tau }_1},F_{{\tau }_2},\dots ,F_{{\tau }_{q_{\tau }}}\}$, which contains q_τ tasks. The error ${\xi }_{\tau }^G$ corresponding to this niche is defined as

$\begin{array}{ll}{\xi }_{\tau }^G=\frac{1}{1+p_{\tau }}\left({\xi }_{\tau }^H+{\displaystyle \sum _{i=1}^{p_{\tau }}\frac{1}{q_{\tau }}}{\displaystyle \sum _{j=1}^{q_{\tau }}{\xi }_{{\tau }_i,{\tau }_j}^M}\right),& (8)\end{array}$

where ${\xi }_{\tau }^H$ and ${\xi }_{{\tau }_i,{\tau }_j}^M$ are defined in Equations (6, 7), respectively (in the latter case the subscripts j and τ change to τ_i and τ_j respectively, since we have to account for the different microbial networks and functions associated to the niche G_τ).

The holobiont error for the specialized interaction case is computed by averaging ${\xi }_{\tau }^G$ over all the niches: ${\xi }^L=\frac{1}{P_G}{\displaystyle {\sum }_{\tau =1}^{P_G}{\xi }_{\tau }^G}$.

Author Contributions

SH, SS-M, AF, and MA: conceived the experiments, analyzed the data, and wrote and revised the manuscript. SH and MA: designed the experiments. SH: designed the software used in the analysis.

Funding

MA thanks the Marcos Moshinsky Foundation for the 2014 Fellowship. SH thanks CONACyT for a Ph.D. scholarship. This work was partly supported by PAPIT-UNAM grant IN226917.

Conflict of Interest Statement

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.

Acknowledgments

We thank Philippe Cluzel, Osbaldo Resendis and Pablo Vinuesa for useful comments and discussions.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2018.01836/full#supplementary-material

References

Albert, R., and Othmer, H. G. (2003). The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in drosophila melanogaster. J. Theor. Biol. 223, 1–18. doi: 10.1016/S0022-5193(03)00035-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Aldana, M. (2003). Boolean dynamics of networks with scale-free topology. Physica D 185, 45–66. doi: 10.1016/S0167-2789(03)00174-X

CrossRef Full Text | Google Scholar

Aldana, M., Balleza, E., Kauffman, S., and Resendiz, O. (2007). Robustness and evolvability in genetic regulatory networks. J. Theor. Biol. 245, 433–448. doi: 10.1016/j.jtbi.2006.10.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Andrew, W., Kohl, K. D., Brucker, R. M., Edward, J., Opstal, V., Phylosymbiosis, S. R. B., et al. (2016). Phylosymbiosis : relationships and functional effects of microbial communities across host evolutionary history. PLoS Biol. 14:e2000225. doi: 10.1371/journal.pbio.2000225

CrossRef Full Text | Google Scholar

Balleza, E., Alvarez-Buylla, E. R., Chaos, A., Kauffman, S., Shmulevich, I., and Aldana, M. (2008). Critical dynamics in genetic regulatory networks: examples from four kingdoms. PLoS ONE 3:e2456. doi: 10.1371/journal.pone.0002456

PubMed Abstract | CrossRef Full Text | Google Scholar

Bates, J. M., Mittge, E., Kuhlman, J., Baden, K. N., Cheesman, S. E., and Guillemin, K. (2006). Distinct signals from the microbiota promote different aspects of zebrafish gut differentiation. Dev. Biol. 297, 374–386. doi: 10.1016/j.ydbio.2006.05.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Blaser, M. J. (2017). The theory of disappearing microbiota and the epidemics of chronic diseases. Nat. Rev. Immunol. 17, 461–463. doi: 10.1038/nri.2017.77

PubMed Abstract | CrossRef Full Text | Google Scholar

Blaser, M. J., and Falkow, S. (2009). What are the consequences of the disappearing human microbiota? Nat. Rev. Microbiol. 7, 887–894. doi: 10.1038/nrmicro2245

PubMed Abstract | CrossRef Full Text | Google Scholar

Bordenstein, S. R., and Theis, K. R. (2015). Host biology in light of the microbiome: ten principles of holobionts and hologenomes. PLoS Biol. 13:e1002226. doi: 10.1371/journal.pbio.1002226

PubMed Abstract | CrossRef Full Text | Google Scholar

Brucker, R. M., and Bordenstein, S. R. (2013). The hologenomic basis of speciation: gut bacteria cause hybrid lethality in the genus nasonia. Science 466, 667–669. doi: 10.1126/science.1240659

CrossRef Full Text | Google Scholar

Camp, J. G., Frank, C. L., Lickwar, C. R., Guturu, H., Rube, T., Wenger, A. M., et al. (2014). Microbiota modulate transcription in the intestinal epithelium without remodeling the accessible chromatin landscape. Genome Res. 24, 1504–1516. doi: 10.1101/gr.165845.113

PubMed Abstract | CrossRef Full Text | Google Scholar

Cho, I., and Blaser, M. J. (2012). The human microbiome: at the interface of health and disease. Nat. Rev. Genet. 13, 260–270. doi: 10.1038/nrg3182

PubMed Abstract | CrossRef Full Text | Google Scholar

Consortium, H. M. P. (2012). Structure, function and diversity of the healthy human microbiome. Nature 486, 207–214. doi: 10.1038/nature11234

CrossRef Full Text | Google Scholar

Contreras, A. V., Cocom-Chan, B., Hernandez-Montes, G., Portillo-Bobadilla, T., and Resendis-Antonio, O. (2016). Host-microbiome interaction and cancer: potential application in precision medicine. Front. Physiol. 7:606. doi: 10.3389/fphys.2016.00606

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Daniels, B. C., Kim, H., Moore, D., Zhou, S., Smith, H. B., Karas, B., et al. (2018). Criticality distinguishes the ensemble of biological regulatory networks. Phys. Rev. Lett. 121, 138102. doi: 10.1103/PhysRevLett.121.138102

PubMed Abstract | CrossRef Full Text | Google Scholar

Davidich, M. I., and Bornholdt, S. (2008). Boolean network model predicts cell cycle sequence of fission yeast. PLoS ONE 3:e1672. doi: 10.1371/journal.pone.0001672

PubMed Abstract | CrossRef Full Text | Google Scholar

Davidson, E. H. (2010). Emerging properties of animal gene regulatory networks. Nature 468, 911. doi: 10.1038/nature09645

PubMed Abstract | CrossRef Full Text | Google Scholar

Davidson, E. H., and Levine, M. S. (2008). Properties of developmental gene regulatory networks. Proc. Natl. Acad. Sci. U.S.A. 105, 20063–20066. doi: 10.1073/pnas.0806007105

PubMed Abstract | CrossRef Full Text | Google Scholar

Deines, P., Lachnit, T., and Bosch, T. C. (2017). Competing forces maintain the hydra metaorganism. Immunol. Rev. 279, 123–136. doi: 10.1111/imr.12564

PubMed Abstract | CrossRef Full Text | Google Scholar

Derrida, B., and Pomeau, Y. (1986). Random networks of automata: a simple annealed approximation. Europhys. Lett. 1, 45–49. doi: 10.1209/0295-5075/1/2/001

CrossRef Full Text | Google Scholar

Donaldson, G. P., Lee, S. M., and Mazmanian, S. K. (2015). Gut biogeography of the bacterial microbiota. Nat. Rev. Microbiol. 14, 20–32. doi: 10.1038/nrmicro3552

PubMed Abstract | CrossRef Full Text | Google Scholar

Doolittle, W. F., and Booth, A. (2017). It's the song, not the singer: an exploration of holobiosis and evolutionary theory. Biol. Philos. 32, 5–24. doi: 10.1007/s10539-016-9542-2

CrossRef Full Text | Google Scholar

Doolittle, W. F., and Inkpen, S. A. (2018). Processes and patterns of interaction as units of selection: an introduction to ITSNTS thinking. Proc. Natl. Acad. Sci. U.S.A. 115, 201722232. doi: 10.1073/pnas.1722232115

PubMed Abstract | CrossRef Full Text | Google Scholar

Douglas, A. E., and Werren, J. H. (2016). Holes in the hologenome: why host-microbe symbioses are not holobionts. MBio 7, e02099–15. doi: 10.1128/mBio.02099-15

CrossRef Full Text | Google Scholar

Ereshefsky, M., and Pedroso, M. (2015). Rethinking evolutionary individuality. Proc. Natl. Acad. Sci. U.S.A. 112, 10126–10132. doi: 10.1073/pnas.1421377112

PubMed Abstract | CrossRef Full Text | Google Scholar

Espinosa-Soto, C., Padilla-Longoria, P., and Alvarez-Buylla, E. R. (2004). A gene regulatory network model for cell-fate determination during arabidopsis thaliana flower development that is robust and recovers experimental gene expression profiles. Plant Cell 16, 2923–2939. doi: 10.1105/tpc.104.021725

PubMed Abstract | CrossRef Full Text | Google Scholar

Farrell, J. J., Zhang, L., Zhou, H., Chia, D., Elashoff, D., Akin, D., et al. (2011). Variations of oral microbiota are associated with pancreatic diseases including pancreatic cancer. Gut 61, 582–588. doi: 10.1136/gutjnl-2011-300784

PubMed Abstract | CrossRef Full Text | Google Scholar

Fernández, L., Langa, S., Martín, V., Maldonado, A., Jiménez, E., Martín, R., et al. (2013). The human milk microbiota: origin and potential roles in health and disease. Pharmacol. Res. 69, 1–10. doi: 10.1016/j.phrs.2012.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Fitzpatrick, B. M. (2014). Symbiote transmission and maintenance of extra-genomic associations. Front. Microbiol. 5, 1–15. doi: 10.3389/fmicb.2014.00046

PubMed Abstract | CrossRef Full Text | Google Scholar

Foster, K. R., and Bell, T. (2012). Competition, not cooperation, dominates interactions among culturable microbial species. Curr. Biol. 22, 1845–1850. doi: 10.1016/j.cub.2012.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Francescone, R., Hou, V., and Grivennikov, S. I. (2014). Microbiome, inflammation and cancer. Cancer J. 20, 181. doi: 10.1097/PPO.0000000000000048

PubMed Abstract | CrossRef Full Text | Google Scholar

Gilbert, S. F., Sapp, J., and Tauber, A. I. (2012). A symbiotic view of life: we have never been individuals. Q. Rev. Biol. 87, 325–341. doi: 10.1086/668166

PubMed Abstract | CrossRef Full Text | Google Scholar

Gonzalez, A., Stombaugh, J., Lozupone, C., Turnbaugh, P. J., Gordon, J. I., and Knight, R. (2011). The mind-body-microbial continuum. Dialogues Clin. Neurosci. 13, 55–62.

PubMed Abstract | Google Scholar

Gordon, J., Knowlton, N., Relman, D. A., Rohwer, F., and Youle, M. (2013). Superorganisms and holobionts. Microbe 8, 152–153. doi: 10.1128/microbe.8.152.1

CrossRef Full Text | Google Scholar

Grice, E. A., and Segre, J. A. (2011). The skin microbiome. Nat. Rev. Microbiol. 9, 244–253. doi: 10.1038/nrmicro2537

PubMed Abstract | CrossRef Full Text | Google Scholar

Guerrero, R., Margulis, L., and Berlanga, M. (2013). Symbiogenesis: the holobiont as a unit of evolution. Int. Microbiol. 16, 133–143. doi: 10.2436/20.1501.01.188

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, L., Zhao, G., Xu, J. R., Kistler, H. C., Gao, L., and Ma, L.-J. (2016). Compartmentalized gene regulatory network of the pathogenic fungus fusarium graminearum. New Phytol. 211, 527–541. doi: 10.1111/nph.13912

PubMed Abstract | CrossRef Full Text | Google Scholar

Halfvarson, J., Brislawn, C. J., Lamendella, R., Vázquez-Baeza, Y., Walters, W. A., Bramer, L. M., et al. (2017). Dynamics of the human gut microbiome in inflammatory bowel disease. Nat. Microbiol. 2, 1–7. doi: 10.1038/nmicrobiol.2017.4

PubMed Abstract | CrossRef Full Text | Google Scholar

Hameya, F. K., Nestorowaa, S., Kinstona, S. J., Kenta, D. G., Wilsona, . N. K., and Göttgens, B. (2017). Reconstructing blood stem cell regulatory network models from single-cell molecular profiles. Proc. Natl. Acad. Sci. U.S.A. 114, 5822–5829. doi: 10.1073/pnas.1610609114

CrossRef Full Text | Google Scholar

Hooper, L. V., Littman, D. R., and Macpherson, A. J. (2012). Interactions between the microbiota and the immune system. Science 336, 1268–1273. doi: 10.1126/science.1223490

PubMed Abstract | CrossRef Full Text | Google Scholar

Hooper, L. V., Wong, M. H., Thelin, A., Hansson, L., Falk, P. G., and Gordon, J. I. (2001). Molecular analysis of commensal host-microbial relationships in the intestine. Science 291, 881–884. doi: 10.1126/science.291.5505.881

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, S., Eichler, G., Bar-Yam, Y., and Ingber, D. E. (2005). Cell fates as high-dimensional attractor states of a complex gene regulatory network. Phys. Rev. Lett. 94, 128701. doi: 10.1103/PhysRevLett.94.128701

PubMed Abstract | CrossRef Full Text | Google Scholar

Kauffman, S. (1969a). Homeostasis and differentiation in random genetic control networks. Nature 224, 177–178. doi: 10.1038/224177a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Kauffman, S. A. (1969b). Metabolic stability and epigenesis in randomly constructed genetic nets. J. Theor. Biol. 22, 437–467. doi: 10.1016/0022-5193(69)90015-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Kinouchi, O., and Copelli, M. (2006). Optimal dynamical range of excitable networks at criticality. Nat. Phys. 2, 348–351. doi: 10.1038/nphys289

CrossRef Full Text | Google Scholar

Krajmalnik-Brown, R., Ilhan, Z.-E., Kang, D.-W., and DiBaise, J. K. (2012). Effects of gut microbes on nutrient absorption and energy regulation. Nutr. Clin. Pract. 27, 201–214. doi: 10.1177/0884533611436116

PubMed Abstract | CrossRef Full Text | Google Scholar

Laland, K., Uller, T., Feldman, M., Sterelny, K., Müller, G. B., Moczek, A., et al. (2014). Does evolutionary theory need a rethink? Nature 514, 161–164. doi: 10.1038/514161a

PubMed Abstract | CrossRef Full Text | Google Scholar

Langton, C. G. (1990). Computation at the edge of chaos: phase transitions and emergent computation. Physica D 42, 12–37. doi: 10.1016/0167-2789(90)90064-V

CrossRef Full Text | Google Scholar

Ley, R. E. (2010). Obesity and the human microbiome. Curr. opin. Gastroenterol. 26, 5–11. doi: 10.1097/MOG.0b013e328333d751

PubMed Abstract | CrossRef Full Text | Google Scholar

Ley, R. E., Hamady, M., Lozupone, C., Turnbaugh, P. J., Ramey, R. R., Bircher, J. S., et al. (2008). Evolution of mammals and Their Gut Microbes. Science 320, 1647–1651. doi: 10.1126/science.1155725

PubMed Abstract | CrossRef Full Text | Google Scholar

Ley, R. E., Peterson, D. A., and Gordon, J. I. (2006a). Ecological and evolutionary forces shaping microbial diversity in the human intestine. Cell 124, 837–848. doi: 10.1016/j.cell.2006.02.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Ley, R. E., Turnbaugh, P. J., Klein, S., and Gordon, J. I. (2006b). Microbial ecology: human gut microbes associated with obesity. Nature 444, 1022. doi: 10.1038/4441022a

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, F., Long, T., Lu, Y., Ouyang, Q., and Tang, C. (2004). The yeast cell-cycle network is robustly designed. Proc. Natl. Acad. Sci. U.S.A. 101, 4781–4786. doi: 10.1073/pnas.0305937101

PubMed Abstract | CrossRef Full Text | Google Scholar

Lloyd-Price, J., Abu-Ali, G., and Huttenhower, C. (2016). The healthy human microbiome. Genome Med. 8, 1–11. doi: 10.1186/s13073-016-0307-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Lynch, M. (2010). Evolution of the mutation rate. Trends Genet. 26, 345–352. doi: 10.1016/j.tig.2010.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Lynch, M., Ackerman, M. S., Gout, J. F., Long, H., Sung, W., Thomas, W. K., et al. (2016). Genetic drift, selection and the evolution of the mutation rate. Nat. Rev. Genet. 17, 704–714. doi: 10.1038/nrg.2016.104

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, S., Snyder, M., and Dinesh-Kumar, S. P. (2017). Discovery of novel human gene regulatory modules from gene co-expression and promoter motif analysis. Sci. Rep. 7, 5557. doi: 10.1038/s41598-017-05705-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Mai, V. (2004). Dietary modification of the intestinal microbiota. Nutr. Rev. 62, 235–242. doi: 10.1111/j.1753-4887.2004.tb00045.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Manor, O., Levy, R., and Borenstein, E. (2014). Mapping the inner workings of the microbiome: genomic- and metagenomic-based study of metabolism and metabolic interactions in the human microbiome. Cell Metab. 20, 742–752. doi: 10.1016/j.cmet.2014.07.021

PubMed Abstract | CrossRef Full Text | Google Scholar

McCabe, L., Britton, R. A., and Parameswaran, N. (2015). Prebiotic and probiotic regulation of bone health: role of the intestine and its microbiome. Curr. Osteoporos. Rep. 13, 363–371. doi: 10.1007/s11914-015-0292-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Moeller, A. H., Li, Y., Ngole, E. M., Ahuka-Mundeke, S., Lonsdorf, E. V., Pusey, A. E., et al. (2014). Rapid changes in the gut microbiome during human evolution. Proc. Natl. Acad. Sci. U.S.A. 111, 16431–16435. doi: 10.1073/pnas.1419136111

PubMed Abstract | CrossRef Full Text | Google Scholar

Moran, N. A., and Sloan, D. B. (2015). The hologenome concept: helpful or hollow? PLoS Biol. 13:e1002311. doi: 10.1371/journal.pbio.1002311

CrossRef Full Text | Google Scholar

Morgan, X. C., Tickle, T. L., Sokol, H., Gevers, D., Devaney, K. L., Ward, D. V., et al. (2012). Dysfunction of the intestinal microbiome in inflammatory bowel disease and treatment. Genome Biol. 13, R79. doi: 10.1186/gb-2012-13-9-r79

PubMed Abstract | CrossRef Full Text | Google Scholar

Morris, J. J., Lenski, R. E., and Zinser, E. R. (2012). The black queen hypothesis: evolution of dependencies through adaptative gene loss. Mbio 3, e00036–12. doi: 10.1128/mBio.00036-12

CrossRef Full Text | Google Scholar

Motta, S. S., Cluzel, P., and Aldana, M. (2015). Adaptive resistance in bacteria requires epigenetic inheritance, genetic noise, and cost of efflux pumps. PLoS ONE 10:e0118464. doi: 10.1371/journal.pone.0118464

PubMed Abstract | CrossRef Full Text | Google Scholar

Nicholson, J. K., Holmes, E., Kinross, J., Burcelin, R., Gibson, G., Jia, W., et al. (2012). Host-gut microbiota metabolic interactions. Science 336, 1262–1268. doi: 10.1126/science.1223813

PubMed Abstract | CrossRef Full Text | Google Scholar

Nykter, M., Price, N. D., Larjo, A., Aho, T., Kauffman, S. A., Yli-Harja, O., et al. (2008). Critical networks exhibit maximal information diversity in structure-dynamics relationships. Phys. Rev. Lett. 100, 058702. doi: 10.1103/PhysRevLett.100.058702

PubMed Abstract | CrossRef Full Text | Google Scholar

Oliveri, P., Tu, Q., and Davidson, E. H. (2008). Global regulatory logic for specification of an embryonic cell lineage. Proc. Natl. Acad. Sci. U.S.A. 105, 5955–5962. doi: 10.1073/pnas.0711220105

PubMed Abstract | CrossRef Full Text | Google Scholar

Pérez-Landero, S., Sandoval-Motta, S., Martínez-Anaya, C., Yang, R., Folch-Mallol, J. L., Martínez, L. M., et al. (2015). Complex regulation of hsf1-skn7 activities by the catalytic subunits of pka in saccharomyces cerevisiae: experimental and computational evidences. BMC Syst. Biol. 9:42. doi: 10.1186/s12918-015-0185-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Queller, D. C., and Strassmann, J. E. (2016). Problems of multi-species organisms: endosymbionts to holobionts. Biol. Philos. 31, 855–873. doi: 10.1007/s10539-016-9547-x

CrossRef Full Text | Google Scholar

Rawls, J. F., Samuel, B. S., and Gordon, J. I. (2004). Gnotobiotic zebrafish reveal evolutionarily conserved responses to the gut microbiota. Proc. Natl. Acad. Sci. U.S.A. 101, 4596–4601. doi: 10.1073/pnas.0400706101

PubMed Abstract | CrossRef Full Text | Google Scholar

Resendis-Antonio, O., Hernández, M., Mora, Y., and Encarnación, S. (2012). Functional modules, structural topology, and optimal activity in metabolic networks. PLoS Comput. Biol. 8:e1002720. doi: 10.1371/journal.pcbi.1002720

PubMed Abstract | CrossRef Full Text | Google Scholar

Rogers, G., Keating, D., Young, R., Wong, M., Licinio, J., and Wesselingh, S. (2016). From gut dysbiosis to altered brain function and mental illness: mechanisms and pathways. Mole. Psychiatry 21, 738. doi: 10.1038/mp.2016.50

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenberg, E., Sharon, G., Atad, I., and Zilber-Rosenberg, I. (2010). The evolution of animals and plants via symbiosis with microorganisms. Environ. Microbiol. Rep. 2, 500–506. doi: 10.1111/j.1758-2229.2010.00177.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenberg, E., and Zilber-Rosenberg, I. (2016). Microbes drive evolution of animals and plants: the hologenome concept. mBio 7, e01395–15. doi: 10.1128/mBio.01395-15

PubMed Abstract | CrossRef Full Text | Google Scholar

Roughgarden, J., Gilbert, S. F., Rosenberg, E., Zilber-Rosenberg, I., and Lloyd, E. A. (2018). Holobionts as units of selection and a model of their population dynamics and evolution. Biol. Theory 13, 44–65. doi: 10.1007/s13752-017-0287-1

CrossRef Full Text | Google Scholar

Rupnik, M., Wilcox, M. H., and Gerding, D. N. (2009). Clostridium difficile infection: New developments in epidemiology and pathogenesis. Nat. Rev. Microbiol. 7, 526–536. doi: 10.1038/nrmicro2164

PubMed Abstract | CrossRef Full Text | Google Scholar

Sachs, J. L., and Hollowell, a. C. (2012). The origins of cooperative bacterial communities. mBio 3, e00099–12. doi: 10.1128/mBio.00099-12

PubMed Abstract | CrossRef Full Text | Google Scholar

Sachs, J. L., Skophammer, R. G., Bansal, N., and Stajich, J. E. (2014). Evolutionary origins and diversification of proteobacterial mutualists. Proc. R. Soc. B Biol. Sci. 281, 20132146. doi: 10.1098/rspb.2013.2146

PubMed Abstract | CrossRef Full Text | Google Scholar

Sagan, L. (1967). On the origin of mitosing cells. J. Theor. Biol. 14, 225IN1–274IN6.

PubMed Abstract | Google Scholar

Sandoval-Motta, S., Aldana, M., and Frank, A. (2017). Evolving ecosystems: Inheritance and selection in the light of the microbiome. Arch. Med. Res. 48, 780–789. doi: 10.1016/j.arcmed.2018.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Sears, C. L., and Garrett, W. S. (2014). Microbes, microbiota, and colon cancer. Cell Host Microbe 15, 317–328. doi: 10.1016/j.chom.2014.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Seekatz, A. M., and Young, V. B. (2014). Clostridium difficile and the microbiota. J. Clin. Invest. 124, 4182–4189. doi: 10.1172/JCI72336

PubMed Abstract | CrossRef Full Text | Google Scholar

Serra, R., Villani, M., Graudenzi, A., and Kauffman, S. (2007). Why a simple model of genetic regulatory networks describes the distribution of avalanches in gene expression data. J. Theor. Biol. 246, 449–460. doi: 10.1016/j.jtbi.2007.01.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharpton, T. J. (2018). Role of the gut microbiome in vertebrate evolution. mSystems 3:e00174–17. doi: 10.1128/mSystems.00174-17

PubMed Abstract | CrossRef Full Text | Google Scholar

Shin, S. C., Kim, S.-H., You, H., Kim, B., Kim, A. C., Lee, K.-A., et al. (2011). Drosophila microbiome modulates host developmental and metabolic homeostasis via insulin signaling. Science 334, 670–674. doi: 10.1126/science.1212782

PubMed Abstract | CrossRef Full Text | Google Scholar

Shmulevich, I., Kauffman, S. A., and Aldana, M. (2005). Eukaryotic cells are dynamically ordered or critical but not chaotic. Proc. Natl. Acad. Sci. U.S.A. 102, 13439–13444. doi: 10.1073/pnas.0506771102

CrossRef Full Text | Google Scholar

Spor, A., Koren, O., and Ley, R. (2011). Unravelling the effects of the environment and host genotype on the gut microbiome. Nat. Rev. Microbiol. 9, 279. doi: 10.1038/nrmicro2540

PubMed Abstract | CrossRef Full Text | Google Scholar

Stern, M. D. (1999). Emergence of homeostasis and “noise imprinting” in an evolution model. Proc. Natl. Acad. Sci. U.S.A. 96, 10746–10751.

PubMed Abstract | Google Scholar

Taxis, T. M., Wolff, S., Gregg, S. J., Minton, N. O., Zhang, C., Dai, J., et al. (2015). The players may change but the game remains: network analyses of ruminal microbiomes suggest taxonomic differences mask functional similarity. Nucleic Acids Res. 43, 9600–9612. doi: 10.1093/nar/gkv973

PubMed Abstract | CrossRef Full Text | Google Scholar

Theis, K. R., Dheilly, N. M., Klassen, J. L., Brucker, R. M., Baines, J. F., Bosch, T. C. G., et al. (2016). Getting the hologenome concept right: an eco-evolutionary framework for hosts and their microbiomes. mSystems 1, e00028–16. doi: 10.1128/mSystems.00028-16

PubMed Abstract | CrossRef Full Text | Google Scholar

Torres-Sosa, C., Huang, S., and Aldana, M. (2012). Criticality is an emergent property of genetic networks that exhibit evolvability. PLoS Comput. Biol. 8:e1002669. doi: 10.1371/journal.pcbi.1002669

PubMed Abstract | CrossRef Full Text | Google Scholar

Tropini, C., Earle, K. A., Huang, K. C., and Sonnenburg, J. L. (2017). The gut microbiome: connecting spatial organization to function. Cell Host Microbe 21, 433–442. doi: 10.1016/j.chom.2017.03.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Turnbaugh, P. J., and Stintzi, A. (2011). Human health and disease in a microbial world. Front. Microbiol. 2:190. doi: 10.3389/fmicb.2011.00190

PubMed Abstract | CrossRef Full Text | Google Scholar

van Nood, E., Vrieze, A., Nieuwdorp, M., Fuentes, S., Zoetendal, E. G., de Vos, W. M., et al. (2013). Duodenal infusion of donor feces for recurrent clostridium difficile. New Engl. J. Med. 368, 407–415. doi: 10.1056/NEJMoa1205037

PubMed Abstract | CrossRef Full Text | Google Scholar

van Opstal, E. J., and Bordenstein, S. R. (2015). Rethinking heritability of the microbiome. Science 349, 1172–1173. doi: 10.1126/science.aab3958

PubMed Abstract | CrossRef Full Text | Google Scholar

Wagner, G. P. (2007). The developmental genetics of homology. Nat. Rev. Genet. 8, 473–479. doi: 10.1038/nrg2099

PubMed Abstract | CrossRef Full Text | Google Scholar

West, S. A., Fisher, R. M., Gardner, A., and Kiers, E. T. (2015). Major evolutionary transitions in individuality. Proc. Natl. Acad. Sci. U.S.A. 112, 10112–10119. doi: 10.1073/pnas.1421402112

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilson, I. D., and Nicholson, J. K. (2017). Gut microbiome interactions with drug metabolism, efficacy, and toxicity. Transl. Res. 179, 204–222. doi: 10.1016/j.trsl.2016.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J., Tan, Q., Fu, Q., Zhou, Y., Hu, Y., Tang, S., et al. (2017). Gastrointestinal microbiome and breast cancer: correlations, mechanisms and potential clinical implications. Breast Cancer 24, 220–228. doi: 10.1007/s12282-016-0734-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Yatsunenko, T., Rey, F. E., Manary, M. J., Trehan, I., Dominguez-Bello, M. G., Contreras, M., et al. (2012). Human gut microbiome viewed across age and geography. Nature 486, 222. doi: 10.1038/nature11053

PubMed Abstract | CrossRef Full Text | Google Scholar

Zackular, J. P., Baxter, N. T., Iverson, K. D., Sadler, W. D., Petrosino, J. F., Chen, G. Y., et al. (2013). The gut microbiome modulates colon tumorigenesis. MBio 4, e00692–13. doi: 10.1128/mBio.00692-13

PubMed Abstract | CrossRef Full Text | Google Scholar

Zilber-Rosenberg, I., and Rosenberg, E. (2008). Role of microorganisms in the evolution of animals and plants: the hologenome theory of evolution. FEMS Microbiol. Revi. 32, 723–735. doi: 10.1111/j.1574-6976.2008.00123.x

PubMed Abstract | CrossRef Full Text | Google Scholar