We developed a kinetic model to study metabolic complementation endosymbiont consortium of two different bacteria, when they are allowed to exchange intermediate metabolites of a hypothetical biosynthetic pathway. The model describes a consortium consisting of two different endosymbiotic bacterial species and their host (Figure 1). Each of the two bacterial populations is represented by a large number (N_i with i = 1, 2) of identical cells, so that in practice it is sufficient to consider two different compartments, modeling a single cell of each bacterial strain, enclosed in a larger compartment, which allows for exchanges of metabolites between cells and the host. Both endosymbionts code for the same hypothetical biosynthetic pathway, which allows the conversion of a precursor S (produced by the bacteria), to an end product P, which is essential for any of the three members of the consortium. The biosynthetic pathway is a linear chain composed by two reactions (see Figure 1), with the precursor S being first converted into an intermediate metabolite X, and then X being converted to the final product P.

The intermediate metabolite X and the final product P can freely diffuse outside the cells (Figure 1). Moreover, since the product P is considered essential for both endosymbionts and its host, three demand fluxes are also introduced to represent the incorporation of P to the respective biomasses of the three member of the consortium. Finally, product inhibition has been suggested to be present in some endosymbionts, e.g., anthranilate synthase of B. aphidicola shows a highly conserved allosteric binding site (Lai et al., 1994), albeit in other enzymes like prephenate dehydratase the inhibition could be mitigated by point mutations in the allosteric site (Jiménez et al., 2000). Thus, the model includes a negative feedback inhibition mechanism that regulates the production of P.

The two reactions in the biosynthetic pathway, S →X and X→P, are modeled using reversible Michaelis–Menten kinetics (Noor et al., 2013); we name the corresponding fluxes V_S and V_X, respectively. Moreover, in order to include the effect of product inhibition on the first step (S→X), a term representing such effect was introduced in the kinetic equation of the first reaction (Cornish-Bowden, 2013). There are diverse mechanisms of enzyme inhibition, including allosteric interactions, but in this approach we have analyzed one of the simplest situations and the kinetics of the inhibition was considered, in a first place, to be uncompetitive, although the case of competitive inhibition is discussed later. Thus, the kinetic law of the first reaction S→X, with flux V_S, is the following: (1)VS = kcat[ES]1 +[P]/KI·([S]−[X]/Keq)/KM+1 +[S]/KM+ + [X]/KM− where k_cat is the turnover number, or catalytic constant, of the reaction, KM+ (KM-) is the forward (backward) Michaelis constant, K_eq is the equilibrium constant, and K_I is the inhibition constant. It is worth to note that in the limit 1/K_I → 0 the reversible Michaelis–Menten kinetic law, without inhibition, is recovered. The kinetic for the second reaction (X→P, with flux V_X) is analogous to Equation (1) but without the inhibition term ([P]/K_i).

In order to allow the flux of metabolites between the members of the consortium, the metabolites X and P are allowed to be exchanged with the surrounding environment. Endosymbiotic bacteria usually lack transport systems (Charles et al., 2011); thus, the exchanges of X and P are modeled using first-order, non-saturating, reversible kinetic laws. The fluxes U_α take the following form: (2)Uα=Dα([α]in−[α]out), α= X,P where D_α are related to membrane permeability and diffusion constants of the metabolites (Laidler and Shuler, 1949). Note that the fluxes are assumed to be positive when the metabolites are excreted by the cells. On the other hand, if there were any kind of protein-mediated transport, the corresponding membrane proteins should also have been taken into account. Finally, the product P is assumed to be consumed during growth by the three members of the consortium, with consumption fluxes in each compartment. For the bacterial cells, we model the consumption of the metabolite P using an irreversible Michaelis–Menten reaction: (3)VP = kcat[EP][P][P] + KM

The fluxes V_{P, i} (the index i = 1, 2 corresponding to the cell type) are constrained to be larger than some fixed flux J_{P, i}. The extracellular product is instead consumed by the host. It is useful to normalize this rate to the size of the bacterial populations N = N₁ + N₂, so that the consumption rate equals d[P]_out/dt = −NJ_{P, 0}. The three parameters J_{P, 0} and J_{P, i} thus set the demand of the metabolite P from the bacteria or the host.

Finally, the size of the populations of the two bacteria naturally influence both the total protein concentration and the magnitude of the transport fluxes between the cells and the host, and cannot be disregarded. It is convenient to introduce the normalized populations n₁ = N₁/N and n₂ = N₂/N, with N₁ + N₂ = N. We will therefore first show results obtained by fixing the cell relative populations to be n₁ = n₂ = 1/2, and then exploring the case in which the two population fractions differ. When the whole model is considered, the balance equations for all the metabolites of the system can be written as follows: (4)d[Xi]dt = VS,i−VX,i−UX,i, i=1,2d[Pi]dt = VX,i−VP,i−UP,i, i=1,2d[X0]dt=N(n1UX,1+n2UX,2)d[P0]dt=N(n1UP,1+n2UP,2−JP,0) where N corresponds to the total number of cells in both populations. Is worth to note that the equations for intracellular metabolites are symmetric for the two species, while the balance equations of the extracellular metabolites depend on the size of each symbiont population and, in the case of P₀, on the host demand flux J_{P, 0}. The substrate S is a precursor produced by each of the two cell types, and its concentration is assumed to be fixed (Heinrich and Schuster, 1996). As a consequence, the balance equation of S is not included in the model. For the sake of simplicity, in the following we will set all kinetic parameters to one, unless stated otherwise. Also, we will set the population fractions to n₁ = n₂ = 1/2 and the concentration of the substrate to [S] = 10.

Endosymbiont bacteria involved in nutritional interactions overproduce some essential products, such as essential amino acids or vitamins, and share them with their host (and in some cases with other endosymbionts; Douglas, 1992). On the other hand, as discussed in the introduction, such over-production can be ensured by tuning the enzyme concentrations along the pathway, even if end-product inhibition may not allow to tune the flux easily (Lai et al., 1994; Rouhbakhsh et al., 1997; Baumann et al., 1999). Accordingly, improvements in the efficiency of over-productions, in terms of the economy on protein synthesis, should provide of some advantage to the whole system (Dekel and Alon, 2005; Ponce de León et al., 2008; Shachrai et al., 2010; Goelzer and Fromion, 2011; Flamholz et al., 2013; Kafri et al., 2015; Mori et al., 2016). Thus, in this work we will focus on the trade-off between the maximization of the flux through a pathway and the efficient allocation of the proteins. Using constraint-based optimization this trade-off can be modeled as two complementary optimization problems: (i) maximization of the flux subject to a constraint on the maximum amount of total enzyme concentration; and (ii) minimization of enzyme concentrations subject to a constraint (lower bound) on the flux of interest, which ensure that a certain demand is satisfied.

Herein, we will adopt the second formulation, i.e., enzyme minimization subject to a certain demand flux of the product P, although one can show that both formulations yield completely analogous results (see Supplementary Text 1). The optimization problem is defined in the following way. Given the demand fluxes J_{P, i}, i = 0, 1, 2: (5)VP,i ≥ JP,i ∀i= 0, 1, 2 we compute the metabolite and enzyme concentrations which minimize the total enzyme concentration (normalized by the total bacterial population), defined as follows: (6)C=∑i = 1,2∑α = S,X,Pni[Eα,i] where the cost or objective function C correspond to the total amount of protein invested by both cell types in order to produce an amount of P enough to satisfy the demand required by the three members of the consortium. All the intermediate metabolites are subject to the steady-state condition, which is imposed by setting to zero the time derivatives in (Equation 4). Accordingly, each solution will consist of a set of metabolites and enzymes concentrations, from which one can compute the optimal fluxes using the kinetic laws. A complementary formulation of the optimization problem is possible, i.e., maximization of the flux subject to a constraint on the maximum amount of total enzyme concentration. This other approach can be shown to be completely equivalent to the “enzyme” minimization approach, as shown in Supplementary Note 1.

In order to highlight the patterns of proteome allocation between the cell populations, we introduce the measure of enzyme asymmetries A_α, defined for each enzyme α as follows: (7)Aα=[Eα,1]−[Eα,2][Eα,1]+[Eα,2],α=S,X,P

This quantity reflects the relative difference in the concentration of each enzyme between the two different endosymbiont. Thus, A_α will be zero if the enzyme concentrations in the two cells are the same, while it will equal either 1 or −1 if one cell type completely lacks the corresponding enzyme. In particular, if the two type of cells have the same sets of kinetic parameters, one would expect that the protein asymmetries would be equal to 0, unless some metabolic interaction emerges.

Optimizations were solved using the fmincon function in Matlab 2015a forLinux. Due to the non-linear formulation of the optimization problem, the cost function exhibits multiple local minimums. Therefore, local optimization strategies such that gradient descent are not always able to find the correct minimum. However, if the parameters of the model are varied smoothly, such as in the case of the permeabilities D_X and D_P in Figure 2, one may rely on the fact that the optimal solution should also vary smoothly (excluding the frontiers between the symmetric and the non-symmetric solutions). Therefore, for each combination of parameters, we optimized the cost function starting from a random seed until we got the correct minimum in all cases.

In the sensitivity analysis, we chose a point in the parameter space, roughly corresponding to the border between the symmetric and asymmetric region Then, we randomized by various amounts the kinetic parameters, and for each set of parameters we computed the optimal concentrations, {x^*}, and the corresponding optimal E_S protein asymmetry, A_S = ([E_{1, S}] − [E_{2, S}]) / ([E_{1, S}] + [E_{2, S}]), where the asterisk indicates the optimal solution. Then, for each parameter k_i, we compute the following correlation: (8)Ci=〈AS(x∗)·ki−AS(x∗)·ki〉 where the angled brackets 〈·〉 represent the average over the N samples. In order to ensure that the optimal concentration corresponded to real absolute minima of the cost function, an iterative approach was developed, which is described in detail in Supplementary Text 2. The results of the sensitivity analysis are instead in Supplementary Table 1.

Amino acid sequences of the anthranilate synthase component I (trpE) from twelve different Buchnera strains available, as well as E. coli (strain K12) and S. marcescens, were downloaded from the UniProtKB database (UniProt-Consortium, 2014). The sequences in fasta format can be found in the Supplementary file Data 1. The multiple sequences alignment was performed using T-cofee package using default parameters (Notredame et al., 2000). The alignment representation was done using Jalview (Waterhouse et al., 2009).

A precise quantification of the diffusion constants and membrane permeability for any intermediate metabolite is not a trivial task. Therefore, we used different parameters as estimators of the permeability or diffusion capacity of the intermediate metabolites in different amino acid biosynthetic pathways we analyzed. For each compound we considered: (i) the molecular weight; (ii) number of hydrogen bond acceptors and donors; (iii) the octanol-water partition coefficient; and (iv) the topological polar surface area (TPSA). The structures of the amino acids biosynthetic pathways were retried from MetaCyc (Caspi et al., 2014). The list of pathways analized includes the super-pathway of aromatic amino acid biosynthesis (L-Trp, L-Tyr, L-Phe) and the super-pathway of branched chain amino acid biosynthesis (L-Leu, L-Ile, L-Val). Furthermore, all the physicochemical parameters listed above were computed using OpenBabel (O'Boyle et al., 2011), [accessed through the Python wrapper PyBel O'Boyle et al., 2008] with the exception of the logarithm of the octanol-water partition coefficient (LogP) for which two different estimators were chosen, named ALogP (Tetko and Poda, 2004) and XLOP3 (Cheng et al., 2007). Results can be found in Supplementary Table 2.

Using the model and the optimization criteria introduced in previous sections, we now turn on the analysis of the optimal proteome allocation obtained, for different scenarios. As a starting point, we will set the host demand flux to zero, J_{P, 0} = 0. As explained in Supplementary Note 1, since the total enzyme concentrations is proportional to the demand fluxes, the sum of the demand fluxes J_{P, 1} + J_{P, 2} is not relevant for the optimal strategy of proteome allocation, and only their ratio J_{P, 1}/J_{P, 2} is (as long as J_{P, 0} = 0). We will therefore choose J_{P, i} = 1 for both the two cell types. (The results for different values of J_{P, 0} = 0 and J_{P, 1}/J_{P, 2} are shown below).

We show in Figure 2 the optimal protein asymmetries obtained by varying the value of the diffusion constants, D_X and D_P, for different values of the inhibition constant K_I. Given the symmetry of the problem, i.e., the reactions in the two cell types have identical kinetic parameters, we will only consider the solutions with A_S ≥ 0. In the extreme case where D_X = D_P = 0, the two cells cannot interact in any way and no complementation can arise. However, when the metabolites are allowed to permeate, the two cell populations are in principle allowed to interact, and an asymmetric solution corresponding to metabolic complementation may emerge. If the effect of the product inhibition is too low (e.g., 1/K_I ≤ 0.1), then the optimal solution is always characterized by a symmetric enzyme allocation, i.e., A_α = 0, for all enzymes involved (Figure 2A), and for all values of the permeability constants D_X and D_P. Conversely, when the inhibition gets larger and larger (e.g., 1/K_I ≥ 1), asymmetric solutions starts to emerge as the optimal strategy, depending on the values of the diffusion constants D_X and D_P (Figures 2B–D). For a given value of 1/K_I, the diffusion constants define a transition line, which divides the plane into two zones, corresponding to the range of parameter for which the symmetric or asymmetric solution is optimal. When the condition favors the asymmetric strategy, in general all the enzymes exhibit certain degree of asymmetry. In particular in the case of the enzyme E_S, the asymmetry only takes the extreme values, A_S = 0 or |A_S| = 1, whereas for the other two enzymes the degree of asymmetry depends on the diffusion constants. Interestingly, we find that large values of D_X promote the asymmetry, while for large values of D_P the symmetry is restored. In other words, asymmetric solutions arise when D_X/D_P is large enough. Summarizing, asymmetric proteome allocation becomes an optimal strategy in the presence of by product inhibition and when the product diffuses slower than the intermediate metabolite X. If P diffuses too fast, the concentrations within both cell types will be roughly the same, and therefore the inhibition cannot be by-passed. Conversely, if P diffuses much slower than X, its concentration on the cell coding for the first enzyme (i.e., E_S) will be low and its contribution as an inhibitor will be negligible.

As shown in Figure 2, in the regions where asymmetry emerges as the optimal solution, we have the following relations for enzyme asymmetries: A_S = 1, A_X < 0, and A_P > 0. Thus, in order to have a deeper insight on the relation between diffusion and the sudden change from symmetric to asymmetric solution, we focus now on single value of the inhibition constant (1/K_I = 10). We can individually study the symmetric and asymmetric solutions by imposing that A_S = 0 (in the symmetric case) or A_S = 1 (in the asymmetric case). Then, for these two cases, the optimal proteome allocation was computed for different values of the diffusion constants, D_X and D_P.

Figure 3A shows that the protein cost of the asymmetric strategy grows when the metabolites cannot easily diffuse across the membranes, i.e., when the values of D_X and D_P decrease (yellow surface/lines). On the other hand, the asymmetric strategy becomes optimal for large values of the diffusion rates. For large values of the diffusion constants we see that the asymmetric solution is optimal, with the A_X enzyme asymmetry approaching the value −1 (Figure 3B). In this conditions, the pathway is completely split between the two cells, with one enzyme (E_S) only present in the first cell (Figure 3B), and the second enzyme (E_X) only present in the other cell (Figure 3C). We see in Figure 3D that the asymmetric condition, A_S = 1, is always accompanied by an excretion of the X metabolite from the first cell, which is gathered by the second cell U_{X, 1} = −U_{X, 2} > 0. Interestingly, the exchange of X (Figure 3D) still takes place even in the case in which P cannot cross the cell membrane, i.e., D_P = 0 (Figure 3E). In this condition, the second population becomes a cheater, and takes advantage of the leakage of X from the first population. Nevertheless, if the diffusion constant D_P is increased, the second cell starts sharing the product P with the first cell, thus establishing a metabolic complementation relationship. Finally, if D_P is too large, the concentrations of the P metabolites are very similar in both cells, so that there is no way to reduce the effect of inhibition by partitioning the metabolic route. Therefore, the cell switches from metabolic complementation to a “neutralist” strategy, with completely symmetric enzyme allocation and no exchange of metabolites between the two cells types.

In this section, we extend our results to more general settings. First, we check that our results do not depend qualitatively on the inhibition kinetics. Instead of a non-competitive mechanism (Equation 1), we use a competitive mechanism (Cornish-Bowden, 2013) so that the concentration of the product metabolites affects the saturation of the enzyme, and not the turnover rate. Although this mechanism is more complicated than the uncompetitive case, due to the interplay between the different concentrations of the S, X, and P metabolites, the results obtained are analogous to the simpler case of uncompetitive inhibition (see Supplementary Figure 1). This suggests that while product inhibition is a necessary condition for the emergence of metabolic complementation, it does not depend on specific inhibition mechanisms, and the metabolic complementation strategy is favored by strong feedback inhibition.

Then, we probed the effect of different demand fluxes, both to the host (J_{P, 0}) and to the bacterial species (J_{P, 1} and J_{P, 2}). As one increases the host demand flux J_{P, 0}, the qualitative pattern of the different optimal strategies is only slightly affected by the presence of a flux of products to the host (see Supplementary Figure 2), even though the total amount of protein to be allocated clearly increases with the demand flux. Similarly, we consider the case where the two bacterial populations have different demands for the product P, for instance because of a different biomass composition. We see that an asymmetry in the demands for P can promote, even more, the asymmetry in the proteome allocation (see Supplementary Figure 3); for large enough 1/K_I and demand flux asymmetry ρ_J = (J_{P, 1} − J_{P, 2}) / (J_{P, 1} + J_{P, 2}), the optimal solution is characterized by a complete splitting of the pathway within the two cells, [E_{S, 1}] = [E_{X, 2}] = 0.

A similar effect is found when the two bacterial species have different populations, that is, n₁ ≠ n₂. As n₁ is moved away from ½(with n₂ = 1 − n₁), the protein asymmetry is generally enhanced, with larger asymmetries A_X (see Supplementary Figure 4). All these results show that the complementation patterns described by the model are robust against variation in the demand fluxes and the relative population sizes.

The number of kinetic parameters in the model is quite large, and it is very difficult to probe the effect of perturbations for each different combination of parameters. Therefore, we performed a sensitivity analysis, on the model with uncompetitive inhibition, by randomly perturbing the parameters, and studying how the perturbation of each parameter correlates with the insurgence of the asymmetric optimal solution (See Section Materials and Methods for details). We see that the parameters which promote the enzyme asymmetry most are, as expected, the inverse inhibition constant 1/K_I (reflecting the amount of inhibition), the concentration of the S metabolite, and the diffusion constant D_X of the X metabolite (Supplementary Table 1). We also find that the equilibrium constant of the S→X reaction also positively affects the asymmetry, albeit slightly; this is because increased concentrations of the X metabolite allows for larger shuttling fluxes across the two cells. On the other hand, the diffusion constant D_P of the P metabolite is found to promote the symmetry between the two cells, decreasing the average value of A_S.

In the previous section, we gave an in depth characterization of a model in which the biosynthetic pathway is composed by two reactions. However, real biosynthetic pathways are usually composed by a larger number of reactions. In order to study how the permeability affects the emergence of a complementation, an extended version of the model was developed, which keeps the fundamental features of its simplest version, while adding more intermediates in the pathway. The extended version of the model was obtained by introducing five intermediate metabolite, X₁, …, X₅, along with five enzymes EX₁, …, EX₅. Each enzyme EX_i, with the index i ranging between 1 and 5, catalyses the reversible conversion of X_i into X_i+1 or, in the case of EX₅, the product metabolite P. This extended model allows us to study how a long biosynthetic pathway is differently partitioned within the two cells for different permeabilities of the intermediate metabolites.

In the extreme case in which only one intermediate metabolite can permeate through the cell membranes, we find that the general behavior of the extended model is totally consistent with the results obtained using the simplified version (See Emergence of Metabolic Complementation Section). As shown in Figure 4, the pathway is split in correspondence of the most permeable intermediate metabolite. The first cell converts the substrate S into the intermediate metabolite with large permeability, which in turn is (mostly) shuttled to the other cell and converted into the product P, which is then shared between the two bacteria. A wide range of behaviors are obtained when different, more complex, scenarios are considered (see Supplementary Figure 5). For example, when all the intermediate metabolites can permeate across the membranes, the optimal solution still exhibits metabolite complementation. However, enzymes levels change smoothly along the pathway with many metabolite being exchanged, and therefore the qualitative pattern lacks a unique breaking point.

In this section, we discuss the well-documented case of the complementation of the tryptophan biosynthetic pathway between the members of an endosymbiotic consortium in relation with the results here obtained (Gosalbes et al., 2008; Sloan and Moran, 2012). This event of metabolic complementation has been reported to be convergent in two different insect systems and evolutionary unrelated. On one hand, the complementation of the tryptophan biosynthesis has been described in the endosymbiotic consortium composed by B. aphidicola and S. symbiotica SCc hosted by the cedar aphid (Gosalbes et al., 2008), and also the phylogenetically related aphid T. salignus (Manzano-Marín et al., 2016). On the other hand, the complementation has also been described in the system of the psyllid H. cubana, where the coding genes for the tryptophan pathway are distributed among the genomes of two endosymbionts, Ca. C. ruddii and a tentatively identified secondary symbiont (Sloan and Moran, 2012), thus, exhibiting the same pattern than the one described for the case of the aphid C. cedri (Martínez-Cano et al., 2015). It is worth to note that although aphids and psyllids belongs to the same suborder Sternorrhyncha, their corresponding endosymbionts belong to different bacteria lineages, indicating the independence of the symbiotic events (Baumann et al., 2006). This convergence suggests the existence of some evolutionary advantage of the complementation strategy. In the following we will focus on the case of C. cedri, because this system has been extensively characterized (Pérez-Brocal et al., 2006; Gosalbes et al., 2008; Burke and Moran, 2011; Lamelas et al., 2011).

The tryptophan biosynthetic pathway involves five enzymes, encoded by seven genes (trpA-F), which allows the conversion of chorismate to tryptophan (Crawford, 1989; see Figure 5). Moreover, as in the case of most amino acid biosynthesis, the first step (catalyzed by anthranilate synthase) is subject to allosteric feedback inhibition by tryptophan the end product of the pathway (Rouhbakhsh et al., 1996, 1997; Kwak et al., 1999; Spraggon et al., 2001). Previous studies on different B. aphidicola strains have showed that all the key residues in the allosteric site are conserved in the trpE sequences (Lai et al., 1994; Rouhbakhsh et al., 1996). However, since the trpE gene of B. aphidicola BCc was not considered in the mentioned study, we have performed a multi-sequence analysis to compare this sequence with its homologs in other eleven B. aphidicola strains, and in Escherichia coli and Serratia marcescens. Accordingly, we have found that the key residues P21 and S40, as well the motif LLES, purportedly involved in allosteric feedback inhibition (Caligiuri and Bauerle, 1991; Kwak et al., 1999; Tang et al., 2001), are even better conserved in B. aphidicola BCc than in the other B. aphidicola strains (See Supplementary Figure 6, and Materials and Methods Sections for more details), indicating that a product inhibition mechanism could be still active.

If there is an evolutionary pressure to increase the production of tryptophan, then the inhibiting effect of tryptophan accumulation can be overcome by producing large amounts of the anthranilate synthase enzyme (Lai et al., 1994). In fact, in several B. aphidicola strains, trpEG genes are amplified, either with multiple copy of the genes in the chromosome, or with multiple-copy plasmids, suggesting an adaptation to overproduce tryptophan for hosts (Rouhbakhsh et al., 1997; Baumann et al., 2006). However, increasing protein levels to overcome the inhibition effect, may represent a large cost in terms of protein synthesis. According to our model, the metabolic complementation strategy can be used overcome feedback inhibition, while reducing the total amount of protein required to supply enough product. Indeed, in the C. cedri consortium the genes coding for the tryptophan biosynthetic pathway are distributed among the two different endosymbionts: the first two genes, trpEG encoding the anthranilate synthase activity, are hosted in a plasmid in the P-symbiont B. aphidicola BCc, whereas the gene products necessary to catalyze the synthesis of tryptophan from anthranilate (trpDCBA) are encoded by the genome of the co-primary S. symbiotica SCc, which in turn, lacks the trpEG genes (Gosalbes et al., 2008; Lamelas et al., 2011). Thus, we argue that the presence of metabolic complementation in these bacterial endosymbionts could be interpreted as an evolutionary innovation to reduce the protein cost required to supply tryptophan to all the members of the consortium.

Furthermore, the model also predicts that when a metabolic complementation is an optimal strategy, the point where the pathway should be split will correspond to the most permeable metabolite. To contrast this prediction to the present case of study, we have collected structural data of all the intermediate metabolite of the tryptophan biosynthesis pathway, and we used different physicochemical properties and rule-based estimators to evaluate membrane permeability (See Material and Methods Section), including the molecular weights and the solubility of the compounds. Using these parameters, we identify plausible candidate metabolites for becoming breaking points in the tryptophan biosynthetic pathway. As seen in Table 1, among all the intermediate metabolites in this pathway, both anthranilate and indole appear as most plausible “breaking points,” while the other intermediates are less likely to cross membranes at an appreciable rate. Indole is an important signaling molecule in organisms such as E. coli (Martino et al., 2003), which is known to diffuse easily through cell membranes (Piñero-Fernandez et al., 2011). This raises the question: why the pathway breaking point is anthranilate instead of indole? An answer may lie in the structure of the tryptophan synthase complex (Dunn, 2012; Miles, 2013). The complex is a tetramer α2-β2. The α subunit, encoded by trpA gene, catalyses the synthesis of indole from indole-3-glycerol-phosphate. Then, indole is directly channeled to the β subunit, encoded by trpB, where it is used to synthesize tryptophan. Interestingly, it has been shown experimentally that if the α subunit is absent, the catalytic activity of tryptophan synthase is reduced by 1–2 orders of magnitude (Kirschner et al., 1991). Therefore, keeping both trpA and trpB in the same compartment appears to be a constraint and thus, indole is not a feasible candidate for the breaking point, since it is not allowed to diffuse freely in the cell. Anthranilate is thus the intermediate most likely to diffuse through cell membranes.

An additional example of a pathway breaking point concurrent with the most permeable metabolite is the biosynthesis of biotin in some aphids. All aphid members of the Lachninae family that have been analyzed up to now through genomic sequencing of their endosymbionts (namely, C. cedri, C. tujafilina, and T. salignus) harbor B. aphidicola and S. symbiotica in their bacteriocytes, and the biotin pathway is fragmented in two parts (Manzano-Marín et al., 2016). Thus, Buchnera cells seem responsible for the synthesis of the intermediate 8-amino-7-oxononanoate (KAPA) whereas Serratia genomes code for the enzymes required for biotin biosynthesis from KAPA. Since all the precursors of this intermediate are covalently bounded to an acyl carrier protein, KAPA is the first metabolite in the pathway showing an adequate permeability to cross membranes and paving the way for a metabolic complementation, in agreement with our model (see Supplementary Table 2).

Finally, we will briefly discuss the case of the transaminase activities losses in different lineages of B. aphidicola (Hansen and Moran, 2011; Poliakov et al., 2011; MacDonald et al., 2012), as well other endosymbionts, such as Tremblaya and Portiera (Zientz et al., 2004; Baumann et al., 2006; McCutcheon and von Dohlen, 2011; Zelezniak et al., 2015). Transaminase activities are a critical step in the biosynthesis of amino acids; transamination corresponds to the last step in the production of branched-chain amino acids (BCAs) and aromatic amino acids. Moreover, all BCAs share the same transaminase enzyme (EC 2.6.1.42), and the same happens for the aromatic amino acids (i.e., phenylalanine and tyrosine; with enzyme EC 2.6.1.57). In several lineages of endosymbionts, such as most strains of Buchnera, these transamination steps are lost, while the rest of the pathway is retained. Transaminases are expressed in the host cells (Hansen and Moran, 2011; Russell et al., 2013) originally with a catabolic function but now recruited to perform the last step in a biosynthetic pathway. As the in the case of tryptophan, the synthesis of these amino amino acids is commonly regulated at different levels, which usually include product feedback inhibition (White et al., 2012).

The above mentioned complementation does not involve two different endosymbionts, but is instead a metabolic relationship between the endosymbiont and its host. Still, this case includes all the main features of the model discussed in this work, i.e., a metabolite demand, feedback inhibition, and the absence of transporters. Thus, in order to test if the breaking point of these pathway corresponds to the most permeable metabolite, we have computed the permeability estimators for the intermediates in the biosynthesis of the branched amino acid as well the aromatic amino acid. The results show that in all the cases, the breaking point correspond to the most permeable metabolite (see Supplementary Table 2), thus suggesting that the complementation may have emerged due to an evolutionary pressure on the consortium.

The evolution of metabolic complementation in bacterial communities is most likely shaped by different forces, many of them may acting simultaneously (Tan et al., 2015). On one side, endosymbionts are subject to a process of irreversible gene loss and genome shrinkage, known as the Muller's ratchet, due to its confinement to intracellular life and the bottleneck caused by its small population size (Moran, 1996). Thus, under this scenario, the gene losses with no detrimental effect on the system could be fixed at the population level by the effect of genetic drift. In this case, metabolic complementation would not necessarily confer a selective advantage, and would be a consequence of random gene losses.

Instead, a selective explanation focuses on the cost of protein synthesis, which affects both the microbial cells and the host (which also has to provide nutrients to the endosymbionts). Therefore, hosts whose endosymbionts are able to sustain a given demand of certain product while reducing the protein cost will exhibit an improvement in their fitness (Shachrai et al., 2010). The same can apply to the fitness of microbial communities as a whole; interestingly, experimental results on synthetic cross-feeding bacteria has showed that the interactions that emerge in this scenarios confer a significant fitness advantage to the genotypes involved in the interaction, and thus, stabilize the coexistence of cross-feeding organisms (Mee et al., 2014; Pande et al., 2014; Germerodt et al., 2016).

Another evolutionary mechanism for the emergence of metabolic complementation is given by the Black Queen Hypothesis (BQH). According to the BQH, whenever a metabolic function is redundant among the member of a microbial community, and the final product of the process leaks or escape from the producer cells, becoming a public good, the rest of the members can take advantage of this “common good” (Morris et al., 2012). Thus, the BQH suggest that the loss of a leaky function will have some selective advantage at the level of each microbial species within the community; eventually, most of the members of the community will loss the genes coding for the function, until the production of public goods is just sufficient to support the community. This is not necessarily in contrast with the previous hypothesized fitness advantage for cross-feeders communities, since the loss of redundant functions may also increase the fitness of the system as a whole; the minimization of the total protein amounts might be a good approximation of the final state of a BQH-like dynamics of gene losses, where “intermediate goods” are available to all community members.

In perspective, understanding the evolutionary and mechanistic properties of natural microbial communities might allow the design of synthetic ones of potential industrial interest (Hosoda and Yomo, 2011; Großkopf and Soyer, 2014). For instance, the common laboratory bacteria E. coli have been used to study artificial communities (Hosoda et al., 2011; Mee et al., 2014; Pande et al., 2014). Interestingly, synthetic communities of cross-feeders based on E. coli mutant strains (with both over-producers and auxotrophs) were found to grow efficiently, even faster than the growth of the wild-type strain (Pande et al., 2014; Germerodt et al., 2016). Our model can be extended to these microbial communities by including the cost of protein synthesis in a genome-scale metabolic network, which has been recently been done in the context of constrained-based modeling for E. coli (Mori et al., 2016). However, the inclusion of allosteric regulation and diffusion still present some challenges (Tepper et al., 2013; MacHado et al., 2015). This approach may eventually lead to automated design of microbial communities, which would be of extreme interest for industrial applications.

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.