Discrimination of thermophilic and mesophilic proteins
© Taylor and Vaisman; licensee BioMed Central Ltd. 2010
Published: 17 May 2010
Skip to main content
Volume 10 Supplement 1
© Taylor and Vaisman; licensee BioMed Central Ltd. 2010
Published: 17 May 2010
There is a considerable literature on the source of the thermostability of proteins from thermophilic organisms. Understanding the mechanisms for this thermostability would provide insights into proteins generally and permit the design of synthetic hyperstable biocatalysts.
We have systematically tested a large number of sequence and structure derived quantities for their ability to discriminate thermostable proteins from their non-thermostable orthologs using sets of mesophile-thermophile ortholog pairs. Most of the quantities tested correspond to properties previously reported to be associated with thermostability. Many of the structure related properties were derived from the Delaunay tessellation of protein structures.
Carefully selected sequence based indices discriminate better than purely structure based indices. Combined sequence and structure based indices improve performance somewhat further. Based on our analysis, the strongest contributors to thermostability are an increase in ion pairs on the protein surface and a more strongly hydrophobic interior.
Organisms that thrive at very high temperatures have been actively studied since the discovery of Thermus aquaticus in the hot springs of Yellowstone in the 1960’s . Heat tolerant organisms are often separated into two classes: thermophiles, which have optimum growth temperatures (OGT) in the range 45-80 °C, and hyperthermophiles with OGTs above 80 °C. Mesophilic organisms are defined as those with OGT’s between 15 °C and 45 °C, while psychrophiles, which we do not address here, have OGT’s no greater than 15 °C. Sometimes the break points between these classes are assigned slightly differently. Hyperthermophiles come mostly from the kingdom Archea, but there are two genera of hyperthermophilic Eubacteria, namely Thermotogales and Aquifex. Thermophiles are more phylogenetically diverse and include Eubacteria, Archea, and some fungi.
In addition to providing insights into the principles of protein folding and stability, understanding what makes some proteins more thermostable than others is of practical interest. Thermophilic proteins are more resistant to proteolysis and chemical denaturation, hence there is interest in engineering hyperstable biocatalysts relying on the same mechanisms that nature uses [2–4]. Thermophilic polymerases, proteases, and xylanases already have industrial applications [4, 5].
The search for the physical basis of thermostability in proteins goes back 30 years to the work of Perutz . Since then, a great many papers have been written on the subject. Some of the proposed mechanisms/indicators of increased thermostability include: a more highly hydrophobic core [7, 8], tighter packing or compactness , deleted or shortened loops [10, 11], greater rigidity [3, 12, 13] (for example through increased Proline content in loops), higher secondary structure content , greater polar surface area , fewer and/or smaller voids [14, 16], smaller surface area to volume ratio , fewer thermolabile residues [16, 18], increased hydrogen bonding , higher isoelectric point , and more salt bridges/ion pairs and networks of salt bridges [6, 20–25].
Statistically significant changes in sequence composition between mesophilic and thermophilic proteins have been reported. The amino acids Asn, Gln, Met, and Cys are thermolabile—they are not stable at high temperatures and tend to undergo deamidation (Asn and Gln) or oxidation (Met and Cys) . These amino acids are less common in thermophilic proteins and the thermolabile residues that do occur are usually buried . Ile is preferred to Leu in hydrophobic regions of the structure because the side chain carbons can exist in all three χ rotameric states compared to only two for Leu which can result in tighter side chain packing . Farias and Bonato  have reported that Gly, Lys, Tyr, and Ile are preferred in thermophilic organisms while Gln, His, Ala, and Cys are preferred in mesophiles. Camillau and Claverie  have reported that thermophilic proteins have less Gln, Ala, and His on their surfaces than mesophilic proteins do and more charged residues on their surfaces, particularly Lys and Glu. Haney et al.  have compared 115 proteins from Methanococcus jannaschii to mesophilic proteins from other Methanococcus species and found that the frequencies of Ile, Glu, Arg, Lys, Pro, and Tyr are significantly greater in the thermophile and the frequencies of Gly, Met, Gln, Thr, Asn, and Ser are smaller.
More ion pairs have been strongly and consistently linked with thermostability in the literature. Water has a dielectric constant of about 80 at 0°C, which drops to 55 at 100°C and is lower still at the extreme pressures near hydrothermal vents in the deep sea where some hyperthermophilic organisms live. A lower dielectric constant makes electrostatic interactions stronger and therefore ion pairs should have a greater stabilizing effect at high temperatures and pressures [21, 29].
Evidence for some of these proposed mechanisms/indicators is equivocal. For instance, Karshikoff and Ladenstein found no significant difference in packing between thermophilic and mesophilic proteins  and salt bridges in a protein core have been reported to be destabilizing [31, 32]. Das and Gerstein  have reported that the lengths of proteins from the eubacterium Aquifex aeolicus are greater than those of archeal hyperthermophilic orthologs and therefore hyperthermophilic proteins may be shorter than their mesophilic counterparts simply because most hyperthermophiles are archeal, not necessarily because shorter loops promote enhanced thermostability. Querol et al.  surveyed 122 references for 195 single point mutants which have been unambiguously linked to greater thermostability and found that greater rigidity, as measured by crystallographic B-factors, is not a good indicator of thermostability.
Liang et al. have studied the proteomes of 15 thermophiles and 74 mesophiles using the tendencies of residue pairs separated by no more than 20 in primary sequence to occur together to discriminate mesophilic from thermophilic proteomes . Farias and Bonato  have devised a sequence composition based index capable of correctly classifying organisms. The index ri is characteristic of a single protein and is defined as ri=(E+K)/(Q+H) where E, K, Q, and H are the percent compositions of these amino acid types in protein i. Those authors took the average of ri over all the proteins in an organism to give an average r that, without exception, fell in different ranges for the mesophiles (r < 2.5), thermophiles (3.2 < r < 4.6), and hyperthermophiles (r > 4.5) in their test set. Further, they showed that r is high in chaperonins (heat shock proteins) in both mesophiles and thermophiles  thereby concluding that the sequence signal is indicative of thermostability and not a phylogenetic artifact. Similarly, Claverie and colleagues [27, 37] have devised the CvP bias (charged versus polar), defined as (D+E+K+R)-(N+Q+S+T), the sum of amino acid compositions, that they have also used to classify an organism as mesophilic, thermophilic, or hyperthermophilic by computing the average CvP for all the proteins from that organism. Zeldovich et al.  have reported a sequence composition based index defined as I+V+Y+W+R+E+L (abbreviated IVYWREL) that was an extremely good predictor of thermostability when averaged over whole proteomes, and even for just the membrane proteins from these proteomes. IVYWREL, again when averaged over proteomes, also correlates very well with OGT. None of these authors claim that their indices work well at discriminating individual pairs of thermophilic and mesophilic orthologs, however it is natural to ask if they can, and we will test this question here. Glyakina et al. have done one of the few large scale structure based discrimination studies .
A Delaunay tessellation derived four body statistical contact pseudo-potential has been reported previously [42, 43] which has been shown to contain more information than pairwise contact potentials . Under this pseudo-potential, the score of some particular amino acid quadruplet (i, j, k, l), which corresponds to the residues at the vertices of a Delaunay simplex, is defined as:
where: fijkl is the observed frequency of simplices with amino acid types i, j, k, and l at their vertices in a large non-redundant training set S; ai, aj, ak, and al are the observed frequencies of the individual amino acid types in S; and c is a combinatorial factor. Variations of this potential have been successfully applied to fold recognition [43, 50] and the analysis of protein stability  and activity . We will denote by Q the sum of the log-likelihoods qijkl from the residue quadruplets corresponding to all Delaunay simplices in the tessellation of a protein structure.
The radius of the sphere circumscribing a Delaunay tetrahedron gives a measure of its eccentricity. The small, nearly equilateral tetrahedra in the interior of the tessellation have small circumsphere radii, on the order of the size of the simplices. Radically skewed, nearly coplanar tetrahedra on the surface of the tessellated protein, on the other hand, can have circumsphere radii orders of magnitude larger than the diameter of the molecule. The tetrahedrality T is another measure of simplex eccentricity. Denote the length of the six edges of a simplex as L1 - L6 . The tetrahedrality is then defined by:
A few protein structure descriptors we will use are based on molecule contact graphs, so a brief introduction to this is also in order. The residues in contact with one another in a protein can be thought of as a graph or network (Fig. 2) and analyzed using techniques from elementary graph theory and the theory of complex networks. In the literature, residues are typically defined to be in contact by a simple proximity cutoff, but in this work, graph nodes correspond to residues and graph edges join nodes when the corresponding residues are joined by a Delaunay simplex, shorter than some fixed cutoff, in the tessellation of the protein structure.
Several contact network derived quantities have been used before to analyze protein structures [61, 62]. The degree k of a node in an undirected graph is the number of edges impinging on it. The average degree over all nodes in the contact graph will be referred to as the coordination number. A minimum path between nodes i and j is one for which the sum of weights of the edges along the path is smallest from among all possible paths. The minimum path length Lij between nodes i and j is the sum of the weights along a minimum path. In our case here, edges have weight one, and a minimum path is one for which the fewest edges are traversed. The characteristic path length L of a network is the average of the minimum paths between all node pairs i, j where i ≠ j. In general, there will be many paths between distinct nodes i and j that have the minimum path length. Some classes of networks have the clustering property, which means that two nodes which are both joined by edges to a third, are more likely to also be joined to each other than are two nodes picked at random . In such networks, there are well defined neighborhoods with subsets of nodes tending to be connected to each other and tending not to be connected to nodes in other neighborhoods. The clustering property is measured by the clustering coefficient of a node Cn is the number of actual edges En between neighbors of node n divided by the number of possible connections between those neighbors: Cn = 2En/(k(k-1)), where k is the degree of node n. The clustering coefficient C for the entire network is the average of all the Cn.
Linear least square best fit lines of number of residues in hyperthermophilic (Nh) or thermophilic proteins (Nt) to the number of residues in their mesophilic counterparts (Nm) are: Nh=0.94 Nm + 5.00 and Nt=0.97 Nm + 5.15 where there are 122 pairs in the hyperthermophile pairs and 127 in thermophile pairs (see Additional files 1 and 2 for lists). Our data show that thermophilic proteins are usually somewhat shorter than their mesophilic counterparts, and hyperthermophilic proteins are shorter still. This observation is in line with the results of Eisenberg .
t-tests of sequence percent compositions of mesophilc, thermophilic, and hyperthermophilic proteins.
Statistical significance of fractions in Table 1.
P for thermophile pairs
P for hyperthermophile pairs
6.85 x 10-4
6.97 x 10-6
1.65 x 10-5
1.94 x 10-8
4.99 x 10-8
3.12 x 10-12
8.34 x 10-12
2.66 x 10-16
6.78 x 10-16
<2.20 x 10-16
<2.20 x 10-16
<2.20 x 10-16
<2.20 x 10-16
Discriminatory power of structure and sequence derived quantities
Thermophile (127 pairs)
Hyperthermophile (122 pairs)
Contact Network Derived Quantities
coordination number (no cutoff)
clustering coefficient (no cutoff)
characteristic path (no cutoff)
Combined Sequence and Structure Including Threading Potentials
total count 400 over-rep quads/residue
4-body potential/residue (20Å cutoff)
4-body potential/residue (no cutoff)
4-body potential/res (hyper only, no cutoff)
4-body potential/res (meso only, no cutoff)
4-body potential/res (thermo only,no cutoff)
ProsaII combined score
Delaunay Simplex Geometry
median circumsphere radius(no cutoff)
mean tetrahedrality (no cutoff)
number simplices/residue (10Å cutoff)
number simplices/residue (no cutoff)
Naccess solvent accessible area
Delaunay surface area (no cutoff)
van der Waals area
Delaunay volume (no cutoff)
Van der Waals volume
Delaunay area/volume (10Å cutoff)
Delaunay area/volume (no cutoff)
van der Waals area/volume
van der Waals volume/Delaunay volume
secondary structure content (H+E 3 state DSSP)
number of residues
total Kyte-Doolittle hydrophobicity
sd Kyte-Doolittle hydrophobicity
The best discriminatory indices we tested were: one version of the 4-body Delaunay threading potential; the count of over-represented quadruplets; the ratio of Delaunay surface area to volume (for hyperthermophiles); the standard deviation of Kyte-Doolittle hydrophobicity (for hyperthermophiles); Delaunay area/volume=and van der Waals volume/Delaunay volume, particularly for hyperthermophiles; the CvP bias, and IVYWREL. Few of the tested indices (when computed for individual proteins not averaged over proteomes) correlate even moderately well with OGT. Those for which the correlation is strongest are the 4-body potentials, overrep400, CvP, IVYWREL, and E+K/Q+H which all have r~0.4-0.7.
The Delaunay simplex-based descriptors (overrep400 and the 4-body potentials) that work best for discrimination use large simplex edge length cutoffs ( >20 Å). This implies that there are important residue contacts on the surface of the proteins because that is invariably where simplices with very long edge lengths reside. This combined with threading potential data presented later leads us to believe that the presence of more charged residues on the protein surface is at least one of the things these descriptors pick up.
Since increased hydrophobicity of the protein core has been proposed as a mechanism for thermostability, one might expect the sum SKD of the Kyte-Doolittle hydrophobicities of all residues would be a good discriminator. SKD is not, but the variance of the Kyte-Doolittle hydrophobicities is (Table 3). Apparently, then, the increase in core hydrophobicity is accompanied by an increase in hydrophilic residues.
We have studied the mean contribution to the 4-body potential score Q under a 7-letter reduced alphabet to each structure in the mesophilic, thermophilic, and hyperthermophilic subsets of 521nonredundant. The reduced alphabet was used to simplify the analysis by bringing down the possible number of residue quadruplets that can reside at the vertices of a Delaunay tetrahedron from 8855 for a 20 amino acid alphabet to 210. The reduced alphabet is D, E, K, R, I=(I,V), A=(A,F,G,L,P), N=(C,H,M,N,Q,S,T,W,Y. The four-body Delaunay threading score Q for a hyperthermophile with the average quadruplet composition is 36.3, for a thermophile Q is 28.2, and for a mesophile Q is 17.6. For hyperthermophiles, the ΔQ (with respect to mesophiles) attributable to quadruplets with at least two hydrophobic residues is +11.07, for quadruplets with at least two charged residues it is +12.78, and for quadruplets with at least two polar residues it is –2.53. For thermophilic quadruplets the ΔQ figures are: at least two hydrophobic residues +5.99, at least two charged residues +4.55, and at least two polar residues –0.34. We see ,therefore, that the 4-body potential goes up for (hyper)thermophiles both due to associations between charged residues but also due to quadruplets of more highly hydrophobic residues. The increase in charged residues produces a larger change than the stronger hydrophobics with hyperthermophiles, but the situation is reversed in thermophiles.
It is possible to accurately discriminate (hyper) thermophilic proteins from their mesophilic counterparts based on sequence and structural properties. Sequence based indices used to discriminate entire proteomes also work well on individual thermophile/mesophile ortholog pairs. Purely structure-based indices are, generally speaking, poor discriminators. Combined sequence structure indices like the threading potential are somewhat better than sequence alone.
The primary factors differentiating thermophilic from mesophilic proteins according to our analysis are surface ion pairs and more strongly hydrophobic core residues. The conclusion of previous authors that the basis for the thermostability of thermophiles and hyper-thermophiles is somewhat different is also borne out here (e.g. the preference of thermophiles for Arg and of hyperthermophiles for Lys).
Extensions of this work currently underway include compiling larger test sets and breaking them down by kingdom of origin as well as OGT. Heat shock proteins should be compared to regular proteins from the same organism and non-thermophilic archeal proteins should be compared to orthologs from thermophiles. Proteins from psychrophiles should be analyzed too. It would be a small step to use more sophisticated pattern recognition methods to discriminate or classify based on multiple indices.
Finally, it may be possible to design a thermostable protein from a non-thermostable one by an adaptive walk in sequence space, threading the altered sequences onto the structure of the non-thermostable protein, such that one or more of the good discriminators described here always increases.
In this paper, we have addressed the discrimination problem where given sequences or structures from a mesophilic protein and a thermophilic or hyperthermophilic counterpart, the objective is to determine which is which. This was done by assembling a large set of thermophilic protein chains from the PDB and their corresponding mesophilic analogs and another large non-redundant set of hyperthermophilic PDB protein chains along with their mesophilic analogs. They will be referred to as the pairs sets: pairs-thermophile and pairs-hyperthermophile. We have computed several structure and sequence based numerical indices, based on the quantities that other authors have reported are associated with thermostability, and tested their ability, individually, to successfully discriminate between the thermophile/mesophile pairs. One could apply more sophisticated classification or regression techniques to a combination of these quantities, but for now we have opted for a very simple test of each quantity in isolation in order to verify if it is indeed consistently associated with increased thermostability.
The pairs sets were constructed to contain pairs of high quality x-ray structures with high structural and functional similarity that differ only in that one is mesophilic and the other thermophilic or hyperthermophilic. They were assembled by compiling all PDB x-ray structures from a large list of organisms categorized as mesophile, thermophile or hyperthermophile using OGT’s obtained from the ATCC website (http://www.atcc.org/common/catalog/bacteria/bacteriaIndex.cfm). All structures with missing Cα coordinates, insertion codes, or alternate atoms were then eliminated. The resulting two lists of thermophilic and hyperthermophilic proteins were submitted separately to PISCES to generate two much smaller non-redundant sets in which all members had crystallographic resolution no greater than 2.2 Å, an R-factor no greater than 0.23, and where no pair of structures had a sequence identity greater than 30%. The members of these non-redundant sets of chains from thermophilic organisms were each then submitted to structure comparison and alignment servers (CE , SSM , DALI, VAST) to obtain mesophilic structure neighbors with rmsd no greater than about 4 Å with respect to the thermophile where the structural alignment included ~80% or more of each structure, and where the thermophile and mesophile had identical or close EC numbers or functional annotation. In some cases, more than one mesophile per thermophilic protein was kept. When multiple mesophilic analogs to a single thermophilic protein were included in the pairs sets, no restriction was placed on their similarity with respect to each other except that the sequences not be identical and that they come from different organisms. Lists of the resulting structure pairs and structure alignment data can be found among the supplementary material. The structural alignment data in these tables were computed using CE .
To train threading potentials and compute sequence composition biases, we compiled a second set of PDB x-ray structures (521nonredundant), larger and more representative than the pairs sets. The set contained 175 mesophilic, 162 thermophilic, and 184 hyperthermophilic protein structures, none of which was in either pairs set. As with the pairs sets, 521nonredundant was assembled with the help of PISCES  All members had crystallographic resolution no greater than 2.2 Å, an r-factor no greater than 0.23, and no pair of structures had a sequence identity greater than 66%. The pairwise similarity threshold was set higher for this set than pairs in order to allow the possibility that it could contain mesophilic and thermophilic orthologs, however a lower similarity threshold would have made little difference—a 30% similarity cutoff would have eliminated only 41 structures.
We have computed several structure and sequence based numerical indices to see if they can successfully discriminate between the related mesophilic and (hyper)thermophilic proteins in the pairs sets (Table 3). Among the quantities tested for discriminatory power were: the three contact-network derived quantities described in the introduction (coordination number, characteristic path length, and clustering coefficien), the ratio (E+K)/(Q+R) of Farias and Bonato , the CvP bias defined earlier , the sum and standard deviation of the Kyte-Doolittle hydrophobicities of all residues in the protein, several four-body Delaunay tessellation based threading Q-scores , the ProsaII two-body distance dependent threading score , the median circumsphere radius, the mean tetrahedrality, the surface area to volume ratio (both Delaunay and van der Waals), and the number of Delaunay simplices per residue. We have also tested the packing density defined as the van der Waals volume of the protein divided by the all-atom Voronoi volume. The algorithm of Gavezotti  was used to calculate van der Waals volumes and the Geometry Code Library 2.0 of Tsai et al.  to calculate Voronoi volumes.
We have defined discriminatory power of a given numerical index as simply the fraction f of (hyper)thermophile-mesophile pairs from a pairs set where that index is systematically larger or smaller for the (hyper)thermophilic protein, by any amount. For example, 86.1% of the time for structure pairs in the pairs hyperthermophile set, the ratio (E+K)/(Q+H) is greater for the hyperthermophilic than for the mesophilic protein, and so the discriminatory power is 0.861.
Highly over-represented and highly under-represented residue quadruplets at the vertices of tessellated thermostable proteins and the factors by which they differ with respect to mesophiles.
The work was supported by NIST and NSF. We thank Lev Krayzman for his help with the compilation of datasets.
This article has been published as part of BMC Structural Biology Volume 10 Supplement 1, 2010: Selected articles from the Computational Structural Bioinformatics Workshop 2009. The full contents of the supplement are available online at http://www.biomedcentral.com/1472-6807/10?issue=S1.
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.