extfieldsofstudy:
Computer Science
Medicine

text:
Understanding the evolution of functional redundancy in metabolic networks

Abstract Motivation Metabolic networks have evolved to reduce the disruption of key metabolic pathways by the establishment of redundant genes/reactions. Synthetic lethals in metabolic networks provide a window to study these functional redundancies. While synthetic lethals have been previously studied in different organisms, there has been no study on how the synthetic lethals are shaped during adaptation/evolution. Results To understand the adaptive functional redundancies that exist in metabolic networks, we here explore a vast space of ‘random’ metabolic networks evolved on a glucose environment. We examine essential and synthetic lethal reactions in these random metabolic networks, evaluating over 39 billion phenotypes using an efficient algorithm previously developed in our lab, Fast-SL. We establish that nature tends to harbour higher levels of functional redundancies compared with random networks. We then examined the propensity for different reactions to compensate for one another and show that certain key metabolic reactions that are necessary for growth in a particular growth medium show much higher redundancies, and can partner with hundreds of different reactions across the metabolic networks that we studied. We also observe that certain redundancies are unique to environments while some others are observed in all environments. Interestingly, we observe that even very diverse reactions, such as those belonging to distant pathways, show synthetic lethality, illustrating the distributed nature of robustness in metabolism. Our study paves the way for understanding the evolution of redundancy in metabolic networks, and sheds light on the varied compensation mechanisms that serve to enhance robustness. Supplementary information Supplementary data are available at Bioinformatics online.

Introduction
Biological systems are highly robust against perturbations, which includes loss of gene function (Wagner, 2005). Robustness against null mutations in living systems is carried out by many different mechanisms, one of them being redundancy (Hartman et al., 2001). Redundant genes or functionally redundant pathways tend to cope up for any perturbations in the system. Although there exist a few genes that are indispensable to the survival of the organism, duplicate genes are often retained during evolution in order to combat mutations (Papp et al., 2004). Some genes might be inactive in certain environmental conditions but become active in a different condition to compensate. Even though duplicate genes exist, some alternate metabolic pathways also play a role in establishing robustness (Mahadevan and Lovley, 2008). These perform distinct functional roles and are not redundant in all environmental conditions (Harrison et al., 2007;Ihmels et al., 2004). Adaptation also facilitates the addition of new functional associations, which do not arise to overcome mutational effects per se but are evolutionary byproducts and later enable robustness against deleterious mutations (Pál et al., 2005;Wang and Zhang, 2010).
Metabolism forms the core of all cellular processes, building up components required for the growth and survival of organisms. Metabolic networks play a key role in synthesizing components necessary for growth from the nutrients available to the organism. Perturbations to key enzymatic pathways often abolish growth, thus causing lethality. A comprehensive understanding of the metabolic reactions in an organism establishes that some fraction of the reaction functions/pathways is redundant, to enable robustness against single mutations (Wang and Zhang, 2009). Metabolic networks are highly robust against perturbations and utilize the backup  This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com Bioinformatics, 34, 2018, i981-i987 doi: 10.1093/bioinformatics/bty604 ECCB 2018 reactions/pathways that exist to compensate any lethal damage to the organism (Behre et al., 2008). Synthetic lethality refers to functional associations that exist between genes/reactions that compensate for one another upon single gene deletions (Hartman et al., 2001). Synthetic lethals refer to a pair of reactions that perform distinct functionalities although they enable organisms to combat lethal perturbations in the system. A synthetic double lethal (DL), therefore, consists of a pair of reactions, which when deleted individually enable the organism to survive, but when deleted together results in lethality. Synthetic DLs enable an organism to survive single mutations/knock-outs, thus enhancing the robustness of the metabolic networks. In a majority of the cases, one reaction involves in a key metabolic pathway and upon deletion of the reaction, fluxes are re-routed through the lethal pair to enable the growth of the organism (Ghim et al., 2005;Gü ell et al., 2014). Synthetic lethals play an important role in providing backup functions, thus, helps the organism sustain perturbations. Experimental identification of these synthetic lethals through highthroughput techniques is still difficult (Kaelin Jr, 2005). Therefore, many computational methods, based on flux balance analysis (FBA) (Varma and Palsson, 1994;Kauffman et al., 2003), have been developed in the past, viz. SL Finder (Suthers et al., 2009), MCSEnumerator (von Kamp andKlamt, 2014) and Fast-SL (Pratapa et al., 2015). Of these, Fast-SL, developed previously in our laboratory, provides a very efficient way to predict synthetic lethals in a large number of networks, under different growth conditions.
With the above perspective, we aimed to address some of the key questions: How do metabolic networks evolve to acquire new functional redundancies? How do reactions compensate for one another under perturbation? Do the redundancies vary based on the environmental conditions? Does evolution retain native functional redundancies? We build on a previously reported methodology to evolve metabolic networks in silico (Barve and Wagner, 2013) to answer these questions. We generated 2000 metabolic networks in minimal glucose environment and identified single lethals and DLs in these random metabolic networks. We show that across networks, many reactions compensate for one another under various environmental conditions. We observed that key functional reactions evolve to acquire new redundancies and that they tend to partner with a higher number of reactions. We also noticed that real networks possess a significantly higher number of essential reactions and redundant reactions when compared with random networks. We further found that reactions which play very different functional roles can form redundancies in different environmental conditions. Our results also suggest that some of the functional redundancies are environment specific and some compensate for one another in all environmental conditions. Overall, our approach builds towards a much better understanding of the evolution of functional redundancies in metabolism.

Universe of possible reactions
The 'universe' of possible reactions is the set of all the biochemical reactions that are known to occur in organisms. The reaction universe was adapted from a previously described, pruned dataset from Szappanos et al. (2016), which is based on the ModelSEED database (Henry et al., 2010). This universal set of reactions consists of 5162 reactions, including exchange and transport reactions, comprising 2805 metabolites.

Flux balance analysis
The cornerstone of our analyses is the ability to predict the growth phenotype of any given metabolic network. For this, we employ FBA, which is a widely used constraint-based approach to identify an optimized flux distribution, for a given metabolic network, based on maximizing/minimizing a given (linear) objective function (Varma and Palsson, 1994;Kauffman et al., 2003). FBA solves a linear programming problem to identify a flux distribution, assuming that the network is in steady state. The formulation of FBA is as below: where c denotes the objective function, v is the vector of fluxes of all the biochemical reactions and S corresponds to the stoichiometric matrix, of dimensions m Â r that encodes the stoichiometry of all m metabolites in r reactions present in the network. FBA is a reliable method to predict the growth rate of organisms in many conditions (Blank et al., 2005;Papp et al., 2004;Harrison et al., 2007;Szappanos et al., 2016) and has been previously used for studying metabolic network robustness and plasticity (Barve et al., 2012;Almaas et al., 2005;Barve and Wagner, 2013;Gü ell et al., 2014).

Generation of random networks
Random networks were generated using a modification of the previously devised method developed by Barve and Wagner (2013). The method uses a Markov-chain Monte Carlo (MCMC) approach to swap a random reaction from a known genome-scale metabolic network with a randomly selected reaction from the universe. At each random swap, an FBA was performed to evaluate the growth of the generated network on a given carbon source. The threshold for identifying growth on a particular carbon source was set to 5% of the wild-type flux, i.e. the biomass flux of E.coli in that particular carbon source. We generated 2000 such random networks starting from the E.coli iAF1260 metabolic network (1983 reactions þ 303 exchange reactions), with each random network produced as a result of 5000 such reaction swaps. For each of the random networks, a genotype vector was obtained, which denotes the presence/absence of reactions from the universe. Each of our random networks includes all 553 exchange reactions from the universe, resulting in metabolic networks each having 2536 reactions.

Growth in different environments
We analysed synthetic lethals and the growth capabilities of the metabolic networks in a set of ten environments, a subset of the fifty environments reported earlier (Barve and Wagner, 2013). These are minimal growth environments that contain oxygen, ammonium, inorganic phosphate, sulphate, sodium, potassium, cobalt, iron, protons, water, molybdate, copper, calcium, chloride, magnesium, manganese and zinc along with a sole carbon source that differs for each environment considered. We considered the 10 different carbon sources based on their biochemical diversity; these included pyruvate (Pyr), D-glucose (Glu), acetate (Ace), D-fructose (Fru), D-galactose (Gal), D-xylose (Xyl), lactose (Lac), inosine (Ino), D-sorbitol (Sor), trehalose (Tre). The random networks generated were analysed for growth in each of the above ten environments using FBA, and a phenotype vector (of length 10) was formed for each of the random networks. The phenotype vector is a bit vector indicating the growth ('1')/no growth ('0') phenotype of each of the random networks in the ten different environments considered.

Identification of synthetic lethals
Synthetic lethals in the metabolic networks were identified using Fast-SL (Pratapa et al., 2015). Fast-SL rapidly and efficiently identifies synthetic lethals by pruning the search space followed by a subsequent exhaustive enumeration of all possible combinations in the pruned space. This enables us to perform synthetic lethality studies at a scale not possible before. Overall, we studied synthetic lethals by efficiently evaluating all possible double deletion phenotypes for 2000 metabolic networks in ten different environments 1983 2 Â 2000 Â 10 % 39 billion double deletions using Fast-SL, identifying hundreds of thousands of synthetic lethals. For the purposes of simulation, we defined lethality as the inability to produce 1% of wild-type growth rate (biomass flux).

Results
We evolved metabolic networks starting from E.coli iAF1260 and performed exhaustive analyses of the single and DLs in the evolved networks. We analysed how the redundancies are comparable with the E.coli network and how they differ in various environmental conditions. Further, we study how different reactions can form functional redundancies. We also compare the lethals in the universe set of reactions with those in the random networks.

E.coli harbours a higher level of redundancy compared with random networks
Identification of lethals in random networks establishes the essentiality of reactions for the growth of organisms in particular growth conditions. Identification of reactions that can form DL sets in random networks helps to identify the different possible compensation mechanisms that can occur. In order to figure out the possible single lethals and DL pairs that could exist in metabolic networks, we generated 2000 random networks with a constraint for growth on glucose. For each of the 2000 random networks, single lethals and DLs were identified using Fast-SL (as described in Section 2). Single and synthetic DL reactions were identified in all 10 environments ( Fig. 1). In glucose, the essential reactions in random networks varied from 229 to 319. The mean number of single lethals that existed in random networks was 262, as against 266 lethals in E.coli iAF1260. Around 20% of the random networks had single lethals between 260 and 265. The E.coli iAF1260 model contains essential reactions significantly higher than the mean (and median) of the essential reactions that existed in random networks. Figure 1A also shows the distribution of single lethals for all the environments in random networks.
We also analysed the synthetic DLs in random networks and E.coli (Fig. 1B). As an illustration, in glucose environment, observations of DLs in random networks showed that there exist networks which have as low as 48 synthetic DLs. Also, the maximum number of DLs that we observed in a random network was 294 (in glucose environment). The mean number of DLs that were identified in the 2000 random networks was 108. Around 17% of the networks had DLs between 90 and 100. However, the real network, E.coli iAF1260, contains 126 DLs, higher than the average number of DLs in random networks. Figure 1B shows the distribution of the DLs found in all the environments compared with those in E.coli.
Also, the fraction of networks that have single lethals and DLs less than the real network was computed as in Table 1. The results show that, in almost all environments, the fraction is greater than 0.5. This suggests that real networks possess a significantly higher fraction of essential reactions and compensation mechanisms than random networks.
To find if the real networks had a higher level of redundancy when compared with random networks, we performed the sign test to compare the number of lethals in E.coli with those that are present in random networks. We found that in the case of single lethals, eight environments had a P-value < 0.05 (Supplementary Table S1), indicating that the number of essential reactions in E.coli is significantly higher than that in random networks, with exceptions in the environments acetate and inosine. As in the case of DLs, the sign test resulted in a P-value < 0.01 (Supplementary Table S1) for nine environments, signifying that the real network has more number of redundant reactions when compared with random networks with the exception being that of lethals in acetate environment. Acetate is an outlier since we could identify only four networks that could show growth out of the 2000 random networks generated on glucose environment.
We further went ahead to identify how many of the DLs in E.coli exist in random networks. Multiple scenarios are possible: first, the DL may exist as is in the random network. In this case,  there is no different compensation mechanism that is present (additionally). Second, one of the reactions in the DL may now be singly lethal, owing to the loss of the partner that could compensate. Third, the DL may not show up at all, since there is probably a higher level of redundancy (e.g. a third compensating reaction) that is now present in the network. Here, we investigated for each of the DLs in E.coli, in which networks are the lethals present together in random networks and if present, are they are found as synthetic lethals. Of the lethals that were present together in random networks, 68% of the lethals existed in random networks. The remaining 32% of the DLs were not found as DLs in the random networks, thus could have been compensated by another reaction forming a higher order lethal. Thus, real networks harbour higher redundancies in comparison with random networks, demonstrating the fact that the real networks are not minimal networks, but evolved to be robust, preserving redundant functions.

Reactions possess very different functional redundancies
An interesting aspect of identifying synthetic lethals in random networks is that we can find (nearly) all possible functionally redundant reactions that could exist in metabolic networks. Understanding the functional redundancies gives insight into the multiple ways organisms exhibit robustness against perturbations. Thus, we aimed to illustrate the existence of multiple synthetic lethal interactions other than those that are known to exist in nature till date. Through our extensive analysis of DLs in random networks, we are able to illustrate a gamut of possible compensations that can potentially occur in metabolic networks. We show how a single reaction can 'partner' with a very large number of different reactions, compensating for single reaction deletions.
We analysed the DLs of the networks that were evolved on glucose, for such partners. Of the 4609 reactions in the universe (excluding the exchange reactions), only 2306 reactions had potential pairs with other reactions. We define a metric for the number of partners that exist for a given reaction as 'Reaction Compensation Index' (RCI). Thus, the RCI for a reaction is the number of partners that this reaction can form a synthetic DL. For instance, an RCI of 200 indicates that the reaction has 200 different partners in DLs, across the 2000 random networks that we sampled. We identified the RCI for all the 2306 reactions; Figure 2 shows the histogram of RCI for all reactions. We find that 82 reactions have an RCI ! 100, 12 of which have an RCI exceeding 200. We also observed that there are 1032 reactions with RCI 10. The reaction with the maximum RCI was found to be the formation of ATP from ADP as below: À cytosol; p À periplasmÞ The above reaction for the formation of ATP, the key energy currency of the cell, from ADP is a reaction that is absolutely necessary for the organism. This reaction is catalysed by ATP synthase, a key enzyme necessary for organisms to synthesize the energy storage molecule of cells.
Further, we analysed the reactions that compensated for ATP synthase. Of the 266 reactions, 35 reactions involved ATP (Supplementary Table S2B), remaining does not involve ATP (Supplementary Table S2A). When identifying the pathways these compensatory reactions were present in, we found that most of the reactions that compensate for ATP synthase belong to the alternate carbon metabolism, glycolysis/gluconeogenesis, and nucleotide salvage pathway. The reactions that compensated predominantly involve pathways whose end products can in a way be converted to ATP.
Further, we analysed the reactions which had RCI of greater than 200 (Supplementary Table S3). Of these reactions, most of them belonged to the central carbon metabolism, and few to the folate metabolism, methylglyoxal metabolism, and anaplerotic reactions. A closer study of these reactions revealed that glyceraldehyde-3-phosphate production produced in two different reactions one from fructose-6-phosphate and the other from fructose-1, 6-bisphosphate were among the twelve reactions with top RCI. The isomerization of glucose-6-phosphate to fructose-6phosphate was one among the twelve reactions that had higher compensation mechanisms. Also, the reaction from the citric acid cycle that involves the inter-conversion of fumarate to malate was among the top reactions that had maximum RCI. We also found the reaction involving the conversion of ribulose-5-phosphate to xylulose-5-phosphate belonging to the pentose phosphate pathway was one among the reactions with the top RCI.
The above findings suggest that compensation occurs in a range of ways, thus enabling the network to sustain growth. There exists a high level of redundancy which further enhances the organisms' robustness to perturbations. It has also been established that compensation can occur in different ways and across different reaction pairs, thus enabling the organism to survive single mutations. Also, a key interesting finding is that most of the reactions that have very  high RCI values are absolutely necessary for growth in glucose medium and metabolic networks tend to evolve in a way that has maximum compensation for such reactions.

How do synthetic lethals vary across environments?
Synthetic lethals across environments has been studied previously in the yeast metabolic network (Harrison et al., 2007). However, a study of synthetic lethals across random networks in different environments helps us to better understand the space of synthetic lethals and how reactions can compensate for one another. We identified DLs in all the 2000 random networks that were evolved on glucose on all 10 different environments under consideration (as described in Section 2). Figure 3 shows the distribution of DLs in a given number of environments; the DLs are dependent on environments in which the network grows on. There were 163 DLs which were identified in all environments thus, establishing the fact that there are some synthetic lethals which occur irrespective of the environments the networks thrive on. A closer look at the DLs in different environments suggests that there exists DLs which are constrained by the environments the networks grow on. Of the unique DLs that were identified on the networks evolved on glucose, 8043 DLs were unique to environments. This implies that these reaction pairs are present only in one environment and absent in the other nine environments.
The environmentally independent DLs can be regarded as 'core' synthetic lethals which tend to compensate for each other under all environmental conditions. On the other hand, there exist lethals that are unique to environments or occur in a few environments that compensate for each other in only a specific growth medium. This is in agreement with the published work on synthetic lethals of yeast on varied environments (Harrison et al., 2007) that establishes two types of synthetic lethals: environmentally independent and environmentally dependent.
On analysis of the 163 DLs present in all environments, seven of these reaction lethals were found as lethals in E.coli in all environments. A surprising finding is that the 163 lethals had reactions nonnative to E.coli acquired from the universe as the networks evolved and still are also common in all environments. Thus, it illustrates that metabolic networks could evolve to acquire reactions that can compensate for the existing reactions in the network under multiple environmental conditions. The number of environmentally dependent DLs was identified as in Table 2. We could observe that inosine medium has the maximum fraction of DLs specific to its environment. The fraction of lethals that are environmentally dependent in acetate environment is the least in comparison to all other environments. This could potentially be since reactions that involve in the growth on inosine environment are different from those networks that assimilate the other carbon sources.
The distribution of synthetic DLs across different environments begs an interesting question: suppose we were to consider the entire universe of reactions as a giant metabolic network, how do the characteristics of the lethals change? Thus, we went ahead to analyse the lethals in the universe of reactions.

Identification of higher order lethals in the universe
An analysis of essential reactions in a universe of reactions has been carried out previously (Barve et al., 2012), which illustrates that single lethals in the universe are also essential in any (random) network they are present in. Such reactions were termed superessential, as they are absolutely required for growth on a particular medium.
We here extended this idea for synthetic lethals to understand if all the synthetic DLs present in the universe occur in random networks. We thus identified the DLs in the universe for glucose environment.
We identified 102 lethal pairs in the universe of reactions, on a glucose environment. These 102 pairs are constituted by 135 reactions. We analysed the number of times each lethal pair in the universe appears in our random networks. We found that 12 lethal pairs were present in more than 400 random networks. Around 14 DLs were present in less than 10 networks. The 102 lethal pairs can be termed 'universal double lethals' since they exist as DLs in the universe, and thus cannot possess further reactions that can compensate for the above reactions.
We then compared the DLs in the universe with those in E.coli. We envisioned three scenarios: first, if both the reactions that form a DL in the universe are present together in E.coli, they would remain a DL. Second, if one of the reactions is blocked, owing to nonavailability of reactants, the other would be a single lethal. Third, if only one of them is present in E.coli, the reaction would end up as a single lethal. Of the 102 DLs in the universe, 16 lethals had both reactions present in E.coli. Of these 16, 12 lethal pairs remained as lethals in E.coli, while in four of them, one of the reactions in the pair was singly lethal. In case of the remaining 86 of the DLs from the universe, only one of the reactions in a lethal pair was present in E.coli, and the reaction was single lethal in each of the 86 cases. It would be very interesting to further extend these analyses to all the random networks studied herein, which might point to other interesting compensatory mechanisms.
We further went ahead to analyse the other DLs that were present in the random networks, beyond these 102 DLs, and if it is possible to use the information on higher order lethals that exist in the universe to identify some of these DLs. To this end, we enumerated the triple lethals (TL) in the universe on glucose environment using Fast-SL, which identified 122 TLs. Given any TL set, the loss of one of the reactions would make it a DL set. Thus, there are a potential 3 2 Â 122 ¼ 366 DL sets that could exist across the random networks. Dropping repetitions, we identify that there are 240/366 unique possible DLs. We analysed the overlap of these enumerated lethals with the DLs in each of the random networks and found that  (8043) were unique to only one environment, while the last bar shows that a small number of lethals (163) were found in each of the ten environments studied here only a fraction of the DLs found in random networks could be found using the TLs in the universe. This could possibly illustrate that the additional DLs that exist in random networks are a subset of further higher order lethals present in the universe.
In summary, the universe tends to have varied compensations, which may exist as such in random networks or predominantly result in lower order lethals. The analyses also reveal the layered compensatory mechanisms that exist in different metabolic networks owing to the different set of reactions present in them. We could find many DLs in random networks that neither exist as DLs in the universe nor are a subset of TLs in the universe, thus pointing towards the existence of higher order redundancies in metabolism.

Random networks illustrate how very diverse reactions can compensate for one another
Compensatory mechanisms are sometimes very straightforward, such as producing the same key metabolite. However, there are instances where the synthetic lethals occur in different pathways, and their compensatory functions are less understood until explored thoroughly. To identify such lethals from the random networks that perform varied functions, we studied some of the lethals identified (in random networks) for their biological significance.
We observed that some widely different reactions, which belong to different pathways, exist as synthetic lethal pairs. For example, the reaction below forms a synthetic lethal pair in random networks while they are absent in E.coli and the universe of reactions. The conversion of chorismate to pyruvate and 4-Hydroxybenzoate, which forms a part of the ubiquinone biosynthesis forms a synthetic lethal pair with the production of 4-Hydroxymandelate from p-hydroxyphenylpyruvate, which is involved in the biosynthesis of antibiotics. When we analysed the enzymes which catalyse these reactions, we found that the former reaction is catalysed by a lyase (EC 4.1.3.40) and the latter by oxidoreductase (EC 1.13.11.46).
Another example of a non-intuitive lethal is the conversion of L-Glutamyl-tRNA from L-Glutamate (ModelSEED: rxn06436) and the reaction which involves in the inter-conversion of succinyl-CoA and glycine to aminolevulinate involved in glycine, serine and threonine metabolism (ModelSEED: rxn00599). The two reactions are catalysed by different types of enzymes, one by a ligase (EC 6.1.1.24) and the other by a transferase (EC 2.3.1.37). A few other examples of such non-intuitive lethals are listed in the Supplementary Table S4.
Thus, we could identify lethals which are absent in E.coli and universe set of reactions but constitute a synthetic lethal pair in random networks. These included many reactions that perform very diverse functions, yet could form a synthetic lethal pair. Our results pave the way to understand how diverse reactions (or pathways) can form functionally redundant pairs.

Discussion
Redundancy is a common mechanism employed in cells to achieve robustness. In the context of metabolism, it is interesting to identify how, under different environmental contexts, very different reactions (or, enzymes) within a cell can compensate for one another.
In this study, we carry out an extensive analysis of metabolic networks, to understand how redundancy can exist and evolve in metabolism. Beginning with the metabolic network of E.coli, we construct 2000 'random' metabolic networks using a previously established approach (Barve and Wagner, 2013). Each of these metabolic networks has the same size as E.coli and exhibits growth on a minimal glucose medium. These random networks essentially mimic the evolution of E.coli on such an environment. Through a comprehensive study of synthetic lethals across these 2000 random networks in 10 different environments, we shed light on how very diverse reactions (or enzymes) can compensate for one another in metabolic networks. Our results also point towards deep redundancy in the organization of cellular metabolism.
Firstly, we see that E.coli harbours a significantly higher number of single lethals and DLs in nearly all of the environments we considered here. Real networks, thus, are not minimal networks, which merely comprise the least number of reactions necessary for survival in the growth medium but have a high level of redundancy. In other words, the real network tends to retain a larger number of functionally redundant reactions in comparison to random networks. This finding demonstrates that nature preserves redundant functions, enabling survival even when they are subjected to single mutations.
Next, we studied the number of possible different compensations that a single reaction can have. We found that around 12 reactions had an RCI of > 200, showing how a very large number of reactions can form a synthetic lethal pair with a given reaction. Our results show that very important essential reactions necessary for growth on a particular growth medium tend to evolve more functional redundancies. We also noticed that the key enzyme ATP synthase exhibited a higher number of lethal pairs (266) indicating the significance of the reaction for the survival of organisms. Here, it must be noted that while the reactions may not necessarily exactly compensate for one another, in all cases, the cell cannot survive unless at least one of the reactions in the pair is active. These reactions may act as important links in their respective pathways, which compensate for one another in the metabolic network.
Another interesting finding is that synthetic lethals do vary across environments, although there are certain reaction redundancies that tend to exist in all environments. This was in accordance with earlier published work of synthetic lethals in the metabolic network of yeast under different environmental conditions (Gü ell et al., 2014), but we additionally observe how the synthetic lethals are shaped, across thousands of random networks that show growth in a given environment.
Our analyses of lethal pairs in the universe of reactions give insight into the layered compensations that occur in different metabolic networks. Also, the enumeration of higher order lethals sheds important light on the different compensatory mechanisms that could exist in any metabolic network. Further, we were able to identify lethal pairs which perform very different distinct functions, yet compensate one another to resist perturbations. These reactions, although they do not have any metabolites in common and involve in very different pathways, form a lethal pair, indicating the essentiality of either of the reactions for survival of the organism.
Although we focus primarily on reactions, the 'unit of deletion' in the model, it is easy to extend this approach to study genes and construct corresponding genetic interaction networks. Our study generally points towards how genes with very different functions (for example, EC 4.1.3.40 and EC 1.13.11.46) can compensate for one another in an organism' s metabolism. The reaction universe is also not necessarily complete; but as and when new reactions are discovered, it will enable us to identify even more novel compensation mechanisms. Although we have listed the possible reaction lethals, a closer look into the pathways would enable us to understand the mechanisms involved in the formation of the lethal pairs. Overall, our work emphasises the adaptive functional redundancies that are acquired, as microbes evolve on a particular environment. We could thus identify novel synthetic lethals and throw light on the wide repertoire of possible functional redundancies, that enhance robustness in metabolic networks.

Funding
Conflict of Interest: none declared.


extfieldsofstudy: Biology Medicine

text:
Comparative Genomic Analysis of Holospora spp., Intranuclear Symbionts of Paramecia

While most endosymbiotic bacteria are transmitted only vertically, Holospora spp., an alphaproteobacterium from the Rickettsiales order, can desert its host and invade a new one. All bacteria from the genus Holospora are intranuclear symbionts of ciliates Paramecium spp. with strict species and nuclear specificity. Comparative metabolic reconstruction based on the newly sequenced genome of Holospora curviuscula, a macronuclear symbiont of Paramecium bursaria, and known genomes of other Holospora species shows that even though all Holospora spp. can persist outside the host, they cannot synthesize most of the essential small molecules, such as amino acids, and lack some central energy metabolic pathways, including glycolysis and the citric acid cycle. As the main energy source, Holospora spp. likely rely on nucleotides pirated from the host. Holospora-specific genes absent from other Rickettsiales are possibly involved in the lifestyle switch from the infectious to the reproductive form and in cell invasion.

INTRODUCTION
Symbiotic associations with bacteria are common to many, if not all, eukaryotes. Endosymbiotic bacteria are well studied in plants, insects, and some other organisms (Moran et al., 2005;Kawaguchi and Minamisawa, 2010;McCutcheon and Moran, 2010;Bennett and Moran, 2013;Becker et al., 2015;Wu et al., 2015;Wexler et al., 2016;Zipfel and Oldroyd, 2017). By the original definition (Oulhen et al., 2016), a symbiotic relationship does not imply benefits for the host. Nevertheless, in symbiotic relationships, symbiont explores its host as a habitat. In turn, host can acquire some necessary products and/or properties which it unable to synthesize, e.g., ammonium produced for plants by nitrogen fixing bacteria, essential amino acids produced by endosymbionts of sap-sucking insects (Moran et al., 2005;Archibald, 2015;Sabater-Muñoz et al., 2017), or vitamins for blood-feeding tsetse fly (Dale and Welburn, 2001). Some symbionts, like some strains of Wolbachia spp., protect the host from pathogenic viruses (Hedges et al., 2008); or are necessary for the maturation of the immune system, as e.g., Wigglesworthia for tsetse flies (Weiss et al., 2011). Caedibacter creates growth advantages for infected paramecia and helps them to out-compete uninfected ones (Kusch et al., 2000). An endosymbiont of ciliates from the genus Euplotes is essential for the host cell proliferation (Vannini et al., 2012. However, most relations are not well understood, and an advantage for the host is not always demonstrated, especially for facultative endosymbionts (Lu et al., 2016;Yañez et al., 2016).
Endosymbiotic bacteria and endonuclear parasites typically have several characteristic features, such as small genome size, low GC content, and short intergenic spacers (reviewed in McCutcheon and Moran, 2012;Batut et al., 2014;Martínez-Cano et al., 2015). Another important trait of these bacteria is an accelerated rate of genome evolution, caused by the minimal gene flow and substantial genetic drift due to the small effective population size of endosymbionts (Marais et al., 2008;Boscaro et al., 2017;Sabater-Muñoz et al., 2017). This leads to formation of pseudogenes and gene loss, and an overall decrease of the metabolic capacity (McCutcheon and Moran, 2012;Boscaro et al., 2017). The loss of genes involved in the reparation drives an even faster evolution of symbiotic species (McCutcheon and Moran, 2012). It also can lead to the decrease in the GC content (Horst et al., 1999;Moran et al., 2005;Long et al., 2018). A stable environment in the host cell does not require complicated regulatory systems, yielding reduction of intergenic regions (McCutcheon and Moran, 2012).
Holospora spp. are endonuclear symbiotic Alphaproteobacteria from the order Rickettsiales that inhabit either macro-or micronucleus of Paramecium spp. It has been suggested that Holospora spp. are parasites, as other Rickettsiales. A recent analysis of 16S rRNA lead some researchers to separate the basal rickettsial lineage of Holospora-like bacteria in a separate order Holosporales (Ferla et al., 2013;Szokoli et al., 2016), but for convenience we discuss both orders here together as Rickettsiales, because they have similar lifestyle and metabolic capacities. The genus Holospora currently is comprised of nine species, each showing clear nuclear and host specificity (Görtz and Schmidt, 2015). Holospora are not able to reproduce outside the host cell but can leave the host in order to invade a new one. Unlike most of the studied symbiotic species, Holospora have a complex life cycle involving two morphologically different forms. The short reproductive form exists only in the host nucleus and multiplies by binary fission, while the long infectious form does not multiply and is capable of infecting new host cells (Gromov and Ossipov, 1981;Fokin and Görtz, 2009). This form shows a unique cytological organization with a pronounced polarity (a large perinuclear space and a special tip), which seems to be of functional significance for the infection process. The infectious form enters new host via Paramecium food vacuole together with food bacteria, but escapes digestion, exits the vacuole and enters the cytoplasm followed by nuclear invasion. Two main mechanisms of this process have been proposed: Holospora either disrupt the digestive vacuole or enter the transport vacuoles (Schweikert et al., 2013). Then Holospora reach the nucleus by utilizing the host's actin cytoskeleton (Fokin and Görtz, 2009;Sabaneyeva et al., 2009). In Holospora obtusa, the 89 kDa periplasmic protein was shown to be associated with cell invasion (Iwatani et al., 2005). The invasion of macronucleus by H. obtusa is associated with release of the 63 kDa periplasmic protein into the macronucleus of the host (Abamo et al., 2008). Additionally, 15 and 39 kDa periplasmic proteins are released from the cell tip during the macronucleus infection . Mechanisms underlying the transition from the reproductive to the infectious form and back are not well understood, but the 5.4 kDa protein with a signal peptide is only detected in the intermediate and infectious forms .
In the presence of Holospora, Paramecium is able to grow, divide and mate (Schweikert et al., 2013). Holospora contribute to the heat-shock resistance in Paramecium caudatum, as cells infected with H. obtusa express high levels of hsp70 mRNA . They may also assist in acquiring high-salt and osmotic-shock resistance to the host (Fujishima, 2009). The presence of H. caryophila in the Paramecium biaurelia nucleus is advantageous in several cell lines during the exponential growth (Bella et al., 2016). Despite these observations of benefits for the host, Holospora have been shown to negatively affect host cells. For example, the presence of Holospora elegans in Paramecium leads to the formation of dysfunctional macronucleus during conjugation (Görtz and Fujishima, 1983). Holospora undulata increases mortality of the host, especially at low-food treatment (Restif and Kaltz, 2006), and high concentrations of the infectious form in the macronucleus inhibit the host cell growth (Fujishima, 2009). The relationship between Holospora spp. and their paramecia hosts seems to be a complex system in which benefits or damages for the host can be highly context-dependent.
The metabolism of Rickettsiales, and in particular of Rickettsia spp., has been studied in detail (Andersson et al., 1998;Hotopp et al., 2006;Fuxelius et al., 2007;Georgiades et al., 2011). Rickettsia export at least 51 metabolites from the host (Driscoll et al., 2017). They are not able to synthesize amino acids, nucleotides, lack glycolysis, and have to import such compounds as coenzyme A, pyruvate, FAD, biotin, etc.
So far, three Holospora genomes have been sequenced, all of which are endosymbionts of Paramecium caudatum: macronucleus-specific H. obtusa and micronucleus-specific H. elegans and H. undulata (Dohra et al., 2013(Dohra et al., , 2014. Analysis of common orthologous genes has yielded 572 single-copy core genes shared by the three genomes, and Holospora have been shown to rely on the host for energy production (Dohra et al., 2014). However, no detailed metabolic pathway reconstruction has been performed.
Here, we report a comparative analysis of four Holospora species, including Holospora curviuscula, a newly sequenced macronuclear endosymbiont of P. bursaria. We propose that Holospora use host-produced nucleotides as its energy source. We also identify the essential compounds that Holospora are likely able to synthesize, and compare their metabolism to that of other Rickettsiales.

Growth Conditions and Genomic DNA Preparation
Holospora curviuscula has an obligate association with its host, Paramecium bursaria, and is therefore uncultivable, so the bacteria have been grown inside host cells. Holospora curviuscula strain NRB217 from Paramecium bursaria isolated in the Novosibirsk Akademgorodok were obtained from the infected clones maintained in the CCCS (Culture Collection of Ciliates and their Symbionts), Research Park, Saint-Petersburg State University. The host cells were cultivated at the room temperature on the lettuce medium inoculated with Enterobacter aerogenes as a food resource for paramecia. The culture of P. bursaria bearing H. curviuscula was concentrated by centrifugation (10 min at 4500 g) and then homogenized using 10% solution of detergent Nonidet P-40 (Sigma-Aldrich Cat No. 21-3277 SAJ). The infectious forms of H. curviuscula were isolated from the homogenate by centrifugation in Percoll density gradient (Sigma-Aldrich Cat No. P1644) as described previously (Rautian and Wackerow-Kouzova, 2013). Genomic DNA was isolated with the DNeasy Blood and Tissue kit (QIAGEN Cat No. 69504) using a modified protocol -the time of incubation of cell homogenate with ATL buffer and proteinase K was increased to 16 h. All the subsequent steps were performed according to the standard Quick-Start Protocol.

Genome Assembly and Annotation
Two libraries were generated, paired-end MiSeq Illumina library (2 × 250) (PE), and mate-pair library with insertion size 4 kbp (MP). The initial MP library size was 2.5 million reads, and the PE library size, 4.5 million. MP reads were filtered with NextClip v1.3 (Leggett et al., 2014) and only category A pairs of reads (both reads in a pair contain adapters) were selected for further analysis. Both MP and PE were filtered by quality with trimmomatic version 0.33 (Leading 15, sliding −15, slidingwindow −4:25), after filtering 1.3 million MP reads and 1.7 million PE reads were retained for assembly. Processed reads were assembled with platanus version 1.2.1. The Whole Genome Shotgun project has been deposited at DDBJ/ENA/GenBank under the accession PHHC00000000. The version described in this paper is version PHHC01000000.
To estimate the completeness of the Holospora assemblies, we used HMMer to search for essential PFAM domains (Rinke et al., 2013). A domain was considered to be present, if it was found by hmmsearch (with default parameters) (Mistry et al., 2013) and the bias was lower than the hit score. For all missing genes we checked whether they were present in 60 genomes of Rickettsiales and 23 other endosymbiotic Proteobacteria with tiny genomes less than 0.8 Mb (the complete list of genomes is given in Supplementary Table 1).

16S rRNA-Based Phylogeny
To reconstruct the phylogenetic tree of Holospora spp., the PhyML software (Guindon et al., 2010) was used. 16S rRNAs from H. curviuscula NRB217, H. elegans E1, H. obtusa, and H. undulata HU1 were extracted from the full genome sequences, while 16S rRNAs from H. acuminata and Candidatus Gortzia infectiva were obtained from the SILVA database (accession numbers: KC164379 and HE797907, respectively). To root the tree, we used 16S rRNA of Rickettsiales (SILVA accession number: CP009217). The phylogenetic tree was constructed with PhyML v. 3.1 with 100 bootstrap replicates, GTR substitution model, and -o tlr parameter.

Analysis of Repeats
To estimate the fraction of repetitive DNA in the Holospora genome assemblies, reads were realigned with Bowtie 2 (-end-toend mode) (Langmead and Salzberg, 2012) back to the respective assemblies, and the coverage for each nucleotide was calculated. In the case of H. curviuscula, we used PE reads generated in this study, reads for H. obtusa and H. undulata were downloaded from DDBJ (accession numbers DRP001203 and DRA001008, respectively). A fragment was presumed to be duplicated if it was longer than 20 bp, and its coverage was at least 1.6-fold higher than the median contig coverage.

Orthologous Groups and Genome Comparisons
To investigate the phylogenetic relations between H. undulata and H. elegans, pairwise genome alignments of all Holospora genomes were constructed with MAUVE (snapshot 2015-02-13) (Darling et al., 2010). Alignments of orthologous genes from the essential list (see above) were extracted from the MAUVE output, and all gaps were removed. The number of substitutions was calculated for each gene for each pair of genomes.
To investigate the gene repertoire of all Holospora and other Rickettsiales, we constructed groups of orthologous genes using OrthoMCL v. 2.0.9 (Li et al., 2003) with the MCL software v. 14-137. The output gene clusters were grouped in five subgroups using an ad hoc perl script which considered orthologs, co-orthologs and inparalogs. All Rickettsiales and other Holospora genomes (for the complete list see Supplementary  Table 2) were downloaded from NCBI Genbank (Benson et al., 2013).

Holospora curviuscula Genome
The H. curviuscula draft genome sequenced here is comprised of 152 scaffolds (210 contigs) with N50 of ∼39 kbp, total length of 1.7 Mb, and GC content of 37.6% ( Table 1). The longest scaffold is 153367 bp and contains 126 CDSs. We predicted 1594 genes, including 1555 protein-coding genes with the average encoded protein length of 218 aa (Supplementary Figure 1). We assigned protein function to 683 genes. The genome assembly contains all rRNA genes and a set of 36 tRNAs necessary to recognize all codons. The 16S rRNA analysis revealed that the sequenced genome indeed belongs to the genus Holospora, and specifically to H. curviuscula, and is closer to Holospora acuminata, another endosymbiont of P. bursaria, than to Holospora from P. caudatum (Figure 1). The H. curviuscula genome is the largest Holospora genome sequenced so far ( Table 1). Its GC content is slightly higher than that of most other non-Rickettsiales obligate endosymbionts (McCutcheon and Moran, 2012), but is typical for Rickettsiales (Supplementary Figure 2).
We estimated the quality of the assembly by searching for nearly universal bacterial genes (see Materials and Methods). Among the 139 universal PFAM domains, 121 are present in the H. curviuscula genome. Nine of the 18 missing domains have been found with a lower similarity threshold, which makes their presence uncertain. The remaining nine of the 18 missing domains are absent in all other Holospora as well. Moreover, these domains are also absent either in some endosymbiotic bacteria with small genomes or in some Rickettsiales genomes (Figure 2), meaning that their absence can be tolerated. We further searched for these 18 missing domains in unassembled sequencing reads of H. curviuscula (14% of all reads), and found no indication of their presence. We manually checked 615 ORFs for frameshifts, and found only nine ORFs with frameshifts. This number may be an underestimate, because H. curviuscula genome contains short ORFs with unknown functions that have no homologs outside Holospora and can be remnants of some ancestral genes. Still, it implies that the H. curviuscula genome contains few pseudogenes.
Rickettsiales have been reported to have an atypical rRNA operon (Andersson et al., 1998). However, all Holospora have a standard rRNA operon, with 23S and 16S genes located close to each other and separated by several tRNA genes.

Reanalysis H. undulata and H. elegans Genomes
Holospora undulata HU1 and H. elegans E1 extracted from the micronucleus of P. caudatum have been recognized as separate species based on differences in the cell morphology (Görtz and Schmidt, 2015). However, our comparison of 16S rRNA genes from the available genome assemblies (Genbank IDs NZ_ARPM00000000 and NZ_BAUP00000000, respectively) of these species revealed just a single-nucleotide difference (Figure 1), yielding > 99.9% sequence identity, which is above the common thresholds of 97 or 98.7% used to define species (Janda and Abbott, 2007). That led us to inquire whether the genomes of H. undulata and H. elegans indeed represented different species or just strains of the same species. We constructed whole-genome pairwise alignments of the Holospora genomes and counted mismatches in the essential genes (Figure 3). The H. elegans-H. undulata pair carried ∼100-fold fewer mismatches than all other genome pairs.

Holospora spp. Have Multiple Repeat Sequences and Contain Prophages
All available Holospora assemblies are comprised of multiple contigs, and reassembly has not led to better assemblies even with added genome coverage. Moreover, the genome length of H. curviuscula has been estimated at ca. 2 Mb (unpublished data), while the total assembly length is around 1.7 Mb. This suggests that Holospora genomes may contain a considerable fraction of repetitive DNA. To test this, we searched for fragments with high read coverage (see Materials and Methods). In H. curviuscula, we found ∼300 repeats with the total length of ∼300 kb ( Figure 4A) and the estimated copy number varying from 2 to 21 ( Figure 4B). In H. obtusa and H. undulata, the similarly estimated lengths of repetitive DNA were 60 and 161 kb, respectively. We were not able to estimate the fraction of repetitive DNA for H. elegans, as no raw sequencing data were available. The calculated repeat lengths could be underestimates, as we have applied strict coverage cutoffs yielding conservative repeat boundaries. For the same reason, we could overestimate the number of repeats by splitting long repeats into individual repeats in regions of low coverage. Still, the match between the difference in the assembly length and the genome length and the total coverage of repetitive DNA suggests that unassembled repeats comprise ∼15% of the genome in H. curviuscula. The high copy-number repetitive regions in H. curviuscula include short ORFs of  Table 1). Color code: dark blue, the domain is present; gray, the domain is absent; cyan represents domains found with decreased similarity threshold (see Materials and Methods). The number in each cell is the number of genomes where the domain has not been found. possible transposases. In addition, the identified repeats include two copies of rRNA operons with the similar structure, and also genes encoding prophage-like proteins.
To characterize prophages further, we searched for complete prophages and, additionally, for phage-related PFAM domains in all Holospora genomes. We found two possible prophage regions in H. elegans, one in H. obtusa, and two in H. undulata. One of the H. undulata prophages was intact, contained all proteins necessary for the phage replication and was surrounded by integration sequences (Supplementary Figure 3); the remaining copies were severely disrupted and lacked insertion sequences. Although no intact prophages were found in H. curviuscula, scaffold565 had a locus with two DDE superfamily endonucleases and two putative transposases, and other scaffolds carried genes encoding phage capsid and phage portal proteins similar to the ones found in other Holospora as well as some full-length and fragmentary transposases (data not shown).

Holospora spp. Cannot Synthesize Amino Acids and Some Other Essential Small Molecules
Holospora have a reduced metabolism even by the standards of the already gene-poor order Rickettsiales (Georgiades et al., 2011). In particular, Holospora lack genes encoding enzymes of the citric acid cycle (CAC) (Supplementary Figure 4). By contrast, Rickettsia and other Rickettsiales have at least some of these enzymes and use them to convert small molecules (Driscoll et al., 2017).
All sequenced Holospora lack amino acid synthesis pathways. For most amino acids, all enzymes of the pathway are missing. The exceptions are the partial pathways for the biosynthesis of L-tryptophane, L-lysine, and L-glutamate and for the conversion of L-alanine to D-alanine, an essential compound for formation of the cell wall ( Figure 5A). This means that all amino acids have to be imported from the host. While Rickettsia also cannot produce amino acids and probably import them from the host (Driscoll et al., 2017), other Rickettsiales are capable of synthesizing some amino acids such as L-lysine and L-glutamine (Fuxelius et al., 2007). Similar to Rickettsia, Holospora cannot produce chorismate or use it as a substrate for downstream reactions, and need to import it or its derivates from the host.
Like Rickettsia, Holospora are not capable of purine and pyrimidine synthesis (Figure 5B), although, also like Rickettsia, they can convert UTP to ÑTP. Therefore, they need to obtain all nucleoside triphosphates from the host.  On the other hand, Holospora carry genes for ribonucleotide reductases and therefore can convert ribonucleotides into deoxyribonucleotides.
The ubiquinone biosynthesis pathway is also partial, as it is comprised of only dimethylallyl diphosphate (DMAPP) and the enzymes downstream of it (Supplementary Figure 6), similarly to what has been reported for Rickettsia (Driscoll et al., 2017). Unlike Rickettsia, Holospora are unable to convert DMAPP to isopentenyl diphosphate, although they carry geranyldiphosphate synthase which is missing in Rickettsia.
Therefore, DMAPP seems to be obligatorily imported from the host.
The only major metabolic pathway that is almost intact is the fatty acids synthesis (Supplementary Figure 5).
In order to determine how the missing compounds are imported from the host, we analyzed predicted transporters. In H. curviuscula, thirty transport-associated genes have been identified. As expected, Holospora genomes encode oligopeptide and amino acid transporters, as well as proteases, some of which are periplasmic. In particular, we have found arginine/glutamine, proline, acetyl-serine/cysteine, choline/glycine/betaine, and putative branched amino acids candidate transporters. In addition, we found putative transport systems for magnesium, ferric ions, ribose, purines, sulfoacetate, L-galactonate or other sugars, and putrescine. However, these transport systems are insufficient to deliver all missing compounds to Holospora; in particular, it is not clear, how the remaining amino acids are delivered, unless as components of oligopeptides.

Holospora spp. Use Nucleotides as the Main Energy Source
All Holospora lack most genes involved in energy production. In particular, all enzymes necessary for the glycolysis except phosphoglycerate mutase ( Figure 5C) and all enzymes of the citric acid cycle except malate dehydrogenase (Supplementary Figure 4) are missing. The F 1 F 0 -ATPase is also missing. In addition, Holospora cannot produce coenzyme A from scratch, although as Rickettsia, they seem to have CoaE (PF01121, Figure 2) and thus are able to synthesize coenzyme A from dephospho-CoA (Driscoll et al., 2017). Of the pentose phosphate pathway, only the non-oxidative branch is present, as well as the downstream enzyme ribose-phosphate mutase that converts phosphorybosil pyrophosphate (PRPP) to sugars ( Figure 5C); the energy-producing oxidative pathway leading to the ribulose-5phosphate is missing. By contrast, all Holospora have the pyruvate dehydrogenase complex and are able to convert pyruvate to acetyl-CoA, and further to acetoacetyl-CoA and acetoacetate, with the production of ATP ( Figure 5D). All Holospora have a set of ribonucleotide reductases, which would allow them to use either nucleotides or ribonucleotides as an energy source. No obvious source of energy other than nucleotides was found.

Secretory Systems and Putative Invasins
Although most Rickettsiales are parasites and their genomes are highly reduced, they retain secretory systems such as the Tat and Sec pathways, a type IV system, and the TolC protein from a type I secretion system . Moreover, the VirB protein from a type IV secretory system is thought to be essential for the host invasion in most Rickettsiales Gillespie et al., 2016). By contrast, all studied Holospora demonstrate a significant decay of secretion pathways. They still possess the complete Sec system, additional systems helping to translocate proteins to the outer membrane (LolA, LolD, and possibly LolE) or transport them outside the cell (BamA, BamB, BamD, and chaperone DegP), and a TolC-like protein. However, we found no proteins similar to components of the Tat-system (TatA, TatB, or TatC). Since Tat-system proteins may be hard to identify , we performed a genome-wide search for proteins with the twin-arginine signal required for the recognition by the Tat-system and found no significant hits. Together, these observations show that the Tat transport system is indeed missing, which means that Holospora are unable to export folded substrates.
In Rickettsia, invasion is associated with proteins RalF, RickA, and Sca  which are missing in Holospora, implying that Holospora use other mechanisms to invade the host cell. There are also no proteins with ankyrin domains assumed to play a role in the pathogenesis in Rickettsia . Although Sec proteins are present in all Holospora, only one protein with a putative autotransporter domain (PF03797) was found in H. curviuscula (HCUR_00103). Even though this gene has homologs in other Holospora, this domain was not predicted in them, implying that it can be a false positive. However, we found multiple (3-9 copies per genome) genes similar to ompA in all Holospora. OmpA has been previously shown to be involved in the pathogenesis by Rickettsiales (reviewed in Palmer and Noh, 2012).

Holospora-Specific Genes
To search for genes that could be essential for the survival in the nucleus, we compared orthologous groups (OGs) for all Holospora with a set of all complete Rickettsiales genomes (Figure 6). We identified 102 OGs that occur in all Holospora and do not occur in any other Rickettsiales species (hereafter, Holospora-specific OGs, HOGs), 97 of HOGs contained no paralogs. HOGs encoded transporters, regulators, and hypothetical proteins. We expect proteins essential for survival and adaptation to the host either to be located on the surface of the bacterial cell or to be secreted. Therefore, we focused on those proteins that contained transmembrane helices or signal peptides. Thirty seven HOGs contained proteins with predicted transmembrane helices; 12 other OGs, proteins with predicted classic signal peptides; eighteen OGs were predicted to be secreted by the non-classical pathway (Supplementary Table 3). One of the latter HOGs contained the 89 kDa periplasmic protein which has been reported to be associated with cell invasion in H. obtusa (Iwatani et al., 2005). One HOG contained proteins similar to OmpA. These proteins are highly conserved among Holospora, and all Holospora except H. curviuscula carry two paralogs of this gene.
FIGURE 6 | Groups of orthologous proteins in Holospora spp. and Rickettsiales. Numbers in cells indicate either the number of orthologous groups, or, in the case of unique genes, the number of singletons. Numbers in bold represent groups of orthologous groups common for either macronuclear (8) or all (102) Holospora. OGs common for the macronuclear endosymbionts are two OGs encoding alpha/beta hydrolases; one, a possible aspartate/glutamate racemase or a malate isomerase; one, nucleotide sugar epimerase; one OG had an unknown function; and three remaining OGs encode short proteins, of which one was similar to the DDE superendonuclease, and another one has a predicted signal peptide. Of 102 Holospora-specific OGs: 30 are secreted, 37 contain putative transmembrane helices; 17 HOGs contain short peptides, of which 5 are potentially secreted.
Holospora genomes encode numerous short (<100 aa) peptides with unknown functions. Short proteins have been shown to be involved in spore formation, regulation of transport, regulation of transcription, and signal transduction, or to possess antimicrobial or other toxic activities (Wang et al., 2008;Storz et al., 2014). Of the 102 HOGs, seventeen contained short proteins, including five HOGs that contained proteins with secretory signal peptide sequences. One of these HOGs has been described earlier as the 5.4 kDa protein involved in the switch between the reproductive and infectious forms .
To determine whether there are proteins that could determine the nuclear specificity of Holospora, we searched for HOGs specific to the two macronuclear species, H. curviuscula and H. obtusa. We found eight HOGs present in both H. curviuscula and H. obtusa and absent in H. undulata, H. elegans and in other studied Rickettsiales. Of these, two HOGs encoded alpha/beta hydrolases; one, a possible aspartate/glutamate racemase or a maleate isomerase; one, nucleotide sugar epimerase; one HOG had an unknown function; and three remaining HOGs encoded short proteins, one of which was similar to DDE superendonuclease, and another one had a predicted signal peptide.

DISCUSSION
Here we report a comparative analysis of the newly sequenced genome of H. curviuscula NRB217 and other available Holospora genomes. The phylogenetic analysis of all available Holospora has confirmed (Rautian and Wackerow-Kouzova, 2013) that they are clustered by the host species rather than by micro-vs. macronuclear specificity, so that H. curviuscula is the outgroup to the previously analyzed genomes. Furthermore, the fact that we found only one substitution in the 16S rRNA gene and few substitutions genome-wide between the sequenced samples of H. elegans and H. undulata suggests that these are in fact strains of the same species.
The previous genomic analysis of Holospora spp. (Dohra et al., 2014) included only a general description of the available COG categories and some reconstruction of the metabolism. It has been observed that many pathways are missing in Holospora. The pyruvate dehydrogenase complex that is present in Holospora has been proposed to be a possible relic of an ancestral pathway. Additionally, it has been suggested that Holospora strongly relies on the host for energy production.
Addition of H. curviuscula allows for a detailed comparative genomic analysis of the metabolism of Holospora spp. We found that all Holospora have reduced metabolic capacities and have to import many metabolites from the host.
Even though the genomes of Holospora are relatively large ( Table 1) in comparison with other symbiotic species such as Buchnera, Candidatus Baumannia, or Candidatus Carsonella, Holospora are unable to produce most of the essential compounds. Furthermore, all available Holospora genomes seem to contain a large fraction of repetitive sequences, which complicates the genome assembly. Arguably, as all available Holospora genomes consist of multiple contigs, some of the enzymes could have been missed. However, quantitative analysis suggests that all universal genes are indeed present and that nearly all non-repetitive regions are contained in contigs; any potential missing part of the genome would have to be small, and cannot account for a large number of missing genes. Furthermore, our analysis is based on independent assemblies of multiple moderately related species, and it is unlikely that the same gene would be missing in several assemblies.
While smaller endosymbionts can produce at least some amino acids or retain some parts of the central metabolism, Holospora are unable to synthesize any amino acid, and have to acquire them from the host, despite the fact that some relics of pathways are seen, such as a partial pathway for the tryptophan synthesis. Holospora lack glycolysis, the Entner-Doudoroff pathway, and the pentose phosphate pathway. Surprisingly, Holospora have no enzymes of the citric acid cycle, which is unusual for Rickettsiales, as even the most reduced Rickettsia retain it.
Although we have not performed a complete metabolic reconstruction, it is evident that a broad range of compounds, including amino acids, DMAPP, and chorismate derivates have to be imported from the host for Holospora to survive. It is unlikely that these compounds are accumulated by Holospora during its brief residence outside the nucleus. Indeed, the process of Holospora infection is a rapid one, whereas Holospora stay in the nucleus for long periods of time, which means that all these compounds have to be acquired from the host. Amino acids can be obtained by degrading host nuclear proteins, and indeed the Holospora infection is associated with increased proteolytic activity in the macronucleus (Freiburg, 1985). The situation with other compounds is more intriguing. All these small molecules can passively diffuse through the nuclear pore complex (Knockenhauer and Schwartz, 2016). Further, Holospora alter the host's gene expression and increase RNA synthesis (Freiburg, 1985;Fujishima, 2009). These alterations may help bacteria to acquire the necessary nutrients. To understand how the missing nutrients are delivered to Holospora, we searched for transport proteins and attempted to predict their specificity. This has explained some, but not all of the missing nutrients, and this topic has to be investigated further.
We propose that Holospora not only relies on the host for the energy production, but specifically use nucleotides or ribonucleotides as the energy source. Indeed, they are able to interconvert them and to convert UTP to CTP, and putative ribose transport proteins, which can also transport nucleotides, are present. We suggest that ribonucleotides are the preferred energy source, as the intracellular abundance of ribonucleotides can be 10-to 100-fold higher than that of dNTPs (Traut, 1994). This seems natural, given the nuclear habitat of Holospora.
The interactions between Holospora and their hosts have been studied in some depth. Holospora infection can either decrease or increase the Paramecium viability under a variety of conditions Hori et al., 2008;Fujishima, 2009;Bella et al., 2016). Addition of Holospora curviuscula to the analysis allowed us not only to investigate proteins specific to Holospora in a particular host, but also to find proteins that may play a role in the macronucleus infection. While the details of the mechanism remain unclear, it is likely that the interaction and stable infection rely on Holospora-specific proteins, in particular, secreted short conserved peptides. The 89 kDa periplasmic protein shown to be involved in the H. obtusa infection (Iwatani et al., 2005) is conserved in all studied Holospora spp., which indicates its importance. Finally, we propose that OmpA-like proteins may be involved in the Holospora invasion.
Overall, the analysis of the Holospora genomes demonstrates that their metabolic capabilities are unusually restricted, especially given the genome size, and provides a list of candidate genes essential for their unique lifestyle.

AUTHOR CONTRIBUTIONS
SG, MR, and MG designed the research. AB and MR isolated and cultivated bacteria. ML prepared the sequencing libraries and sequenced the genome. ML and SG assembled the genome. SG, AB, DM, and MG annotated the genome and performed the comparative analysis. SG, AB, MR, and MG wrote the paper.

FUNDING
This study was supported by the Russian Science Foundation under grant 14-50-00150.