Abstract

Acidithiobacillus ferrooxidans obtains its metabolic energy by reducing extracellular ferrous iron with either downhill or uphill electron transfer pathway. The downhill electron transfer pathway has been substantially explored in recent years to underpin the mechanism of iron respiration but, there exists a wide gap in our present understanding on how these proteins are organized as a supercomplex and what sort of atomic level interactions governs their stability in the iron respiratory chain. In the present study, we aimed at unraveling the structural basis of supermolecular association of respirasomes using protein threading, protein-protein docking, and molecular dynamics (MD) simulation protocols. Our results revealed that Phe312 of outer membrane cytochrome c plays a crucial role in diffusing electrons from heme C group to Asp73 of rusticyanin. In line with the previous experimental results, His143 of rusticyanin was found to have a stable interaction with Glu121 of periplasmic cytochrome c4. Cytochrome c4 interacts with subunit B of cytochrome c oxidase through Lys146 and Thr148 of the conserved hydrophobic/aromatic motif 145-WKWTFSY-151 to attain stability during simulation. Phe468 of cytochrome c oxidase was found indispensable for stabilizing heme aa3 during MD simulation. Taken together, we conclude that the molecular interactions of charged and hydrophobic amino acids present on the surface of each respirasome form a hypothetical electron wire in the iron respiratory supercomplex of A. ferrooxidans.

1. Introduction

Acidithiobacillus ferrooxidans is an extensively studied Gram-negative acidophilic chemolithoautotrophic bacterium adored for its remarkable ability to oxidize insoluble ferrous iron (Fe2+) to soluble ferric iron (Fe3+). This obligate prokaryote derives its metabolic energy by oxidizing reduced sulfur compounds or ores containing Fe2+ under acidic condition using intracellular O2 as oxidant [1]. The electrons released during oxidation of Fe2+ are transported across the cytoplasmic membrane either through a thermodynamically favorable “downhill” pathway for the reduction of O2 to H2O or pushed against the redox potential gradient along an “uphill” pathway to reduce NAD(P)+ to NAD(P)(H) [24].

Unlike many other anaerobic respiratory chains, the bioenergetic metabolism of this organism involves several proteins with highest redox potential and is encoded by rus and petI operons for the downhill and uphill pathways, respectively [4, 5]. The rus operon constitutes a set of eight genes encoding an outer membrane (OM) cytochrome c (Cyc2), a periplasmic blue copper (Cu) rusticyanin (RcY), a periplasmic cytochrome c4 (Cyc1), a putative OM protein of unknown function (ORF1), and an inner membrane (IM) aa3 Cox [58]. It has been demonstrated that expression of the rus operon at transcriptional and posttranscriptional level is specifically triggered by the presence of Fe2+ at the extracellular environment. Further, biochemical studies have established that Fe2+ is oxidized by high molecular Cyc2, present at the OM, and electrons are transferred downward to O2 through periplasmic RcY and Cyc1, finally reaching inner membrane aa3 Cox [6, 912] for energy transduction.

Recently, it has been proposed that the proteins required for downhill electron transfer are organized as a supercomplex spanning OM, periplasm, and IM [13]. Amongst the four respirasomes of downhill pathway, crystal structures of periplasmic RcY [14, 15] and Cyc1 [16] have been determined, but experimental information on the structure of OM Cyc2 and IM Cox is yet to be elucidated. This limitation is hindering postulation of a valid model for the specific binding surfaces of respirasomes. However, a working model of binding interaction between RcY and Cyc1 has been proposed by Mukhopadhyay et al., who showed that water molecules are implicated in transferring electrons from RcY to Cyc1 [1719]. It is intriguing to note that, despite more than 40 years of research, the structural features of Cyc2 and its interactions with RcY are poorly understood. In addition, the interaction between Cox and Cyc1 is still elusive. In this scenario, theoretical models for the unsolved structures of respirasomes would be highly beneficial to establish a crude hypothesis on the probable electron channel in A. ferrooxidans.

In the present context, we employed rigid-body protein-protein docking and explicit membrane water simulation protocols to construct abstract models of the respirasomes (Cyc2, RcY, Cyc1, and Cox) for identifying the domains likely to persuade transferring of electrons. In order to achieve this, first we constructed the three-dimensional (3D) model of Cyc2 using threading algorithms followed by docking with the crystal structure of RcY. Second, the X-ray structure of Cyc1 was docked with RcY that is bound to Cyc2. Third, the docked homology models of CoxA and CoxB subunits were tied up against Cyc1 of the Cyc2-RcY-Cyc1 complex. Subsequently, the whole respirasome complex was subjected to molecular dynamics (MD) simulation for obtaining completely relaxed molecular conformations. The results thus obtained provided us a convincing picture of popular electron-wire model postulated to occur within downhill electron transfer pathway [6, 16].

2. Materials and Methods

2.1. Protein Threading of Cyc2

The primary sequence of Cyc2 protein (accession number: O33823) was obtained from the UniProtKB database. Template search using BLAST and PSI-BLAST web servers [20] against Protein Data Bank (PDB) revealed that Cyc2 sequence had extremely low (10% to 14%) identity with all the solved structures available in PDB. At this low sequence identity, structure prediction using standard homology modeling protocol is prone to numerous structural errors. Therefore, the 3D model of Cyc2 was built with commonly and successfully employed fold recognition or threading approaches [21]. The target sequence was subjected to several fold recognition servers such as PHYRE [22], FUGUE [23], SAMT02 [24], pGenTHREADER [25] and two metaservers such as Genesilico [26] and Pcons [27]. The best hit of each server was compared against the secondary structure of the target protein predicted by PSIPREDv3.0 [28]. Our analysis revealed that the second top scored template (PDB ID: 2VSE [29]) suggested by pGenTHREADER had comparable secondary structure as that of Cyc2. The target-template pair-wise sequence alignment was performed using Discovery Studio v3.1 software (DS 3.1; Accelrys Software Inc.). The alignment was carefully edited to comply with the sequence alignment proposed by pGenTHREADER. The N-terminal residues of Cyc2 corresponding to the standard signal sequence were truncated before constructing the model. Ten models of Cyc2 protein were generated using DS 3.1, and the best model was selected according to minimum discrete optimized potential energy (DOPE) score. The resulting Cyc2 3D model was optimized and refined by energy minimization with steepest descent integrator using GROMOS96-43a1 force field of GROMACS software [30]. Stereochemical evaluation of the model was done in Auto Dep Input Tool (ADIT) validation server (http://validate.rcsb.org/), SAVeS server (http://nihserver.mbi.ucla.edu/SAVES/), and Z-score was calculated by ProSA-web [31].

2.2. Docking of heme C with Cyc2

In A. ferrooxidans, heme C group is attached to the CXXCH signature motif present at the N terminus of Cyc2 protein [32]. The heme C group was obtained from NCBI’s PubChem database [33], using the PubChem ID: 444125. The docking of heme C with Cyc2 was performed by PatchDock server [34] using a clustering RMSD value of 1.5 Å, which is recommended for clustering protein-small molecule docking solutions. Top ten docking solutions were downloaded from the server, and the single best solution was selected according to highest geometric shape complementarity score [35]. The interaction of heme C with its binding residues and correctness of docking were ascertained by visual inspection using graphical interface of DS 3.1.

2.3. Docking of Cyc2 with RcY

The X-ray structure of RcY was obtained from PDB using PDB ID: 2CAK [14]. The surface amino acids of Cyc2 and RcY were analyzed for charged residues that might form potential interacting interfaces. It was reported that the Cu atom of RcY interacts with heme A of Cyc1 to transfer the electron [16]. Hence, it is reasonable to postulate that opposite surface residues of the Cu binding site must interact with Cyc2. Therefore, residues near the heme C group of Cyc2 and residues opposite to Cu binding site of RcY were chosen as probable binding interface. The structure of RcY was corrected for errors like nonstandard names, alternate conformations, missing N and C termini, and incomplete atom sets in amino acids. Missing hydrogen atoms were added using DS 3.1 before subjecting to docking analysis. The docking of Cyc2 with RcY was performed with PatchDock. In order to get relevant solutions, potential binding site residues of Cyc2 and RcY were included in a text file and subjected to PatchDock for docking search stage. Residues surrounding the surface region near heme C were considered as receptor binding site, while residues opposite to Cu site of RcY were specified as ligand binding site. The final clustering of docked complexes was done by specifying a clustering RMSD value of 4 Å, which is recommended for protein-protein docking. PatchDock implements different sets of parameters optimized for different types of complexes. We specified the “default complex” type, which is recommended for general purpose protein-protein docking experiments. The best docked conformations were ranked based on geometric shape complementarity score.

2.4. Docking of Cyc1 with Cyc2-RcY Complex

The crystal structure of Cyc1 was downloaded from PDB using the PDB ID: 1H1O [16]. The structure preparation of Cyc1 was similar to RcY as described in the previous section. All hydrogen atoms were added before performing docking experiment using PatchDock. In this docking experiment, the Cyc2-RcY complex was taken as receptor, and the Cyc1 structure was considered ligand. Cyc2 and RcY were represented by chains A and B, respectively. The electrons are transferred from the Cu center of RcY to the heme A group of Cyc1 through the interaction of charged residues on RcY and Cyc1 [16]. Therefore, the oppositely charged surfaces of both the proteins were supplied as binding sites to PatchDock server. Other docking parameters are same as those used for docking of Cyc2 and RcY.

2.5. Modeling and Docking of CoxA and CoxB

The 3D model for CoxA and CoxB subunits was constructed using respective atomic coordinates of the crystal structure of Paracoccus denitrificans (PDB ID: 1AR1 [36]). Sequence alignment and model building were performed with DS 3.1 software. Accuracy of the built models was confirmed by ADIT, SAVeS, and ProSA-web servers. Final docking of CoxAB dimer with Cyc1 of Cyc2-RcY-Cyc1 complex was carried out in PatchDock using same parameters as described for protein-protein docking above. The complete molecular association of the respirasomes (Cyc2-RcY-Cyc1-CoxAB) was then energy minimized and optimized by MD simulation in explicit membrane aqueous environment before continuing further interaction studies.

2.6. MD Simulation

The proteins responsible for electron transfer in A. ferrooxidans are arranged in the form of an electron wire spanning OM, periplasm, and IM or cytoplasmic membrane [13]. Therefore, structural refinement of the docked complex was performed by placing the proteins in a two-membrane compartment system in such a manner that Cyc2 was partially integrated in OM, while CoxAB dimer was embedded in IM. For membrane construction, preequilibrated palmitoyl-oleoyl-phosphoethanolamine (POPE) bilayer was obtained from the D. Peter Tieleman website [37]. InflateGRO methodology [38] was implemented for the insertion of the docked respirasome complex within the membrane bilayer. MD simulation was carried out with a hybrid force field containing Gromos96-53a6 [39] and Berger-lipid [40] parameters. Water layers were added to the simulation system using a van der Walls radius of 0.375 Å for carbon atoms of lipid acyl chains, thus avoiding water molecules within hydrophobic core of the bilayers. A physiological ionic strength (0.15 M) of counter ions was added to attain electroneutrality. The resultant system was energy minimized in GROMACS 4.5 simulation package for 10,000 steps constraining backbone atoms of the proteins with steepest descent integrator until convergence point at maximum force of 1000 kJmol−1 nm−1 was observed. This was followed by 5,000 steps of conjugate gradient minimization without backbone constraints to obtain strain-free molecular conformations. The energy minimized structures were then divided into separate index groups for calculating energy values and other dynamic properties for individual proteins after MD simulation. Position restraints dynamics was performed by constraining backbone atoms of the proteins first in NVT ensemble for 100 ps followed by NPT ensemble for 1 ns. Finally, production run of 5 ns was performed without backbone constraints. The simulations were performed at 300 K using a time step of 0.002 ps. A distance cutoff of 12 Å with switching at 10 Å was used for van der Waals interactions. The particle mesh Ewald (PME) method was applied to treat electrostatic interactions. The trajectory data was saved at every 2,500 steps (corresponding to 5 ps), generating 1,000 structures over 5 ns of simulation for analysis. Trajectory data was analyzed using VMD 1.8.7 [41] and gnuplot 4.2 (http://www.gnuplot.vt.edu/) programs. Molecular graphics and visualization were prepared using DS 3.1 and VMD software.

3. Results and Discussion

3.1. Construction of Cyc2 3D Model

In the downhill iron respiratory chain of A. ferrooxidans, electrons are transferred along a series of membrane-located and periplasmic soluble proteins which are arranged in the form of a supercomplex through direct protein-protein interactions spanning both OM and IM [13]. Cyc2 protein is the first electron carrier in the Fe2+ respiratory pathway [6] and is one of the most diverse kinds of protein molecules having extremely low sequence identity with the homologous proteins of known structures. Our sequence analysis of Cyc2 in PDB using BlastP tool revealed that all the resultant hits have ~14% sequence identities, which is far below the minimum identity of 30% required for homology modeling [42]. Therefore, the task of template selection and target-template alignment was carried out by employing fold recognition or threading method. Previously, Schoonman et al. have suggested that the threading methods are able to predict more accurate models than comparative modeling in cases of low sequence identity (<30%) [43]. Thus, based on the recent studies of successful implementation of threading approach [4449], we carried out template search against the fold databases using several threading algorithms including PHYRE, FUGUE, pGenTHREADER, SAMT02, Genesilico, and Pcons metaservers. Amongst all the results of threading algorithms, only the mosquitocidal holotoxin from Bacillus sphaericus (2VSE) proposed by pGenTHREADER had similar secondary structure as that of Cyc2 with a value of 0.0006, SolvE score of −4.2, and net score of 48.983 amid a high confidence prediction range. Pairwise sequence alignment of 2VSE and Cyc2 is shown in Figure 1(a). Accuracy of the constructed model was validated by stereochemical analysis (dihedral angles; Phi/Psi) using Ramachandran plot [50], which indicated that 86.9% of residues fall in the most favored region, 12.2% of residues in the additional allowed region, 0.9% of residues in the generously allowed region, and no residues in the disallowed region (Figure S1 of the Supplementary Material available online at http://dx.doi.org/10.1155/2013/295718). The resulting Z-score of Cyc2 was found to be −7.25 (Figure S2), which is comparable to that of the known X-ray structures of similar size available in PDB. For additional verification of the compatibility of the generated model with its primary sequence, Verify 3D was implemented, which revealed that 91.84% of the residues had an average 3D-1D score >0.2, indicating a good quality model. Secondary structure prediction of the modeled Cyc2 protein revealed predominant presence of beta sheets and random coils with few alpha-helical segments (Figure S3). Further, the modeled protein was structurally analyzed by superimposing on the template using Cα-atoms. Although the sequence identity between Cyc2 and 2VSE was only 14%, the low RMSD value of 0.9Å indicated that the quality of the generated model is reasonably good having entire beta sheets and coil regions with similar spatial orientations (Figure 1(b)). With this well-validated geometry, the modeled Cyc2 structure was considered for further analysis. Hydropathy plot analysis suggests that Cyc2 contains no putative transmembrane (TM) segments (Figure S4). Hence, it is plausible that Cyc2 is partially integrated in the OM and the heme C containing portions; that is, CXXCH signature motif at the N terminus is likely to accept electrons released during the oxidation of Fe2+. The OM localization of Cyc2 finds support from some earlier study showing presence of cytochrome c in the OM region of some neutrophilic microorganisms deriving metabolic energy from insoluble metals and minerals [5154].

3.2. Interaction between heme C and Cyc2

The heme-containing proteins carry out numerous biological functions due to the different interactions between heme group and surrounding amino acids. In A. ferrooxidans, the heme C group of Cyc2 performs one of the most important biological functions by binding and transporting electrons released during oxidation of Fe2+. Hence, we studied the interaction of heme C with neighboring residues. In order to accomplish this, docking of heme C was carried out in the binding site of Cyc2. The heme C was docked close to the conserved heme binding motif, 43-CXXCH-47 that is present in all heme-containing proteins, and His residue of this motif is found covalently attached to the central iron atom of the heme group [32]. A close inspection of docked complex indicated that His47 of 43-CXXCH-47 is present at a distance of 2.5 Å from the central Fe atom of heme C (Figure 2(a)), and in biological conditions a covalent bond is expected between the Fe and imidazole ring of His47 [55]. In addition, heme C was observed to be surrounded by a large number of hydrophobic residues, namely, Phe227, Tyr441, Phe312, Tyr50, Trp436, Leu135, Phe136, and Trp445. The two carboxyl groups (COO) of heme C form strong hydrogen bonds (H-bonds) with the side chains of Glu438 and Asn358 (Figure S5). The stability of heme C was further confirmed after observing a strong H-bond network around its 10 Å radius (Figure S5). We found a number of intermolecular H-bonds between the residues of Cyc2 and heme C that has been shown in Table S1. The heme C contains two sulfur atoms, S1 and S2, of which S1 interacts with Thr48 through H-bonds (Figure 3). The heme C is found entangled in a hydrophobic pocket at the surface of the protein. Because the hydrophobic region should always be excluded from solvent environment, we can speculate that the heme binding pocket could be the most probable binding surface for RcY, thus indicating a probable path of electron transfer to RcY (Figure 2(b)).

3.3. Molecular Interaction between Cyc2 and RcY

RcY is the intermediate electron carrier in the downhill Fe2+ pathway of A. ferrooxidans. In the initial phase of energy transduction, Cyc2 spontaneously transfers electrons to RcY. Therefore, we examined the intermolecular interactions of Cyc2 and RcY to gain more insight about their interacting surfaces. First, we ascertained their conformational stability in bound condition during MD simulation. As can be seen from Figure 3, Cyc2 and RcY displayed backbone (C-N-CA) RMSDs of 0.44 nm and 0.12 nm, respectively. The low RMSD value of RcY as compared to Cyc2 indicates that RcY is structurally more rigid than Cyc2, probably because it is bound on either side by Cyc2 and Cyc1 in the respiratory supercomplex. On the other hand, Cyc2 being freed at extracellular end or possibly due to its larger size (449 amino acids) took longer equilibration time, but it eventually attained stability for rest of the simulation.

Average coordinates Cyc2 and RcY complex were extracted from the most stable portion of their respective dynamics trajectories. A close analysis revealed that Cyc2 interacts with RcY mainly through hydrophobic and partly with electrostatic interactions. A single H-bond (interatomic distance of 2.7 Å) was observed between main chain CO of Gly310 from Cyc2 and side chain OH of Thr126 from RcY (Figure 2(b)). By carefully observing the interacting surface of RcY, we can say that Asp73, the only negatively charged residue properly situated on the opposite surface of Cu center, could be the most probable component to facilitate electron transfer from heme C group of Cyc2 to the Cu atom of RcY (Figures 2(b) and 4(a)). Abergel et al. [16] have previously reported that in RcY-Cyc1 complex, the negatively charged Glu121 of Cyc1 is crucial for the stability during the electron transfer. Similarly, in our Cyc2-RcY interaction model, the negatively charged Asp73 of RcY might be critical for the stability of the complex. Further, previous reports have indicated that Asp73 might serve as one of the Cu coordinating ligands [56]. In contrast, subsequent site-specific mutants and gene expression results revealed that His85 acts as the Cu coordinating ligand [57]. In our present docking model, we observed that Asp73 is actually situated at the opposite surface of Cu site in RcY and thus is unlikely a coordinate Cu atom. From these observations, it can be deduced that Asp73 might be transferring electrons to the Cu center through a passage made by the core hydrophobic residues. Further, the presence of Phe312 of Cyc2 at intermediate position between heme C and Asp73 could be presumed an important factor influencing electron transfer from the donor site of heme C to the acceptor site of Asp73 (Figure 2(b)). Previous studies on the cytochromes had clearly demonstrated the significance of phenylalanine/tyrosine residues in mediating electron transfer [58, 59]. A close visual inspection of the core region of RcY structure after MD simulation allowed us to observe that several hydrophobic interactions comprising the residues Phe125, Val74, Phe76, Phe54, Phe87, and Trp7 form a channel required for electron transfer to the Cu center (Figure 4(a)). These observations suggest that hydrophobic interactions potentially play crucial roles in upholding the structure and function of Cyc2-RcY complex. The stable association between Cyc2 and RcY observed by protein-protein docking is in congruence with recently reported far-Western blotting experiments showing direct interaction between the two proteins in the respiratory supercomplex [13].

3.4. Interaction between RcY and Cyc1

Cyc1 is the third electron acceptor in the downhill electron pathway of A. ferrooxidans, which receives the electron from RcY. After determining the stable interaction between Cyc2 and RcY, we analyzed the binding interface between RcY and Cyc1. To understand these interactions, first we investigated the interaction of amino acids around the Cu center of RcY during MD simulation. The crystal structure of RcY was found to have ten pairs of electronegative and electropositive atoms within a distance of 3.5 Å around the Cu atom. We monitored breaking and formation of this H-bond network throughout the 5 ns MD trajectories. As can be seen from Table 1, seven H-bonds persisted for longer than 3 ns, and interestingly, all Cu coordinating ligands except Cyc138 continued to form H-bonds for longer than 4 ns. Similar to the calculated RMSD of RcY, the crystal structure of Cyc1 showed a low RMSD profile (0.19 nm) during MD simulation (Figure 3). Since Cyc1 is bound on one side by RcY and on the other side by Cox, it remained reasonably stable. Average coordinates were obtained from the stable dynamics trajectories, and H-bonds were analyzed. The side chain amino group (NH2) of Lys166 from Cyc1 forms two H-bonds with CO groups of Gly48 and Phe49 from RcY. The NH2 group of Gln163 from Cyc1 interacts with CO group of Phe49 from RcY. CO group of Gln115 from Cyc1 interacts with NH2 group of Lys60 and OH of Ser53 from RcY. Main chain CO of Thr112 from Cyc1 interacts with main chain NH group of Gly147 from RcY. Main chain CO of Ala145 from RcY forms two H-bonds with main chain NH of Thr112 and main chain NH of Val111 from Cyc1 (Figure 4(b) and Table S2). Our docking result identified a strong H-bond between His143 of RcY and Glu121 of Cyc1, which is in agreement with the X-ray crystallography study of Abergel et al. [16] (Figure 4(c)). Cyc1 contains two types of heme groups: heme A and heme B. The former interacts with the Cu atom of RcY, and the latter receives electron form heme A and transfers it to the Cox. Analysis of the heme groups within the Cyc1 structure revealed that carboxyl groups of heme A and heme B interact directly with each other indicating a possible path of electron transfer to the next recipient protein (Figure 4(d)). The heme A and heme B exhibit high redox potentials of +365 mV and +480 mV, respectively [60], which was found to be largely attributed to the environment of these two molecules. The CO groups of both the hemes were found to be involved in numerous H-bond interactions with surrounding amino acids of Cyc1 (Table S3). Apart from H-bonds, the two hemes are stabilized by several nonbond interactions comprising residues Cys16, Cys19, His20, Pro33, Leu35, and Met64 for heme A and Cys119, Cys122, His123, Pro134, Leu136, Leu147, and Met161 for heme B. A closer look on the number of H-bonds and the residues that participate in H-bonds revealed that heme A is stabilized by five H-bonds and six direct nonbond interactions, whereas for heme B seven H-bonds and seven nonbond interaction were observed. Several past studies have suggested that RcY is present at the intermediate position between Cyc2 and Cyc1 in A. ferrooxidans, as is shown in Figure 5. However, this model is probably specific to A. ferrooxidans, because of the fact that some acidophilic bacteria involved in Fe2+ oxidation reportedly do not possess the gene encoding RcY and that OM c type cytochrome is always located next to periplasmic c4 type cytochrome [61]. In such a situation, the OM cytochrome c and the periplasmic cytochrome c could be directly interacting in vivo while transferring electrons for energy metabolism [62].

3.5. Construction of CoxA and CoxB 3D Models

Cox is the terminal electron acceptor in the electron pathway of A. ferrooxidans. The diheme Cyc1 transfers electrons from OM Cyc2 to Cox’s catalytic subunit for the reduction of O2 to form H2O [5, 6]. Cox is comprised of four subunits, namely, CoxA, CoxB, CoxC, and CoxD. CoxB is implicated in electron transfer from Cyc1 to CoxA subunit, where O2 is converted into a molecule of H2O, whereas CoxC and CoxD mainly have a structural role [5, 36, 63]. Therefore, we constructed homology models for CoxA and CoxB subunits based on the respective subunits of Paracoccus denitrificans. BLAST search revealed that the target protein had 34% identity, 51% positives, and 73% query coverage with the template. The structure evaluation of the modeled protein using Ramachandran plot showed that 89.6% of residues fall within the most favored region and no residue was found in the disallowed region of the plot indicating conformational accuracy (Figure S6). The calculated Z-score of the modeled protein was −5.2, which falls within the Z-score range of X-ray determined structures of similar size (Figure S7). Structure evaluation using Verify 3D showed that 80.09% of the residues had an average 3D–1D score >0.2. Normally, a score >80% is a sign of satisfactory quality of the predicted model. The cofactors in the core region of the twelve TM domains, namely, one low spin heme a, one high spin heme a3, and one Cu atom, as reported in the experimental structure of the template, were obtained to the same coordinate space in the target model using structure alignment in DS 3.1. Similarly, homology model for CoxB was obtained from CoxB of Paracoccus denitrificans. The structure evaluation programs including Ramachandran plot indicated that 88.9% of residues were in most favored region with 0.9% residues in disallowed region. Verify 3D analysis revealed that 80.2% of the residues had averaged 3D-1D score >0.2, and ProSA-web calculation showed that the modeled structure had a Z-score that is within the range characteristic of native proteins of similar size. CoxB consists of a total of 254 amino acids with a calculated molecular weight of 28, 240  dalton [5]. The first 51 amino acids of this protein correspond to the standard signal sequence, and the mature protein is comprised of two N terminal TM domains and a C terminal periplasmic domain which are believed to transfer electrons to CoxA subunit [36, 63].

3.6. Interaction between Cyc1 and CoxAB

Our study revealed that the N terminal TM helices of CoxB interact with 8th and 9th helices of CoxA subunit mainly through hydrophobic interactions involving side chains of both the subunits. The periplasmic domain of CoxB interacts with the periplasmic loops of CoxA through numerous H-bonds. The H-bond interaction between CoxB and CoxA subunits is shown in Table S4. As Cyc1 transfers electrons to the catalytic site of CoxA by interacting with the Cu center of CoxB [5], we carried out a detail analysis of the interaction between Cyc1 and CoxB. During MD simulation, the coordinates of CoxA and CoxB were grouped into a single index so that dynamic properties and energy profiles could be collectively obtained for comprehensive analysis. Stability of the CoxAB dimer was ascertained by Cα RMSD from the starting structure, which was calculated to be 0.31 nm. Here, it may be noted that like Cyc2, CoxAB is freed from the cytoplasmic end and hence required longer simulation time to reach equilibration. The overall compactness of the individual proteins with respect to their center of masses in the complex was analyzed by plotting radius of gyration (Rg) values during MD simulation. Interestingly, we did not find an appreciable difference between the Rg values for Cyc2 and CoxAB (Figure 6). Altogether, as exhibited in Figure 6, Cyc2 and CoxAB structures showed Rg values of ~2.4 nm. In contrast, both RcY and Cyc1 structures showed Rg values of ~1.5 nm. Since Cyc2 and CoxAB are comparative models, structural artifacts are inevitable to occur in the backbone or side chain conformations, and therefore, it might have resulted in higher Rg values during equilibration simulation. On the other hand, since RcY-Cyc1 complex is packed against both Cyc2 and CoxAB in the bacterial periplasm, it experienced greater compactness as indicated by lower Rg values.

Post-MD analysis for the average coordinates of Cyc1-CoxAB revealed that Cyc1 forms a tight complex with the periplasmic domain of CoxB subunit involving numerous H-bonds and electrostatic interactions. The H-bond network between these two proteins is as follows: the C terminal CO of Phe254 and OH of Thr157 from CoxB form H-bonds with NH2 of Gln172 and epsilon nitrogen of His152 from Cyc1, respectively. The main chain NH of Gly155 from CoxB interacts with the side chain CO of Asn153 from Cyc1. One Zn2+ atom was found at the junction between CoxB and Cyc1 coordinated by His152 at a distance of 2.0 Å. The Zn2+ atom was found to be surrounded by Val159, Thr157, Gly155 of CoxB, and His152 of Cyc1 with a distance of 3.9, 2.1, 1.9, and 2.0 Å, respectively. The side chain NH2 of Lys146 from CoxB interacts with the main chain CO of Gly142 from Cyc1 through H-bonds with a bond distance of 2.7 Å. The side chain NH of Gln58 from Cyc1 forms a comparatively loose H-bond with the main chain CO of Asn228 from CoxB at a distance of 3.1 Å (Figure 7(a) and Table S5). Apart from these H-bonds, the interaction between Cyc1 and CoxB structures was found to be stabilized by several hydrophobic residues comprising Val170, Tyr63, Tyr29, and Ile62 of Cyc1 and Val159, Trp145, and Val179 of CoxB.

A closer look into the binding interface of Cyc1 and CoxB revealed that CoxB possessed a stretch of seven residues, mostly hydrophobic in nature (145-WKWTFSY-151), which forms boundary line between surfaces of the two proteins. Lys146 and Thr148 of this motif were found to be responsible for stabilizing Cyc1-CoxB complex by forming two strong H-bonds with Gly142 and Gly155 of Cyc1, respectively (Figure 7(a)). These observations agree well the experimental data suggesting that the binuclear Cu center of CoxB subunit accepts electrons from heme B of Cyc1 channeled through the WKWTFSY hydrophobic motif [5, 36, 63]. The CO group of heme B might be the electron donor in this case. The two Cu atoms at CoxB periplasmic domain are separated from each other by a distance of 2.6 Å, and a disulfide bridge is formed between Cyc222 and Cyc226. CuA is coordinated by His181, Cyc226, and Met233, whereas CuB is coordinated by Cyc222 and His230 (Figure 7(b)). From this Cu site, electron is passed to the CoxA catalytic site where CO of heme a, situated at a distance of ~7 Å from the Cu center of CoxB, acts as the electron acceptor.

3.7. Stereochemistry of heme aa3 in the Catalytic Site of CoxA

CoxA of A. ferrooxidans contains heme a at two different sites, each with a different function. The heme a acts as an electron carrier and transfers electrons to heme a3, whereas the latter binds molecular O2 and reduces it into H2O with the help of an electron received from extracellular environment (Figure 7(c)). While Fe of heme a is hexacoordinated, four from the porphyrin ring of heme a and the remaining two from His159 and His469 of CoxA, Fe of heme a3 is coordinated by five different nitrogen atoms, four from the nitrogen atoms of the porphyrin ring and one from His467 of CoxA, thus leaving the sixth site available to bind O2 [64]. The heme a3, in turn, is involved in reduction of O2 to form H2O with the help of a Cu cofactor. The Cu atom is coordinated by His333, His382, and His383 (Figure 7(b)). One of the Cu coordinating histidine amino acids might play electrophilic role required for the catalysis. The details of stabilizing H-bonds between heme groups and amino acids of CoxA are presented in Table S6. In addition to this, it was proposed that the conserved Phe468 found in between His467 and His469 of CoxA might be important for the transfer of electrons from heme a to heme a3 [65, 66], and our simulation results support this hypothesis. During the 5 ns MD simulation, two heme groups were found to be quite stable, probably due to the hydrophobic effect put forth by Phe468 as indicated by the RMSD graph (Figure 8). The results indicated that heme a experienced less fluctuation during MD simulation because it is connected to His159 and His469 from both sides, whereas heme a3 experienced comparatively larger fluctuation due to the absence of O2 at one end. The stability of heme aa3 binding region was ascertained by observing root mean square fluctuation of backbone residues during MD simulation (Figure 9). The residues that connect heme a to heme a3, that is, His467, Phe468, and His469, were found to be completely stable during simulation. Our modeling and simulation results suggest that the residues of Cox involved in stabilizing heme aa3 are found at the upper side of cytoplasmic membrane near the periplasm. From this observation, we can postulate that O2 reduction takes place at a pH nearly equivalent to periplasmic pH of 3.0 [67, 68].

4. Conclusions

The interdependence among the docked structures was first observed when OM Cyc2 was shown to interact with RcY through heme C and Phe312. At the interacting surface, Asp73 of RcY was considered to be involved in direct electron transfer between two complexes, which might be facilitated by physical interactions of Phe312 and Asp73. A detailed analysis revealed that the hydrophobic residues enclosed within core region of RcY form a passage for electrons to reach the Cu center. Subsequent docking of Cyc1 with Cyc2-RcY complex indicated that His143 of RcY associates with Glu121 of Cyc1, suggesting a probable channel of electron transfer between the two proteins. The heme A and heme B of Cyc1 play crucial roles in transferring electron received from RcY to the binuclear Cu center of CoxB through the aromatic/hydrophobic sequence motif, that is, 145-WKWTFSY-151. Phe468 was found to stabilize heme a and heme a3 byhydrophobic interaction. Our results provide a comprehensive insight on the structural organization of the respiratory supercomplex and on the probable path of electrons within downhill pathway of A. ferrooxidans. Altogether, this computational model can be used to understand the popular electron wire concept of metalloproteins.

Conflict of Interests

The authors declare that no competing interests exist with any commercial bodies.

Acknowledgments

The authors are thankful to Department of Biotechnology, Government of India, for providing computational facility to carry out the present research. The authors also declare that valid license has been acquired for the commercial software used in this study.

Supplementary Materials

The supplementary material contains hydrogen bond interactions between heme C and the residues of modelled Cyc2 within 10 Å distance from heme C; molecular interaction between RcY and Cyc1; interaction between heme groups and the amino acid ligands of Cyc1; interaction between CoxA and CoxB; interaction between Cyc1 and CoxB; and interaction between heme aa3 groups and CoxA subunit. Ramachandran plots of Cyc2 and CoxA models has been provided depicting more than 90% of residues of both proteins within allowed region of the plots. The predicted Z-scores of modeled Cyc2 and CoxA subunits has been provided as two dimensional plots with central black spots indicating Z-scores of the query protein with respect to that of the solved proteins available in PDB. The predicted secondary structure of Cyc2 protein has been shown. Hydropathy analysis is shown indicating Cyc2 protein does not contain potential transmembrane regions. The graphical representation of intermolecular interaction between heme C and the surrounding residues has been illustrated.

  1. Supplementary Materials