extfieldsofstudy: Biology Medicine

text:
Genome-scale Co-evolutionary Inference Identifies Functions and Clients of Bacterial Hsp90

The molecular chaperone Hsp90 is essential in eukaryotes, in which it facilitates the folding of developmental regulators and signal transduction proteins known as Hsp90 clients. In contrast, Hsp90 is not essential in bacteria, and a broad characterization of its molecular and organismal function is lacking. To enable such characterization, we used a genome-scale phylogenetic analysis to identify genes that co-evolve with bacterial Hsp90. We find that genes whose gain and loss were coordinated with Hsp90 throughout bacterial evolution tended to function in flagellar assembly, chemotaxis, and bacterial secretion, suggesting that Hsp90 may aid assembly of protein complexes. To add to the limited set of known bacterial Hsp90 clients, we further developed a statistical method to predict putative clients. We validated our predictions by demonstrating that the flagellar protein FliN and the chemotaxis kinase CheA behaved as Hsp90 clients in Escherichia coli, confirming the predicted role of Hsp90 in chemotaxis and flagellar assembly. Furthermore, normal Hsp90 function is important for wild-type motility and/or chemotaxis in E. coli. This novel function of bacterial Hsp90 agreed with our subsequent finding that Hsp90 is associated with a preference for multiple habitats and may therefore face a complex selection regime. Taken together, our results reveal previously unknown functions of bacterial Hsp90 and open avenues for future experimental exploration by implicating Hsp90 in the assembly of membrane protein complexes and adaptation to novel environments.

Introduction
In eukaryotes, the universally conserved and essential chaperone Hsp90 aids the folding of key proteins in development and responses to environmental stimuli [1][2][3]. In yeast, up to 10% of all proteins are estimated to be Hsp90 clients under standard culture conditions [4]. Hsp90 function is even more important under stressful conditions that challenge protein folding, such as increased temperature [5]. The activity of eukaryotic Hsp90 is further modulated by various co-chaperones, which confer substrate specificity and alter protein folding kinetics [2,5]. Depletion of eukaryotic Hsp90 in vivo increases phenotypic variation, reveals 'cryptic' heritable variation, and increases penetrance of mutations [6][7][8][9]. Accordingly, eukaryotic Hsp90 enables organisms to maintain a stable phenotype in the face of environmental and genetic perturbation and to correctly interpret environmental stimuli.
In stark contrast, in prokarya, Hsp90 is not essential [10] and many bacterial genomes lack Hsp90 altogether [11]. Among Archaea, only very few species contain Hsp90, and those are thought to have gained Hsp90 horizontally from bacteria [11,12]. This fragmented phylogenetic pattern likely results from multiple independent gains and losses, though phylogenetic reconstructions are confused by ancient Hsp90 paralogy [11,12]. At the amino acid level, the Escherichia coli Hsp90 (High-temperature protein G or HtpG) is 42% identical to its human homolog, suggesting strong stabilizing selection consistent with functional conservation [13]. Indeed, E. coli Hsp90 appears to retain generic protein chaperone activity [14] and homologous Hsp90 mutations cause chaperone defects in both the prokaryotic E. coli and eukaryotic yeast [15]. However, there are no identified obligate Hsp90 co-chaperones in bacteria, adding to the uncertainty regarding the extent of its client spectrum and specificity.
To date, only three proteins have been implicated as Hsp90 clients in bacteria, with non-overlapping functions in ribosome assembly, the assembly of light-harvesting complexes, and the CRISPR/Cas immunity system [16][17][18]. Several other proteins have been shown to physically interact with the chaperone [19,20]. Together with our knowledge of eukaryotic Hsp90 function, these data have given rise to the speculation that Hsp90 may facilitate the assembly of oligomeric protein complexes in bacteria, much like it does in eukaryotes [21]. Unlike in eukaryotes, however, further exploration of Hsp90's functional role in bacteria has proven challenging because there are no pleiotropic Hsp90-dependent phenotypes.
To address this challenge, we used a genome-scale coevolutionary 'guilt-by-association' approach [22,23] to explore the spectrum of conserved Hsp90-associated genes, functions, and organismal traits. Hsp90-associated genes tended to function in flagellar assembly, chemotaxis, and secretion. Consistent with these functions, Hsp90-associated organismal traits included the ability to inhabit multiple environments. To add to the sparse list of known bacterial Hsp90 clients, we further developed a statistical method to predict putative Hsp90 clients, which included flagellar, ribosomal, and chaperone proteins. We validated our predictions experimentally, focusing on two candidates functioning in motility and chemotaxis. Indeed, both the flagellar protein FliN and the kinase CheA were found to be Hsp90 clients in vivo. Our findings demonstrate the power of co-evolutionary inference to correctly identify substrates and functions of conserved genes like bacterial Hsp90.

Hsp90 paralogs in bacteria
Our method for inferring the function of bacterial Hsp90 is based on the analysis of its distribution across the bacterial phylogeny. However, this analysis is complicated by the existence of multiple ancient Hsp90 paralogs in bacteria. These paralogs may be older than existing phyla in bacteria [11,12], and may have evolved distinct functions on this enormous time scale. To address this issue and to identify each paralog, we first clustered bacterial Hsp90s by sequence identity. We identified 897 bacterial Hsp90 protein sequences in the KEGG database [24] and built a neighbor-joining gene tree of bacterial Hsp90s ( Figure S1A-B). We observed two well-supported long-branching clades as well as several less confident divisions in the tree ( Figure S1B). These two long-branching clades contain sequences corresponding to the 'hsp90B' and 'hsp90C' paralogs that were described previously [11,12]. All other branches correspond to 'hsp90A' [11], which is the largest of the Hsp90 families in bacteria ( Figure S1C, Text S1). Notably, hsp90A is the lineage out of which all eukaryotic Hsp90s (excluding mitochondrial and chloroplast Hsp90s) are derived. Moreover, the E. coli gene htpG belongs to the hsp90A family, and its gene product is the best-studied bacterial Hsp90 protein. For these reasons, we restricted our analysis to hsp90A.
Genome-wide detection of genes co-evolving with hsp90A We set out to identify orthologous groups whose presence and absence profiles across bacterial species are associated with the presence and absence profile of hsp90A. To avoid spurious associations, any such comparative analysis must go beyond a naïve comparison of presence/absence patterns across genomes and incorporate phylogenetic information [25]. To this end, we used BayesTraits [26][27][28], a computational framework for phylogenetic analysis of character evolution. Given the states (e.g., presence/absence) of two characters across some set of species and a phylogenetic tree relating these species, BayesTraits evaluates the likelihood of various evolutionary models throughout the tree. This approach can be utilized, for example, to determine whether these two characters evolve in a mutually dependent vs. an independent fashion.
We used BayesTraits to detect associations between hsp90A and 4646 other orthologous groups in bacteria (which hereafter we shall refer to as 'genes' for simplicity). We used the tree constructed by Ciccarelli et al. [29] as a model phylogeny ( Figure 1). In this initial analysis, we tested for any kind of dependency between hsp90A and other genes, and did not make specific assumptions about the nature of the relationship between hsp90A and the genes in question [28]. Specifically, we compared a model in which the

Author Summary
Hsp90 is a chaperone protein that aids the folding of many other proteins (clients), which tend to be signal transduction proteins. Hsp90 is particularly important when organisms are under environmental or mutational stress (e.g. in cancerous cells). Although Hsp90 is well-studied in eukaryotic species from yeast to humans, little is known about its counterpart in bacteria. To address this challenge, we analyzed the presence and absence of thousands of genes across numerous bacterial species and identified genes that co-evolved with Hsp90. These genes provide insights into potential functions of Hsp90 in bacteria. We found that Hsp90 co-evolves with membraneassociated protein complexes such as the flagellum and that Hsp90 is associated with a preference for inhabiting multiple habitats. We extended our analysis to identify genes that exhibit evolutionary dynamics characteristic of Hsp90 clients. Many of the putative clients were involved in flagellar assembly, suggesting a crucial role of Hsp90 in the regulation of bacterial motility. We experimentally confirmed that E. coli Hsp90 interacts with selected candidates and demonstrated Hsp90's role in flagellar motility and chemotaxis. The computational approach described here, identifying novel functions and specific clients of bacterial Hsp90, further provides exciting starting points for research in bacterial chaperone biology.
rate of gain and loss of a given gene is independent of the rate of gain and loss of hsp90A (independent evolution) vs. a model in which the rate of gain and loss of this gene is affected by the presence or absence of hsp90A or vice-versa (co-evolution).
In total, we found 327 genes that co-evolve with hsp90A (Dataset S1). We will refer to this set as hsp90A-associated genes. These hsp90A-associated genes were significantly enriched for annotations related to the flagellum and to bacterial secretion systems (Table 1). Moreover, out of the 16 hsp90A-associated bacterial secretion genes, 10 were part of the non-flagellar Type III secretion system, suggesting that hsp90A is associated specifically with this system rather than with secretion systems in general. Using a different and markedly more extensive phylogeny [30] provided similar results (see Text S1, Table S1), as did a pruned Ciccarelli tree without the species containing the hsp90B or hsp90C (see Text S1).

Characterization of co-evolutionary dynamics
The associations of hsp90A with other genes identified above are agnostic to the specific nature of the dependency between hsp90A and the gene in question. For example, our initial analysis could not distinguish between a positive association (i.e. genes tend to be gained and lost together) and a negative association (i.e. genes tend not to co-occur in genomes). Similarly, this analysis did not distinguish between genes whose gains and losses are affected by the presence of hsp90A (but that do not themselves affect hsp90A evolution) and genes that exhibit mutually dependent dynamics with hsp90A. Without a quantitative estimate of the effects that hsp90A and its co-evolving partners have upon one another, inference of Hsp90A function and its relationship with other genes is challenging.
To characterize the specific nature of the dependency between hsp90A and hsp90A-associated genes, we therefore examined rates of gain and loss inferred by BayesTraits. We focused on the two major non-overlapping hsp90A-associated functional categories, flagellar assembly and bacterial secretion. Considering, for example, fliI, a representative flagellar gene, we found that its gain and loss was strongly affected by the presence of hsp90A.
Specifically, in the presence of hsp90A, fliI was often gained and rarely lost, whereas it was rarely gained and often lost when hsp90A is absent (Figure 2A). This pattern was common to all hsp90Aassociated flagellar genes ( Figures 2C, S2), suggesting a positive association between hsp90A and flagellar genes throughout evolution. In contrast, the co-evolutionary relationship between hsp90A and yscN, a representative nonflagellar type III secretion system gene, was markedly different, with yscN presence strongly affecting the gain and loss of hsp90A ( Figure 2B). Specifically, the presence of yscN was associated with a large increase in the rates of gain and (even more dramatically) loss of hsp90A relative to these rates in its absence. Again, this pattern was common to all hsp90Aassociated bacterial secretion genes ( Figures 2D, S3, S4), suggesting a negative association between hsp90A and nonflagellar secretion genes throughout evolution.
To further validate the fundamentally distinct co-evolutionary dynamics of these two groups of genes, we considered four different co-evolutionary models: (1) hsp90A and the gene in question are independent (null); (2) hsp90A and the gene in question are mutually dependent; (3) hsp90A is dependent on the gene in question but not vice versa, and (4) the gene in question is dependent upon hsp90A but not vice versa (Methods). We used the Akaike Information Criterion (AIC [31]) to determine which of these 4 models best fit the co-evolutionary dynamics of each hsp90A-associated gene. As expected, none of the hsp90A-associated genes fit the independent model. Of the 27 hsp90A-associated flagellar genes, 25 were classified as being dependent on hsp90A but not vice-versa (model 4). Of the 16 hsp90A-associated secretion system genes, 10 genes were classified as mutually dependent with hsp90A (model 2; 6 of which were Type III secretion system genes), whereas 6 were classified as affecting the evolution of hsp90A (model 3). Furthermore, considering all hsp90A-associated genes, we found that genes that best fit each of the evolutionary dependency models above (models 2, 3, and 4) were enriched for different functions (Table 1). Specifically, among genes dependent on hsp90A, flagellar motility was strongly enriched, whereas among genes mutually dependent on hsp90A, secretion system components were enriched. Taken together, these patterns suggest that flagellar genes and secretion system genes had markedly different regimes of co-evolution with hsp90A.

Prediction of Hsp90A clients
Although many genes exhibited distinct patterns of co-evolution with hsp90A, these patterns could be the result of indirect evolutionary relationships rather than the outcome of a direct interaction with Hsp90A. We therefore aimed to predict specific genes that encode putative hsp90A clients. Our method is based on the assumption that strong, conserved clients should be heavily dependent on Hsp90A, and thus should be found only rarely in the absence of hsp90A throughout evolution. To estimate the expected frequency of each hsp90A-associated gene with and without hsp90A, we used the inferred BayesTraits rates to calculate the steady-state probabilities of each of the 4 possible two-gene presence/absence states (Methods). These probabilities represent the proportion of the time that some arbitrary bacterial lineage will spend in each of the presence/absence states throughout evolution. From these probabilities we calculated a Putative Client Index (PCI) for each hsp90A-associated gene to evaluate how often it was present without hsp90A throughout evolution, compared to a null expectation (see Methods). This index is close to zero for genes Arrow widths in each diagram are scaled to represent the rate of each transition. The average transition rate and standard deviation across multiple BayesTraits runs are displayed (see Methods). Box plots of the rates of gain and loss of all hsp90A-associated flagellar genes (n = 27; C) and all hsp90Aassociated Type III secretion genes (n = 10; D) further demonstrate consistent co-evolutionary dynamics of genes in these categories. A box plot of all hsp90A-associated secretion genes (including all types) is provided as Figure S4. doi:10.1371/journal.pgen.1003631.g002 that were infrequently present without hsp90A and were hence likely to be Hsp90A clients. We defined the genes with the lowest PCI values as putative clients (Table 2; see also Text S1).

Novel and known functions of putative Hsp90 clients
Consistent with our prior analysis, several flagellar genes behaved as potential clients (Table 2). In particular, our set of putative clients included several genes (fliH, fliI, fliN) whose products had been previously shown to physically interact with Hsp90A in E. coli [19]. The products of these genes are cytoplasmic components of the flagellar rotor and export apparatuses. In contrast, nonflagellar type III secretion genes were all absent from the list of potential clients. In fact, nonflagellar type III secretion system components were rated as some of the least likely clients by our index (Figure 3). This disparity in predicted client status mirrors the different evolutionary relationships of these complexes with hsp90A ( Figure 2). Chaperone/proteases (e.g. ClpA and PpiD) also ranked high in our list of potential clients. Hsp90A is known to collaborate with other chaperone systems such as DnaK [14,32] but to date no obligate co-chaperones have been described. The identified chaperone/proteases may represent such co-chaperones or collaborating chaperone systems, since our index cannot discriminate between Hsp90 clients and Hsp90 cochaperones (or other collaborating proteins). Alternatively, these observed associations could simply indicate that components of the cytoplasmic stress response are dependent upon Hsp90A.
We also found several unexpected putative clients, such as the 3hydroxybutyryl-CoA dehydrogenase PaaH and the transcription termination factor Rho, which we predict to be the two strongest clients. Further study will be necessary to understand these associations and the underlying cause of the co-evolutionary association between these genes and hsp90A.

Swimming motility and chemotaxis assays of Hsp90Adefective E. coli
Our putative clients and the predicted chaperone role of Hsp90A in flagellar assembly are consistent with previous observations. Specifically, the deletion of E. coli hsp90A, also known as htpG, resulted in reduced surface swarming movement [33]. We also previously observed physical interactions between the HtpG protein and certain flagellar proteins [19]. Yet, these observations lacked a clear demonstration of client status or mechanism, and E. coli swarming is a complex behavior that depends on numerous factors in addition to flagellar function [34]. We therefore set out to test our hypothesis that Hsp90A is physiologically important for flagellar assembly and function and that flagellar components are indeed Hsp90A clients.
We examined the swimming motility phenotype of DhtpG E. coli strains on soft-agar plates (Methods). In contrast to surface swarming, swimming is a less complex behavior, in which bacteria use functional flagella and chemotaxis components to swim from an inoculation point through agar pores, following nutrient gradients that are created by nutrient depletion within the colony. The soft-agar assay is routinely used to assay bacterial swimming motility and chemotaxis. To enhance our ability to detect differences between wild-type and DhtpG cells, the assays were performed competitively. Competitive assays emphasize small differences between strains and reduce experimental error, thereby increasing the sensitivity of the assay. After mixing equal amounts of YFP-labeled WT and CFP-labeled DhtpG strains, this mixture was inoculated in the center of a soft-agar plate and incubated at 34uC for 8 hrs. We then counted cells of each strain in the plate center vs. the outer edge using fluorescence microscopy ( Figure 4A). DhtpG mutants migrated less efficiently to the plate's outer edge relative to WT, confirming that they are partially deficient in their motility and/or chemotaxis ( Figure 4B). This defect is apparently subtle, since little difference between WT and DhtpG cells was observed in a non-competitive assay ( Figure S5), but it could be revealed due to strong selection for cells with optimal motility and chemotaxis at the outer edge of the spreading bacterial population. We also tested the phenotype of the HtpG(E34A) mutant, which has reduced rates of ATP hydrolysis and is deficient in substrate refolding [14,35]. Since HtpG ATPase activity is necessary for release of clients, HtpG(E34A) is less efficient at releasing clients [36][37][38]. Indeed, this mutant showed stronger motility/chemotaxis defects than the DhtpG strain ( Figure S5), presumably due to sequestration of its client proteins. We therefore employed the HtpG(E34A) mutant in all subsequent assays as a more sensitive test of HtpG involvement. Taken together, our observations suggest that the motility defect may be due to the improper function or sequestration of HtpG clients.

FRET observation of HtpG interactions with flagellar motor components
To further investigate the in vivo interaction of HtpG with flagellar components, we used htpG-yfp and htpG(E34A)-yfp constructs expressed in WT cells to perform acceptor photobleaching FRET between HtpG and FliN-CFP over an E. coli growth curve. Motility of E. coli is known to increase at the transition from the early exponential to post-exponential phase of growth [39], and this experimental design enabled us to examine the HtpG-FliN interaction in the context of the flagellar assembly process. If HtpG is indeed involved in the assembly process of these structures, the interaction of HtpG with FliN should correspond temporally to the timing of flagellar assembly. Indeed, we found that the interaction with FliN peaked at OD 600 = 0.2 ( Figure 5A) and correlated well with the onset of cell motility in wild-type cells ( Figure 5B). Moreover, the interaction of HtpG(E34A) with FliN was stronger and delayed compared to the binding of wild-type HtpG. Correspondingly, the onset of motility was delayed in cells expressing HtpG(E34A) ( Figure 5B). This is consistent with the delayed release of clients by HtpG(E34A), suggesting that HtpG's role in motility derives from a direct involvement in flagellar complex assembly.
Given that both bacterial and eukaryotic Hsp90s are known to collaborate with Hsp70 in refolding proteins [14,[40][41][42], we considered the possibility that this was also the case for bacterial flagellar assembly. We previously showed that some flagellar motor components interact with DnaK, the E. coli Hsp70 homolog [19]. Therefore, we repeated the FRET experiments testing for interactions between HtpG or HtpG(E34A) and FliN in a DcbpADdnaJ background. CbpA and DnaJ are DnaK co-chaperones and are essential for DnaK-dependent refolding activity [14]. DnaK should not be able to pass substrates to HtpG in this mutant background. Indeed, we found that FRET interactions with FliN disappear for both HtpG proteins in this background ( Figure S6A), suggesting that DnaK-dependent remodeling precedes HtpG action in flagellar complex assembly.

FRET observation of HtpG interactions with chemoreceptor components
Since a recent high-throughput assay showed kinases to be overrepresented among eukaryotic Hsp90 clients [43,44], we next examined whether the HtpG-dependent defects in chemotaxis may also be due to defective chemoreceptor kinase activity. Although no chemotaxis proteins were found in our list of the strongest putative clients, we did observe a significant enrichment of these components in the hsp90A-associated set (Table 1). We thus tested interactions between six chemoreceptor cluster components and HtpG(E34A) using, as before, acceptor photobleaching FRET (Table S4). We observed a strong interaction of HtpG(E34A) with the chemoreceptor kinase CheA. Our results suggest that the FliN/HtpG and CheA/HtpG interactions are direct and do not depend on other flagellar or chemotaxis proteins, since these interactions are robust to deletion of flhC, which ablates expression of all endogenous flagellar and chemotaxis genes (Table S4) [19]. Moreover, the CheA dimerization domain was required for association with HtpG, supporting the hypothesis that HtpG aids oligomerization of its clients [17,45]. Testing HtpG interactions with other chemotaxis proteins of E. coli revealed an additional strong interaction with the dimeric phosphatase CheZ but not with other proteins (Table S4).
We again examined the temporal dynamics of these interactions. Due to the hierarchical order of flagellar and chemotaxis gene expression [39,46], the assembly of chemoreceptor clusters is delayed compared to the assembly of flagellar motors as nonmotile cells transition into motile cells. Indeed, the interaction of HtpG with CheA peaked at OD 600 = 0.3, after the FliN peak ( Figure 5A). Just as for FliN, the interaction of HtpG(E34A) with CheA was stronger and delayed compared to wild-type HtpG, and the HtpG-CheA interaction disappeared in a DcbpADdnaJ background ( Figure S6B). Collectively, these findings suggest that HtpG plays an important role in the assembly of both the flagellar motor and chemoreceptor clusters through separate client interactions.

Association of hsp90A with life history traits in bacteria
Given the role of HtpG in chaperoning proteins that mediate interactions with the environment, and the known role of eukaryotic Hsp90 in phenotypic robustness, we finally examined whether hsp90A directly co-evolved with certain bacterial organismal traits. We considered several organismal traits, including aerobism, thermophilicity, halophilicity, the ability to form endospores, pathogenicity, motility, and habitat preferences (see Methods). We used BayesTraits and the Ciccarelli tree to identify traits that co-evolve with hsp90A. Out of the 11 analyzed traits, 4 exhibited significant associations with hsp90A (p,0.05; Table S5), with the strongest association observed between hsp90A and the  capacity to inhabit multiple habitats. Moreover, examining the gain and loss rates obtained, we found that hsp90A is gained and lost at significantly higher rates in organisms that inhabit multiple habitats (with no gains inferred in single habitat organisms), suggesting that a preference for multiple habitats imposes a different selection regime on hsp90A ( Figure 6). We also tested whether the co-evolutionary dependency between hsp90A and multiple-habitat preferences was unidirectional, as we observed for some hsp90A-associated genes. Comparing the four co-evolutionary models described above and applying AIC to identify the bestfitting model, we found that hsp90A gain and loss depended on habitat preference, but not vice versa. This observation suggests that in organisms inhabiting multiple environments hsp90A is subjected to dynamically shifting selective pressures, potentially alternating between selection for and against hsp90A.

Discussion
We set out to discover Hsp90 functions conserved throughout the bacterial tree of life. We found that hsp90A, the most common paralog of bacterial Hsp90, bore strong signatures of co-evolution with several hundred genes and with specific life history traits, shedding light on its function and impact on evolutionary history. Most notably, we found that hsp90A co-evolved with membrane protein complexes such as flagella and other Type III secretion (T3S) systems. Our results suggest that Hsp90's role in sensing and responding to environmental stimuli is conserved between bacteria and eukaryotes.
Similar to verified eukaryotic Hsp90 clients [5], our predicted putative Hsp90A clients were a diverse group of proteins (e.g. the flagella protein FliN, the chaperone ClpA, and the ribosomal protein RluB; see Table 2) that tended to belong to specific functional categories (e.g. flagellar proteins, chaperones, and ribosomal components). As our methods can only infer associations between genes that are frequently gained and lost, we may substantially underestimate the number of hsp90A-associated genes and clients. However, the non-essentiality and frequent loss of hsp90A throughout bacterial diversity argues that genes not captured in our analysis (since they are not frequently gained Error bars indicate standard errors from three replicates. For these assays, a truncated form of CheA lacking the first 97 amino acids (CheA s ) was used because this fusion was more stable against spontaneous proteolysis than the fusion to full-length CheA, but showed similar interaction with HtpG (Table S4)   and lost) are unlikely to be strongly dependent on the chaperone throughout bacteria.
The subtlety of the bacterial Hsp90 mutant phenotypes that we (and others) report implies that Hsp90's role in cellular physiology has diverged between eukaryotes and prokaryotes [17,45,47]. In other words, either essential pieces of cellular physiology changed, or Hsp90 function changed. We favor the first hypothesis, because Hsp90 is well-conserved among bacteria, archaea, and humans at the sequence level [13], and retains a similar quaternary structure [48] and biochemical activity [15,37,44]. In contrast, bacterial and archaeal cells differ significantly from eukaryotic cells. Eukaryotic cells have higher cell compartmentalization, longer and multifunctional proteins with multiple domains [49], and increased protein interactome complexity [50]. Together with the existence of many eukaryotic Hsp90 co-chaperones, all these features may contribute to the greater essentiality of Hsp90 in eukaryotes.
The dependence of HtpG-client interactions upon the DnaK chaperone system, as observed by us and by others [14,15], argues that Hsp90A is well-integrated with other chaperone systems. Our putative clients included ClpA, the substrate adaptor for the ClpAP/ClpAXP chaperone/protease complexes, and PpiD, a periplasmic chaperone [51]. Like HtpG, PpiD is necessary for optimal swarming motility [33], suggesting that it may participate in flagellar assembly. We speculate that these proteins act as Hsp90A co-chaperones in some bacteria; alternatively, their dependence on Hsp90A may represent an example of collaborating chaperone systems.
The best-characterized Hsp90 client in bacteria is the structural ribosomal protein L2 [15,18], which is near-universally conserved throughout life (and hence not detectable by our method). In addition to L2, other ribosomal proteins were found to interact with HtpG in large-scale proteomics analyses. In agreement with these observations, we found the ribosomal proteins RlmE and RluB among the predicted hsp90A clients.
Although these chaperone and ribosomal proteins were predicted to be stronger clients than flagellar proteins, our experimental validation focused on the latter as their client status was suggested by previous observations [19,33]. We present four lines of evidence for HtpG client status for the flagellar protein FliN and the chemoreceptor kinase CheA, including direct interactions with HtpG, physiologically relevant timing of HtpG-FliN/CheA interactions, phenotypic consequences of reduced HtpG function in CheA/FliN-dependent traits, and dependence of CheA/FliN interactions with HtpG upon the Hsp40-Hsp70 pathway. The identification of FliN and CheA as HtpG clients is consistent with the hypothesis that bacterial Hsp90 facilitates the assembly of large membrane-associated protein complexes [17,45].
Curiously, whereas the flagellar T3S system contained Hsp90A clients, the nonflagellar T3S system is predicted to have an antagonistic relationship with Hsp90A. Nonflagellar T3S systems and the flagellar T3S systems are closely related (NF-T3SS and F-T3SS) [52,53]. 9 NF-T3SS components are directly homologous to flagellar components, of which 8 were found to co-evolve with hsp90A in our analysis. Yet, these 8 genes are predicted to coevolve antagonistically with hsp90A (Figure 3), whereas their flagellar homologs are mostly predicted to be clients (for instance, the fliI and yscN genes shown in Figure 2 are homologous). This result suggests that some relationship with Hsp90A is conserved between the two T3S systems, but with apparently opposite effects in each system. This result may reflect the fact that each of these systems is an adaptation to different ecological challenges. Specifically, we have shown that Hsp90A is important for flagellum-enabled motility and chemotaxis in E. coli. This mode of motility is strongly adaptive in certain physical environments [34,54,55], and thus Hsp90A is likely to be associated with fitness in these environments through flagellar assembly. The presence of NF-T3SS is likewise an adaptation to certain biotic environments [55,56]. Our observation that organisms inhabiting multiple habitats experience fluctuating selection for hsp90A is also consistent with competing selection pressures. Representative genes of these homologous T3S families were not significantly associated with habitat preferences, arguing that hsp90A's association with habitat preferences is not a byproduct of associations with T3S systems. Nonetheless, we suggest that these two T3S systems constitute a link between Hsp90A and phenotypic robustness across different environments.
Inferring function from evolutionary associations has some caveats. For instance, F-T3S systems can be found in genomes that lack hsp90A. If F-T3S systems include Hsp90A clients, then what may render Hsp90A-dependent stabilization dispensable in some bacteria? Experimental validation will be necessary to answer such questions, and to distinguish true client relationships from indirect co-evolutionary associations. As discussed before, our method is subject to gene set bias, in that only genes that are gained and/or lost frequently will have enough statistical power to reject the null hypothesis. Similarly, as our method assumes that relationships are maintained throughout the analyzed phylogeny, we cannot reliably detect genes that are associated with hsp90A in some organisms but not in others.
Although much work remains to articulate the precise mechanistic relationships between hsp90A and its co-evolving genes, our results highlight the tremendous potential of evolutionary inference for guiding experimental research. More generally, our study provides a successful example of how evolutionary perspectives and phylogenetic analyses can inform and advance the study of complex biological systems and the inference of elusive biological functions.

Prokaryotic Hsp90 paralogs
We downloaded all Hsp90 amino acid sequences (including all paralogs) for bacteria with full KEGG genome annotations from the KEGG database [24,57]. We aligned these sequences using ClustalO [58], and used the PHYLIP package [59] to construct neighbor-joining trees and assess their phylogenetic support through bootstrapping. We assigned Hsp90 families to branches according to bootstrap support for the branch and previous classifications [11,12].

Genome data
We acquired presence/absence patterns of genes across organisms from the KEGG database release 60.0 (in the form of KEGG Orthology/KO profiles) [57], and functional annotations from KEGG Class. Genes that were either present in fewer than five species or absent in fewer than five species in the tree of interest were dropped from our analysis, as these genes are unlikely to show meaningful signatures of co-evolution by this method.

Phylogenetic trees
We obtained the tree constructed by Ciccarelli et al. (Ciccarelli tree) [29] and pruned it to 148 bacterial species for which KEGG genome data was available. We also obtained the LTP104 version of the 16S/23S rRNA tree from the All-Species Living Tree Project (Yarza tree) [30,60]. We used ARB [61] to prune this tree to bacterial species for which KEGG genome data was available.
We further pruned this tree to omit clades placed paraphyletically at the taxonomic levels of phylum, class, order, and family. This filtered tree included 797 bacterial species. As BayesTraits cannot process trees with zero-length branches, all branch lengths equal to zero were replaced with a negligible branch length (0.00001, approximately an order of magnitude smaller than the next smallest branch length in each tree).

Organismal trait data
We acquired organismal trait data from the NCBI Entrez genome project, November 2011 [62]. We recoded all traits into presence/absence patterns for the trait in question. For instance, an organism found to be pathogenic towards any other organism was coded as '1' for the trait of pathogenicity, whereas an annotated organism that was never found to be pathogenic was coded as '0'. Similarly, we coded both thermophilic and hyperthermophilic organisms as '1' for the trait of thermophilicity, whereas all other annotated organisms were coded as '0'; anaerobic organisms were coded as '0' for the trait of aerobicity, whereas all other annotated organisms were recoded as '1'. We define as inhabiting multiple habitats any organism that inhabits more than one of NCBI's habitat categories. For BayesTraits analysis, the tree was pruned to include only species annotated for the trait in question (each trait analysis was accordingly performed on a slightly different set of species; see Table S5 for details on species number for each analysis).

Detecting evolutionary associations with BayesTraits
A complete description of the BayesTraits (v1.0) framework can be found elsewhere [26]. Briefly, consider a character with 2 states, 0 and 1. If a species has 2 such distinct characters, it can occupy 4 possible states: 1:(0,0), 2:(0,1), 3:(1,0), and 4: (1,1). Specifically, if these 2 characters represent the presence or absence of two genes, hsp90A and gene X, these four states correspond to (hsp90A2, X2), (hsp90A+, X2), (hsp90A2, X+), and (hsp90A+, X+). Evolution is then the process by which these genes are gained and lost over time. Consider accordingly an evolutionary process where only one character can change state at a time. Such a process can then be described by 8 parameters for the rates of transition per unit time between these 4 states: Q = [q12, q13, q21, q31, q24, q34, q42, q43], where qxy is the rate of transition from state x to state y. BayesTraits implements this model of evolution as a continuous-time Markov process and estimates each of these rate parameters by maximum-likelihood (ML). We further validated that these ML-based rates are consistent with reversible-jump Markov chain Monte Carloderived estimates (Methods; Text S1). This estimation is based on a phylogeny and on the states of the two characters at the tips of the phylogeny. Having estimated these rates, BayesTraits additionally calculates the likelihood of the model based on the character states at the tips of the phylogeny.
We can further compare different models of evolution by forcing certain parameters to be equal. We specifically considered the following 4 models: 1) hsp90A and X are independent (Q: q12 = q34, q21 = q43, q13 = q24, q31 = q42; 4 parameters total) 2) hsp90A and X are mutually dependent (No parameter restrictions; 8 parameters total) 3) X depends on hsp90A but not vice versa (Q: q12 = q34, q21 = q43; 6 parameters total) 4) hsp90A depends on X but not vice versa (Q: q13 = q24, q31 = q42; 6 parameters total) Identifying hsp90A-associated genes We used discrete from the BayesTraits package [26][27][28] to infer associations between hsp90A and other bacterial genes and between hsp90A and various organismal traits. We first tested for an evolutionary association with hsp90A by comparing model 1 to model 2 above with a likelihood ratio test (LRT), as previously described [28]. In our likelihood-ratio tests, the 2Log(LR) approximates a x 2 test statistic for rejecting the independent model as a null hypothesis, and is calculated as twice the difference of the log-likelihoods of a co-evolutionary model and a model of evolutionary independence. The set of genes for which model 2 is preferred (i.e., model 1 is rejected as a null hypothesis) have an evolutionary association with hsp90A. Since different runs of the BayesTraits maximum likelihood method can potentially produce different parameter values, we repeated this procedure 100 times, each potentially resulting in a different gene set. We validated that these sets are similar and the choice of gene set does not substantially affect downstream analysis (Text S1). Any gene that was found to be associated with hsp90A in at least 90 runs was defined as hsp90A-associated gene. See Text S1 for more details.

Reversible-jump Markov chain Monte Carlo analysis
We selected 10 genes at random from the hsp90A-associated set and used the BayesTraits implementation of reversible-jump Markov chain Monte Carlo to estimate the rate parameters for their gain and loss in concert with hsp90A [63]. For each of these 10 genes, we used an exponential rate prior with mean and variance equal to 30, and ran the chain for 150 million iterations while sampling every 100 iterations. We discarded the first 75 million iterations as burn-in and used the remaining iterations as a posterior distribution of rate parameter estimates. We used Tracer v1.5 [64] and previously described criteria to evaluate chain convergence in this remaining sample [65]. For each rate, we used the median of its posterior distribution in this sample as a point estimate.

Co-evolutionary model selection
To provide an accurate description of the co-evolutionary dynamics of hsp90A-associated genes, we further applied Bayes-Traits to these genes, estimating the likelihood of each of the four models described above. We identified the best fit model for each gene using the Akaike Information Criterion (AIC) [31], taking into account both the likelihood score and the number of parameters in each model. We again repeated this procedure 100 times and classified a gene into a specific co-evolutionary model only if it fit this same model in at least 90 runs (see Text S1 for more details). This two stage scheme, first identifying associated genes and then selecting a model that best describes their evolutionary relationship with hsp90A, provides a more stringent test of co-evolution and supports a simple approach for multiple testing correction.

Prediction of Hsp90A clients in bacteria
We used BayesTraits-derived evolutionary transition rates under the fully unrestricted model to estimate residence times in specific states (for instance, the proportion of time spent by bacteria in a state where both hsp90A and some other gene are present, vs. the time when only the other gene is present) under steady state dynamics. For a given gene, the probability of being in one of the four states, A: (hsp90A absent, Gene absent), B: (hsp90A present, Gene absent), C: (hsp90A absent, Gene present), D: (hsp90A present, Gene present) at a very small increment of time Dt after time t is given by: We can differentiate this to obtain the instantaneous change in each probability: At steady state dA/dt = 0, dB/dt = 0, etc., and therefore: This set of linear equations can be solved for A, B, C, and D, with the requirement that A+B+C+D = 1. We replaced 0 rates with the smallest nonzero rate in the model multiplied by 0.001 to allow transitions between all states. The positive nonzero solution for A, B, C, and D can then be conceived as the expected residence times along some arbitrary bacterial lineage. We used these residence times to estimate a Putative Client Index, PCI, denoting the normalized residence time in state C: Notably, if Hsp90A and the gene's product have no client relationship, the proportion of time spent in state C is expected to be equal to (C+D)(A+C), so a PCI close to 1 indicates that the observation does not differ from the expectation. Smaller values of PCI therefore indicate that a gene is observed less frequently than expected without hsp90A, and is thus more likely to be a client. Since no obvious threshold value can be defined, we considered the 20 genes with the lowest PCI values as putative clients ( Figure 3 and Table 2; Methods). To account for variation in rates between BayesTraits runs we repeated this procedure 100 times and defined as putative clients those that were identified as clients in at least 90 of these runs (see Text S1). PCI scores shown in Table 2 and Figure 3 are averages across all runs.

Functional enrichment analysis
We used a hypergeometric test to assess whether each KEGG Class functional annotation is overrepresented in the various Hsp90-associated gene sets. As a background set in each case we used the entire set of genes analyzed. Any annotation present in less than 4 copies in the background set was not considered. We accepted enrichments at a 5% FDR.
TB soft agar plates were prepared by supplementing TB with 0.3% agar (Applichem) and when necessary with 100 g/mL ampicillin and 1 mM IPTG. Equal amounts of cells from different overnight cultures, adjusted depending on their optical density to the equivalent of 2.5 ml of culture with OD 600 of 2.0, were inoculated and allowed to spread at indicated temperatures for indicated times. Following incubation, photographs of plates were taken with a Canon EOS 300D (DS6041) camera. Images were analyzed with ImageJ (Wayne Rasband, NIH, http://rsb.info.nih.gov/ij/) to determine the diameter of the rings of spreading colonies. For analysis of motility at different growth stages (indicated by OD 600 value), percentages of motile cells were estimated from the microscopy movies of swimming cells. The experiment was performed with the RP437 strain, which is non-motile above 37uC. Cells were grown overnight in TB medium at 37uC to completely inhibit their motility. After dilution in fresh TB medium to OD 600 0.01, cells were grown at 34uC for measurements.

Fluorescence imaging
For microscopy, cells were taken from the soft-agar plates and applied to a thin agarose pad (1% agarose in tethering buffer). Fluorescence imaging was performed on a Zeiss AxioImager microscope equipped with an ORCA AG CCD camera (Hamamatsu), a 1006NA 1.45 objective, and HE YFP (Excitation BP 500/25; Emission BP 535/30) and HE CFP (Excitation BP 436/25; Emission BP 480/40) filter sets. Each imaging experiment was performed in duplicate on independent cultures. All images were acquired under identical conditions. Images were subsequently analysed using ImageJ software.

Acceptor photobleaching FRET measurement
FRET measurements by acceptor photobleaching were performed on a custom-modified Zeiss Axiovert 200 microscope as described before [66]. Briefly, cells expressing YFP and CFP fusions of interest were concentrated about tenfold by centrifugation, resuspended in tethering buffer and applied to a thin agarose pad (1% agarose in tethering buffer). Excitation light from a 75 XBO lamp, attenuated by a ND60 (0.2) neutral-density filter, passed through a band-pass (BP) 436/20 filter and a 495DCSP dichroic mirror and was reflected on the specimen by a Z440/532 dual-band beamsplitter (transmission 465-500 and 550-640 nm; reflection 425-445 and 532 nm). Bleaching of YFP was accomplished by a 20 sec illumination with a 532 nm diode laser (Rapp OptoElectronic), reflected by the 495DCSP dichroic mirror into the light path. Emission from the field of view, which was narrowed with a diaphragm to the area bleached by the laser, passed through a BP 485/40 filter onto a H7421-40 photon counter (Hamamatsu). For each measurement point, photons were counted over 0.5 s using a counter function of the PCI-6034E board, controlled by a custom-written LabView 7.1 program (both from National Instruments). CFP emission was recorded before and after bleaching of YFP, and FRET was calculated as the CFP signal increase divided by the total signal after bleaching. DflhC strains were used to define direct interactions between HtpG and flagellar and chemotaxis components. In this background expression of endogenous flagellar and chemotaxis genes is inhibited, thus eliminating indirect interactions that may result from concomitant binding of HtpG and tested protein to a third flagellar or chemotaxis protein.

Supporting Information
Dataset S1 hsp90A-associated genes and relevant information. (XLS) Figure S1 Phylogenetic clustering of bacterial hsp90 paralogs. (A) Neighbor-joining phylogeny of 897 bacterial Hsp90 amino acid sequences. Groups Hsp90A, Hsp90B, and Hsp90C as defined by Chen et al. [11] are illustrated. (B) Consensus neighbor-joining tree for 100 bootstraps with clades collapsed to highlight deep branch structure. Bootstrap support for each branch is displayed and is also reflected by the branch lengths. One species (ZIN, representing Hsp90 from the organism Candidatus Zinderia insecticola CARI), never grouped within the other divisions shown, and was excluded from our analysis. The branch separating Hsp90B and Hsp90C from the Hsp90A clades is present in 99/100 bootstrap trees. (C) Hsp90A, B, and C presence/absence patterns mapped onto a 16S/23S rRNA phylogeny of 797 bacterial species [30] (see Text S1). Branch lengths are ignored for ease of display. (TIF) Figure S2 Co-evolutionary gain and loss rates of all hsp90Aassociated flagellar genes. The layout of each diagram is similar to that used in Figure 2. (TIF) Figure S3 Co-evolutionary gain and loss rates of all hsp90Aassociated secretion genes. The layout of each diagram is similar to that used in Figure 2. (TIF) Figure S4 Box plots of the rates of gain and loss of all hsp90Aassociated secretion genes (n = 16). See also Figure 2D     Text S1 Additional details on Hsp90 paralog distribution, consistency of BayesTraits runs, and robustness of co-evolutionary associations to choice of phylogeny. (DOC)


text:
Characterization of Enterococcus faecalis Phage IME-EF1 and Its Endolysin

Enterococcus faecalis is increasingly becoming an important nosocomial infection opportunistic pathogen. E. faecalis can easily obtain drug resistance, making it difficult to be controlled in clinical settings. Using bacteriophage as an alternative treatment to drug-resistant bacteria has been revitalized recently, especially for fighting drug-resistant bacteria. In this research, an E. faecalis bacteriophage named IME-EF1 was isolated from hospital sewage. Whole genomic sequence analysis demonstrated that the isolated IME-EF1 belong to the Siphoviridae family, and has a linear double-stranded DNA genome consisting of 57,081 nucleotides. The IME-EF1 genome has a 40.04% G+C content and contains 98 putative coding sequences. In addition, IME-EF1 has an isometric head with a width of 35 nm to 60 nm and length of 75 nm to 90 nm, as well as morphology resembling a tadpole. IME-EF1 can adsorb to its host cells within 9 min, with an absorbance rate more than 99% and a latent period time of 25 min. The endolysin of IME-EF1 contains a CHAP domain in its N-terminal and has a wider bactericidal spectrum than its parental bacteriophage, including 2 strains of vancomycin-resistant E. faecalis. When administrated intraperitoneally, one dose of IME-EF1 or its endolysin can reduce bacterial count in the blood and protected the mice from a lethal challenge of E. faecalis, with a survival rate of 60% or 80%, respectively. Although bacteriophage could rescue mice from bacterial challenge, to the best of our knowledge, this study further supports the potential function of bacteriophage in dealing with E. faecalis infection in vivo. The results also indicated that the newly isolated bacteriophage IME-EF1 enriched the arsenal library of lytic E. faecalis bacteriophages and presented another choice for phage therapy in the future.

Introduction
Enterococcus faecalis and Enterococcus facium belong to Gram-positive bacteria that commonly inhabit the lower intestinal tract, oral cavity, and vaginal tract of humans or animals. E. faecalis and E. facium are opportunistic pathogens. In a healthy individual, the two bacteria have no adverse effect on the host; however, the bacterium can cause life-threatening infection in immune-compromised patients [1][2][3]. Although E. faecalis is used as a starter for fermented food and probiotics [4], the bacterium can cause endocarditis and bacteremia, urinary tract infections, meningitis, and other infections in humans. In addition, Enterococcus has been found to be associated with diabetic foot infection and root canal treatment failure [5][6][7]. Several virulence factors have been thought to contribute to E. faecalis infections [3,8]. A plasmid-encoded hemolysin, called the cytolysin, is important for pathogenesis in animal models of infection, whereas cytolysin combined with high-level gentamicin resistance is associated with a five-fold increase in risk of death in human bacteremia patients [9][10][11]. Enterococci survive in very harsh environments including extreme alkaline pH (9.6) and salt concentrations. Enterococci resist bile salt detergents, heavy metals, ethanol, azide, and desiccation, and can grow at 10 °C to 45 °C and survive at 60 °C for 30 min [7]. Nosocomial infection caused by E. faecalis or E. facium is becoming a major concern in hospital settings because acquisition of virulence or antibiotic resistance makes this bacterium more difficult to be controlled.
To date, even the so called "last defense line" of vancomycin has no inhibiting effect on some vancomycin-resistant strains of E. faecalis and E. facium [12,13]. Thus, finding an alternative agent or a strategy for treating antibiotic-resistant E. faecalis or E. facium infection has become increasingly indispensable [14].
Lytic bacteriophage or its endolysin has been reconsidered as a potential agent for treating clinical infection under the condition of multi-drug resistant bacteria emerging recently [15]. However, the characteristics of species-specific infection, or the narrow host spectrum of bacteriophage and inadequate methodologies in applying bacteriophage in animal assay or clinical research, greatly limit its further application in therapy or decontamination [16]. Thus, isolating new bacteriophage to enrich its therapeutic arsenal and assess its function in animal model will compensate the above shortages or limitations, and accelerate the application of bacteriophage in future therapy or decontamination. In this study, we isolated and characterized a lytic E. faecalis bacteriophage named IME-EF1 and its endolysin. The results indicated that bacteriophage IME-EF1 or its endolysin has potential in treating E. faecalis infection or contamination.

Phage isolation and characterization
When incubated with E. faecalis 002 (Table 1), the treated sewage sample resulted in the formation of a clear plaque. Picking up the plaque and inoculating it into the E. faecalis 002 culture at the mid-exponential phase led to the lysis of the bacteria. The isolated bacteriophage was named IME-EF1. After purification and concentration, the electron microscopy result showed that phage IME-EF1 has an icosahedral head and a non-contractile tail. The width of the head is about 35 nm to 60 nm, the length is about 75 nm to 90 nm, and the tail is about 130 nm to 220 nm ( Figure 1A). Thus, the IME-EF1 was classified into the Siphoviridae family morphology. The optimal multiplicity of infection (M.O.I.) of the phage IME-EF1 was determined, as shown in Table 2. The results indicated that the same amount of phage as that of the bacteria (namely M.O.I=1) led to the highest production of progeny phage when inoculating; thus, the optimal M.O.I. of phage IME-EF1 is 1. The one-step growth kinetics is shown in Figure 1B. The results showed that the phage IME-EF1 had a latent period time of 25 min and a burst size of 60 PFU/infected bacteria when infecting its host bacteria E. faecalis 002.

Whole genome sequencing and annotation
High-throughput sequencing of the phage genomic DNA generated 311,503 valid reads, with which the complete sequence of the genome was assembled using both Velvet and CLC Genomic Workbench, with an average coverage of 75× [10]. The complete genome of phage IME-EF1 consists of 57,081 bp, with an average GC content of 40.04% (GenBank Accession number: KF192053). Running BLASTN with whole genome sequence showed that the isolated IME-EF1 was highly homologous to Enterococcus phage SAP6 (GenBank Accession Number JF731128), with a coverage of 87%, and Enterococcus phage BC-611 (DDBJ Accession Number AB712291), with a coverage of 81%, both of which also belong to Siphoviridae family and were isolated in South Korea and Japan, respectively. Numerous studies on several other E. faecalis bacteriophage have been reported, such as EFAP-1, phiFL1-3, EFRM31, EfaCPT1 (Beheshti et al., GenBank Accession Number JX193904), EF24c, EF11, SAP6, and BC-611, which have been thoroughly genetically characterized [17][18][19][20][21][22][23]. However, IME-EF1 has high similarity to SAP6 and BC-611 ( Figure 2A).
The IME-EF1 genome contains 98 putative CDSs, but no tRNA was found in the IME-EF1 genome using online program tRNAscan-SE 1.21 [24,25]. The culture results in the bacterium and the absence of tRNA or integrase gene in the IME-EF1 genome suggest that IME-EF1 is a lytic bacteriophage. A total of 71 CDS of IME-EF1 encode proteins of unknown function. The other 27 CDS were predicted to encode morphogenesisrelated proteins (major tail protein, major capsid/head protein, phage head morphogenesis protein, phage portal protein, phage terminase large subunit, phage terminase small subunit, phage minor structure protein, tail fiber protein, minor capsid protein, and tail protein), replication-associated proteins, and DNA manipulation proteins (methyltransferase, phage DNA binding protein, cytidine deaminase, metalloprotease, DNA polymerase, LPS glycosyltransferase, adenylate kinase, endodeoxyribonuclease, HNH homing endonuclease/phage intron, replicative DNA helicase, DNA replication protein, transcription regulator, and DNA primase) ( Figure 2B).

Endolysin identification
Annotation found a 714 nt length sequence in the negative strand of the IME-EF1 genome, encoding a putative 237 aa endolysin protein. The conserved domain prediction indicated that the putative endolysin is an amidase containing a conserved CHAP domain in its N-terminal end and shares highly similar amino acid sequence to the corresponding putative lysins in SAP6 and BC-611, both of which also encode a CHAP domain endolysin. When phylogenized with seven putative or confirmed endolysin from other E. faecalis phage, the putative endolysins from IME-EF1, SAP6, and BC-611 are located in the same sub-clade and the other containing similar domains in N-terminal sequences were also clustered ( Figure  3A). The putative endolysin gene of IME-EF1 was cloned and fused with His tag sequence into vector pColdI. Through the induction of 1 mM isopropyl beta-D-thiogalactoside (IPTG), the transformed E. coli M15 expressed 30 kDa soluble protein, and the band can be probed using anti-His Tag antibody ( Figure  3B).

Bactericidal spectrum determination and Endolysin binding assay
Among the 10 clinical strains of E. faecalis, 10 clinical strains of E. faecium, and 4 clinical strains of S. aureus, only 3 strains of E. faecalis and 1 strain of E. faecium can be lysed by IME-EF1. However, all strains of E. faecalis tested in this assay and the same strain of E. faecium could be lysed by purified IME-EF1 endolysin ( Table 3). The results further demonstrated that endolysin has wider bactericidal spectrum than its original phage. To further evaluate the bactericidal activity of IME-EF1 Table 1. Bacterial strains employed in this study and their antibiotic-resistant spectrum. Characterization of Phage IME-EF1 PLOS ONE | www.plosone.org endolysin, bacterial culture in the mid-exponential phase was collected, suspended in PBS, and then mixed with purified endolysin. The descent of OD600 was measured and used to assess the bactericidal ability. The results showed that all the strains of E. faecalis used in this assay could be lysed by the purified endolysin, and resulted in the decreasing in the OD600 value mainly within 1.5 h after endolysin addition ( Figure 4A). Endolysin binding assay was conducted in this study to evaluate whether binding to bacteria is required for IME-EF1 endolysin functioning during lysis. When incubated with purified endolysin, the susceptible bacterial strains can bind to endolysin. After pelleting the bacteria, the adsorbed endolysin can be pulled down accordingly, resulting in endolysin reduction in the supernatant. The bacterial strains that cannot be lysed by the endolysin, such as E. faecium 81 and S. aureus, cannot significantly reduce the supernatant endolysin compared with the sensitive bacteria or input control ( Figure  4B). The trivial reduction in E. faecium 81 bar or S. aureus bar may be due to the unspecific absorption. The results showed that the lysis by IME-EF1 endolysin requires the binding of the endolysin to the host cells.

Animal assay in vivo
When the mice were inoculated intraperitoneally with different amount of E. faecalis 002, the minimum amount that resulted to death in all mice within 56 h was 2×10 9 CFU. Thus, 2×10 9 CFU of E. faecalis 002 was used to construct the sepsis model in mice through intraperitoneal injection. The results indicated that one dose of bacteriophage IME-EF1 rescued the mice to some extent. Administering bacteriophage at 30 min post inoculation led to a therapeutic efficacy with a survival rate of 60%, whereas administering at 4 h provided a 40% protection rate ( Figure 5A). Similar rescue pattern was observed when the mice were administrated intraperitoneally with purified endolysin. Administering endolysin at 30 min post inoculation protected 80% of the mice from lethal challenge, but the protection efficacy was poor when endolysin was administered at 4 h. The results indicated that both phage IME-EF1 and its endolysin could make an effective protection against E. faecalis 002 challenge in mice, whereas the early administration resulted in better therapeutic efficacy than the late administration. At 4, 8, and 12 h post infection with E. faecalis 002, blood samples were collected from the mice tail vein, and bacterial count was figured out on BHI culture plate. Intraperitoneal administration with phage endolysin significantly reduced the number of bacteria in the blood, and early treatment at 30 min post infection had better bactericidal effect than at 4 h ( Figure 5B). No bacteria were detected on the plate of the phage treatment group (data not shown). In this assay, the administration with phage IME-EF1 and endolysin without E. faecalis 002 inoculation had no effect on mouse survival.

Discussion
E. faecalis is a Gram-positive, non-motile, facultative anaerobic microbe, formerly classified as part of the Group D Streptococcus system, and commensally inhabits the gastrointestinal tracts of humans and other mammals. E. faecalis has been found to be responsible for some lifethreatening infections in humans, especially in nosocomial environment, where the naturally high levels of antibiotic resistance found in E. faecalis contribute to its pathogenicity.
With the emergence of antibiotic-resistant (especially vancomycin-resistant) E. faecalis, alternative intervention tools or strategies have been pursued for treating these drugresistant Enterococci infection in future. Under these circumstances, a therapy based on phage or its lytic component acquired reconsideration. In this research, we isolated a lytic bacteriophage IME-EF1 from hospital sewage using a strain of multi-drug resistant E. faecalis, and explored its potential application in fighting E. faecalis infection.
Although genetic characterization of numerous E. faecalis bacteriophage genome have been determined in previous reports, the genomic sequence of the newly isolated IME-EF1 only shared highly homologous sequence with the recently isolated E. faecalis phage BC-611 and SAP6. The three strains of E. faecalis phage belong to the Siphoviridae family and own a similar length of genomic sequence (57,081, 53996, and 58619 nt) and G+C content (40.04%, 40.45%, and 40.00%) [22,23]. The three strains of E. faecalis phage can be clustered in the same sub-clade based on endolysin gene ( Figure 3A) or whole genome sequence (data not shown). The mean identity of the whole genome between IME-EF1 and SAP6 is 81.46%. The minor structure protein gene, tail fiber protein gene, and DNA replicative gene are more heterologous than the whole genome sequence, and share an identity between 72% and 77%. The results demonstrated that host determining genes had a vast diversity even in highly homologous phage. Although basic biological characters have been determined, the bactericidal activity of the E. faecalis phage of BC-611 and SAP6 had not been evaluated yet. In this study, we characterized IME-EF1 and explored its potential application in controlling E. faecalis infection in vitro and in vivo. Considering their homology, the results obtained from this research could predict the potential bactericidal functions of BC-611 and SAP6. The sepsis mouse model was constructed in BALB/c mice through injecting E. faecalis intraperitoneally. One dose administration with bacteriophage IME-EF1 or its endolysin can rescue the mice from the lethal challenges of E. faecalis. Under the same amount of bacteriophage or endolysin, early administration at 30 min post infection presented better protection than at 4 h post infection. Early administration with bacteriophage or endolysin has been reasonably speculated to effectively control the dissemination of bacteria. The bacterial count result in blood supported the aforementioned speculation and proved that early administration of endolysin significantly reduces the amount of E. faecalis in the blood. Administering bacteriophage intraperitoneally also caused the significant reduction of bacteria in the blood (data not shown). However, considering the fact that bacteriophage in the blood sample Table 3. Lytic spectrum of IME-EF1 and its endolysin. possibly adsorbed to bacteria during the preparation and disturbed the subsequent clonal growth, the bacterial count of the groups rescued with bacteriophage was not presented in this study. The adsorption rate, latent time, and burst size are used to evaluate the phage lytic ability. The average latent time and burst size in the tailed phage are typically 40 min to 50 min and 50 PFU/infected to 100 PFU/infected bacteria, respectively [20]. The newly isolated phage IME-EF1 can form an approximately 1 mm diameter of plaque in the sensitive strains of bacteria, and has a latent time of 25 min and a burst size of 60 PFU/infected bacteria. Thus, although IME-EF1 has short latent time in host bacteria E. faecalis, it should be recognized as a moderately lytic phage. A broad lytic spectrum was considered as a prerequisite character of therapeutic candidate phage. In this study, although a clinical strain of E. faecium can be lysed by IME-EF1, the phage showed relatively narrow lytic spectrum on E. faecalis (3/10). E. faecalis phage SAP6 has been found to lyse its host bacteria, but could not lyse the other 10 strains of E. faecalis [23]. The results indicated that phages in the same sub-clade having homologous genomic components and arrangement may have similar bactericidal characteristics. Phage endolysin generally had broader lytic spectrum than its original phage, which was also observed in this research. The endolysin identified from IME-EF1 can lyse 10 strains of E. faecalis and 1 strain of E. faecium, which can also be infected and lysed by the original phage. Fortunately, the two strains of vancomycin-resistant E. faecalis that cannot be lysed by IME-EF1 could be lysed by its endolysin.
Bacteria have evolved numerous mechanisms to counteract phage infection [26]. Endolysin has been proven to target specific molecules in the host peptidoglycan that are essential for bacterial viability [27,28]. Thus, the occurrence of endolysinresistant bacteria is unlikely or seldom because phage have naturally evolved with their bacterial hosts over millions of years to produce the enzymes essential for the release of progeny phage. Repeated exposure of Streptococcus pneumoniae and Bacillus cereus to low concentrations of endolysin Pal and PlyG on agar plates and in liquid culture, respectively, did not result in the emergence of resistant mutants even after numerous cycles [29,30], making endolysin a hopeful candidate in treating multi-drug resistant bacterial infection in future.
Endolysin has been generally recognized to contain a twodomain structure [31,32]. The N-terminal catalyzing domain is responsible for cleaving peptidoglycan in the host membrane, whereas the C-terminal binding domain serves to recognize the host cell and bind to it. The catalytic domain is categorized into four different groups upon the 14 cleavage sites, which include N-acetylmuramidases, N-acetyl-β-D-glucosamidases, Nacetylmuramoyl-L-alanine amidases, and L-alanoyl-Dglutamate endopeptidases [28,33]. In this research, predicting the endolysin domain was performed with NCBI specialized BLAST tool conserved domains [34]. The CHAP (cysteine, histidine-dependent aminohydrolases) domain was found in the confirmed or putative endolysins from highly homologous E. faecalis bacteriophage IME-EF1, SAP6, and BC-611. This domain has been previously proven to correspond to an amidase function [35]. To date, all four kinds of N-terminal catalytic domains have been found in these submitted sequences of E. faecalis bacteriophage putative endolysin in GenBank. The N-acetylmuramoyl-L-alanine amidase domain is found in the N-terminal of PlyV12 [36], phiEF24C ORF9 [37], and EfaCPT1 putative endolysin (GenBank Accession Number AFO10817.1). Both EF62phi (GenBank Accession Number YP_006218727) and phiFL1A putative endolysin (GenBank Accession Number YP_003347517) gained an N-terminal GH25 (glycosylhydrolase family 25) muramidase [18,38,39].

Figure 4. Endolysin lytic activity and binding assay.
A. The OD600 descent is used to evaluate the lytic activity of IME-EF1 endolysin on the isolated E. faecalis, E. faecium, and S. aureus. B. Assay of IME-EF1 endolysin binding to sensitive bacterial strains. The expressed IME-EF1 endolysin was incubated with sensitive bacteria E. faecalis 002 and E. faecium 010, as well as the insensitive bacteria E. faecium 81 and S. aureus for 10 min, and the bacteria was then pelleted and the remaining supernatant endolysin was measured with ELISA at OD490. "Input" means no bacteria were added. Significance was determined using independent t test (* P<0.05). However, in contrast to EF62phi, phiFL1A putative endolysin contained a C-terminal LysM (Lysine Motif) domain, a small globular domain with approximately 40 amino acids and widespread protein module involved in binding peptidoglycan in bacteria and chitin in eukaryotes [40]. The putative amidase encoded by Enterococcus phage phiEf11 contains both peptidase family 23 and glucosamidase domains in its Nterminal (GenBank Accession Number YP_003358818) [21].
No domain or motif was found to be essential in C-terminal of endolysin ΦEF24C ORF9, according to the bioinformatics analysis, but a few amino acid deletion in the C-terminal resulted in the loss of function [37]. Bioinformatics prediction did not find any conserved domain in the C-terminal of IME-EF1 endolysin, but the 50 aa deletion in this region caused the mutant protein to be insoluble in E. coli (data not shown). Online prediction using ProtScale on the ExPASy server revealed that the C-terminal of IME-EF1 endolysin was full of hydrophilic amino acids and the deletion may destroy its hydrophilicity and stability in the solution [41,42]. The results demonstrated that, although some conserved domains had not been identified yet in the C-terminal of IME-EF1 endolysin or ΦEF24C ORF9, its role in endolysin functioning was essential.

Ethics statement
All animal studies were approved by the Institutional Animal Care and Use Committee at Beijing Institute of Microbiology and Epidemiology, Beijing, China, and performed according to institutional guidelines for animal welfare. Mice were housed (8 per cage) and maintained on a 12 h light: dark cycle (lights on at 7:00 am), with continuous access to food and water. At the end of the animal study, the mice were sacrificed by CO 2 asphyxiation in accordance with the guideline. Bacterial isolation and sewage collection were performed in Affiliated Hospital, Academy of Military Medical Sciences by Xiuyun Yin, who is a staff of the hospital. Animal studies were performed in the Experimental Animal Center at Beijing Institute of Microbiology and Epidemiology, and did not involve any endangered or protected species. No specific permissions were required for these locations and activities.

Phage isolation
The bacterial strains employed in this study were isolated in the Affiliated Hospital, Academy of Military Medical Sciences, Beijing, China, and shown in Table 1. Brain Heart Infusion (BHI) was purchased from Becton Dickinson and Company (Sparks, MD). The sewage water sample was collected from the same hospital, centrifuged at 10,000 × g for 10 min, and then filtered through a 0.45 µm pore membrane. Approximately 100 µl of the filtered water sample was added to 5 ml of logarithmic phase E. faecalis 002 in the BHI culture medium, and incubated overnight at 37 °C. After centrifugation (10,000× g, 4 °C, 10 min), the supernatant was collected and filtered through a filter with 0.45 µm pore membrane. The filtrate (0.1 ml) was mixed with 0.5 ml of E. faecalis 002 in BHI culture medium (OD600=0.6) and 5 ml of molten top soft nutrient agar (0.75% agar), which was then overlaid onto solidified base nutrient agar (1.5% agar) [43]. Following incubation for 6 h at 37 °C, the clear plaques by phage were picked from the plate. The phage titer was determined using the double-layered method previously described by Adams [44] .

Purification and characterization of bacteriophage IME-EF1
A single plaque was picked into 5 ml of BHI medium containing E. faecalis 002 (OD600 = 0.6) and cultured at 37 °C for 6 h. The entire 5 ml of suspension was transferred into 500 ml of BHI medium for culture at 37 °C overnight. Chloroform was then added into the 500 ml of culture to a final concentration of 0.1% before being mixing gently and allowed to stand at room temperature for about 30 min. Solid sodium chloride was added to the culture to a final concentration of 1 M, which was then incubated in an ice water bath for 1 h. The culture was centrifuged at 11,000 × g for 10 min to remove the cell debris. Subsequently, polyethylene glycol 8000 (PEG8000) was added into the supernatant to a final concentration of 10% (w/v) while slowly stirring with a magnetic stirrer at room temperature. The solution was then transferred to a polypropylene centrifuge tube in an ice water bath and incubated for at least 1 h to precipitate the phage particles. Following centrifugation (11,000 × g for 10 min at 4 °C), the phage-containing precipitate was suspended in 5 ml of SM buffer (5.8 g of NaCl, 1.2 g of MgSO4, 50 ml of 1 M Tris-HCl (pH7.5), 0.1 g of gelatin in 1000 ml of deionized water) [45]. An equal volume of chloroform was then added to separate the phage particles from PEG8000. After centrifugation at 3,000 × g for 10 min, the aqueous phase was recovered and filtered through a 0.22 µm pore-size membrane filter to remove debris. The purified phage was then stored at 4 °C. Bacteriophage was directly stained with 2% phosphotungstic acid on the carboncoated grid and examined with a transmission electron microscope (Philips TECNAI-10).

Determination of one-step growth curve
The one-step growth kinetics curve of IME-EF1 was measured as previously reported [46]. Briefly, a midexponential-phase culture of E. faecalis 002 (OD600=0.4) was harvested and infected with bacteriophage IME-EF1 at an M.O.I. of 0.001, and then allowed to adsorb for 10 min at room temperature. The mixture was centrifuged, suspended in 5 ml of BHI medium, and cultured at 37 °C. Samples were obtained at 5 min or 10 min intervals, and immediately diluted and plated for phage titration.

Extraction of genomic nucleic acid of IME-EF1
Approximately 200 µl of purified phage suspension was added to the lysis solution (final concentration, 0.5% SDS, 50 µg/ml proteinase K, 20 mM EDTA), and the mixture was incubated for 1 h at 56 °C. An equal volume of saturated phenol was added to the mixture, gently vortexed, and the culture medium was centrifuged (10,000×g, 4 °C, 5 min). The upper aqueous was collected, and an equal volume of a solution (phenol: chloroform: isoamyl alcohol=25:24:1) was added to the upper aqueous. After gently vortexing, the mixture was centrifuged (10,000×g, 4 °C, 5 min). The upper aqueous was collected and added with an equal volume of isopropanol. The mixture was allowed to stand at −20 °C for more than 30 min, and then centrifuged (10,000×g, 4 °C, 20 min). The precipitation was collected and washed twice with 75% ethanol.
The nucleic acids were suspended in deionized water and stored at −20 °C.

Whole genome sequencing and bioinformatics analysis
The genomic DNA of the isolated bacteriophage IME-EF1 was subjected to high-throughput sequencing using a Life Technologies Ion Personal Genome Machine Ion Torrent sequencer (San Francisco, CA) according to the instructions of the manufacturer. The complete genome sequence of the phage was assembled using Velvet [47] and CLC Bio (Aarhus, Denmark), and then annotated using the online tool RAST [48]. Sequence similarity analysis was performed using the NCBI BLAST program (National Library of Medicine). Finding domains in endolysin protein sequence was conducted with NCBI Specialized BLAST in CCD (conserved domain database). The alignment of putative endolysin proteins was carried out with software CLC Bio in the parameter of less accurate, Gap open cost 10, and Gap extension cost 1.

Identification and characterization of the phage endolysin
The genome of the phage was annotated with the online tool RAST. The coding sequence of the putative endolysin of the phage was amplified with primers (Forward: 5'gcggatccatggttaaattaaacgatgtac-3'; Reverse: 5'gctctagactatactttaacttgtggtaaagc-3'). The PCR fragment was cloned into the expression vector pColdI using restriction endonuclease XbaI and BamHI (NEB). The plasmid was sequenced and transformed into E. coli strain M15 for expression. When the transformed E. coli grew to an optical density 0.6 at 600 nm under the selection of 100 µg/ml ampicillin, the suspension was allowed to stand for 30 min at 16 °C. The growth medium was supplemented with IPTG to a final concentration of 1 mM, and cultured aerobically for 24 h at 16 °C. After centrifugation (10,000×g, 10 min, 4 °C), the cell pellets were collected, and then stored at −70 °C. To check the expression of the putative endolysin gene in E. coli, the cell pellet was sonicated and centrifuged (10,000×g, 5 min). The supernatant was dissolved in loading buffer (50 mM Tris-HCl, 2% sodium dodecyl sulfate-polyacrylamide, 0.1% bromophenol blue, 10% glycerol, and 1% β-mercaptoethanol), and then boiled for 5 min. The samples were separated on a 12% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE), and stained with Coomassie brilliant blue, and then transferred to polyvinylidene difluoride membranes. The membrane was blocked overnight at 4 °C by incubating with 5% skim milk, followed by incubation for 2 h at 37 °C with the mouse anti-His monoclonal antibody (Sigma Aldrich). Bound antibody was detected with an HRP-conjugated goat antimouse IgG (Sigma Aldrich) combined with enhanced chemiluminescent reagents (Millipore, Bedford, MA, USA). For endolysin purification, the transformed E. coli cells were sonicated at 200 W to 300 W for 5 seconds with an 8 second interval. After being centrifuged (10,000×g, 20 min, 4 °C), the supernatant was added to 5 volumes of Ni-Native-0 buffer (50 mM NaH 2 PO4, 300 mM NaCl, pH 8.0). The mixture was loaded on the chromatography column and flowed through the column with a rate of 0.5 to 1 ml/min. The column was washed with five volumes of Ni-Native-500 (50 mM NaH 2 PO4, 300 mM NaCl, 500 mM imidazole, pH 8.0) buffer with a flow rate of 1 ml/ min. The elution solution was then dialyzed against PBS buffer.

Lytic spectrum determination
A spot test was adopted to determine the lytic spectrum of the bacteriophage IME-EF1 and its endolysin. Briefly, 10 strains of E. faecalis, 10 strains of E. faecium, and 4 strains of S. aureus were inoculated respectively on BHI semi-solid agar medium (0.8%). After solidification, each culture plate was divided into three grids. Approximately 2 µl of IME-EF1 and 2 µl of endolysin were dropped onto the center of a grid respectively, and the culture plate was then incubated at 37 °C for 12 h. The same amount of bacteriophage suspension buffer or endolysin elution buffer was spotted on a grid as the negative control. A clear and transparent plaque indicates that the host bacteria can be lysed by IME-EF1 or its endolysin.
Lytic activity and binding assay of IME-EF1 endolysin E. faecalis, E. faecium, and S. aureus were grown to midexponential growth phase, washed three times with PBS, and then suspended in 3 ml of PBS. Twenty-five microliter of the purified endolysin protein (1 µg/ml) was added into the bacterial suspension, and the mixture was cultured on a shaking platform at 37 °C. The value of OD600 was detected at 1.5 h and 4 h, respectively. To evaluate whether binding to host cell is necessary for endolysin functioning, E. faecalis 002, E. faecium 010, and a strain of S. aureus suspended in 200 µl of PBS in 1.5 ml Eppendorf tube were incubated with 25 µl of the purified endolysin protein (1 µg/ml) for 10 min, and then pelleted (10,000×g, 1 min). The endolysin in the supernatant was measured with Enzyme-linked immunosorbent assay (ELISA). Briefly, 25 µl of the supernatant in 75 µl of coating buffer (0.2 M sodium carbonate/bicarbonate, pH 9.4) was immobilized on a polystyrene microtiter plate (Thermo). The primary antibody is mouse anti-His monoclonal antibody, whereas the secondary is an HRP-conjugated goat anti mouse IgG antibody. TMB (tetramethyl benzene) was used as substrate, and optical absorbance is measured at 490 nm.

Mouse model infected with E. faecalis 002 and therapeutic assay
Female (6 to 8 weeks old) BALB/c mice were purchased from the Laboratory Animal Center of the Academy of Military Medical Sciences, Beijing, China. E. faecalis 002 was grown to OD600=1.6, and then centrifuged (10,000×g, 10 min, 4 °C). The bacterial pellet was suspended in saline and counted on the plate of the solid BHI culture medium. To determine the minimum lethal dosage, 5×10 7 , 1×10 8 , 2×10 8 , 5×10 8 , 1×10 9 , 2×10 9 , and 5×10 9 CFU (Clony Forming Unit) bacteria were injected intraperitoneally in 100 µl of saline. Each concentration group contains 8 mice. The mouse sepsis model with E. faecalis 002 was constructed by injecting the minimum lethal dosage of E. faecalis 002 intraperitoneally. The 10 M.O.I of bacteriophage IME-EF1 or 0.2 mg of expressed endolysin in 100 µl of saline was administrated intraperitoneally on the other side of abdomen at 30min, 4 h after bacterial inoculation respectively. Then, 100 µl of saline was used as non-therapeutic control. Each group contained 8 mice. The 10 M.O.I. of the bacteriophage IME-EF1 and 0.2 mg of the expressed endolysin in 100 µl of saline was administrated intraperitoneally as negative control. All the survival mice were sacrificed by CO2 asphyxiation at 96 h post inoculation.
The value of the statistical analysis of OD600 value descent, and endolysin binding capacity data was performed using an independent t test. Survival data were analyzed using a Kaplan-Meier survival analysis with a log-rank method of statistics. P<0.05 were considered significant.


extfieldsofstudy:
Biology
Computer Science
Medicine

text:
Integration-free Methods for Generating Induced Pluripotent Stem Cells

Induced pluripotent stem (iPS) cells can be generated from mouse or human fibroblasts by exogenous expression of four factors, Oct4, Sox2, Klf4 and c-Myc, and hold great potential for transplantation therapies and regenerative medicine. However, use of retroviral vectors during iPS cell generation has limited the technique’s clinical application due to the potential risks resulting from genome integration of transgenes, including insertional mutations and altered differentiation potentials of the target cells, which may lead to pathologies such as tumorigenesis. Here we review recent progress in generating safer transgene-free or integration-free iPS cells, including the use of non-integrating vectors, excision of vectors after integration, DNA-free delivery of factors and chemical induction of pluripotency.

Introduction
In 2006, induced pluripotent stem (iPS) cells were first generated from mouse fibroblasts by forced expression of four protein factors, Oct4, Klf4, Sox2 and c-Myc, mediated by retrovirus [1]. Human cells can also be reprogrammed using the same set of four factors [2] or Oct4, Sox2, Lin28 and Nanog [3]. iPS cells were proven to exhibit pluripotency at levels similar to embryonic stem (ES) cells. Thus patient-specific pluripotent stem cells can be produced without the problems of immune rejection for disease models, transplantation therapies and regenerative medicine.
However, iPS technology is complicated by the potential risks posed by genome-integrating viruses. The transgenes are randomly but permanently integrated into the host genome at multiple sites together with viral vector backbone. The expression of exogenous transgenes is silenced after reprogramming is achieved, and this silencing is indispensable for maintaining pluripotency. Such genome-integrating viral vectors can produce insertional mutations which may influence differentiation potential, or even result in tumorigenesis due to reactivation of the c-Myc oncogene [4]. Reprogramming methods that utilize viral vectors are therefore judged too risky to be used in clinical therapies. Recently, several approaches have been developed to generate safer transgene-free or integration-free iPS cells, as listed in Table 1.

Non-integrating vectors
A series of non-integrating vectors, including viral vectors and non-viral vectors, have been used successfully for generating iPS cells.

Adenovirus
Replication-incompetent adenoviral vectors were used to generate iPS cells from mouse [5] and human [6] somatic cells. These vectors allow for transient expression of exogenous genes without integration into the host genome, although with a low efficiency compared to integrating vectors [5].

Sendai virus
Sendai virus (SeV) is thought to be able to replace retrovirus, since it is also an RNA virus and can be easily removed by antibody-mediated negative selection. SeV replicates in the form of single-stranded RNA in the cytoplasm of infected cells, therefore DNA would be neither produced nor integrated into the host genome during transduction. Thus it has been used as a safe vector in the field of gene therapy. SeV-derived transgenes were expressed three days after transduction. The resulting iPS cells were generated with a high efficiency [7], and the SeV vectors were then diluted by cell passages [7].

Expression plasmids
Several classes of vectors have been used as viral substitutes to produce transgene-free iPS cells without integration, and expression plasmids were among the first to be investigated. Repeated transfection of the plasmids containing the four Yamanaka factors into mouse embryonic fibroblasts (MEF) resulted in iPS cells without integration, although with a much lower efficiency than viral vectors [8].

Episomal vectors
Other types of non-integrating vectors were explored to introduce reprogramming factors into mouse or human somatic cells. For example, portions of the Epstein-Barr human herpesvirus were used to create a vector that can be transfected without viral packaging. The resulting extrachromosomal oriP/EBNA1 (Epstein-Barr nuclear antigen-1) episomal vector also features a drug selection mechanism for removing the vector from cells. The oriP/EBNA1 vector undergoes stable extrachromosomal replication only once per cell cycle, without integrating into the host genome. However, efficiency for iPS generation using episomal vector was very low (3-6 · 10 À6 ) [9]. Okita et al. [10] further improved the episomal vector by using p53 suppression and non-transforming L-Myc instead of c-Myc, resulting in markedly enhanced iPS generation (1 · 10 À5 -3 · 10 À4 ).

Minicircle vector
Minicircle expression vectors (i.e., supercoiled DNA molecules that lack a bacterial origin of replication) have higher transfection efficiencies and longer-term transgene expression compared to plasmids due to reduced silencing for exogenous genes. A single minicircle vector containing four reprogramming factors was introduced into human ES cells by nucleofection, and the minicircle plasmid backbone was then excised and degraded using the PhiC31-based intramolecular recombination system [11].

Liposomal magnetofection
For liposomal magnetofection (LMF), cationic lipids mediate the self-assembly of complexes containing plasmids and nanoparticles of superparamagnetic iron. These ''CombiMag-DNA'' ternary complexes can be concentrated at cell surfaces using a strong magnetic field to transfect vectors into the targeted cells. An optimized LMF protocol for generating iPS cells achieved short reprogramming times of 8 days or less. In addition, out of seven independently-generated iPS lines, two iPS lines that had no integrated vector were produced [12].

Deletion after integration
Two systems have been applied to remove integrated transgenes from mouse or human iPS cells.

Cre-loxP system
Cre/loxP recombination was used to excise integrated transgenes from iPS cells. Human iPS cells were produced using a single retroviral vector carrying OCT4, SOX2, KLF4 and c-MYC linked via picornaviral 2A plasmids. This cassette was removed by transfection of Cre recombinase after reprogramming was achieved. However, residual vector sequences were left behind, therefore insertional mutations remain a risk [13]. After complete reprogramming, excision of the piggyBac vector and its induction factors from iPS cell lines was achieved with a pulse of PB transposase expression. In 10 out of 11 subclone lines, the sequence at the transposon integration site reverted to wild type after vector excision [14,15].

DNA-free delivery Proteins
Alternative methods to avoid introducing genetic modifications include delivery of the reprogramming proteins or mRNA directly into cells, rather than their expression from DNA. These methods have been successfully demonstrated but can be much more complicated to perform.
To generate recombinant proteins that can penetrate and cross the plasma membrane of somatic cells, the carboxy termini of four reprogramming factors were fused to a poly-arginine protein transduction domain. After being added to the cell culture medium, the recombinant transcription factors readily entered cells within 6 h and could translocate to the nucleus. iPS cells were obtained after four repeated protein transductions at 8 mg/ml with 1 mM valproic acid (VPA), a histone deacetylase inhibitor that can significantly improve reprogramming efficiency [16].
Similar protein transduction was also accomplished by Kim et al. [17]. Human iPS cell-like colonies were established within 8 weeks after six rounds of transduction using four reprogramming factors fused with 9-arginine and myc tags. The resulting induction efficiency was approximately twice as high as that using viral transduction.

mRNAs
The use of synthetic RNAs for reprogramming was reported by Warren et al. with surprisingly high reprogramming efficiencies [18]. mRNAs including 3 0 and 5 0 untranslated regions were synthesized via in vitro transcription. RNA transfection with a cationic vehicle into MEF and human somatic cells revealed a high, dose-dependent cytotoxicity. Modified ribonucleotides and phosphatase treatment, in conjunction with media supplemented with the interferon inhibitor B18R, were used to optimize the synthetic mRNA approach and improve cell viability. Maximal protein expression was observed 12-18 h after transfection, followed by a rapid turnover that necessitates subsequent daily transfections. A five-factor cocktail including a modified LIN28 RNA produced reprogramming efficiencies of 1.4%, much higher than using the virus-based approach (10 À4 ), and with a shorter time period (around 18 days) [18].

Chemical induction
Small molecules have advantages because they are non-immunogenic, and can be more easily administered and standardized. Moreover, their effects on inhibiting or activating specific proteins are often reversible and dose-dependent. Small molecule libraries and combinations of compounds have been screened to identify substitutes for transgenic DNA delivery. Initially, iPS cells were generated from mouse fibroblasts using a single gene, Oct4, with a small molecule combination termed ''VC6T'' [VPA, CHIR99021 (CHIR), 616452 and tranylcypromine] [19]. Next, chemical substitutes for Oct4 were screened as supplements to VC6T to obtain transgene-free iPS cells. In the Oct4 promoterdriven GFP expression (OG) system, some GFP-positive clusters were induced by VC6T plus forskolin (VC6TF), however they lacked expression of Oct4 and Nanog, indicating incomplete reprogramming [20]. To identify the chemicals critical for late reprogramming, a doxycycline (DOX)-inducible Oct4 expression system was used in which DOX was added only in the first 4-8 days. Several cAMP agonists including forskolin, prostaglandin E2 and rolipram, and epigenetic modulators such as 3-deazaneplanocin A (DZNep), 5-azacytidine, sodium butyrate and RG108 were identified [20]. Among all the candidates, DZNep was proven to work well. DZNep was added 16 days after treatment with VC6TF, followed by a day 28 switch to 2i-medium for dual inhibition (2i) of glycogen synthase kinase-3 (GSK-3) and mitogen-activated protein kinase (MAPK) signaling. The resulting chemically induced pluripotent stem cells (CiPSCs) developed into colonies with ES-like properties. Efficiency can be boosted by a synthetic retinoic acid receptor ligand, TTNPB, for up to 0.2% cells induced [20].
DZNep is critical for activating endogenous Oct4. As an S-adenosyl homocysteine (SAH) hydrolase inhibitor, DZNep may repress the SAM-dependent cellular methylation process and significantly decrease DNA and H3K9 methylation at the Oct4 promoter, which may contribute to Oct4 activation. These findings indicate that appropriate combinations of inducing chemicals and previously characterized drugs can provide attractive mechanisms and applications for potential clinical utility, bypassing nucleic acid manipulations and the need to deploy protein factors.

Perspective
The efficiency and safety of generating and using iPS cells seem to show a negative correlation, and thus clinical applications of iPS technology await the validation of a mature protocol that better balances these two important elements. Generation of iPS cells by recombinant proteins is still worth considering for clinical use due to its high safety, despite being quite expensive and having very low efficiency. Recent reports point out that activation of toll-like receptor-3 (TLR3) leads to epigenetic modification and facilitates nuclear reprogramming by the four inducing factors [21]. Therefore, TLR3 manipulation may support nuclear reprogramming with a much higher efficiency compared to using retrovirally encoded proteins alone. Additionally, mRNA transduction was reported to have high efficiency, but the rapid turnover of mRNA complicates this approach and necessitates repeated, daily transfections. The newly reported chemical induction method is promising for clinical use, owing to its high safety and feasibility. Moreover, the chimeric mice generated from CiPSCs were 100% viable and apparently healthy, unlike those generated from retrovirus-induced iPS cells. Chemical library screening for TLR3 activators could lead to refinements in the CiPSC induction cocktail and protocol that further boost efficiency. Additional verification of CiPSCs including tetraploid complementation is needed to confirm stem cell function before its extensive application. Nonetheless, CiP-SCs provide a new tool to study the mechanisms of reprogramming that may lead to more improvements in iPS cell production and safety, and eventually powerful iPS cell therapeutics.