Proteome-wide structural analysis quantifies structural conservation across distant species

Global analysis of the AlphaFold2-predicted structures reveals distinct statistical characteristics between H. sapiens and E. coli or M. jannaschii proteomes

We obtained a total of 36,764 structures from the PDB (Berman et al. 2000; Burley et al. 2021) and 29,527 modeled structures for H. sapiens, E. coli, and M. jannaschii from the AlphaFold2 database (Tunyasuvunakool et al. 2021) (for a schematic representation for the workflow and methods used in this study, see Supplemental Fig. S1). Although the number of the structures from the PDB is greater than that from the AlphaFold2 database, the experimental structures cover no more than one-third of the proteomes for any of the three species because there are a lot of redundant structures obtained under various experimental conditions for a same protein. By comparison, the AlphaFold2 database almost covers the entire proteome, except for the noncanonical isoforms for some protein-coding genes.

Before the global analysis of all the structures, we used a graph-based community clustering approach to trim the unfolded regions and separate multidomain proteins into single domains based on the predicted alignment error matrix generated by AlphaFold2. An example of the structure trimming is depicted in Figure 1, A and B, in which human HMGB1, an abundant architectural chromosomal protein, was separated into two HMG box domains. The predicted aligned error (PAE) matrix is a heatmap image generated by AlphaFold2 that displays the distance error for every pair of residues (Fig. 1A). The PAE value for any of the two residues in a given protein indicates how confident the relative three-dimensional position between the two residues can be determined. Therefore, if two individual domains are connected by a flexible linker, as in the case of the two HMG box domains in HMGB1 (Fig. 1A), then the relative position between inter-domain residues cannot be confidently determined, but the relative position between two residues within a single domain can be confidently determined with a high PAE value. In addition, the PAE matrix can suggest nonstructured regions in a protein. In the case of HMGB1, the nonstructured C-terminal tail shows low PAE values between most intra-domain residues (Fig. 1B). To define individual domains and trim nonstructured regions without human intervention for hundreds of thousands of protein structures, we used a graph-based clustering algorithm (Leiden) to detect clusters of residues based on the graph connectivity in the heatmap image (Fig. 1B).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 1.

Analysis of the AlphaFold2-predicted structures reveals distinct statistical characteristics between the proteomes of the three model organisms. (A,B) An example illustrating how AlphaFold modeled structures were trimmed into single domains. The orange box indicates the range for each domain. The boundaries of the orange boxes were determined using the Leiden algorithm. (A) The alignment error matrix of the protein HMGB1 (UniProt ID P09429). The axes represent the positions of residues along the sequence of the protein. The x-axis corresponds to the reference position, and the y-axis corresponds to the predicted position. The color at the coordinate (x,y) indicates the AlphaFold2's expected position error of residue y if residue x serves as the reference and the predicted and true positions of residue y are compared. Darker color indicates lower errors or higher accuracy, whereas lighter color indicates higher error or lower accuracy. Thus, the matrix is an indication of the confidence of the relative position between any of the two residues. (B) The modeled tertiary structure of HMGB1 and its relationship with the two individual domains. Blue color indicates high model confidence, and orange color indicates low model confidence. (C) Distribution of the number of domains in each protein structure for all the modeled protein structures from each of the three species in the AlphaFold2 database. Unfolded domains and domains smaller than 50 aa were excluded from the analysis. (D) Distribution of the number of amino acids in each domain for all the modeled protein structures in the AlphaFold2 database. Unfolded domains and domains smaller than 50 aa were excluded from the analysis

Using the trimmed single-domain structures from the AlphaFold2 database, we examined the statistical characteristics of the H. sapiens, E. coli, and M. jannaschii proteomes (Fig. 1C,D). The results show that most proteins in the proteomes of the three species contain only one single domain. A small fraction of genes cannot encode any folded domain >50 aa (11.5% in H. sapiens, 7.6% in E. coli, and 4.2% in M. jannaschii). In general, multidomain proteins are more frequently observed in the human proteome than the other two prokaryote proteomes (average number: 1.15 for M. jannaschii, 1.16 for E. coli, and 1.58 for H. sapiens) (Fig. 1C). This result agrees with the previous finding that the fraction of proteins with two or more domains is greater in eukaryote and multicellular organisms (Björklund et al. 2005; Han et al. 2007). Although the average number of domains (excluding unfolded domain and domain size <50 aa) is greater, the average size of the domains is smaller in the human proteome than in the other two prokaryote proteomes (average size: 225 aa for M.jannaschii, 246 aa for E. coli, and 204 aa for H. sapiens, excluding domains <50 aa) (Fig. 1D). We reckoned that these observations are reasonable, as it has been established that metazoans tend to have more multidomain proteins than fungi and protozoa, which, in turn, have more than bacteria and archaea (Tordai et al. 2005). Additionally, domain size determines the versatility and mobility of the domain building blocks, and smaller domains are more likely to be used in multidomain proteins (Altschul et al. 2005; Pozzoli et al. 2008). Thus, as the complexity of the organism increases, the average size of the domains decreases.

Protein structure comparison can reveal evolutionary relationships between distantly related species

In most cases, apart from the functional sites on the protein surface and the hydrophobic core maintaining the structural integrity, a large portion of amino acid residues in a given structure is not evolutionarily important. Therefore, tertiary structures of proteins are often more conserved than their amino acid sequences.

An example to illustrate this point was shown in Figure 2, A and B, in which structural and sequence comparisons of H. sapiens glutathione S-transferase mu 4 and E. coli glutathione S-transferase yfcF were made. Here, we used a sequence-independent algorithm, TM-align, to measure the similarity between the structures. We used EMBOSS Stretcher (Myers and Miller 1988) to calculate the similarity between the amino acid sequences. As is shown, the structures of H. sapiens and E. coli glutathione S-transferase are nearly superimposable, whereas the amino acid sequences of the two proteins only share 35% similarity and 19% identity.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 2.

Comparison between structural similarity and sequence similarity or identity. (A,B) An example showing that proteins with nearly superimposable tertiary structures may have a low level of sequence similarity or identity. (A) AlphaFold2 modeled protein structures for E. coli and human glutathione S-transferase and alignments of the two structures. (B) Sequence alignment of E. coli and human glutathione S-transferase. The color indicates residue groups (blue, hydrophobic; red, positive charge; magenta, negative charge; green, polar). Sequence identity and similarity are indicated on top of the amino acid letters. An asterisk suggests a fully conserved position; a colon suggests strong similarity, and a period suggests weak similarity. (C,D) Comparison between structural similarity (TM-score) and sequence similarity (C) or identity (D) for 1000 pairs of selected experimental structures. The samples were obtained through stratified sampling, in which random sampling was performed at each stratum level (structural similarity ranging from x to x + 0.1, where x counts from 0.0 to 0.9). In each structure pair, one structure is from H. sapiens, and the other is from E. coli. Note that these two figures illustrate the correlation between structural similarity and sequence similarity or identity. The stratified sampling approach used here means that the figures do not represent the actual distribution of the structural similarity and sequence similarity or identity. For the list of selected experimental structures, see Supplemental Table S1. For the actual distribution at a genomic level, see Figure 3, D through F.

To further elucidate how structural similarity is related to sequence similarity in distantly related species, we randomly selected 1000 unique pairs of experimental structures with similarity scores roughly equally distributed in the range from zero to one (we randomly picked 100 pairs of experimental structures with structural similarity ranging from x to x + 0.1, where x counted from zero to 0.9 with an increment of 0.1). The comparison between the structural and sequence similarities in the selected structure pairs suggests that high sequence similarity is always accompanied by high structural similarity. However, the conclusion is not necessarily vice versa because when the sequence similarity is low, the proteins may still share a high structural similarity (Fig. 2C; Supplemental Table S1). This trend persists when sequence identity, rather than sequence similarity, is examined (Fig. 2D; Supplemental Table S1). Thus, sequence similarity and identity are not good indications of the relationship between two protein structures when the two are only moderately alike.

An enormous number of structural homologs can be discovered through global structural alignment, amounting to hundreds of thousands of pairs

When studying the evolutionary relationships between distantly related species, relying solely on sequence similarity and identity can be challenging due to the limited and ambiguous information they provide. Proteins often diverge to the extent that sequence-based methods alone cannot detect the similarities between them. Therefore, we wanted to investigate how much we could discover by using the structural information of the entire proteome. To accomplish this, we compared all possible combinations of individual domain structures between three species in a pairwise manner. To gain an estimation of the distribution of structural similarities between protein pairs and prevent duplicated counting of numbers, we exclusively used the AlphaFold2 modeled structures.

The global pairwise structural alignment revealed hundreds of thousands of between-species protein pairs with good structural similarity (Fig. 3A–C). For instance, when we compared the proteomes of E. coli and H. sapiens, we found 634,284 protein pairs with a structural similarity above 0.4, 122,102 pairs with a similarity above 0.5, and 33,323 pairs with a similarity above 0.6. Similarly, when we compared the proteomes of M. jannaschii and H. sapiens, we found 296,285 protein pairs with a structural similarity above 0.4, 49,986 pairs with a similarity above 0.5, and 8736 pairs with a similarity above 0.6. As the protein similarities cutoff increases, the number of protein pairs decreases sharply. It should be noted that the proteome sizes of M. jannaschii, E. coli, and H. sapiens are different, with H. sapiens being the largest, E. coli being the next largest, and M. jannaschii being the smallest. Although the absolute number of pairs at each similarity level differs, the percentage numbers are actually similar. Using a cut-off score of 0.5 (this is a widely cited and empirically determined cutoff for the TM-score [Zhang and Skolnick 2005; Xu and Zhang 2010]), we found that 0.09% of protein pairs had similar structures when comparing H. sapiens and E. coli, 0.09% when comparing H. sapiens and M. jannaschii, and 0.19% when comparing E. coli and M. jannaschii. Around 9500 protein-coding genes from human proteomes are involved in the protein pairs between E. coli and H. sapiens with structural similarity above the cutoff, and around 7200 human protein-coding genes from human proteome are involved in protein pairs between M. jannaschii with structural similarity above the cutoff. (Note that one protein may share a similar structure with multiple proteins from another proteome, so the number of protein pairs is greater than the number of protein-coding genes.)

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 3.

Distribution of structural similarity and sequence identity when comparing proteins at the proteomic scale. (A–C) Distribution of the numbers of protein pairs with structural similarity above 0.4 when comparing AlphaFold2 modeled proteins from one species with those from another. The subfigures compare the proteomes of E. coli versus H. sapiens (A), M. jannaschii versus H. sapiens (B), and M. jannaschii versus E. coli (C). Note that the number of protein-coding genes in each species differs, so a greater absolute number does not necessarily indicate a greater proportion. The left y-axis indicates the absolute number, and the right y-axis indicates the proportion in all possible protein combinations. (D–F) Overlay of 2D scatter plot and 2D contour plot depicting the distribution of structural similarity and sequence identity of the proteins. The subfigures compare the proteomes of E. coli versus H. sapiens (D), M. jannaschii versus H. sapiens (E), and M. jannaschii versus E. coli (F). The vertical dashed lines indicate the sequence identity of 0.2 and 0.25, and the horizontal dashed line indicates the structural similarity of 0.5. (G–I) The proportion of protein pairs of all possible combinations in each zone. The subfigures compare the proteomes of E. coli versus H. sapiens (G), M. jannaschii versus H. sapiens (H), and M. jannaschii versus E. coli (I). The twilight zone is defined as a sequence identity between 0.2 and 0.25. The safe zone is defined as sequence identity greater than 0.25, and the midnight zone is defined as sequence identity less than 0.2.

The AlphaFold2 modeled structures enabled traversal of the “twilight zone” at the boundary of the limitation of sequence alignment

The limit for confidently detecting evolutionary relationships among proteins in the sequence-alignment analysis is often referred to as the “twilight zone,” with a threshold of ∼20%–25% or 20%–30% sequence identity (Doolittle 1986; Chung and Subbiah 1996; Rost 1999). Esser and coworkers determined that there are only 850 proteins that share sequence identity >25% between yeast and 60 prokaryotic genomes using a local-alignment approach (Esser et al. 2004). Apparently, our findings reveal a larger number of structurally similar protein-coding genes compared with the number of homologs determined by sequence identity a decade ago, despite the greater evolutionary distance between human and prokaryotes than between yeast and prokaryotes.

So where are these structurally similar proteins located? How many of them are in the twilight zone (sequence identity 20%–25%), and how many are hiding in the “midnight zone” (sequence identity <20%)? What is the size of these zones? To answer these questions, we performed global sequence alignments between the individual domain structures and plotted the distributions of the structural similarity and sequence identity for all protein pairs with a structural similarity above 0.35 (Fig. 3D–F). As we were performing pairwise alignment of all possible combinations, most data points were concentrated in the region below the structural similarity and sequence identity thresholds, which is expected as they represent irrelevant proteins. Most protein pairs with structural similarity above 0.5 have a sequence identity below the 25% threshold. A significant number of proteins with high structural similarity are found in the midnight zone, where sequence identity is almost undetectable.

The global alignment algorithm is typically slow. Because of the limitation in computational resources, we only calculated the sequence identity for the protein pairs with a structural similarity above 0.35. Nevertheless, we used random sampling to pick a fraction of the total and obtain a representative subset to approximate the distribution of the structural similarity and sequence identity below the structural similarity 0.35 (for the distribution of random sampling, see Supplemental Fig. S2). Combining the distribution derived from the random sampling (Supplemental Fig. S2) and the complete data for proteins with structural similarity above 0.35 (Fig. 3D–F), we estimated the sizes of the safe zone (sequence identity > 25%), twilight zone, and midnight zone (Fig. 3G–I). As expected, the midnight zone is much larger than the twilight, which is in turn larger than the safe zone. The size of the safe zone (relative to the number of total possible combinations) is larger in the comparison of E. coli and M. jannaschii proteomes than in the comparisons of H. sapiens with either E. coli or M. jannaschii proteomes. For proteins with evident structural homology (>0.5), the twilight zone (with sequence identities 20%–25%) is two to five times larger than the safe zone, and the midnight zone is even larger, spanning 10–20 times the size of the twilight zone. These results indicate that in large-scale phylogenetic surveys encompassing all alignments with sequence identities >25%, we are merely sampling a tiny fraction of structurally homologous proteins. Furthermore, the proportions of numbers with structural similarity above and below 0.5 suggest that even the safe zone is not entirely safe. However, this largely depends on how we set the boundary between the twilight zone and the safe zone. If the boundary was set to a sequence identity of 0.3, then >98% of the proteins in the safe zone should share more than 0.5 structural similarities (Fig. 4A–C). In fact, the distribution histogram for structural similarities within the twilight zone shows two clear clusters: one centered at around 0.7 structural similarity, representing the true homologous protein pairs; and the other right-skewed cluster represents groups of data points at the extreme tail of the population, corresponding to nonhomologous protein pairs (Fig. 4D–F). Given a sufficiently large number of nonhomologous protein pairs, it is inevitable that the sequence identity between some irrelevant proteins will exceed the cutoff of 0.2 by chance alone, leading to uncertainty in the determination of the homologous relationship solely based on the sequence.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 4.

Examination of the distribution of structural similarity across different ranges of sequence identity. (A–C) Numbers of protein pairs (left, y-axis) and proportion in all possible combinations (right, y-axis) across different sequence identity ranges. The inserted pictures at the top right corner depict zoomed-in figures for protein pairs with sequence identity above 0.3. The subfigures compare the proteomes of E. coli versus H. sapiens (A), M. jannaschii versus H. sapiens (B), and M. jannaschii versus E. coli (C). The light gray color bar indicates those with structural similarity less than 0.5, and the dark gray color bar indicates those with structural similarity equal to or greater than 0.5. The light gray bar is invisible in the zoomed-in pictures at the corners, because very few protein pairs are identified with structural similarity less than 0.5 within the indicated ranges. (D–F) Distribution of structural similarity within the twilight zone, defined as sequence identity between 0.2 and 0.25. The subfigures compare the proteomes of E. coli versus H. sapiens (D), M. jannaschii versus H. sapiens (E), and M. jannaschii versus E. coli (F). The bin size for the distribution plot is 0.25. (G–I) The percentage fraction of protein pairs with structural similarity greater or less than 0.5 across different sequence identity ranges. The subfigures compare the proteomes of E. coli versus H. sapiens (G), M. jannaschii versus H. sapiens (H), and M. jannaschii versus E. coli (I). The x-axis ends at a sequence identity range of 0.55–0.70 because no protein pairs are identified above this value range.

Figure 3, G through I, shows that if the boundary between the safe zone and twilight zone is set at a sequence identity of 0.25, then the safe zone would make up 0.0014% of the total possible protein combinations, whereas the twilight zone would constitute 0.0052% of the total possible protein combinations for E. coli versus H. sapiens. Similar numbers were obtained when we compared the H. sapiens proteome with the M. jannaschii proteome. The fraction numbers increased to 0.009% and 0.021% when we compared the M. jannaschii proteome with the E. coli proteome, suggesting that the selection of species may have a significant impact on the size of the twilight zone. Above the structural similarity of 0.5, the twilight zone for comparing the H. sapiens and E. coli proteomes involves around 2854 (16%) human genes and 1279 (32%) E. coli genes, and the twilight zone for comparing the H. sapiens and M. jannaschii proteomes involves around 1359 (7.5%) human genes and 497 (29%) M. jannaschii genes. In contrast, 9172 (51%) and 6846 (38%) human genes are involved in the midnight zone above a structural similarity of 0.5 when comparing the H. sapiens proteome with the E. coli and M. jannaschii proteomes, respectively. If a structural similarity above 0.5 is used as the criterion for evaluating the false-positive rates of global sequence alignment, then most of the aligned protein detected in the range of sequence identity of 0.2–0.3 should be considered as false-positive hits. The false-positive rate drops sharply when the sequence identity increases to above 0.3, which agrees with the empirical determination for the range of the twilight zone (Fig. 4G–I).

Next, we inspected the identities of the proteins found within the twilight zone. Tables 1 and 2 display the lists of the most structurally similar protein pairs in the twilight zones when comparing the H. sapiens proteomes with either the E. coli proteome or M. jannaschii proteome (Tables 1, 2; Supplemental Tables S2, S3). Membrane transporters are prominent in the most structurally conserved protein list of the twilight zone when comparing H. sapiens to E. coli (Table 2; Supplemental Table S2). Incidentally, it is known that membrane proteins are less conserved than water-soluble proteins across the tree of life in terms of sequence, presumably due to a stronger adaptive selection in changing environments (Sojo et al. 2016). Apart from the membrane proteins, other proteins in the two tables participate in various biological processes, such as small-molecule metabolism, protein synthesis, and protein degradation. No other generalizable features are readily apparent from direct observation of the tables.

Global comparison of the structures between H. sapiens and E. coli or between H. sapiens and M. jannaschii indicates that structurally conserved genes are enriched in specific pathways

Next, we asked which protein is the most structurally conserved when comparing all human proteins with those from E. coli or M. jannaschii. Because we are not examining the distribution of structural similarity and size of the twilight zone, duplicate counting is no longer an issue. We included all the experimental structures along with the computed structures in the analysis. The best aligned score among all the structures was used for the evaluation of the level of structural conservation. Shown in Tables 3 and 4 are lists of proteins (domains) sharing the greatest structural similarities upon comparing the structures of the two species. As is evident from the tables, all the structurally conserved proteins between H. sapiens and E. coli are metabolic enzymes involved in carbohydrate, amino acid, and nucleotide metabolism (Table 3), whereas almost all the most conserved proteins between H. sapiens and M. jannaschii are involved in protein synthesis (Table 4). These results imply that eukaryotic proteins that evolved from bacteria and archaea are enriched in different biological processes, consistent with the findings of Esser and colleagues (Esser et al. 2004).

It is important to note that single gene comparison may not always provide reliable information on the evolutionary history of an organism. For instance, as shown in Table 3, the most structurally conserved gene MTR is associated with methionine synthesis. However, humans have just one enzyme related to methionine synthesis, although E. coli has all the genes required for this pathway. Humans may have lost the other genes in this pathway during evolution, but horizontal gene transfer could also be a reasonable explanation for the high similarity. The use of multiple genes is necessary to gain a more comprehensive understanding of an organism's evolutionary history.

Therefore, to further investigate the enrichment preference of the structurally conserved genes, we ranked all the human protein-coding genes based on the structural conservation scores for each gene (Supplemental Table S4) and performed gene set enrichment analysis (Mootha et al. 2003; Subramanian et al. 2005) with KEGG pathway or module gene sets (Supplemental Table S5). In line with the observation of the most conserved proteins, human proteins bearing a strong resemblance to E. coli proteins were selectively enriched in various small-molecule metabolic pathways (Fig. 5A). The glycolytic pathways and the subsequent catabolism of pyruvate were ranked at the top of the enriched pathway list, which is in line with the observation that the rate of evolution of the glycolytic enzymes is exceptionally slow (Fothergill-Gilmore 1986). Fatty acid metabolism was ranked immediately after the glycolytic pathway, consistent with the fact that eukaryotes harbor the same lipid type as bacteria (Lykidis 2007). Sugar and fat are the primary energy sources for cells. Thus, all three top enriched pathways were linked to the energy supply. In contrast, human proteins resembling those from M. jannaschii were preferentially shown as subunit compositions of cytosolic ribosome and proteasome, two complexes in charge of protein synthesis and degradation, respectively (Fig. 5B). Concordantly, proteasome and cytosolic ribosome proteins are also prominent among the 50 most conserved yeast proteins that are specific to archaea in a previous study, which identifies 850 yeast proteins with homologs in prokaryotes (Esser et al. 2004). Unlike the top enriched pathways for conserved genes to E. coli proteins, these two complexes deal with macromolecule metabolism instead of small-molecule metabolism. To confirm the differential enrichment pattern, we directly compared the structural conservation scores to E. coli and M. jannaschii. Indeed, human proteins related to energy supply are more likely to resemble E. coli structures than M. jannaschii structures (Fig. 5C), whereas the human proteins constituting cytosolic ribosome or proteasome are more similar to M. jannaschii homologs than E. coli homologs (if any) (Fig. 5D).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 5.

Global comparison of the structures between distantly related species indicates that structurally conserved genes are enriched in specific pathways. (A,B) Gene set enrichment analysis for all protein-coding genes in the human genome with structural conservation scores for H. sapiens to E. coli (A) and H. sapiens to M. jannaschii (B) used as ranking metrics. The built-in KEGG pathway gene sets were used to evaluate the enriched gene set. The bar charts show significantly enriched gene sets with normalized enrichment scores over 3.5. The enrichment analysis reveals the likelihood of the genes in the indicated gene set being ranked at the top of the ranked list based on structural similarity. (C,D) Comparison of the structural conservation scores for H. sapiens to E. coli and H. sapiens to M. jannaschii of protein-coding genes in selected gene sets. Points above the diagonal indicate that the conservation score on the y-axis is better than that on the x-axis. Part C corresponds to the gene sets enriched in the conserved genes for H. sapiens to E. coli. Part D corresponds to gene sets enriched in the conserved genes for H. sapiens to M. jannaschii. For the lists of genes in the gene sets, see Supplemental Table S5.

To further investigate the structural similarity between the three species, we compared the similarity scores for all the ribosome subunits from M. jannaschii to their annotated orthologs in the E. coli ribosome and human ribosome (Supplemental Table S6). It is worth noting that humans have two types of ribosomes, those in the cytosol, which synthesize all the protein-coding genes in the nuclear genome, and those in mitochondria, which are responsible for synthesizing 13 protein-coding genes in mitochondrial DNA. Consistent with the analysis results above, the paired comparison revealed the H. sapiens cytosolic ribosome resembles the M. jannaschii ribosome more than the H. sapiens mitochondrial ribosome and E. coli ribosome do (Supplemental Fig. S3). The human mitochondrial ribosome subunits bear the least structural similarity to the M. jannaschii ribosome subunits.

Global structural comparison provides insight into the origin of G protein-coupled receptor (GPCR) signaling pathway

Not every human protein can be assigned with a structural homolog from the E. coli or M. jannaschii proteomes. Some proteins may have homolog(s) in only one of the two prokaryotic organisms, whereas others may not have a homolog in any of the two organisms. Then what is the distribution of these proteins? If the cut-off structural similarity score is set at 0.5, then 42% of the human proteins cannot be assigned with a structural homolog from either the E. coli or M. jannaschii proteomes, 18% of the proteins can be matched to structural homolog(s) in the E. coli proteome only, and 4% of the proteins can be matched to structural homolog(s) in the M. jannaschii proteome only. The remaining 35% of the proteins have homologs in both prokaryotic organisms (Fig. 6A; Supplemental Table S7). In general, the proteins with structural homologs in E. coli outnumbered those with homologs in M. jannaschii.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 6.

Global structural comparison sheds insight into the origin of the GPCR signaling pathway. (A) Distribution of structurally conserved and nonconserved proteins in the human proteome. Fractions of human proteins with structural homologs in the E. coli or M. jannaschii proteomes are shown. A sub-pie chart is presented for proteins with structural homologs found in both the E. coli and M. jannaschii proteomes. “Similar structural conservation scores” indicates a difference smaller than 0.1 between the scores, whereas “closer to E. coli homolog” or “closer to M. jannaschii homolog” indicates that the difference in the structural conservation score is greater than 0.1. (B) Gene overrepresentation analysis for proteins with a structural homolog found in E. coli only. For the list of proteins used for the overrepresentation analysis, see Supplemental Table S7. (C) The tertiary structures of the three GPCR-like proteins found in the E. coli proteome. The orange box indicates the seven-transmembrane domain. (D) Structural analysis for E. coli dgcT. The protein consists of two domains, one seven-transmembrane domain, and one nucleotide cyclase. The most structurally similar homolog for each domain is shown.

Consistent with the gene set enrichment analysis (Fig. 5A), overrepresentation analysis showed that genes involved in fatty acid metabolism are overrepresented in the subset that has homolog(s) in the E. coli proteome only (Fig. 6B). Three seemingly irrelevant pathways (olfactory transduction, renin–angiotensin system, and neuroactive ligand–receptor interaction) are also overrepresented in the subset. In fact, a common feature between the three pathways is that all three pertain to the communication between the extracellular environment and the GPCR-mediated signaling cascade. In particular, the olfactory transduction gene set contains more than 300 GPCR proteins, and all these proteins were present in the subset with structural homologs in the E. coli proteome only.

The GPCR superfamily is the largest superfamily in the human genome, with the seven-transmembrane (7TM) domain being its most prominent feature. Although proteins with the 7TM topology have been found in prokaryotes such as rhodopsin, the phylogenetic relationships between the prokaryotic 7TM proteins and eukaryotic GPCR proteins cannot be inferred unambiguously from sequence similarity alone (Strotmann et al. 2011). However, by referencing protein structures, a previous study has provided evidence that eukaryotic GPCR may share common ancestors with prokaryote rhodopsins (Shalaeva et al. 2015).

Herein, upon close inspection of the overrepresented gene sets, we found that all the overrepresented human GPCR proteins (Fig. 6B) share a moderate structural resemblance to three E. coli proteins: dgcT, cdgI, and trhA (also known as yqfA) (Fig. 6C). What is interesting here is that both dgcT and cdgI contain two domains: one 7TM domain and one nucleotide cyclase domain (Fig. 6D). The cyclase domains from dgcT and cdgI share similar structures with an experimental protein structure for human guanylate cyclase (PDB 2WZ1). As its name indicates, GPCR interacts with G-protein, and G-protein binds to adenylate cyclase to signal the downstream events. The core tertiary structures for adenylate cyclase and guanylate cyclase are essentially the same (Roelofs et al. 2001; Winger et al. 2008). Additionally, structure-guided sequence alignment showed that the 7TM domains of dgcT and cdgI share some sequence similarity with known human GPCR proteins (Supplemental Fig. S4). Therefore, the fusion protein composed of a 7TM domain and a cyclase domain can be considered as a primitive template for the GPCR signaling pathway, in which G-protein is omitted and GPCR directly interacts with the nucleotide cyclase. Protein with such a domain architecture is absent in the M. jannaschii proteome. These results suggest that the link between GPCR-like transmembrane protein and cyclic nucleotide secondary messenger may have been established before the emergence of eukaryotes.

The global structural comparison suggests an archaeal origin of core proteins relating to information storage and processing

Next, we investigated the overrepresented protein groups in the protein subset, of which structurally similar homologs can be found in the M. jannaschii but not the E. coli proteome. The overrepresentation analysis revealed that proteins relating to systemic lupus erythematosus and alcoholism were overrepresented in this protein subset (Fig. 7A). However, further evaluation of the enriched genes suggests that the computational algorithm determined that these two pathways were overrepresented because a good number of histones were included in the gene sets. (Somehow, the KEGG database does not include a gene set dedicated to histones.) Histones have multiple copies in the human genome, and each copy was assigned a unique gene entry, leading to a disproportionated number in the two gene sets. Thus, the overrepresentation analysis actually indicates that histones are among the overrepresented proteins with homologs in M. jannaschii only. In addition to histones, the analysis reveals that proteins involved in RNA transcription (general transcription factors and RNA polymerase) and DNA replication are overrepresented in the protein subset. Given that ribosomal proteins are more conserved to M. jannaschii than E. coli (Fig. 3D), we wondered if most of the proteins relating to the central dogma of molecular biology are more similar to their archaeal homologs.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 7.

The global structural comparison suggests an archaeal origin of core proteins relating to information storage and processing. (A) Gene overrepresentation analysis for proteins with a structural homolog found in M. jannaschii only. For the list of proteins used for the overrepresentation analysis, see Supplemental Table S7. (B–G) Distribution of the structural conservation scores for protein-coding genes in the selected gene sets upon comparing H. sapiens to E. coli, Caulobacter vibrioides NA1000, M. jannaschii, Candidatus Prometheoarchaeum syntrophicum, or Candidatus Lokiarchaeum sp. GC14_75. Violin plots were used to represent the distribution of the structural similarities, with the name of the gene set marked at the top of each panel. The three dashed lines represent 25% quartile, median value, and 75% quartile. For the lists of genes in the gene sets, see Supplemental Table S5.

Herein, we expanded our analysis to include two additional archaeal species, Lokiarchaeum sp. GC14_75 and Candidatus Prometheoarchaeum syntrophicum, and one bacterial species Caulobacter vibrioides. We compared the distribution of structural similarities for genes involved in gene sets relating to DNA storage (histones), DNA replication, RNA transcription, and protein translation (ribosome), as well as the distribution of structural similarities for genes relating to fatty acid metabolism and sugar metabolism (Fig. 7B–F). Lokiarchaeum sp. GC14_75 and C. Prometheoarchaeum syntrophicum are two recently discovered archaeal species that are closely related to eukaryotes (Zaremba-Niedzwiedzka et al. 2017; Imachi et al. 2020). These species possess genes involved in functions typically associated with eukaryotes, such as membrane trafficking and cytoskeletal organization. These features suggest that the two species may be closely related to the ancestral archaea involved in the symbiotic interactions that led to the emergence of the first eukaryotic cells. In contrast, although it has been established that mitochondria originated from an alphaproteobacterial-like ancestor (Martijn et al. 2018), identifying the closest species to mitochondria poses a challenging task. This difficulty arises because mitochondria have lost most of their genes over the course of evolution, transferring them to the nucleus. Although endosymbiotic alphaproteobacteria (e.g., Rickettsia prowazekii) may resemble mitochondria, they have lost many modules and pathways, making them less suitable for evaluating the level of structural similarities of pathways. Therefore, we opted to include the free-living alphaproteobacterium, C. vibrioides, in our analysis to complement the comparison between E. coli and H. sapiens. By including these three additional species, we hoped to gain further insight into the evolution of the key cellular functions.

Indeed, proteins relating to central dogma are more likely to show stronger structural resemblance to homologs in archaeal species than in bacterial species (Fig. 7B–E). Conversely, proteins relating to energy supply are more likely to show a stronger resemblance to homologs in bacterial species than in archaeal species (Fig. 7F,G). Overall, the proteins relating to energy supply tend to be more structurally conserved than those involved in central dogma machineries in the species examined here. Human genes involved in fatty acid metabolism display significantly higher structural similarities to their homologous proteins in the two Asgard archaeal species compared with M. jannaschii, although bacterial proteins still show greater overall structural resemblance for these human genes. In summary, our analysis of structural similarities supports the idea that the central dogma machineries in eukaryotes evolved from those in an ancestral archaeal species, which is consistent with the findings of a recent study that showed the genetic information processing genes of eukaryotes have evolved from archaeal species (Brueckner and Martin 2020).

SHF: an online resource for structural homolog identification

As shown above, we have gained a lot of interesting information on the molecular evolution of proteins by comparing the protein structures and identifying structural homologs at the proteomic scale. We believe what has been discovered here is only the tip of the iceberg. Some more helpful information could be extracted if we analyze the data from other perspectives. Therefore, we have made these data available to the public as a searchable online resource (the Structural Homolog Finder [SHF]). We have additionally included some structural analog information for Saccharomyces cerevisiae, the unicellular eukaryotic model organism, in the online resource. The webpage can be accessed at http://micc.tmu.edu.cn/shf. On the web interface, the user may enter either the gene symbol or UniProt ID of the gene of interest for querying structural homologs in another species. Then, the user should specify the species names for the source and target organisms (Supplemental Fig. S5A). The backend program will search for structural homologs in the collection of protein structures from the target organism for the specified protein from the source organism. A typical run of the query takes <10 sec. An example of the output page is shown in Supplemental Figure S5, B and C. A ranked list based on structural similarity and a visualization of the structure overlay are provided to the user. For multidomain proteins, homologs to all the individual domains were listed and ranked together. If no structurally similar homologs were found, the query will return “No similar structure found” but will still return one most similar protein structure (of course, with a poor similarity score) from the target organism.