Ancestral sequence reconstruction dissects structural and functional differences among eosinophil ribonucleases

Evolutionarily conserved structural folds can give rise to diverse biological functions, yet predicting atomic-scale interactions that contribute to the emergence of novel activities within such folds remains challenging. Pancreatic-type ribonucleases illustrate this complexity, sharing a core structure that has evolved to accommodate varied functions. In this study, we used ancestral sequence reconstruction to probe evolutionary and molecular determinants that distinguish biological activities within eosinophil members of the RNase 2/3 subfamily. Our investigation unveils functional, structural, and dynamical behaviors that differentiate the evolved ancestral ribonuclease (AncRNase) from its contemporary eosinophil RNase orthologs. Leveraging the potential of ancestral reconstruction for protein engineering, we used AncRNase predictions to design a minimal 4-residue variant that transforms human RNase 2 into a chimeric enzyme endowed with the antimicrobial and cytotoxic activities of RNase 3 members. This work provides unique insights into mutational and evolutionary pathways governing structure, function, and conformational states within the eosinophil RNase subfamily, offering potential for targeted modulation of RNase-associated functions.

Bacterial vampirism mediated through taxis to serum

Bacteria of the family Enterobacteriaceae are associated with gastrointestinal (GI) bleeding and bacteremia and are a leading cause of death, from sepsis, for individuals with inflammatory bowel diseases. The bacterial behaviors and mechanisms underlying why these bacteria are prone to bloodstream entry remains poorly understood. Herein, we report that clinical isolates of non-typhoidal Salmonella enterica serovars, Escherichia coli, and Citrobacter koseri are rapidly attracted toward sources of human serum. To simulate GI bleeding, we utilized a custom injection-based microfluidics device and found that femtoliter volumes of human serum are sufficient to induce the bacterial population to swim toward and aggregate at the serum source. This response is orchestrated through chemotaxis, and a major chemical cue driving chemoattraction is L-serine, an amino acid abundant in serum that is recognized through direct binding by the chemoreceptor Tsr. We report the first crystal structures of Salmonella Typhimurium Tsr in complex with L-serine and identify a conserved amino acid recognition motif for L-serine shared among Tsr orthologues. By mapping the phylogenetic distribution of this chemoreceptor we found Tsr to be widely conserved among Enterobacteriaceae and numerous World Health Organization priority pathogens associated with bloodstream infections. Lastly, we find that Enterobacteriaceae use human serum as a source of nutrients for growth and that chemotaxis and the chemoreceptor Tsr provides a competitive advantage for migration into enterohaemorrhagic lesions. We term this bacterial behavior of taxis toward serum, colonization of hemorrhagic lesions, and the consumption of serum nutrients, as 'bacterial vampirism' which may relate to the proclivity of Enterobacteriaceae for bloodstream infections.

Navigating contradictions: Salmonella Typhimurium chemotactic responses to conflicting chemoeffector signals show parity with bacterial growth benefits

Many bacteria that colonize the guts of animals use chemotaxis to direct swimming motility and select sites for colonization based on sources of effectors derived from the host, diet, and microbial competitors of the local environ. The complex ecosystem of the gastrointestinal tract contains mixtures of chemoattractants and chemorepellents, but it remains poorly understood how swimming bacteria navigate conflicting signals. The enteric pathogen Salmonella Typhimurium possesses Tsr, a chemoreceptor protein that directs both chemoattraction and chemorepulsion responses, which we employed as a model to study chemotaxis in the presence of conflicting effector stimuli. We investigated how S. Typhimurium responds to human fecal matter, an effector source in the enteric lumen that contains high concentrations of indole, a bacteriostatic chemorepellent produced by the native commensals of the microbiota, and also nutrients such as l-serine, a chemoattractant. The indole concentration in human feces is more than 12-fold the concentration required for half-maximal chemorepulsion, however, we find S. Typhimurium, and various clinical isolates of non-typhoidal S. enterica serovars, are strongly attracted to liquid fecal matter. We further investigated the chemotactic responses of S. Typhimurium to titrations of indole and l-serine and revealed that chemorepulsion to indole is overridden in the presence of excess l-serine. We capture the inversion of these two opposing taxis behaviors in a phenomenon we define as "chemohalation" in which the bacteria organize into a halo around the treatment source with an interior zone of avoidance, which represents a compromise between chemoattraction and chemorepulsion. Growth analyses reveal that the chemotactic responses to these opposing effectors align chemoattraction and chemorepulsion with the relative growth of the bacteria in culture. Hence, our study supports the view that evolution has finely tuned chemotaxis to assess environmental habitability by evaluating the tradeoffs in bacterial growth based on the local combination of effectors.

Ancestral Reconstruction and the Evolution of Protein Energy Landscapes

A protein's sequence determines its conformational energy landscape. This, in turn, determines the protein's function. Understanding the evolution of new protein functions therefore requires understanding how mutations alter the protein energy landscape. Ancestral sequence reconstruction (ASR) has proven a valuable tool for tackling this problem. In ASR, one phylogenetically infers the sequences of ancient proteins, allowing characterization of their properties. When coupled to biophysical, biochemical, and functional characterization, ASR can reveal how historical mutations altered the energy landscape of ancient proteins, allowing the evolution of enzyme activity, altered conformations, binding specificity, oligomerization, and many other protein features. In this article, we review how ASR studies have been used to dissect the evolution of energy landscapes. We also discuss ASR studies that reveal how energy landscapes have shaped protein evolution. Finally, we propose that thinking about evolution from the perspective of an energy landscape can improve how we approach and interpret ASR studies. Expected final online publication date for the Annual Review of Biophysics, Volume 53 is May 2024. Please see http://www.annualreviews.org/page/journal/pubdates for revised estimates.

Simultaneous Mapping of DNA Binding and Nucleosome Positioning with SpLiT-ChEC

The organization of chromatin - including the positions of nucleosomes and the binding of other proteins to DNA - helps define transcriptional profiles in eukaryotic organisms. While techniques like ChIP-Seq and MNase-Seq can map protein-DNA and nucleosome localization separately, assays designed to simultaneously capture nucleosome positions and protein-DNA interactions can produce a detailed picture of the chromatin landscape. Most assays that monitor chromatin organization and protein binding rely on antibodies, which often exhibit nonspecific binding, and/or the addition of bulky adducts to the DNA-binding protein being studied, which can affect their expression and activity. Here, we describe SpyCatcher Linked Targeting of Chromatin Endogenous Cleavage (SpLiT-ChEC), where a 13-amino acid SpyTag peptide, appended to a protein of interest, serves as a highly-specific targeting moiety for in situ enzymatic digestion. The SpyTag/SpyCatcher system forms a covalent bond, linking the target protein and a co-expressed MNase-SpyCatcher fusion construct. SpyTagged proteins are expressed from endogenous loci, whereas MNase-SpyCatcher expression is induced immediately before harvesting cultures. MNase is activated with high concentrations of calcium, which primarily digests DNA near target protein binding sites. By sequencing the DNA fragments released by targeted MNase digestion, we found that this method recovers information on protein binding and proximal nucleosome positioning. SpLiT-ChEC provides precise temporal control that we anticipate can be used to monitor chromatin under various conditions and at distinct points in the cell cycle.

Disentangling contact and ensemble epistasis in a riboswitch

Mutations introduced into macromolecules often exhibit epistasis, where the effect of one mutation alters the effect of another. Knowing the mechanisms that lead to epistasis is important for understanding how macromolecules work and evolve, as well as for effective macromolecular engineering. Here, we investigate the interplay between "contact epistasis" (epistasis arising from physical interactions between mutated residues) and "ensemble epistasis" (epistasis that occurs when a mutation redistributes the conformational ensemble of a macromolecule, thus changing the effect of the second mutation). We argue that the two mechanisms can be distinguished in allosteric macromolecules by measuring epistasis at differing allosteric effector concentrations. Contact epistasis manifests as nonadditivity in the microscopic equilibrium constants describing the conformational ensemble. This epistatic effect is independent of allosteric effector concentration. Ensemble epistasis manifests as nonadditivity in thermodynamic observables-such as ligand binding-that are determined by the distribution of ensemble conformations. This epistatic effect strongly depends on allosteric effector concentration. Using this framework, we experimentally investigated the origins of epistasis in three pairwise mutant cycles introduced into the adenine riboswitch aptamer domain by measuring ligand binding as a function of allosteric effector concentration. We found evidence for both contact and ensemble epistasis in all cycles. Furthermore, we found that the two mechanisms of epistasis could interact with each other. For example, in one mutant cycle we observed 6 kcal/mol of contact epistasis in a microscopic equilibrium constant. In that same cycle, the maximum epistasis in ligand binding was only 1.5 kcal/mol: shifts in the ensemble masked the contribution of contact epistasis. Finally, our work yields simple heuristics for identifying contact and ensemble epistasis based on measurements of a biochemical observable as a function of allosteric effector concentration.

Topiary: Pruning the manual labor from ancestral sequence reconstruction

Ancestral sequence reconstruction (ASR) is a powerful tool to study the evolution of proteins and thus gain deep insight into the relationships among protein sequence, structure, and function. A major barrier to its broad use is the complexity of the task: it requires multiple software packages, complex file manipulations, and expert phylogenetic knowledge. Here we introduce topiary, a software pipeline that aims to overcome this barrier. To use topiary, users prepare a spreadsheet with a handful of sequences. Topiary then: (1) Infers the taxonomic scope for the ASR study and finds relevant sequences by BLAST; (2) Does taxonomically informed sequence quality control and redundancy reduction; (3) Constructs a multiple sequence alignment; (4) Generates a maximum-likelihood gene tree; (5) Reconciles the gene tree to the species tree; (6) Reconstructs ancestral amino acid sequences; and (7) Determines branch supports. The pipeline returns annotated evolutionary trees, spreadsheets with sequences, and graphical summaries of ancestor quality. This is achieved by integrating modern phylogenetics software (Muscle5, RAxML-NG, GeneRax, and PastML) with online databases (NCBI and the Open Tree of Life). In this paper, we introduce non-expert readers to the steps required for ASR, describe the specific design choices made in topiary, provide a detailed protocol for users, and then validate the pipeline using datasets from a broad collection of protein families. Topiary is freely available for download: https://github.com/harmslab/topiary.

Higher-order interactions shape microbial interactions as microbial community complexity increases

Non-pairwise interactions, or higher-order interactions (HOIs), in microbial communities have been described as significant drivers of emergent features in microbiomes. Yet, the re-organization of microbial interactions between pairwise cultures and larger communities remains largely unexplored from a molecular perspective but is central to our understanding and further manipulation of microbial communities. Here, we used a bottom-up approach to investigate microbial interaction mechanisms from pairwise cultures up to 4-species communities from a simple microbiome (Hafnia alvei, Geotrichum candidum, Pencillium camemberti and Escherichia coli). Specifically, we characterized the interaction landscape for each species combination involving E. coli by identifying E. coli's interaction-associated mutants using an RB-TnSeq-based interaction assay. We observed a deep reorganization of the interaction-associated mutants, with very few 2-species interactions conserved all the way up to a 4-species community and the emergence of multiple HOIs. We further used a quantitative genetics strategy to decipher how 2-species interactions were quantitatively conserved in higher community compositions. Epistasis-based analysis revealed that, of the interactions that are conserved at all levels of complexity, 82% follow an additive pattern. Altogether, we demonstrate the complex architecture of microbial interactions even within a simple microbiome, and provide a mechanistic and molecular explanation of HOIs.

Probing the origin of prion protein misfolding via reconstruction of ancestral proteins

Prion diseases are fatal neurodegenerative diseases caused by pathogenic misfolding of the prion protein, PrP. They are transmissible between hosts, and sometimes between different species, as with transmission of bovine spongiform encephalopathy to humans. Although PrP is found in a wide range of vertebrates, prion diseases are seen only in certain mammals, suggesting that infectious misfolding was a recent evolutionary development. To explore when PrP acquired the ability to misfold infectiously, we reconstructed the sequences of ancestral versions of PrP from the last common primate, primate-rodent, artiodactyl, placental, bird, and amniote. Recombinant ancestral PrPs were then tested for their ability to form β-sheet aggregates, either spontaneously or when seeded with infectious prion strains from human, cervid, or rodent species. The ability to aggregate developed after the oldest ancestor (last common amniote), and aggregation capabilities diverged along evolutionary pathways consistent with modern-day susceptibilities. Ancestral bird PrP could not be seeded with modern-day prions, just as modern-day birds are resistant to prion disease. Computational modeling of structures suggested that differences in helix 2 could account for the resistance of ancestral bird PrP to seeding. Interestingly, ancestral primate PrP could be converted by all prion seeds, including both human and cervid prions, raising the possibility that species descended from an ancestral primate have retained the susceptibility to conversion by cervid prions. More generally, the results suggest that susceptibility to prion disease emerged prior to ~100 million years ago, with placental mammals possibly being generally susceptible to disease.

Evolution avoids a pathological stabilizing interaction in the immune protein S100A9

Stability constrains evolution. While much is known about constraints on destabilizing mutations, less is known about the constraints on stabilizing mutations. We recently identified a mutation in the innate immune protein S100A9 that provides insight into such constraints. When introduced into human S100A9, M63F simultaneously increases the stability of the protein and disrupts its natural ability to activate Toll-like receptor 4. Using chemical denaturation, we found that M63F stabilizes a calcium-bound conformation of hS100A9. We then used NMR to solve the structure of the mutant protein, revealing that the mutation distorts the hydrophobic binding surface of hS100A9, explaining its deleterious effect on function. Hydrogen-deuterium exchange (HDX) experiments revealed stabilization of the region around M63F in the structure, notably Phe37. In the structure of the M63F mutant, the Phe37 and Phe63 sidechains are in contact, plausibly forming an edge-face π-stack. Mutating Phe37 to Leu abolished the stabilizing effect of M63F as probed by both chemical denaturation and HDX. It also restored the biological activity of S100A9 disrupted by M63F. These findings reveal that Phe63 creates a molecular staple with Phe37 that stabilizes a nonfunctional conformation of the protein, thus disrupting function. Using a bioinformatic analysis, we found that S100A9 proteins from different organisms rarely have Phe at both positions 37 and 63, suggesting that avoiding a pathological stabilizing interaction indeed constrains S100A9 evolution. This work highlights an important evolutionary constraint on stabilizing mutations, namely, that they must avoid inappropriately stabilizing nonfunctional protein conformations.

Ensemble epistasis: thermodynamic origins of nonadditivity between mutations

Epistasis-when mutations combine nonadditively-is a profoundly important aspect of biology. It is often difficult to understand its mechanistic origins. Here, we show that epistasis can arise from the thermodynamic ensemble, or the set of interchanging conformations a protein adopts. Ensemble epistasis occurs because mutations can have different effects on different conformations of the same protein, leading to nonadditive effects on its average, observable properties. Using a simple analytical model, we found that ensemble epistasis arises when two conditions are met: (1) a protein populates at least three conformations and (2) mutations have differential effects on at least two conformations. To explore the relative magnitude of ensemble epistasis, we performed a virtual deep-mutational scan of the allosteric Ca2+ signaling protein S100A4. We found that 47% of mutation pairs exhibited ensemble epistasis with a magnitude on the order of thermal fluctuations. We observed many forms of epistasis: magnitude, sign, and reciprocal sign epistasis. The same mutation pair could even exhibit different forms of epistasis under different environmental conditions. The ubiquity of thermodynamic ensembles in biology and the pervasiveness of ensemble epistasis in our dataset suggests that it may be a common mechanism of epistasis in proteins and other macromolecules.

Were Ancestral Proteins Less Specific?

Some have hypothesized that ancestral proteins were, on average, less specific than their descendants. If true, this would provide a universal axis along which to organize protein evolution and suggests that reconstructed ancestral proteins may be uniquely powerful tools for protein engineering. Ancestral sequence reconstruction studies are one line of evidence used to support this hypothesis. Previously, we performed such a study, investigating the evolution of peptide-binding specificity for the paralogs S100A5 and S100A6. The modern proteins appeared more specific than their last common ancestor (ancA5/A6), as each paralog bound a subset of the peptides bound by ancA5/A6. In this study, we revisit this transition, using quantitative phage display to measure the interactions of 30,533 random peptides with human S100A5, S100A6, and ancA5/A6. This unbiased screen reveals a different picture. While S100A5 and S100A6 do indeed bind to a subset of the peptides recognized by ancA5/A6, they also acquired new peptide partners outside of the set recognized by ancA5/A6. Our previous work showed that ancA5/A6 had lower specificity than its descendants when measured against biological targets; our new work shows that ancA5/A6 has similar specificity to the modern proteins when measured against a random set of peptide targets. This demonstrates that altered biological specificity does not necessarily indicate altered intrinsic specificity, and sounds a cautionary note for using ancestral reconstruction studies with biological targets as a means to infer global evolutionary trends in specificity.

Evolutionary Conservation of Structural and Functional Coupling between the BRM AT-Hook and Bromodomain

The BAF chromatin remodeling complex is critical for genome regulation. The central ATPase of BAF is either BRM or BRG1, both of which contain a C-terminal bromodomain, known to associate with acetylated lysines. We have recently demonstrated that in addition to acetyl-lysine binding, the BRG1/BRM bromodomain can associate with DNA through a lysine/arginine rich patch that is adjacent to the acetyl-lysine binding pocket. Flanking the bromodomain is an AT-hook separated by a short, proline-rich linker. We previously found that the AT-hook and bromodomain can associate with DNA in a multivalent manner. Here, we investigate the conservation of this composite module and find that the AT-hook, linker, and lysine/arginine rich bromodomain patch are ancient, conserved over ~1 billion years. We utilize extensive mutagenesis, NMR spectroscopy, and fluorescence anisotropy to dissect the contribution of each of these conserved elements in association of this module with DNA. Our results reveal a structural and functional coupling of the AT-hook and bromodomain mediated by the linker. The lysine/arginine rich patch on the bromodomain and the conserved elements of the AT-hook are critical for robust affinity for DNA, while the conserved elements of the linker are dispensable for overall DNA affinity but critical for maintaining the relative conformation of the AT-hook and bromodomain in binding to DNA. This supports that the coupled action of the AT-hook and bromodomain are important for BAF activity.

Identification and Characterization of Zebrafish Tlr4 Coreceptor Md-2

The zebrafish (Danio rerio) is a powerful model organism for studies of the innate immune system. One apparent difference between human and zebrafish innate immunity is the cellular machinery for LPS sensing. In amniotes, the protein complex formed by TLR4 and myeloid differentiation factor 2 (Tlr4/Md-2) recognizes the bacterial molecule LPS and triggers an inflammatory response. It is believed that zebrafish have neither Md-2 nor Tlr4; Md-2 has not been identified outside of amniotes, whereas the zebrafish tlr4 genes appear to be paralogs, not orthologs, of amniote TLR4s We revisited these conclusions. We identified a zebrafish gene encoding Md-2, ly96 Using single-cell RNA sequencing, we found that ly96 is transcribed in cells that also transcribe genes diagnostic for innate immune cells, including the zebrafish tlr4-like genes. In larval zebrafish, ly96 is expressed in a small number of macrophage-like cells. In a functional assay, zebrafish Md-2 and Tlr4ba form a complex that activates NF-κB signaling in response to LPS. In larval zebrafish ly96 loss-of-function mutations perturbed LPS-induced cytokine production but gave little protection against LPS toxicity. Finally, by analyzing the genomic context of tlr4 genes in 11 jawed vertebrates, we found that tlr4 arose prior to the divergence of teleosts and tetrapods. Thus, an LPS-sensitive Tlr4/Md-2 complex is likely an ancestral feature shared by mammals and zebrafish, rather than a de novo invention on the tetrapod lineage. We hypothesize that zebrafish retain an ancestral, low-sensitivity Tlr4/Md-2 complex that confers LPS responsiveness to a specific subset of innate immune cells.

Exploring the Evolutionary History of Kinetic Stability in the α-Lytic Protease Family

In addition to encoding the tertiary fold and stability, the primary sequence of a protein encodes the folding trajectory and kinetic barriers that determine the speed of folding. How these kinetic barriers are encoded is not well understood. Here, we use evolutionary sequence variation in the α-lytic protease (αLP) protein family to probe the relationship between sequence and energy landscape. αLP has an unusual energy landscape: the native state of αLP is not the most thermodynamically favored conformation and, instead, remains folded due to a large kinetic barrier preventing unfolding. To fold, αLP utilizes an N-terminal pro region similar in size to the protease itself that functions as a folding catalyst. Once folded, the pro region is removed, and the native state does not unfold on a biologically relevant time scale. Without the pro region, αLP folds on the order of millennia. A phylogenetic search uncovers αLP homologs with a wide range of pro region sizes, including some with no pro region at all. In the resulting phylogenetic tree, these homologs cluster by pro region size. By studying homologs naturally lacking a pro region, we demonstrate they can be thermodynamically stable, fold much faster than αLP, yet retain the same fold as αLP. Key amino acids thought to contribute to αLP's extreme kinetic stability are lost in these homologs, supporting their role in kinetic stability. This study highlights how the entire energy landscape plays an important role in determining the evolutionary pressures on the protein sequence.

Learning peptide recognition rules for a low-specificity protein

Many proteins interact with short linear regions of target proteins. For some proteins, however, it is difficult to identify a well-defined sequence motif that defines its target peptides. To overcome this difficulty, we used supervised machine learning to train a model that treats each peptide as a collection of easily-calculated biochemical features rather than as an amino acid sequence. As a test case, we dissected the peptide-recognition rules for human S100A5 (hA5), a low-specificity calcium binding protein. We trained a Random Forest model against a recently released, high-throughput phage display dataset collected for hA5. The model identifies hydrophobicity and shape complementarity, rather than polar contacts, as the primary determinants of peptide binding specificity in hA5. We tested this hypothesis by solving a crystal structure of hA5 and through computational docking studies of diverse peptides onto hA5. These structural studies revealed that peptides exhibit multiple binding modes at the hA5 peptide interface-all of which have few polar contacts with hA5. Finally, we used our trained model to predict new, plausible binding targets in the human proteome. This revealed a fragment of the protein α-1-syntrophin that binds to hA5. Our work helps better understand the biochemistry and biology of hA5, as well as demonstrating how high-throughput experiments coupled with machine learning of biochemical features can reveal the determinants of binding specificity in low-specificity proteins.

Inferring a complete genotype-phenotype map from a small number of measured phenotypes

Understanding evolution requires detailed knowledge of genotype-phenotype maps; however, it can be a herculean task to measure every phenotype in a combinatorial map. We have developed a computational strategy to predict the missing phenotypes from an incomplete, combinatorial genotype-phenotype map. As a test case, we used an incomplete genotype-phenotype dataset previously generated for the malaria parasite’s ‘chloroquine resistance transporter’ (PfCRT). Wild-type PfCRT (PfCRT^3D7) lacks significant chloroquine (CQ) transport activity, but the introduction of the eight mutations present in the ‘Dd2’ isoform of PfCRT (PfCRT^Dd2) enables the protein to transport CQ away from its site of antimalarial action. This gain of a transport function imparts CQ resistance to the parasite. A combinatorial map between PfCRT^3D7 and PfCRT^Dd2 consists of 256 genotypes, of which only 52 have had their CQ transport activities measured through expression in the Xenopus laevis oocyte. We trained a statistical model with these 52 measurements to infer the CQ transport activity for the remaining 204 combinatorial genotypes between PfCRT^3D7 and PfCRT^Dd2. Our best-performing model incorporated a binary classifier, a nonlinear scale, and additive effects for each mutation. The addition of specific pairwise- and high-order-epistatic coefficients decreased the predictive power of the model. We evaluated our predictions by experimentally measuring the CQ transport activities of 24 additional PfCRT genotypes. The R² value between our predicted and newly-measured phenotypes was 0.90. We then used the model to probe the accessibility of evolutionary trajectories through the map. Approximately 1% of the possible trajectories between PfCRT^3D7 and PfCRT^Dd2 are accessible; however, none of the trajectories entailed eight successive increases in CQ transport activity. These results demonstrate that phenotypes can be inferred with known uncertainty from a partial genotype-phenotype dataset. We also validated our approach against a collection of previously published genotype-phenotype maps. The model therefore appears general and should be applicable to a large number of genotype-phenotype maps.

Biological macromolecules are built from chains of building blocks. The function of a macromolecule depends on the specific chemical properties of the building blocks that make it up. Macromolecules evolve through mutations that swap one building block for another. Understanding how biomolecules work and evolve therefore requires knowledge of the effects of mutations. The effects of mutations can be measured experimentally; however, because there are a vast number of possible combinations of mutations, it is often difficult to make enough measurements to understand biomolecular function and evolution. In this paper, we describe a simple method to predict the effects of mutations on biomolecules from a small number of measurements. This method works by appropriately averaging the effects of mutations seen in different contexts. We test the method by predicting the effects of mutations on a PfCRT—a macromolecule from the malarial parasite that confers drug resistance. We find that our method is fast and effective. Using a small number of measurements, we were able to gain insight into the evolutionary steps by which this macromolecule conferred drug resistance. To make this method accessible to other researchers, we have released it as an open-source software package: https://gpseer.readthedocs.io.

Evolution of multifunctionality through a pleiotropic substitution in the innate immune protein S100A9

Multifunctional proteins are evolutionary puzzles: how do proteins evolve to satisfy multiple functional constraints? S100A9 is one such multifunctional protein. It potently amplifies inflammation via Toll-like receptor four and is antimicrobial as part of a heterocomplex with S100A8. These two functions are seemingly regulated by proteolysis: S100A9 is readily degraded, while S100A8/S100A9 is resistant. We take an evolutionary biochemical approach to show that S100A9 evolved both functions and lost proteolytic resistance from a weakly proinflammatory, proteolytically resistant amniote ancestor. We identify a historical substitution that has pleiotropic effects on S100A9 proinflammatory activity and proteolytic resistance but has little effect on S100A8/S100A9 antimicrobial activity. We thus propose that mammals evolved S100A8/S100A9 antimicrobial and S100A9 proinflammatory activities concomitantly with a proteolytic 'timer' to selectively regulate S100A9. This highlights how the same mutation can have pleiotropic effects on one functional state of a protein but not another, thus facilitating the evolution of multifunctionality.

Tracing the evolution of novel features of human Toll-like receptor 4

Toll-like receptor 4 (TLR4) is a critical innate immune protein that activates inflammation in response to extracellular cues. Much of the work to understand how the protein works in humans has been done using mouse models. Although human and mouse TLR4 have many shared features, they have also diverged significantly since their last common ancestor, acquiring 277 sequence differences. Functional differences include the extent of ligand-independent activation, whether lipid IVa acts as an antagonist or agonist, and the relative species cross-compatibility of their MD-2 cofactor. We set out to understand the evolutionary origins for these functional differences between human and mouse TLR4. Using a combination of phylogenetics, ancestral sequence reconstruction, and functional characterization, we found that evolutionary changes to the human TLR4, rather than changes to the mouse TLR4, were largely responsible for these functional changes. Human TLR4 repressed ancestral ligand-independent activity and gained antagonism to lipid IVa. Additionally, mutations to the human TLR4 cofactor MD-2 led to lineage-specific incompatibility between human and opossum TLR4 complex members. These results were surprising, as mouse TLR4 has acquired many more mutations than human TLR4 since their last common ancestor. Our work has polarized this set of transitions and sets up work to study the mechanistic underpinnings for the evolution of new functions in TLR4.

pytc: Open-Source Python Software for Global Analyses of Isothermal Titration Calorimetry Data

Here we describe pytc, an open-source Python package for global fits of thermodynamic models to multiple isothermal titration calorimetry experiments. Key features include simplicity, the ability to implement new thermodynamic models, a robust maximum likelihood fitter, a fast Bayesian Markov-Chain Monte Carlo sampler, rigorous implementation, extensive documentation, and full cross-platform compatibility. pytc fitting can be done using an application program interface or via a graphical user interface. It is available for download at https://github.com/harmslab/pytc .

Coevolution of the Toll-Like Receptor 4 Complex with Calgranulins and Lipopolysaccharide

Toll-like receptor 4 (TLR4) induces inflammation in response to both pathogen- and host-derived molecules. Lipopolysaccharide (LPS) recognition by TLR4 has been shown to occur across the amniotes, but endogenous signaling through TLR4 has not been validated outside of placental mammals. To determine whether endogenous danger signaling is also shared across amniotes, we studied the evolution of TLR4-activation by the calgranulin proteins (S100A8, S100A9, and S100A12), a clade of host molecules that potently activate TLR4 in placental mammals. We performed phylogenetic and syntenic analysis and found MRP-126—a gene in birds and reptiles—is likely orthologous to the mammalian calgranulins. We then used an ex vivo TLR4 activation assay to establish that calgranulin pro-inflammatory activity is not specific to placental mammals, but is also exhibited by representative marsupial and sauropsid species. This activity is strongly dependent on the cofactors CD14 and MD-2 for all species studied, suggesting a conserved mode of activation across the amniotes. Ortholog complementation experiments between the calgranulins, TLR4, CD14, and MD-2 revealed extensive lineage specific-coevolution and multi-way interactions between components that are necessary for the activation of NF-κB signaling by calgranulins and LPS. Our work demonstrates that calgranulin activation of TLR4 evolved at least ~320 million years ago and has been conserved in the amniote innate immune system.

Conservation of Specificity in Two Low-Specificity Proteins

Many regulatory proteins bind peptide regions of target proteins and modulate their activity. Such regulatory proteins can often interact with highly diverse target peptides. In many instances, it is not known if the peptide-binding interface discriminates targets in a biological context, or whether biological specificity is achieved exclusively through external factors such as subcellular localization. We used an evolutionary biochemical approach to distinguish these possibilities for two such low-specificity proteins: S100A5 and S100A6. We used isothermal titration calorimetry to study the binding of peptides with diverse sequence and biochemistry to human S100A5 and S100A6. These proteins bound distinct, but overlapping, sets of peptide targets. We then studied the peptide binding properties of orthologs sampled from across five amniote species. Binding specificity was conserved along all lineages, for the last 320 million years, despite the low specificity of each protein. We used ancestral sequence reconstruction to determine the binding specificity of the last common ancestor of the paralogs. The ancestor bound the entire set of peptides bound by modern S100A5 and S100A6 proteins, suggesting that paralog specificity evolved via subfunctionalization. To rule out the possibility that specificity is conserved because it is difficult to modify, we identified a single historical mutation that, when reverted in human S100A5, gave it the ability to bind an S100A6-specific peptide. These results reveal strong evolutionary constraints on peptide binding specificity. Despite being able to bind a large number of targets, the specificity of S100 peptide interfaces is likely important for the biology of these proteins.

Human S100A5 binds Ca2+ and Cu2+ independently

Background

S100A5 is a calcium binding protein found in a small subset of amniote tissues. Little is known about the biological roles of S100A5, but it may be involved in inflammation and olfactory signaling. Previous work indicated that S100A5 displays antagonism between binding of Ca²⁺ and Cu²⁺ ions-one of the most commonly cited features of the protein. We set out to characterize the interplay between Ca²⁺ and Cu²⁺ binding by S100A5 using isothermal titration calorimetry (ITC), circular dichroism spectroscopy (CD), and analytical ultracentrifugation (AUC).

Results

We found that human S100A5 is capable of binding both Cu²⁺ and Ca²⁺ ions simultaneously. The wildtype protein was extremely aggregation-prone in the presence of Cu²⁺ and Ca²⁺. A Cys-free version of S100A5, however, was not prone to precipitation or oligomerization. Mutation of the cysteines does not disrupt the binding of either Ca²⁺ or Cu²⁺ to S100A5. In the Cys-free background, we measured Ca²⁺ and Cu²⁺ binding in the presence and absence of the other metal using ITC. Saturating concentrations of Ca²⁺ or Cu²⁺ do not disrupt the binding of one another. Ca²⁺ and Cu²⁺ binding induce structural changes in S100A5, which are measurable using CD spectroscopy. We show via sedimentation velocity AUC that the wildtype protein is prone to the formation of soluble oligomers, which are not present in Cys-free samples.

Conclusions

S100A5 can bind Ca²⁺ and Cu²⁺ ions simultaneously and independently. This observation is in direct contrast to previously-reported antagonism between binding of Cu²⁺ and Ca²⁺ ions. The previous result is likely due to metal-dependent aggregation. Little is known about the biology of S100A5, so an accurate understanding of the biochemistry is necessary to make informed biological hypotheses. Our observations suggest the possibility of independent biological functions for Cu²⁺ and Ca²⁺ binding by S100A5.

High-order epistasis shapes evolutionary trajectories

High-order epistasis—where the effect of a mutation is determined by interactions with two or more other mutations—makes small, but detectable, contributions to genotype-fitness maps. While epistasis between pairs of mutations is known to be an important determinant of evolutionary trajectories, the evolutionary consequences of high-order epistasis remain poorly understood. To determine the effect of high-order epistasis on evolutionary trajectories, we computationally removed high-order epistasis from experimental genotype-fitness maps containing all binary combinations of five mutations. We then compared trajectories through maps both with and without high-order epistasis. We found that high-order epistasis strongly shapes the accessibility and probability of evolutionary trajectories. A closer analysis revealed that the magnitude of epistasis, not its order, predicts is effects on evolutionary trajectories. We further find that high-order epistasis makes it impossible to predict evolutionary trajectories from the individual and paired effects of mutations. We therefore conclude that high-order epistasis profoundly shapes evolutionary trajectories through genotype-fitness maps.

A key goal for evolutionary biologists is understanding why one evolutionary trajectory is taken rather than others. This requires understanding how individual mutations, as well as interactions between them, determine the accessibility of evolutionary pathways. We used a robust statistical analysis to reveal interactions between up to five mutations in published datasets, meaning that the effect of a mutation can depend on the presence or absence of four other mutations. Simulations reveal that these interactions strongly shape evolutionary trajectories. These interactions lead to profound unpredictability in evolution, as one cannot use the effect of a mutation in the ancestor to predict its effect later in the trajectory.

Robustness of Reconstructed Ancestral Protein Functions to Statistical Uncertainty

Hypotheses about the functions of ancient proteins and the effects of historical mutations on them are often tested using ancestral protein reconstruction (APR)-phylogenetic inference of ancestral sequences followed by synthesis and experimental characterization. Usually, some sequence sites are ambiguously reconstructed, with two or more statistically plausible states. The extent to which the inferred functions and mutational effects are robust to uncertainty about the ancestral sequence has not been studied systematically. To address this issue, we reconstructed ancestral proteins in three domain families that have different functions, architectures, and degrees of uncertainty; we then experimentally characterized the functional robustness of these proteins when uncertainty was incorporated using several approaches, including sampling amino acid states from the posterior distribution at each site and incorporating the alternative amino acid state at every ambiguous site in the sequence into a single "worst plausible case" protein. In every case, qualitative conclusions about the ancestral proteins' functions and the effects of key historical mutations were robust to sequence uncertainty, with similar functions observed even when scores of alternate amino acids were incorporated. There was some variation in quantitative descriptors of function among plausible sequences, suggesting that experimentally characterizing robustness is particularly important when quantitative estimates of ancient biochemical parameters are desired. The worst plausible case method appears to provide an efficient strategy for characterizing the functional robustness of ancestral proteins to large amounts of sequence uncertainty. Sampling from the posterior distribution sometimes produced artifactually nonfunctional proteins for sequences reconstructed with substantial ambiguity.

Multiple Evolutionary Origins of Ubiquitous Cu2+ and Zn2+ Binding in the S100 Protein Family

The S100 proteins are a large family of signaling proteins that play critical roles in biology and disease. Many S100 proteins bind Zn2+, Cu2+, and/or Mn2+ as part of their biological functions; however, the evolutionary origins of binding remain obscure. One key question is whether divalent transition metal binding is ancestral, or instead arose independently on multiple lineages. To tackle this question, we combined phylogenetics with biophysical characterization of modern S100 proteins. We demonstrate an earlier origin for established S100 subfamilies than previously believed, and reveal that transition metal binding is widely distributed across the tree. Using isothermal titration calorimetry, we found that Cu2+ and Zn2+ binding are common features of the family: the full breadth of human S100 paralogs-as well as two early-branching S100 proteins found in the tunicate Oikopleura dioica-bind these metals with μM affinity and stoichiometries ranging from 1:1 to 3:1 (metal:protein). While binding is consistent across the tree, structural responses to binding are quite variable. Further, mutational analysis and structural modeling revealed that transition metal binding occurs at different sites in different S100 proteins. This is consistent with multiple origins of transition metal binding over the evolution of this protein family. Our work reveals an evolutionary pattern in which the overall phenotype of binding is a constant feature of S100 proteins, even while the site and mechanism of binding is evolutionarily labile.

The thermostability and specificity of ancient proteins

Were ancient proteins systematically different than modern proteins? The answer to this question is profoundly important, shaping how we understand the origins of protein biochemical, biophysical, and functional properties. Ancestral sequence reconstruction (ASR), a phylogenetic approach to infer the sequences of ancestral proteins, may reveal such trends. We discuss two proposed trends: a transition from higher to lower thermostability and a tendency for proteins to acquire higher specificity over time. We review the evidence for elevated ancestral thermostability and discuss its possible origins in a changing environmental temperature and/or reconstruction bias. We also conclude that there is, as yet, insufficient data to support a trend from promiscuity to specificity. Finally, we propose future work to understand these proposed evolutionary trends.

Historical contingency and its biophysical basis in glucocorticoid receptor evolution

Understanding how chance historical events shape evolutionary processes is a central goal of evolutionary biology^1–7. Direct insights into the extent and causes of evolutionary contingency have been limited to experimental systems,^7–9 because it is difficult to know what happened in the deep past and to characterize other paths that evolution could have followed. Here we combine ancestral protein reconstruction, directed evolution, and biophysical analysis to explore alternate “might-have-been” trajectories during the ancient evolution of a novel protein function. We previously found that the evolution of cortisol specificity in the ancestral glucocorticoid receptor (GR) was contingent on permissive substitutions, which had no apparent effect on receptor function but were necessary for GR to tolerate the large-effect mutations that caused the shift in specificity.⁶ Here we show that alternative mutations that could have permitted the historical function-switching substitutions are extremely rare in the ensemble of genotypes accessible to the ancestral GR. In a library of thousands of variants of the ancestral protein, we recovered historical permissive substitutions, but no alternate permissive genotypes. Using biophysical analysis, we found that permissive mutations must satisfy at least three physical requirements—they must stabilize specific local elements of the protein structure, maintain the correct energetic balance between functional conformations, and be compatible with the ancestral and derived structures—thus revealing why permissive mutations are rare. These findings demonstrate that GR evolution depended strongly on improbable, nondeterministic events, and this contingency arose from intrinsic biophysical properties of the protein.

Evolutionary biochemistry: revealing the historical and physical causes of protein properties

The repertoire of proteins and nucleic acids in the living world is determined by evolution; their properties are determined by the laws of physics and chemistry. Explanations of these two kinds of causality - the purviews of evolutionary biology and biochemistry, respectively - are typically pursued in isolation, but many fundamental questions fall squarely at the interface of fields. Here we articulate the paradigm of evolutionary biochemistry, which aims to dissect the physical mechanisms and evolutionary processes by which biological molecules diversified and to reveal how their physical architecture facilitates and constrains their evolution. We show how an integration of evolution with biochemistry moves us towards a more complete understanding of why biological molecules have the properties that they do.

Resurrection of an Urbilaterian U1A/U2B″/SNF protein

The U1A/U2B″/SNF family of proteins found in the U1 and U2 spliceosomal small nuclear ribonucleoproteins is highly conserved. In spite of the high degree of sequence and structural conservation, modern members of this protein family have unique RNA binding properties. These differences have necessarily resulted from evolutionary processes, and therefore, we reconstructed the protein phylogeny in order to understand how and when divergence occurred and how protein function has been modulated. Contrary to the conventional understanding of an ancient human U1A/U2B″ gene duplication, we show that the last common ancestor of bilaterians contained a single ancestral protein (URB). The gene for URB was synthesized, the protein was overexpressed and purified, and we assessed RNA binding to modern snRNA sequences. We find that URB binds human and Drosophila U1 snRNA SLII and U2 snRNA SLIV with higher affinity than do modern homologs, suggesting that both Drosophila SNF and human U1A/U2B″ have evolved into weaker binders of one RNA or both RNAs.

Evolution of minimal specificity and promiscuity in steroid hormone receptors

Most proteins are regulated by physical interactions with other molecules; some are highly specific, but others interact with many partners. Despite much speculation, we know little about how and why specificity/promiscuity evolves in natural proteins. It is widely assumed that specific proteins evolved from more promiscuous ancient forms and that most proteins' specificity has been tuned to an optimal state by selection. Here we use ancestral protein reconstruction to trace the evolutionary history of ligand recognition in the steroid hormone receptors (SRs), a family of hormone-regulated animal transcription factors. We resurrected the deepest ancestral proteins in the SR family and characterized the structure-activity relationships by which they distinguished among ligands. We found that that the most ancient split in SR evolution involved a discrete switch from an ancient receptor for aromatized estrogens—including xenobiotics—to a derived receptor that recognized non-aromatized progestagens and corticosteroids. The family's history, viewed in relation to the evolution of their ligands, suggests that SRs evolved according to a principle of minimal specificity: at each point in time, receptors evolved ligand recognition criteria that were just specific enough to parse the set of endogenous substances to which they were exposed. By studying the atomic structures of resurrected SR proteins, we found that their promiscuity evolved because the ancestral binding cavity was larger than the primary ligand and contained excess hydrogen bonding capacity, allowing adventitious recognition of larger molecules with additional functional groups. Our findings provide an historical explanation for the sensitivity of modern SRs to natural and synthetic ligands—including endocrine-disrupting drugs and pollutants—and show that knowledge of history can contribute to ligand prediction. They suggest that SR promiscuity may reflect the limited power of selection within real biological systems to discriminate between perfect and “good enough.”

The functions of most proteins are defined by their interactions with other biological substances, such as DNA, nutrients, hormones, or other proteins. Some proteins are highly specific, but others are more promiscuous and can interact with a variety of natural substances, as well as drugs and pollutants. Understanding molecular interactions is a key goal in pharmacology and toxicology, but there are few general principles to help explain or predict protein specificity. Because every biological entity is the result of evolution, understanding a protein's history might help explain why it interacts with the substances to which it is sensitive. In this paper, we used ancestral protein reconstruction to experimentally trace how specificity evolved in an ancient group of proteins, the steroid hormone receptors (SRs), a family of proteins that regulate reproduction and other biological processes in animals. We show that SRs evolved according to a principle of minimal specificity: at each point in time, these proteins evolved to be specific enough to distinguish among the substances to which they were naturally exposed, but not more so. Our findings provide an historical explanation for modern SRs' diverse sensitivities to natural and man-made substances; they show that knowledge of history can contribute to predicting the ligands to which a modern protein will respond and indicate that promiscuity reflects the limited power of natural selection to discriminate between perfect and “good enough.”

Analyzing protein structure and function using ancestral gene reconstruction

Protein families with functionally diverse members can illuminate the structural determinants of protein function and the process by which protein structure and function evolve. To identify the key amino acid changes that differentiate one family member from another, most studies have taken a horizontal approach, swapping candidate residues between present-day family members. This approach has often been stymied, however, by the fact that shifts in function often require multiple interacting mutations; chimeric proteins are often nonfunctional, either because one lineage has amassed mutations that are incompatible with key residues that conferred a new function on other lineages, or because it lacks mutations required to support those key residues. These difficulties can be overcome by using a vertical strategy, which reconstructs ancestral genes and uses them as the appropriate background in which to study the effects of historical mutations on functional diversification. In this review, we discuss the advantages of the vertical strategy and highlight several exemplary studies that have used ancestral gene reconstruction to reveal the molecular underpinnings of protein structure, function, and evolution.

The pK(a) values of acidic and basic residues buried at the same internal location in a protein are governed by different factors

The pK(a) values of internal ionizable groups are usually very different from the normal pK(a) values of ionizable groups in water. To examine the molecular determinants of pK(a) values of internal groups, we compared the properties of Lys, Asp, and Glu at internal position 38 in staphylococcal nuclease. Lys38 titrates with a normal or elevated pK(a), whereas Asp38 and Glu38 titrate with elevated pK(a) values of 7.0 and 7.2, respectively. In the structure of the L38K variant, the buried amino group of the Lys38 side chain makes an ion pair with Glu122, whereas in the structure of the L38E variant, the buried carboxyl group of Glu38 interacts with two backbone amides and has several nearby carboxyl oxygen atoms. Previously, we showed that the pK(a) of Lys38 is normal owing to structural reorganization and water penetration concomitant with ionization of the Lys side chain. In contrast, the pK(a) values of Asp38 and Glu38 are perturbed significantly owing to an imbalance between favorable polar interactions and unfavorable contributions from dehydration and from Coulomb interactions with surface carboxylic groups. Their ionization is also coupled to subtle structural reorganization. These results illustrate the complex interplay between local polarity, Coulomb interactions, and structural reorganization as determinants of pK(a) values of internal groups in proteins. This study suggests that improvements to computational methods for pK(a) calculations will require explicit treatment of the conformational reorganization that can occur when internal groups ionize.

A buried lysine that titrates with a normal pKa: role of conformational flexibility at the protein-water interface as a determinant of pKa values

Previously we reported that Lys, Asp, and Glu residues at positions 66 and 92 in staphylococcal nuclease (SNase) titrate with pK(a) values shifted by up to 5 pK(a) units in the direction that promotes the neutral state. In contrast, the internal Lys-38 in SNase titrates with a normal pK(a). The crystal structure of the L38K variant shows that the side chain of Lys-38 is buried. The ionizable moiety is approximately 7 A from solvent and ion paired with Glu-122. This suggests that the pK(a) value of Lys-38 is normal because the energetic penalty for dehydration is offset by a favorable Coulomb interaction. However, the pK(a) of Lys-38 was also normal when Glu-122 was replaced with Gln or with Ala. Continuum electrostatics calculations were unable to reproduce the pK(a) of Lys-38 unless the protein was treated with an artificially high dielectric constant, consistent with structural reorganization being responsible for the normal pK(a) value of Lys-38. This reorganization must be local because circular dichroism and NMR spectroscopy indicate that the L38K protein is native-like under all conditions studied. In molecular dynamics simulations, the ion pair between Lys-38 and Glu-122 is unstable. The simulations show that a minor rearrangement of a loop is sufficient to allow penetration of water to the amino moiety of Lys-38. This illustrates both the important roles of local flexibility and water penetration as determinants of pK(a) values of ionizable groups buried near the protein-water interface, and the challenges faced by structure-based pK(a) calculations in reproducing these effects.

Laser light-scattering evidence for an altered association of beta B1-crystallin deamidated in the connecting peptide

Deamidation is a prevalent modification of crystallin proteins in the vertebrate lens. The effect of specific sites of deamidation on crystallin stability in vivo is not known. Using mass spectrometry, a previously unreported deamidation in beta B1-crystallin was identified at Gln146. Another deamidation was investigated at Asn157. It was determined that whole soluble beta B1 contained 13%-17% deamidation at Gln146 and Asn157. Static and quasi-elastic laser light scattering, circular dichroism, and heat aggregation studies were used to explore the structure and associative properties of recombinantly expressed wild-type (wt) beta B1 and the deamidated beta B1 mutants, Q146E and N157D. Dimer formation occurred for wt beta B1, Q146E, and N157D in a concentration-dependent manner, but only Q146E showed formation of higher ordered oligomers at the concentrations studied. Deamidation at Gln146, but not Asn157, led to an increased tendency of beta B1 to aggregate upon heating. We conclude that deamidation creates unique effects depending upon where the deamidation is introduced in the crystallin structure.