A three-dimensional topology of complex I inferred from evolutionary correlations
© Kensche et al.; licensee BioMed Central Ltd. 2012
Received: 13 February 2012
Accepted: 28 June 2012
Published: 3 August 2012
Skip to main content
© Kensche et al.; licensee BioMed Central Ltd. 2012
Received: 13 February 2012
Accepted: 28 June 2012
Published: 3 August 2012
The quaternary structure of eukaryotic NADH:ubiquinone oxidoreductase (complex I), the largest complex of the oxidative phosphorylation, is still mostly unresolved. Furthermore, it is unknown where transiently bound assembly factors interact with complex I. We therefore asked whether the evolution of complex I contains information about its 3D topology and the binding positions of its assembly factors. We approached these questions by correlating the evolutionary rates of eukaryotic complex I subunits using the mirror-tree method and mapping the results into a 3D representation by multidimensional scaling.
More than 60% of the evolutionary correlation among the conserved seven subunits of the complex I matrix arm can be explained by the physical distance between the subunits. The three-dimensional evolutionary model of the eukaryotic conserved matrix arm has a striking similarity to the matrix arm quaternary structure in the bacterium Thermus thermophilus (rmsd=19 Å) and supports the previous finding that in eukaryotes the N-module is turned relative to the Q-module when compared to bacteria. By contrast, the evolutionary rates contained little information about the structure of the membrane arm. A large evolutionary model of 45 subunits and assembly factors allows to predict subunit positions and interactions (rmsd = 52.6 Å). The model supports an interaction of NDUFAF3, C8orf38 and C2orf56 during the assembly of the proximal matrix arm and the membrane arm. The model further suggests a tight relationship between the assembly factor NUBPL and NDUFA2, which both have been linked to iron-sulfur cluster assembly, as well as between NDUFA12 and its paralog, the assembly factor NDUFAF2.
The physical distance between subunits of complex I is a major correlate of the rate of protein evolution in the complex I matrix arm and is sufficient to infer parts of the complex’s structure with high accuracy. The resulting evolutionary model predicts the positions of a number of subunits and assembly factors.
NADH:ubiquinone oxidoreductase (complex I) is with about 1000 kDa [1, 2] the largest of the five complexes of the oxidative phosphorylation (OXPHOS) and a major contributor to the proton motive force that drives the ATP production by ATP-synthase . Complex I has an L-shape with a hydrophilic matrix arm that protrudes into the cytoplasm in bacteria or the mitochondrial matrix in eukaryotes and a hydrophobic membrane arm. The canonical “core” of complex I consists of 14 subunits that originate from three pre-existing evolutionary modules . The N-module at the distal end of the matrix arm contains flavin-mononucleotide (FMN) that accepts electrons from a donor, usually NADH. The electrons are transported through a chain of iron-sulfur (FeS) clusters along the matrix arm towards the joint of the two arms at the membrane. This membrane-proximal part of the matrix arm represents the Q-module in which the electrons are transferred to ubiquinone (Q). The energy freed by the electron-transfer is transmitted along the P-module (NADH1-6/4L) that uses the energy to pump protons across the membrane [5–7].
In diverse taxa, the canonical core of complex I has been extended by further subunits. For instance, complex I in Thermus thermophilus contains an additional subunit located at the interface of the N- and Q-modules  and a recent analysis of complex I in the α-proteobacterium Paracoccus denitrificans identified three additional subunits . Eukaryotes obtained complex I with the endosymbiotic uptake of an α-proteobacterium that gave rise to the mitochondria. Following the endosymbiosis, the mitochondrial genome was reduced and the genes encoding matrix arm subunits of complex I were transferred to the nucleus. Additionally, complex I was extended to up to 45 subunits by so-called “accessory” or “supernumerary” subunits [1, 10]. This set of permanent subunits is further extended by a number of assembly factors absent from the mature complex [11–19].
Up to now, the structures of the complete complex in the eubacterium Thermus thermophilus and the eukaryote Yarrowia lipolytica have been published. However, the latter structure is of a too low resolution to allow identification of the positions of the supernumerary subunits. Approximate subunit positions within the eukaryotic complex are hinted at by various types of experiments, mostly from sub-complexes observed by fractionation or as assembly intermediates (e.g. [1, 20]). For instance, the application of chaotropic detergents to the bovine complex produces the three sub-complexes Iα, Iλ, and Iγ. Because these sub-complexes are large, they provide only rough information about subunit positions. For instance, Iα represents an extended Iλ sub-complex and the additional subunits could in principle be located anywhere on the surface of the Iλ sub-complex. Only limited data are available from yeast-two-hybrid [12, 21], co-immunoprecipitation , or cross-linking  experiments. The identification of the positions of the assembly factors is hampered by the temporariness of the assembly intermediates and our incomplete understanding of the assembly process.
The increasing number of genome sequences allows making predictions of physical interactions by evolutionary correlation methods, including the co-occurrence of genes or phylogenetic profiling [23–26], the mirror-tree approach , and residue correlation [28–31], which have successfully identified new complex I subunits and assembly factors [21, 25] (reviewed in ) and predicted relations between the five OXPHOS complexes . Of these methods, residue correlation is based on the most direct evidence of physical interaction, namely the compensatory mutations at sites of interacting residues to maintain the structure of a protein or complex. By contrast, the mirror-tree method detects co-evolution more indirectly by correlating sequence similarity matrices between orthologous groups . The similarities between protein sequences depend both on species divergence times and on rates of evolution. By removing the similarity due to the species divergence times [34, 35], one obtains similarities that are more related to evolutionary rates. A high correlation in evolutionary rates between protein families can be evidence of a direct physical relation between proteins. For instance, if there is selection to maintain the interaction of two proteins then disrupting mutations have to be compensated for at the rate that they occur. Therefore, to maintain the interaction, an increased rate of change in one protein needs to be compensated for at a similar rate in the other protein. Note that the pairwise correlations in evolutionary rates between proteins can also be due to indirect interactions .
Here we ask whether we can use evolutionary rate correlations to predict the three-dimensional (3D) conformation of complex I. To this aim, we analyzed the evolutionary rate correlation between 38 subunits and 7 assembly factors of human complex I using the mirror-tree method  and find that subunits that are known to be physically close in complex I tend to show a higher degree of correlation in evolutionary rates than those that are physically distant. In the conserved core of the matrix arm, this correlation is strong enough to construct a 3D model with striking similarity to the bacterial reference structure. In a second evolutionary model that includes the 14 canonical core subunits, the membrane and matrix arms appear as clearly separated groups. Finally, we calculated a third evolutionary model including 38 subunits and seven assembly factors. This last model retains some features of the physical structure, including the separation of the matrix and membrane arms and a proximodistal axis in the matrix arm. We discuss the positions of the seven assembly factors in this model and make specific predictions about the association of some assembly factors with each other and with the permanent subunits.
Because our study aims at complex I of human we selected 38 subunits and seven assembly factors of the human complex that have a sufficient number of orthologs for the application of the mirror-tree method (Additional file 1: Table S1). Shortly, we collected homologs of the 45 proteins by querying the nr database. All orthologous sequence sets were aligned and highly variable alignment columns were filtered out using BMGE . Note that NADH3, 4L and 6 could not unambiguously be located in our reference – the structure of the complete complex I of T. thermophilus (PDB:3M9S). However, these subunits are known to be direct neighbors and we decided to treat them as a single unit, termed NADH34L6. We calculated maximum-likelihood trees from the alignments  and obtained a distance matrix for each protein family from which we removed the common signal of the phylogeny . The phylogeny-corrected matrices were correlated and the resulting correlation matrix was transformed into the distance matrix (see Methods). For the purpose of this article we call these distances “co-evolutionary” but note that the signal measured by the mirror-tree method is also determined by expression and general functional relatedness [39, 40].
After the prediction of the topology of the seven matrix arm core subunits, we predicted the topology of the complete core of 14 matrix- and membrane arm subunits that are conserved among bacteria and eukaryotes. Again, the co-evolutionary distances of the subunits can be well embedded in three dimensions (P 3 =0.84) and result in a configuration with separate membrane and matrix arms (Figure 2c). The rmsd of this evolutionary model and bacterial structure is 47.6 Å, which corresponds to about 25% of the length of the membrane arm . Also in the complete-core model, the proximodistal axis of the matrix arm is recovered (r2=0.96, p=9.3×10-5) although the rmsd of the matrix arm core subunits to the bacterial structure is lower (34.0 Å) than that for the configuration that only contained matrix arm subunits (18.7 Å). As expected from the correlation between the pairwise physical distances and their evolutionary correlation, the accuracy of the positioning in the membrane arm is poor (rmsd=61.8 Å). A closer inspection of the evolutionary model suggests that the membrane arm subunits are located as a cluster more sideways of the matrix arm than in the physical structure. In the superimposition of the two structures, this manifests in a tendency of the Q-module away from the membrane and a slight tendency of the distal N-module towards the membrane.
Our results show that the evolutionary rates of complex I subunits contain a significant amount of information about the complex’s quaternary structure. For the matrix arm we found that about 61% of the correlation in evolutionary rate could be explained by the distances of the subunit centers. This is even more striking if we consider that the evolutionary model was derived from eukaryotic sequences and thus should reflect the matrix arm structure in eukaryotes, while our reference structure is from a bacterium. Indeed, the evolutionary 3D model revealed a twist between the N-module and the Q-module when compared to the bacterial structure, a finding that is supported by experimental data (, personal communication).
In the two models that include the membrane arm, mitochondria-encoded subunits were predicted to be separate from nucleus-encoded subunits, which is in line with previous results. The independent variation of evolutionary rates in the nuclear and mitochondrial genomes  may have contributed to the isolation of the mitochondria encoded subunits in our models. Nevertheless, we stress that the position of the membrane arm core subunits, specifically at the proximal end of the matrix arm in both models, indicates a signal of the physical structure in the evolutionary correlation data. Furthermore, the strict separation of the nucleus-encoded Iβ subunits and their mitochondria-encoded counterparts NADH4 and 5 may be explained by other factors, as these two groups also behave differently in experiments [1, 52, 53]. Interestingly, despite its nuclear encoding, the membrane-integral subunit NDUFA1 of the Iα sub-complex [1, 2, 54] is positioned close to the membrane arm core, in particular close to NADH2 (Figures 3a and 3b). In T. thermophilus, NADH2 is located between the two subunits NADH1 and NADH4  both of which are known physical interactors of NDUFA1 . A direct physical interaction of NDUFA1 with NADH2 is therefore likely.
The evolutionary correlation failed to identify the correct topology of the membrane arm core. A number of biological reasons could explain such a lack of signal. First, long-range structural constraints [5, 6] may interfere with the distance-dependent structural constraints that are necessary for a distance-dependent strength of evolutionary rate correlation. Second, the formation of OXPHOS super-complexes with complex I dimers may result in correlations between distant subunits. Indeed, despite their positions at opposite ends of the membrane arm, NADH1 and 5 show a high correlation in evolutionary rates with each other and with subunit CYTb of complex III  consistent with their proximity in OXPHOS complexes organized into respiratory strings . Third, the lack of correlation with physical distance may result from non-adaptive variation in the mitochondria-encoded genes caused by variable and, at least in some eukaryotic taxa, heterogeneous mutation-pressure . Indeed, in a number of animal taxa changes in gene order or mutation-pressure led to non-adaptive changes in mitochondrial genes [57–59]. The mitochondrial genomes of some taxa in our study, such as plants, are clearly different from those in animals (reviewed in ) and their genes are likely under different mutation-pressures . Fourth, the embedding of the membrane proteins in two dimensions might reduce the evolutionary constraints to maintain interactions in comparison to proteins that are embedded in three dimensions.
The integration of multiple proteins in a single model assumes that the interactions are permanent and non-competitive. This is clearly not the case for the model of 45 proteins because it includes assembly factors. This model can therefore not exactly represent a physical structure. According to current models, complex I assembles from independent subcomplexes . Of the assembly factors required for this process and included in our study, only NDUFAF1 (AF1) is required for the assembly of the distal membrane arm sub-complex [13, 63, 64]. In our model, AF1 is located close to the matrix arm, which supports an indirect rather than a direct involvement of AF1 in membrane arm assembly . The distal membrane arm further combines with a pre-formed membrane-anchored proximal matrix/membrane arm that contains the subunits NDUFS2 (S2) and NADH1 (1) and possibly NDUFS3 (S3) and NDUFS7 (S7) [62, 64] and whose assembly involves NDUFAF3 (AF3) and possibly C8orf38 (O38) [17, 66]. Although the membrane-association of AF3 and O38 is not reflected in our data, they form a tightly co-evolving triple with C2orf56 (O56), which is known to bind the proximal matrix arm subunit S2 . The high correlation in evolutionary rates between AF3, O38, and O56 suggest strong selective constraints on their cooperation during the assembly of the proximal matrix/membrane arm sub-complex. The fourth assembly factor that has been experimentally linked to the proximal membrane arm, C20orf7 (O7) [18, 64], is indeed placed close to the proximal matrix arm subunits S2, S3 (Figure 3b, right bottom), and the proximal membrane arm subunits A3, A6, and A8 [2, 47, 48].
After the joining of the two membrane arm intermediates, the proximal matrix arm is further extended. This step involves the NUBPL-mediated assembly of at least one FeS-cluster into the distal matrix arm [11, 67]. In the evolutionary configuration the assembly factor NUBPL is positioned side by side with the permanent subunit NDUFA2 (A2; Figure 3b, right top). Like NUBPL, A2 is associated with the distal matrix arm . The highly conserved A2 subunit is structurally similar to thioredoxin-like proteins with a loop-region of probably variable conformation that contains two cysteines in human (C24 and C58) . These cysteines can form a revertible disulfide bridge with an in-vitro redox-potential in the range of the large majority of isopotential FeS-clusters of complex I [69, 70]. Although the cysteines are not fully conserved, occasionally FeS-clusters are bound by serine, histidine, or aspartate . Indeed, the human serine 30 in NDUFA2 is a good candidate for FeS-cluster binding because it is perfectly conserved in all species, with the notable exceptions of Trypanosoma and Leishmania, in which it is substituted by cysteine. Together these observations and the very strong evolutionary rate correlation of A2 and NUBPL support an involvement of A2 in complex I associated FeS-cluster assembly or maintenance. The peripheral position of A2 and NUBPL in the model could be a consequence of other strong evolutionary constraints not directly related to complex I.
Also NDUFAF2 (AF2, B17.2L) has been linked to the assembly of the distal matrix arm [14, 64]. Interestingly, the evolutionary data position AF2 directly besides its paralog NDUFA12 (A12, B17.2) [10, 14]. Like AF2, A12 is known to be associated to the distal-matrix arm to which it is directly recruited from the mitochondrial matrix . The correlation in evolutionary rates and the independent co-loss in multiple complex I lacking taxa  support an evolutionarily conserved functional relationship of AF2 and A12. It is tempting to speculate that AF2 temporarily binds at the binding site of A12, e.g. to stabilize the local structural context, and is later substituted by its paralog. Such close positioning and physical interaction of homologous proteins within the same protein complex is one of the prevailing trends in the “fate” of duplicated proteins in complexes . Complex I appears to add another twist to this pattern in the sense that the predicted interaction is only temporary.
The rate of protein evolution is influenced by diverse factors , in particular expression and general functional relatedness [39, 40, 74]. It is therefore even more remarkable that we found physical distance to be the major determinant of the evolutionary rate correlation for the complex I matrix arm. However, this result does not apply to the whole complex. Thus, to establish whether the mirror-tree/MDS combination is a good general method to predict quaternary structures, other complexes need to be analyzed. Furthermore, instead of using the mirror-tree method one could use residue correlation to measure the co-evolution of subunits more directly. Residue correlation has been used to predict contact interfaces for protein pairs [30, 75] and to investigate a rotation-symmetric homo-multimeric complex . A simple implementation would be to integrate pairwise residue correlations  or correlations that account for indirect correlations [30, 31, 77] or phylogenetic dependency [76, 78] by in-silico two-hybrid  into subunit distances and map these into three dimensions by multidimensional scaling.
The correlations of evolutionary rates between subunits of the eukaryotic complex I contain detailed information about the structural arrangement of the matrix arm subunits. This allowed us to make specific predictions about the positions of supernumerary subunits and assembly factors of the matrix arm, which may guide further experimental investigations. Multidimensional scaling could not reconstruct the structure of the membrane arm core. A future analysis will have to investigate what may cause this lack a spatial signal along the membrane arm and thus clarify in particular the relevance of conformational dynamics and super-complex arrangement into a respiratory string for the sequence evolution of complex I.
We included 38 permanent subunits and seven assembly factors of human complex I for which a sufficient number of sequences were available. We collected homologous sequences from the nr database  using PSI-BLAST (default parameters). Multiple queries from different species were used whenever PSI-BLAST failed to find known homologs [see Additional file 1. For A6, B9, A12, and AF2, orthologous groups were manually identified in neighbor-joining trees constructed with identity matrices and correcting for multiple substitutions. Species overlap between the partitions was used to divide the trees into separate orthologous groups. The remaining subunits were treated by a different protocol. First, to ensure a separation of the paralogs NADH2, 4, and 5, we built a set of trusted orthologs of NADH2, 4, and 5 from those sequences that had the best bi-directional hit with the human query using PSI-BLAST. From these seed sequence sets we computed three HMM profiles and sorted the remaining homologs into the orthologous group to which they showed the best profile-alignment . For all sequence sets we selected as single ortholog per species the sequence with the highest NEEDLE score in a pairwise alignment to the human query  (default parameters) and/or manual selection based on multiple alignments (MAFFT , CLUSTALW [84, 85], HMMER , HHSEARCH ). The kinetoplastida were excluded from our analysis due to their high level of sequence divergence. To gain high quality alignments, we aligned all sequence sets with CLUSTALW and manually fixed misalignments. The manually curated alignments are provided in Additional file 3. Next, we filtered alignment columns with BMGE  (−m BLOSUM30 -g 0.50 -b 4), removed sequences that had more than 33% gaps, and restricted the alignments to those species for which we found at least eight subunits of the complex. Of the 43 alignments, 39 had more than 75 sequences and there was no alignment with less than 44 sequences. Finally, we calculated phylogenetic trees using RAXML  (Version 7.2.6, PROTGAMMAMTREV for NADH1/2/3/4/4L/5/6, otherwise PROTGAMMAJTT; 4 rate categories) [see Additional file 4. A single tree was calculated for the concatenated alignment of NADH3, 4L, and 6.
with the row vector p T . We derived the reference vector p directly from the subunits’ distance vectors, as suggested by Kann et al.. Specifically, the reference distance between a pair of species was calculated as the average of the distances between these species in the trees of the complex I subunits. Note that this choice of reference as an average of the analyzed vectors also removes the specific pattern of co-variation in evolutionary rates that reflects selective pressure on the complex as a whole. It thus focuses the results on distances between the subunits rather than their distances to unrelated proteins. Finally, the corrected distance vectors were correlated with each other by Spearman rank correlation to yield the subunits' co-evolutionary similarity . We required that the species pairs were present in at least five of our 43 distance vectors. Species pairs occurring in fewer than five vectors were ignored in the correlations. Note that our choice of the set of subunits included the requirement that there are at least 15 species in all pairs of alignments. Only for 17 out of 903 subunit pairs, the correlation values were based on less than 30 species. The mirror-tree method has the advantages of being easily implemented and it requires low computational resources, even with correction for the basic correlation due to the shared phylogeny.
where λ i is the i-th largest eigenvalue of B Δ . Note that the relation between the co-evolutionary dissimilarity and the distance in the cMDS configuration (Shepard diagram, Additional file 1) indicates that the 3D configuration reflects the raw co-evolutionary distances over its whole range.
where |t i - c i | is the distance between the bacterial and predicted center of the i- th subunit.
We thank Richard Notebaart and Fiona Nielsen for valuable discussions. PRK was funded by the BioRange program of the Netherlands Bioinformatics Centre (NBIC) supported by the Netherlands Genomics Initiative (NGI) and the European Union’s Sixth Framework Program EPISTEM (CT-2005-019067). ID was funded by the Portuguese Foundation for Science and Technology (SFRH/32966/2006) and by Bolsas Rui Tavares 2010.
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.