Macroscopic View: Structural Characterization of the Protein’s Global Dynamics
Overall Three-Dimensional Structure. A series of MD simulations up to 0.6 µs of the human WT and F583A mutant in a hydrated lipid bilayer was performed. To ensure the selection of a representative structure, a cluster analysis was performed to select a conformation for further analyses. Representative structure from the WT and F583A mutant, superimposed on to bovine MRP1, is presented in Fig. 2A. We found that our proposed model of the human MRP1 displays the typical topological architecture of eukaryotic ABC transporters. It is composed of 12 TM helices distributed among two domains: TMD1 (helices 6–11) and TMD2 (helices 12–17). The 12 helices from TMD1 and TMD2 are arranged into two pseudo-symmetric bundles: bundle 1, comprising helices 6, 7, 8, 11, 15, and 16; and bundle 2, comprising helices 9, 10, 12, 13, 14, and 17. The two helical bundles are opened to the cytoplasm and comprise a space spanning halfway into the lipid bilayer. The two NBDs, which are located about 30 Å down into the cytoplasm, are completely separated from each other to about 34 Å. Overall, the suggested model highly resembles the bovine MRP1 and other eukaryotic ABC transporters and retains its distinctive arrangement and characteristics.
It is possible to identify the various conformations that the protein assumes during the simulations through cluster analysis of the structures. Accordingly, our 3 Å cut-off cluster analysis results show that the F583A mutant protein acquires more conformations during the simulation, in comparison with the WT protein. 94 % of the WT protein’s conformations belong to cluster #1, presenting one main stable conformation. In contrast, only 32 % of the mutant protein conformations belong to cluster #1, which implies its instability in relation to the WT (Fig. 2B). The WT protein appears to be more rigid than the F583A mutant that is more flexible and less stable. Also, The WT protein converged spontaneously, and relatively quickly into a stable conformation after ~ 30 ns and it maintains that form to the end of the simulation. On the contrary, the mutant protein fluctuates between different conformations and diverges during the entire simulation (Fig. 2C). These findings support that the WT protein being more stable than the mutant.
Evolutionary Conservation. We found that the most conserved residue in the extracellular region of the protein is F583 (Fig. 3A). It was rated similarly to highly conserved sequence elements in the protein, such as Walker A, B, and the signature motifs, from both NBDs, that are crucial for the ATPase function of MRP1. Since conservation is highly predictive in identifying catalytic sites and residues near bound ligands 12, F583’s high conservation score implies its importance for the function of MRP1. Substitution of conserved amino acids was found to lead more frequently to loss of function than substitutions of non-conserved residues 13.
The level of conservation is based on the evolutionary relatedness between MRP1 and its homologs. The human ABCC subfamily of transporters contains 13 members from the ABC superfamily. The ABCC subfamily includes the Cystic Fibrosis Transmembrane Conductance Regulator (CFTR/ABCC7), two Sulfonylurea receptors SUR1/ABCC9 and SUR2/ABCC9, and nine MRPs 14. Using CLUSTALW, sequence alignment of the ECL5 loop (connecting helices 10–11) of the ABCC homologs revealed the invariant conservation of F583, in contrast to other residues in the ECL5 loop (Fig. 3B). Given the high evolutionary conservation of F583 compared to the rest of the ECL5 loop, it is reasonable to assume that F583 may play a major role not only in the protein’s function but in maintaining its structure as well.
Collective Motions of MRP1. The correlated modes of the protein’s motion during the simulations were calculated by Principal Component Analysis (PCA). This analysis filters global, slow, and probable concerted motions, that contribute most to the total atomic displacement. The first five eigenvectors were considered for this analysis, as they are the most important to describe the significant motions of the protein. The covariance matrix of atomic fluctuations was diagonalized for predicting the eigenvalues. The first five eigenvalues in descending order versus the corresponding eigenvector for WT and F583A mutant proteins are shown in Fig. 4A 15. During the simulations, the WT and the mutant had similar eigenvalues for eigenvectors two to five. However, the F583A mutant showed higher correlated motions (161 ± 8 nm2) in the first eigenvector compared with the WT protein (116 ± 7 nm2), as presented in the covariance matrix. The first eigenvalue corresponds to the maximum variance of the projected points. Therefore, a dynamic comparison of the first eigenvector of the WT and the F583A mutant was conducted.
Visualization projection of the first eigenvector, presenting the greatest difference between the WT (blue) and the F583A mutant (red) in the NBDs area is shown in Fig. 4B by a porcupine-like representation. The lengths of the arrows indicate the degree of the significance of the motion; the longer the arrow, the motion is more significant, and vice versa. In the WT protein, the NBDs moving away from each other, creating a putative path for ligand entrance, as expected compared to the transport mechanism (see Introduction). By contrast, the NBDs of the F583A mutant are moving in the opposite direction, towards each other, and by that deviate from the transport mechanism and may prevent proper binding of a ligand. NBD1 and NBD2, which are connected by a flexible hinge, move in opposing directions at the WT protein and the F583 mutant. This striking difference between the correlated significant motions of the protein exemplifies how the single mutation F583A induces changes in the dynamic properties of domains that are 34 Å apart from each other, and 75 Å away from the actual location of the mutation.
Microscopic view: Structural characterization of helices 10–11 and NBD2
To account for the F583A mutation, we focused on the structural alterations of helices 10–11 since ECL5 is located between them. Then, we analyzed the conformational changes of NBD2, followed by the interaction between helix 10 and NBD2, because of F583A’s effect on ATPase activity in NBD2.
Mutual Progression of Helices 10–11 and the Protein. Dynamic comparison of the two-dimensional RMSD (2D-RMSD) of the WT and F583A mutant trajectories is shown in Fig. 5. This two-dimensional RMSD matrix analysis reflects the structural changes of the proteins (Fig. 5A) and of helices 10–11 solely (Fig. 5B) along the simulation time and assesses their stability. The backbone atoms RMSD of the F583A mutant trajectory is presented in relation to the backbone atoms RMSD of the WT protein trajectory. A color code is used for visualizing the RMSD values, colored from blue to red covering the range from 0.3 to 1.22 nm (Fig. 5). The color pattern exhibits progressive yet reversible conformational changes, indicating that both structures are fluctuating and the similarity between them varies.
The mutual evolution of the conformational changes is best manifested when following the diagonal line divided into three phases: a relaxation phase (from ~ 20–30 ns and ~ 20–40 ns for the WT and the mutant, respectively), a progression phase (~ 30–80 ns and ~ 40–80 ns for the WT and the mutant, respectively), and a quiescent phase (from ~ 80 ns until the end of the simulations). The relaxation phase is represented by the dark blue, pale blue, and greenish colors stretching from the bottom of the figure. At this phase, the MRP1 protein at both simulations responds to the absence of the packing forces present at the crystal lattice. Both WT and mutant proteins relax at native physiological membrane conditions. At the progression phase, the F583A mutant undergoes conformational changes (shown as ~ 40–80 ns at the ordinate), rendering it to be less similar to the evolving WT protein which is relatively stable. Consequently, both conformations become more and more distinct one from the other, as indicated by the yellowish-reddish hues. Finally, at the quiescent phase, the mutant MRP1 continues to evolve, interestingly becoming more resembled the WT protein that keeps stable. The higher degree of similarity between the WT and mutant conformations, as reflected by the smaller RMSD values one with respect to the other, is observed at this phase (represented by the greenish hue).
Notably, as the simulations progress, the WT protein experiences limited changes and is more stable, whereas the mutant MRP1 fluctuates between conformations, as demonstrated by our further analyses. The noticeable similarity between the two matrices, the whole protein (Fig. 5A) and helices 10–11 solely (Fig. 5B), suggests that the progression of the protein is mainly shaped by the progression of helices 10–11. It is of interest to point out that the structures of the protein at both simulations do not co-evolve in parallel simultaneously at the same time and hence the two halves of the matrix are not identical albeit similar.
A Change in the Position and Orientation of Helix 10. A close inspection of how the WT protein and the F583A mutant conformations evolve throughout the simulations reveals a high similarity; in both cases, the transmembrane α-helices and the angles between are retained and kept, except from the angle between helices 10 and 11 and the angularity of helix 10. Hence, to account for the difference between the WT protein and the F583A mutant we measured the interhelical angle between helices 10 and 11. We found that for the F583A mutant, the average angle (82.2 ± 2 º) is smaller than the average angle of the WT (85.4 ± 2 º), as measured during the trajectories of the three independent simulations (Fig. 6A). Notably, the distribution of the angle values of the WT is more flattened and is located in close vicinity to the right of the mutant distribution, meaning the WT is more opened to the cytoplasm than the mutant (Fig. 6B). The two distributions showed significant difference (p-value is < .00001). When looking at the most common cluster, one can detect the difference in the degree of openness of the angle between helices 10–11 showing a 5 º difference favoring the WT (Fig. 6C). To reveal the cause for the change in the interhelical angle we performed structural analyses of helices 10 and 11. Since helix 10 is composed of two rod-like structure elements, connected by W553 and V554, we followed the angularity of the helix. In the WT protein helix 10 bends at residues W553-V554 and provides a more open conformation, in comparison to the F583A mutant’s helix (Fig. 6D). The bending angle in helix 10 is 19 ± 3 º for W553 and 18 ± 2 º for V554, while in the F583A mutant the angles are only 11 ± 1 º and 13 ± 0.4 º, respectively (Fig. 6E). The difference in the bending of helix 10 propagated to the edge of helix 10, changing the docking orientation of the helix on NBD2. Both the changes in the opening angle and the bending of helix 10 support our macroscopic results, showing relatively closure of the transport path of the mutant when compared to the WT.
Secondary Structure Motifs in NBD2. To account for the secondary structure differences in NBD2 between the WT and the F583A mutant DSSP analysis for this domain was performed (Fig. 7A). In general, the main differences between the NBD of the WT and the mutant were in the three common secondary structures, α-helices, β-sheets, and turns. The number of residues forming the α-helices and β-sheets is larger in the WT structure when compared to the F583A mutant. There is a difference of five residues and one residue for the formation of α-helices and β-sheets, respectively. By contrast, three other residues in the mutant take part in the formation of turns, causing changes in the direction of the polypeptide chain.
NBD2 consists of conserved motifs that are essential for nucleotide binding, hydrolysis, and NBD-NBD as well as NBD-TMD communication. Zooming in, secondary structural differences in four of those motifs were detected in the F583A mutant. (1) Walker A – unraveling of a helix in residues K1333-S1334, (2) s5/h2 loop – unraveling of an α-helix in residues L1363-K1369, (3) Q-loop – unraveling of a β-strand in residues I1372-P1374 that are located in the starting point of the loop, (4) Walker B - unraveling of a β-strand in residues L1453-D1454 (Fig. 7B). Walker A is responsible for ATP binding (phosphates binding). The conserved K1333 forms H-bonds with the oxygen atoms of the α- and β-phosphates, thereby holding both phosphates in a defined orientation. S1334 forms an H-bond with a D1454 residue of Walker B and interacts with the magnesium ion and the β-phosphate of the bound ATP 16. The loss of the α-helix may interrupt the formation of H-bonds with the incoming ATP phosphates. Walker B takes part in ATP hydrolysis. The changes in the secondary structure, the unraveling of the β-strand, deviating from the conserved secondary structure element, may interrupt the hydrolysis mechanism. The F583A mutation also unravels the β-strand element of the Q-loop, that couples TMD-NBD communication 5,6. Moreover, the most significant change was in the s5/h2 loop which is part of the TMD-NBD2 interface, where the cytoplasmic segment of helix 10 from the TMD1 settles into a socket on the NBD2 surface. It was previously shown that a main secondary structural difference between the degenerate ATPase site (NBD1) to the consensus ATPase site (NBD2) is the lack of s5/h2 element in the degenerated site 3. Moreover, the absence of s5/h2 in the TMD-NBD1 interface renders it to be weaker than the TMD-NBD2 interface. Here, we show the absence of the α-helix h2 and suggest that similarly to NBD1, the lack of this helix impairs the strength of the interaction between TMD-NBD2. Overall, we demonstrate that the F583A mutation causes major changes in critical elements for proper ATPase activity. F583 mediates long-range signaling that ensures interactions between the functional parts of the protein 17,18.
The Interaction between Helix 10 and NBD2. Differences between the WT and the F583A mutant were located both at the macroscopic level and the microscopic level, specifically in helices 10–11 and NBD2. The missing part in the puzzle is how the mutation propagated from the transmembrane helices to the nucleotide-binding site. To answer this question, we scanned for key residues in the communication between the two parts of the protein. We located two residues in the s5/h2 loop, H1364 and R1367, and one residue in helix 10, E521, that together form the communication path (see methods, energy calculations). In the WT protein, the distance between H1364 to E521 stabilizes on an average value of 2 Å while in the mutant it is shakier around an average value of 6 Å. The distance between R1367 and E521 stabilizes on an average value of 3 Å in the WT and 4 Å in the mutant (Fig. 8A).
Next, we performed an interaction energy calculation to analyze the nature of the NBD2-TMD communication (Fig. 8B). For E521-R1367 a salt bridge bond exists through the entire duration of the WT’s and the F583A mutant’s simulations. However, there is a difference in the strength of the bond, presenting an advantage to the WT over the mutant, with an average value of -17 ± 2 kJ/mol and − 12 ± 2 kJ/mol, respectively. For H1364 a bond formed in the WT on an average value of -8 ± 2 kJ/mol while in the F583A mutant no bond is established. The energies reflect the nature of the bonds forming between the residues. In the WT protein, the ionizable side chain of the positively charged H1364 and the guanidinium of the positively charged R1367 from s5/h2 are forming a salt bridge with the anionic carboxylate of the negatively charged E521 from the intracellular section of helix 10. As mentioned above, in the WT protein the average distance of H1364 and R1367 from E521 is 2.5 Å and 3 Å, respectively, suitable for the formation of two salt bridges. However, in the F583A mutant, the distance of H1364 and R1367 from E521 is on average 4 Å, enough for preventing the formation of a stable salt bridge as amino acids in a distance greater than 4 Å do not qualify to form a salt bridge 19.
Having established the putative allosteric effect of the mutation, we further investigated the dihedrals of the key residues. Dihedrals of H1364, R1367 and E521 were calculated, both for the WT and the mutant protein. H1364 and E521 showed a significant difference in the distribution of the dihedral angle of the WT and the F583A mutant protein (p-value = .00374 and .00001, respectively). R1367 showed no difference (p-value = .4413) (Fig. 8C).
Moreover, a comparison of the most common conformation during the simulation (cluster #1) of the WT and the mutant reveals profound differences. H1364 and R1367 rearrange differently in the protein’s conformational space in the mutant relative to the WT, thus become further away from E521 (Fig. 8D). In The WT, the distances between the groups forming the salt bridge of H1364-E521 and R1367-E521 are 2.4 Å and 3.2 Å, respectively, while in the mutant the distances increase to 5.9 Å and 3.6 Å, respectively, and therefore only one salt bridge is formed. Our findings are supported by an additional methodology, an elastic network model (ENM), presented in a recent study showing that E521 together with H1364 and R1367 are involved in preserving the completeness of the TMD-NBD2 interface and are important for the hMRP1 transport 20. Considering the differences in the distances and dihedral angles of H1364-E521 and R1367-E521 it is possible to conclude that the optimal terms for the formation of a communication between the transmembrane region and the ATPase region of the protein are not taking place in the mutant protein.