In order to reconstruct the genome-scale metabolic models (GEM) of S. symbiotica SCc and B. aphidicola BCc, we retrieved the corresponding genomes and the semi-automatically reconstructed Pathway-Genome Databases (PGDB) available in the BioCyc collection in Caspi et al. (2014). The used PGDBs' versions were SSYM568817 19.0 for S. symbiotica SCc, and BAPH37261 19.0 for B. aphidicola BCc (both available in Tier 3 at BioCyc 19.0) (see Supplementary Text S2). The public version of BAPH372461 does not include the information encoded in the pTpr-BCc plasmid (accession number EU660486.1, Gosalbes et al., 2008). This plasmid contains two genes (trpE and trpG) (Table 1), and was manually added to a local version of the B. aphidicola BCc PGDB, created through Pathway Tools v. 19.0 (Karp et al., 2015). Once added, the PathoLogic algorithm was run again in order to update the pathway prediction.

The reconstruction and refinement of the metabolic models was performed following the protocol described by Thiele and Palsson (2010). In order to reconstruct the metabolism of S. symbiotica SCc, the GEM of E. coli K12 MG1655 iJO1366 (Orth et al., 2011) was used as a reference, since this is the phylogenetically closest free-living organism for which a highly refined and validated model exists. First, orthologous genes were identified between the E. coli K12 MG1655 genome (available in EcoCyc 19.0) and the S. symbiotica SCc genome included in the PGDB SSYM568817 19.0. After the gene sequences were extracted and translated, they were compared using Blastp (Altschul et al., 1990) with an e-value maximum of 10⁻¹⁰ and an identity minimum of 75%. Subsequently, this proteome-level comparison was used to infer orthology using OrthoMCL (Li et al., 2003). Combining these results with the set of reactions and pathways present in the SSYM568817 PGDB, a gene-protein-reaction (GPR) table was constructed. Using this GPR table, together with the iJO1366 model and the BiGG database (Schellenberger et al., 2010) a first version of the S. symbiotica SCc model was reconstructed. Moreover, the biomass equation introduced in the model is a modified version of the one present in iJO1366, from which we removed membrane components and cofactors absent in the S. symbiotica SCc network. The coefficients were corrected using the methodology suggested by Henry et al. (2010). The first draft obtained was revised combining the unconnected module (UM) approach (Ponce-de-Leon et al., 2013) together with the previously published genome analyses (Lamelas et al., 2011a). Regarding the metabolism of B. aphidicola BCc, a previously published GEM, named BCc (Belda et al., 2012), was used as a reference. However, this GEM contained several blocked reactions and dead-end metabolites. The resolution of the different gaps was performed by solving the set of UMs combining the pathway inferences present in the BAPH372461 PGDB and the GEM of B. aphidicola Bap (MacDonald et al., 2011), available in BioModels (MODEL1012300000).

The aphid biomass equation was defined by including the various biomass components in corresponding stoichiometric proportions. Since the consortium model does not include the metabolism of the host, the aphid biomass equation only includes the compounds that are provided by the endosymbiotic consortium (i.e., the essential amino acids, vitamins, and cofactors) as previously done (Calle-Espinosa et al., 2016). The stoichiometric coefficients for the set of amino acids provided by the endosymbionts were estimated from the analysis of composition of the Aphis fabae and Acyrtosiphon pisum proteomes (Douglas et al., 2001; Russell et al., 2014). The compositions measured in two different aphid species showed good agreement, a fact that suggests that the values can be extrapolated to C. cedri (see Supplementary Text S1). We used the values from the experimental measurements of the A. pisum amino acid composition obtained by Russell et al. (2014) because this dataset includes the measurements of certain amino acids not included in Douglas et al. (2001). The values were normalized to represent composition in 1 g of dry weight (DW) in the aphid biomass. Additionally, the cofactors provided by the endosymbiotic consortium to the host were added. Since there is no estimation of their proportion in the aphid biomass, their values were set to be several orders of magnitude below the amino acids demand, but reflecting their essentiality to the host (Supplementary Table S5).

The different constraint-based methods used in the present work correspond to the current implementation found in COBRApy toolbox 0.6 (Ebrahim et al., 2013). The boundary conditions used in the individual analyses of the endosymbionts models, iSCc226 and iBCc98, can be found in Supplementary Tables S1, S2, respectively. However, in the case of the compartmentalized model of the consortium, a restriction that limits the total number of carbon atoms diffusing through the membrane of each endosymbiont was introduced (Burgard et al., 2001). By doing so, we are able to separately constrain the exchange of metabolites of each endosymbiont with its surrounding environment without the need of establishing arbitrary boundary for each exchange flux. Thus, instead of having a large set of parameters, one for each exchange flux, we only used two parameters, one for each endosymbiont. The restriction can be expressed as follows:

Where T_k ⊂ J is the set of indexes of the transporters present in the compartment k (endosymbiont) and J the set of all fluxes; v_i is the flux through the transporter i; nc_i is the number of carbon atoms in the transported molecule; and UB is the upper bound, i.e., the parameter establishing the maximum number of carbon atoms that can be exchanged by the α compartment (one of the endosymbionts) and the extracellular compartment (the host cell). This parameter limits the amount of matter flowing through the membrane in terms of the total carbon atoms and was fixed to 100 carbon atoms. In order to guarantee that the transport fluxes are non-negative, the reversible transporters were split into two irreversible transport reactions with opposite direction (Burgard et al., 2001). The biosynthetic capabilities of the endosymbionts for each biomass component, as well as for energy production, were assessed using Flux Balance Analysis (FBA) (Orth et al., 2010). Furthermore, network fragility was predicted through in-silico knockout experiments conducted using FBA as well as the Minimization of the Metabolic Adjustment (MOMA) (Segrè et al., 2002). Details on each method can be found in the extended materials and methods (see Supplementary Text S2).

The metabolic models of S. symbiotica SCc and B. aphidicola BCc were reconstructed individually, and were named iSCc236 and iBCc98, respectively (see Supplementary Text S1 for further details on the reconstruction, and Supplementary Tables S1, S2 and Supplementary File S1 for the complete models). The former consisted of 267 intracellular metabolites and 209 reactions catalyzed by the products of 236 genes, plus 11 orphan reactions (i.e., reactions with unknown coding genes). It also includes 30 transporters associated with a gene, and 49 orphan transport reactions. On the other hand, iBCc98 yielded a smaller network, containing 155 intracellular metabolites and 95 reactions catalyzed by the products of 98 genes, and 8 additional orphan reactions. Additionally, it includes only 5 transporters associated with a gene, and 58 orphan transporters. In both cases, all orphan reactions are required by the model in order to predict biomass formation. The genomes of these endosymbionts contain only a small set of genes coding for substrate-specific transport systems (Wernegreen, 2002; Charles et al., 2011). However, the corresponding metabolic models predict the necessity of metabolite transit through the endosymbiont membrane, albeit the transport mechanisms have not been elucidated in many cases, in which simple diffusion has been proposed as a plausible mechanism (Mori et al., 2016). The metabolic requirements and the biosynthetic capabilities for each of these two models were congruent with those inferred from genomic analyses (Lamelas et al., 2011a), with the exception of the biosynthesis of asparagine, which is predicted by the iSCc236 model to be required as nutritional input provided by the host instead of being synthesized by S. symbiotica SCc (Supplementary Figure S1). Finally, when the energetic capabilities, i.e., the synthesis of ATP, were analyzed for both symbionts, we found that iBCc98 predicts a very low yield of ATP, a limitation that is a direct consequence of the absence of ATP synthase. In turn, this lack of a proton pump mechanism poses a constraint on the regeneration of NADH, through the NADH dehydrogenase complex. The model suggests that part of the NADH may be driven through the conversion of the pair homocysteine and serine, into glycine and methionine (for further details, see Supplementary Text S1).

Since we were interested in modeling the whole consortium, the two metabolic models previously introduced (iBCc98 and iSCc236), were combined to create a single model named iBSCc. Although it is known that B. aphidicola BCc population and S. symbiotica SCc are hosted in different bacteriocytes, a more simplistic representation was chosen where both endosymbiont models are embedded in a single compartment, in a similar way as it has been previously done (Ankrah et al., 2017). Therefore, our model included: (i) a compartment representing the B. aphidicola BCc population; (ii) a compartment representing the S. symbiotica SCc population; and (iii) a single extracellular compartment representing the host cells, where both symbionts are embedded (Figure 1). Consequently, the boundary of the system was defined by the compartment representing the host, and exchange fluxes across the boundary represented the metabolites supplied and consumed by the aphid, as well as the excretion products. The 52 metabolites found in the extracellular compartment include those obtained from the host, common excretion products and metabolites exchanged by the two symbionts. Additionally, five reactions were added to the extracellular compartment, as they have been suggested to be performed by the host and to play a relevant role in the metabolic complementation between host and symbionts (Hansen and Moran, 2011; Poliakov et al., 2011). These reactions include the conversion of phenylalanine into tyrosine (1 reaction), the production of homocysteine and adenosyl-methionine from cysteine (3 reactions), and the assimilation of hydrogen sulfide produced by S. symbiotica SCc, for the production of lipoate (1 reaction) (see Supplementary Table S3).

In order to simulate the growth of the system as a whole, we introduced a combined meta-biomass equation, where each member contributes to the growth of the system with a fixed stoichiometry (Figure 1A). Our model assumes the coupling between the growths of all the members of the system, which would especially apply during the development of the host (Figure 1B). This assumption is justified by different sources of experimental data as well certain theoretical results (see Supplementary Text S1 for further details). Since the exact contribution of the symbionts to the cedar aphid biomass is unknown, we used data obtained from Schizaphis graminum and A. pisum indicating that their symbionts represent 5–15% of the system's biomass (Whitehead and Douglas, 1993; Baumann et al., 2006). Thus, we modeled the proportion between biomass of the symbionts and the host to be 1:9. Finally, since imaging data from the cedar aphid bacteriome indicates similar proportions between the two symbiont species (Gomez-Valero et al., 2004; Pérez-Brocal et al., 2006), we modeled that each bacterial member represents 5% of the system's biomass at every time. It is important to emphasize that this assumption is key for the direct utilization of FBA and other related techniques in the study of the iBSCc model. The reason is that, if both symbionts are in proportion 1:1, the fluxes are normalized by the same quantity (e.g., dry weight of bacteria) and they grow necessarily at the same growth rate. This means that they behave essentially as one “big” bacteria. However, if this assumption does not hold, and the growth rate of the two bacteria are the same (i.e., they are coupled), new (probably non-linear) constraints should be incorporated to the model (see Kerner et al., 2012 as an example). Once the meta-biomass equation was introduced, we use FBA to verify that the iBSCc model is consistent, allowing for a positive meta-biomass flux and the growth of each member of the consortium (see Supplementary Table S3). This optimal flux distribution predicted by FBA (maximization of the meta-biomass, i.e., the growth rate of the iBSCc consortium) is analyzed in the next section.

The overall structure of the models corroborates that the metabolic network of B. aphidicola BCc is specialized in the production of amino acids, while S. symbiotica SCc produces nucleotides and a large number of cofactors (Supplementary Figure S1). Moreover, the model predicts that S. symbiotica SCc can synthesize tryptophan from anthranilate, which has been shown to be provided by B. aphidicola BCc in a paradigmatic case of metabolic complementation (Gosalbes et al., 2008; Lamelas et al., 2011a; Martínez-Cano et al., 2015; Manzano-Marín et al., 2016). Another case of metabolic complementation between the two symbionts predicted by the model occurs in the biotin synthesis pathway, which takes place via the import of the precursor 8-amino-7-oxononanoate (8AONN), produced by B. aphidicola BCc, as recently suggested from genomic data (Manzano-Marín et al., 2016). Lysine biosynthesis also represents a case of complementation, based on the fact that the genome of S. symbiotica SCc encodes all activities of the pathway except the last one, which converts meso-diaminopimelate into lysine (Lamelas et al., 2011a). However, our model indicates that meso-diaminopimelate is used for the synthesis of peptidoglycan, while the complete lysine synthesis pathway is conserved in B. aphidicola BCc, suggesting that this complementation does not occur.

The FBA predictions showed that B. aphidicola BCc synthesizes and provides the host and S. symbiotica SCc with ten amino acids. B. aphidicola BCc also synthesizes anthranilate, 8AONN and shikimate and releases them to the host compartment. In turn, S. symbiotica SCc imports anthranilate, 8AONN and shikimate, and use them as precursors for the biosynthesis of tryptophan, biotin and tetrahydrofolate (THF), respectively (see Figure 2 and Supplementary Table S3). Therefore, three events of metabolic complementation between the two symbionts are predicted, including the biosynthesis of tryptophan, biotin and THF; whereas the former two have been inferred previously from genomic analyses (Gosalbes et al., 2008; Manzano-Marín et al., 2016), the latter is a redefined complementation case (Lamelas et al., 2011a, see Figure 2). Folate cross-feeding has been described between Serratia grimesii and Treponema primitia, both members of the termite Zootermopsis angusticollis gut microbiome (Graber and Breznak, 2005). Since C. cedri, as animals in general, is not able to synthesize folate, and this is not included in the aphid's diet, Serratia probably provides this essential cofactor to both the host and Buchnera.

Regarding S. symbiotica SCc, the iBSCc model predicts that it can synthesize, without need for complementation, cysteine, the four deoxynucleotides, the four triphosphate nucleosides, and thirteen cofactors and coenzymes. Several of these compounds may also be synthesized by the host (e.g., nucleotides, deoxynucleotides and NAD+), who might provide them to B. aphidicola BCc. Indeed, the aphid A. pisum, as most eukaryotes, is able to synthesize nitrogenous bases (Richards et al., 2010; Vellozo et al., 2011), which indicates that C. cedri also should. In this study, however, we will assume that it is S. symbiotica SCc who provides Buchnera with nucleotides and cofactors, and the host with cofactors such as biotin, riboflavin and THF.

The construction of metabolic models of two different endosymbionts with a notable difference in size allowed us to study whether the reductive evolutionary trends indeed generate smaller, more fragile networks. The robustness of the metabolic networks of the cedar aphid endosymbionts was assessed through in-silico knockout analyses conducted with two alternative approaches: FBA and MOMA. In the first place, the fragility of the whole consortium was considered by using the meta-biomass flux as an indicator of the viability. The results show that around 85% (~88% using MOMA) of the genes coded by the endosymbionts are essential in order to sustain the growth of the whole system (see Supplementary Table S4). Then we focussed on the fragility of the individual endosymbionts networks. For the case of iBCc98 FBA predicted that 72 out of the 98 metabolic genes (~74%) are essential, while the MOMA analyses identify 76 (~78%) as essential. The dispensable genes are mostly involved in catabolism, affecting the phosphate pentose pathway, glycolysis, respiratory chain and pyruvate fermentation (see Supplementary Table S4). When performing FBA robustness analyses on iSCc236, 209 genes (~88%) are predicted to be essential, whereas MOMA predicted 5 additional genes as essentials. The 28 dispensable genes predicted by both methods code for 36 enzymatic reactions involved mostly in biosynthetic pathways (e.g., nucleotides and cofactors), but also in the central carbon metabolism (e.g., the pentose phosphate pathway and glycolysis). If we consider the cell wall genes to be dispensable (since they have been repeatedly lost in endosymbiotic bacteria), the number of essential genes drops to 189 (~80%). Additionally, if it is also assumed that the host is who provides the nucleotides and deoxynucleotides, the percentage of essential genes drops to ~70% (data not shown). Although it might seem surprising that these estimates are higher or comparable to those obtained by using the smaller iBCc98 network, the iSCc236 model requires 22 organic compounds to be imported, while iBCc98 requires 29, among which there are nucleotides and cofactors such as NAD⁺ and coenzyme A. Finally, a study of distributed robustness was performed through the analysis of the synthetic lethal (or double lethal) genes, i.e., pairs of non-essential genes whose simultaneous inactivation yields lethality (Wagner, 2005). In the case of iBCc98, ca. 15% of the pairwise combinations between the 26 non-essential genes predicted with FBA are predicted as synthetic lethal. These combinations include 20 of the non-essential genes, indicating that most of the dispensable genes have only a shallow degree of redundancy. On the other hand, iSCc236 predicts only ca. 5% of the possible combinations between the 34 non-essential genes to be lethal. Moreover, by disregarding the cell wall biosynthesis genes as explained above, this number drops to ca. 2%.

Although it is generally accepted that these patterns of complementation are the outcome of the process of genome reduction, whether such organization of the metabolic networks confer a selective advantage to the whole system or not, remains an open question. Previous work has focused on the kinetic aspect of the problem, in particular the role of product inhibition as a plausible condition that may drive the emergence of metabolic complementation (Mori et al., 2016). Herein, we try to approach this problem from a structural point of view, by comparing the metabolic capabilities in alternative scenarios of gene losses and retentions within the three-shared pathways. For this purpose, iBSCc was extended to represent a putative ancestral-state model of the consortium, named iBSCc_Ancest, where the three-shared pathways are still complete in both endosymbionts (see Supplementary Text S2). In this way, it is possible to compare how the different scenarios of gene loss and retention perform with respect to the putative ancestor as well as to the pattern of complementation exhibited by the cedar aphid consortium.

Using iBSCc_Ancest, the space of all viable and scenarios of gene loss and retention patterns (GLRPs) were generated by removing from this model every possible combination of genes, from single genes to the most reduced cases where only one copy of each gene remains present (at least one of the two copies for each gene needs to be functional). A GLRP is considered viable if, after removing the corresponding reactions, FBA predicts a meta-biomass flux greater than zero. Furthermore, in order to reduce the number of combinations of possible GLRPs, the enzyme subsets (ES), i.e., groups of enzymes that always work together under steady state, were first computed for iBSCc_Ancest (see the extended Materials and Methods section in the Supplementary Text S2). Since removing an enzyme from an ES is equivalent to removing all the enzymes in the ES, each ES can be treated as a functional unit. Then, those genes coding for enzymes in the same ES were grouped together. Table 1 shows the structure of the ES for the three pathways considered in this study. Notably, the 32 enzyme activities are grouped in 7 ESs. Given that in iBSCc_Ancest each endosymbiont includes the 7 ESs, the enumeration process yielded a total of 2188 viable GLRPs, which represent alternative consortium models. From this set, 128 are minimal GLRPs, i.e., only one copy of each gene remains active. It is worth noting that the pattern exhibited by the C. cedri consortium is not minimal since the genes coding for the three activities which allow the conversion of shikimate into chorismate are present (see Table 1). In order to simplify the notation each GLRP was coded by a sequence of symbols Δnm, where the subscript n (from 1 to 7) denotes the loss of the enzyme subset by one of the endosymbionts, indicated by the superscript m, which can be either B or S (B. aphidicola or S. symbiotica). Accordingly, Δ1B represents a simple GLRP where B. aphidicola has lost the ES1. A more complex example would be the case of iBSCc (which is the actual pattern exhibited by B. aphidicola and S. symbiotica in the cedar aphid). In this case, the GLRP is coded as follows: Δ1SΔ3SΔ4BΔ5BΔ6SΔ7B. Note that in this case the ES2 is omitted because the genes involved in this ES are present in both endosymbionts, as mentioned above (see Table 1).

The 2188 consortium models, which represent alternative evolutionary scenarios, were evaluated using FBA in two ways by considering: (i) the individual maximal production rate of tryptophan (vTrp*), THF (vTHF*) and biotin (vBtn*); and (ii) the whole system performance, calculated by optimizing the meta-biomass production rate. In all the cases, the optimal values are normalized with respect to the optimal value showed by iBSCc_Anc. Figure 3 summarizes the results of the reduction experiment in terms of the production capabilities of tryptophan, THF and biotin for each GLRP. Firstly, the results show that for any of the three objectives, the production rates of the different GLRPs exhibit great variability (for more details see Supplementary Table S5). This clearly shows that the way in which a pathway is distributed in a complementation event has a profound impact in the pathway capabilities. In the case of the tryptophan production, the different GLRPs can be divided into two main groups: (i) a group of GLRPs (which includes iBSCc) with a tryptophan production rate almost equal to the one exhibited by the putative ancestor (i.e., vTrp*~100%); and ii) a larger group of GLRPs with production rate vTrp*< 90% (Figure 3a). Furthermore, when considering only the minimal GLRPs (i.e., cases with 50% of genes lost) the gap is even larger, and there is only one GLRP with vTrp* ~ 100% whereas the other patterns exhibit a vTrp* < 76%. This minimal GLRP with vTrp* ~100% corresponds to the pattern exhibited by iBSCc with the additional loss of the only set of redundant genes that remain present, which form the ES2 (see Table 1).

On the other hand, when considering the biosynthesis of THF and biotin, the results also show a wide dispersion for the normalized production rate values exhibited by the different GLRPs. However, unlike the case of the tryptophan, in these two cases the optimal production value exhibited by iBSCc decreases considerably with respect to the value exhibited by the putative ancestor (Figures 3b,c). For the case of the production of THF, all the GLRPs with the highest normalized production rate (vTHF*~ 93%) and with more than 40% of the genes lost, shared the following two losses: Δ1S and Δ5S, which correspond to the scenario in where S. symbiotica losses the biosynthetic pathways of shikimate and THF. Whereas Δ1S is consistent with the GLRP shown by iBSCc, Δ5S is the opposite, i.e., in the cedar aphid consortium B. aphidicola has lost the THF biosynthetic pathway. Furthermore, all those GLRPs that involve Δ5B imply an important drop in the normalized production rate vTHF* < 72% (Figure 3b). In the case of the biotin biosynthesis where the GLRP of iBSCc implies vBtn* ~ 64% something similar was found (Figure 3c). According to these results, when the system has lost more than 40% of the genes involved in the analyzed pathways, those GLRPs with a biotin production rate closest to the one exhibited by the putative ancestor, share the same pattern (Δ6BΔ7S). However, this pattern is opposite to the one exhibited by iBSCc, i.e., B. aphidicola losses the capability to synthesize 8AONN and retains the capability to produce biotin from this precursor, whereas S. symbiotica exhibits the complementary pattern. Similar results were found when considering the reduction of each pathway individually, i.e., when only a single pathway is considered redundant, the GLRPs include only the genes involved in the pathway (see Supplementary Figures S2–S4 and Supplementary Table S6).

After the analysis of how the different GLRPs perform over the production rates of individual biomass components, the same study was conducted but using the meta-biomass production rate as a proxy to study the fitness of the alternative evolutionary scenarios. The simulation results, summarized in Figure 4, show that the value of the normalized meta-biomass production (vmeta-biomassn*) varies between 95 and 100%, compared to the putative ancestor (Figure 4a). Although this range is quite narrow (Figure 4b), it may still play a role in a selective process, since herein only the stoichiometric rate is considered, and no other factor, such as the cost of protein synthesis, is considered. Moreover, in a first look the results indicate that the iBSCc GLRP vmeta-biomass* is quite close to the one of the putative ancestor.

A deeper insight into the results shows that among the minimal GLRPs, 48 out of the total 128 (~37%) also exhibit a vmeta-biomass*> 99%. Furthermore, within this set it is possible to find opposite patterns, for example cases such as Δ3SΔ4S and Δ3BΔ4B imply almost identical production rates. Nevertheless, there are some particular patterns that are consistent with the results obtained analyzing the individual pathways. For example, all the minimal GLRPs which include the loss Δ3BΔ4S (i.e., the opposite pattern than the exhibited by iBSCc) imply vmeta-biomassn*< 99%. On the other hand, there are many GLRPs that could represent disadvantages in terms of the system's growth rate. Such would be outcompeted by other more efficient organization of the complementation. Thus, considering only the structure of the networks of the endosymbionts, the simulations indicate that although many GLRPs are possible, the pattern observed in the cedar aphid consortium has retained a large growth rate while undergoing a large degree of gene loss.

The reconstruction and metabolic analysis of B. aphidicola BCc and S. symbiotica SCc, co-primary endosymbionts of the cedar aphid C. cedri, has allowed, in first place, the revision of the annotation of these organisms' genomes. Indeed, through the manual curation of these networks using the UM approach, it was possible to correct annotation errors in different enzymatic activities, and to identify previously unannotated metabolic genes. Moreover, the simulations performed here allowed the refinement of growth conditions and metabolic capabilities of these endosymbionts as compared to previous inferences from genomic analyses. Our study also indicates that GEMs are a possible way to approach the otherwise inaccessible phenotypic features of uncultivable endosymbionts. On the other hand, the analysis of highly reduced metabolic networks, such as the case of iBCc98, bring into consideration methodological issues such as the case of the proton balance. Previously, a sensitivity analysis by Reed et al. (2003) in the genome-scale metabolic model of E. coli K12 iJR904 revealed that the net proton balance could be positive or negative depending on the carbon source used in the growth simulations, which would acidify or basify the environment, respectively. Alas, these predictions remain empirically untested.

The effect caused by proton balance in GEMs is generally low due to the size of these networks. However, in endosymbionts and other small networks, it may generate notorious consequences. In iBCc98, for instance, it considerably reduces the versatility of the metabolism of B. aphidicola BCc by coupling presumably independent processes, such as ATP synthesis and the folate cycle, and over-producing amino acids. Although this metabolic organization would be clearly disadvantageous for a free-living organism, for an endosymbiont member of a nutritional symbiosis it may be selected for at the level of the host. Indeed, a similar behavior has been described recently as applied to the whitefly endosymbiont P. aleyrodidarum (Calle-Espinosa et al., 2016), where the growth of the organism is coupled with the overproduction of amino acids and carotenoids as a consequence of its low energetic capabilities. This phenomenon might play a relevant role in the evolution of nutritional endosymbiosis but it may also represent a methodological artifact as the consequence of a lack of knowledge on how to formulate in such a model the transport of protons through the membrane. One possibility would be the use of protons for the transport of compounds against their gradients. Although the scarcity of annotated transporters in the B. aphidicola BCc genome (Charles et al., 2011) does not seem to support this scenario, this problem falls within a more general umbrella, whereby the nature of the cell envelope (including both the membrane composition and the transport systems) of endosymbionts is largely unknown and might rely on contributions from the host (McCutcheon, 2016).

On the other hand, simulations with iBCc98 and iSCc236 indicate resemblance with previous metabolic analyses from endosymbionts and other bacteria with reduced genomes. We found that these two networks contain very few dispensable genes, with essentiality estimates being around 88% of genes for iSCc236, and 73% for iBCc98 (91 and 78%, respectively, according to MOMA). Although, these results may contradict the idea that the smaller the network, the higher the essentiality of its components, the percentage of essential genes predicted by iSCc236 drastically drops to ~70%, when the cell wall genes are considered dispensable, and that the nucleotides and deoxynucleotides are provided by the host. Moreover, a further 15 and 5% of the genes, respectively, are genes that become essential after the deletion of another non-essential gene, thus displaying only a shallow degree of redundancy. Altogether, the amount of non-essential genes seems to positively correlate with the size of the network (Gil and Peretó, 2015). Moreover, a previous metabolic analysis of B. aphidicola APS showed that the network of this endosymbiont is also highly non-redundant, with 84% (94% according to MOMA) of genes being essential (Thomas et al., 2009). Although differences in the estimation algorithm and the selected criteria make these numbers not directly comparable, it does come to show that the Buchnera lineage evolved nearly minimal networks before the divergence of the Aphidinae and Lachninae aphid subfamilies, about 90 Mya. Moreover, although there are no available metabolic reconstructions for others S. symbiotica strains, the essentiality in the metabolic network of S. symbiotica SCc probably represents the high degree of genome reduction occurring in its symbiotic lineage, at 1.76 Mb and only 672 CDS (Lamelas et al., 2011b). This is likely the result of more recent evolutionary processes, since S. symbiotica from hosts within the Aphidinae subfamily display genomes larger than 2.5 Mb and contain over 2000 CDSs (Burke and Moran, 2011; Foray et al., 2014; Manzano-Marín and Latorre, 2016). This might indicate that the genome reduction process in S. symbiotica SCc occurred after the divergence of the Lachninae subfamily, ca. 55 Mya. Two other genomes from obligate S. symbiotica strains have been described within this subfamily, obtained from Tuberolachnus salignus (tribe Tuberolachnini) and Cinara tujafilina (Eulachnini). The former contains a genome of only 650 kb and 495 CDSs (Manzano-Marín et al., 2016), while the genome of the latter is in an early stage of reduction at 2.5 Mb and 1602 CDSs (Manzano-Marín and Latorre, 2014). If the transmission of Serratia was vertical and no replacement occurred in the C. tujafilina clade, as it has been suggested (Manzano-Marín and Latorre, 2016), extreme reductive processes observed in the iSCc236 network may have been even more recent, possibly occurring no longer than 40 Mya.

Altogether, the two metabolic networks involved in the C. cedri consortium are highly constrained and fragile. This is also shown in the list of metabolic requirements that these organisms exhibit, which is increased by the high number of full and partial pathways that have been lost in both members of the consortium. Moreover, they show a high degree of integration, where both members have suffered massive losses, presumably due to division of labor with the bacterial partner. Cases of such losses are, for instance, the B. aphidicola BCc loss of the ability to produce cofactors like siroheme, biotin or THF, and, in S. symbiotica SCc, of the ability to produce several amino acids such as phenylalanine, threonine and branched amino acids. The model iBSCc establishes three cases of pathway sharing between the two symbionts, namely the biosynthesis of tryptophan, biotin and THF. These three pathways are partitioned between the two bacteria at the level of specific metabolites: anthranilate, 8AONN and shikimate, respectively (see Figure 2). Those three exchanged metabolites are among the most permeable ones of the participant intermediates, as it was predicted by our previous chemoinformatic analysis on metabolic complementation (Mori et al., 2016). Finally, the endosymbiont metabolic networks also show the need for metabolic complementation from compounds synthesized by the host. This is reflected by the requirement for the incorporation of metabolic intermediaries, not just end-products, from the external compartment. Simulations with the metabolic network of S. symbiotica SCc show that two biosynthetic pathways are completed by the host. The loss of the first activities in the biosynthesis of siroheme and coenzyme A generate the need for the import of metabolic intermediaries. A similar case has been previously observed in Blattabacterium, symbiont of the cockroach Blattella germanica, where the initial steps in the biosynthesis of terpenes have been lost (Ponce-de-Leon et al., 2013). These events are likely to commonly evolve in organisms under genome reduction processes, enabled by key factors such as the redundancy of pathways and the feasibility for transport due to permeability of the compound or to the existence or exaptation of generalist transporters.

The simulations performed in this study with the reconstructed model of the consortium, iBSCc, reproduced the cases of metabolic complementation between the two bacterial partners that have been described in the literature, and predicted an additional one, the one involving the synthesis of THF in S. symbiotica SCc from the shikimate provided by B. aphidicola BCc. The emergence of metabolic complementation is a complex phenomenon and is not well understood from a theoretical perspective. Tryptophan biosynthesis is highly regulated at the transcriptional, translational and posttranslational levels. Notably, this pathway is negatively regulated through attenuation of the transcription of anthranilate synthase by tryptophan (Crawford, 1989), a process that has been suggested to facilitate the emergence of metabolic complementation (Mori et al., 2016). Genes for tryptophan synthesis from chorismate are present in most of the sequenced B. aphidicola genomes. In the cedar aphid, the genes encoding this pathway are divided in a B. aphidicola BCc plasmid containing trpEG, and a S. symbiotica SCc operon containing trpABCD. A strikingly identical case has been found to have occurred convergently in separate lineages of Buchnera and Serratia in the aphid T. salignus (Manzano-Marín et al., 2016). Moreover, the same complementation has been described in a completely different system, the one formed by ca. Carsonella ruddii and secondary symbionts related to Sodalis and Moranella, in psyllids (Sloan and Moran, 2012). More complex cases of complementation in tryptophan biosynthesis also exist, such as the case of the mealybug symbionts, where this pathway seems to require the transport of multiple intermediaries, possibly including anthranilate (McCutcheon and von Dohlen, 2011).

Under this perspective, our simulations of all alternative GLRPs of metabolic complementations applied to the three shared pathways yielded novel results. We constructed models representing alternative GLRPs from a hypothetical ancestor containing intact pathways in both symbionts, and compared how well they performed in maximizing the production of tryptophan, THF, biotin and meta-biomass, assuming that BCc and SCc are in proportion 1:1 in all the explored scenarios. Surprisingly, we found that, in the case of the tryptophan, the GLRP exhibited in the cedar aphid consortium behaved nearly optimal, and represented a quasi-minimal design in doing so. These results show that, from a structural point of view, the actual distribution or division of a metabolic pathway between two organisms can perform almost as well as their ancestor while using a smaller gene repertoire. Although the cost of protein synthesis and genome replication cannot be easily integrated in this kind of analysis, it is expected that a reduction in such costs will improve the growth efficiency of an organism (Mori et al., 2016). On the other hand, the results indicate that the GLRP exhibited in the biotin and THF pathways is suboptimal, and there are others GLRPs that allow a greater flux using the same number of genes. However, it is worth noting that the demand for these cofactors is probably much smaller than the demand for tryptophan, and thus the selective pressure for an efficient production of THF and biotin may be less stringent than the case of the amino acid. This situation may also reflect the limits of the simultaneous optimization of the diverse metabolic performances of a given network. In any case, the GLRP of the tryptophan biosynthesis exhibited by the consortium of the cedar aphid is convergent with the metabolic solution observed in the symbiotic consortium of the psyllid Heteropsylla cubana (Sloan and Moran, 2012).

On the other hand, the reduction experiment performed evaluating the meta-biomass equation indicate that the global GLRP exhibited by the cedar aphid consortium, is quasi-optimal in terms of yield, and nearly minimal in terms of gene number. Moreover, additional factors can influence the given structure, and confer more benefits. For instance, the fact that the first and the last steps of the pathway are performed in different compartments reduces the possibility that an accumulation of the end-product tryptophan would inhibit the first biosynthetic step by attenuation. This is despite the fact that the inhibition binding site for tryptophan is highly conserved (Mori et al., 2016), which might be due to constraints in the enzyme functional architecture. Besides the structure of the complementation, the pathway kinetics is also likely to be involved in the function and evolution of the metabolic complementation. However, the complexity associated to a kinetic model and the lack of experimentally-based parameters makes such a model implausible at a genome scale as of today. Future efforts to model complete metabolic systems at a genome scale beyond stoichiometric constraints, by adding reaction kinetics and higher-level processes, such as the cost of protein production and turnover, will shed more light into the structure, function and evolution of metabolism, and in the emergence of metabolic complementation.

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.