Extracellular Mutation Induces an Allosteric Effect Across the Membrane and Hampers the Activity of MRP1 (ABCC1)

doi:10.21203/rs.3.rs-289074/v1

PDF

Research Article

Extracellular Mutation Induces an Allosteric Effect Across the Membrane and Hampers the Activity of MRP1 (ABCC1)

https://doi.org/10.21203/rs.3.rs-289074/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Dynamic conformational changes play a major role in the function of proteins, including the ATP-Binding Cassette (ABC) transporters. Multidrug Resistance Protein 1 (MRP1) is an ABC exporter that protects tissues from toxic molecules. Overexpression of MRP1 has been shown to confer Multidrug Resistance (MDR), a phenomenon in which cancer cells are capable to defend themselves against a broad variety of drugs. Despite an increasing number of structures of MRP1, including a relatively new bovine cryo-EM structure, the accurate molecular details for its drug extrusion mechanism remain vague. In this study, we used varied computational techniques to explore the unique F583A mutation that is known to essentially lock the transporter in a low-affinity solute binding state. We demonstrate how macro-scale conformational changes affect MRP1’s stability and dynamics, and how these changes correspond to micro-scale structural perturbations in helices 10–11 and the nucleotide-binding domains (NBDs) of the protein in regions known to be crucial for its ATPase activity. We demonstrate how a single substitution of an outward-facing aromatic amino acid causes a long-range allosteric effect that propagates across the membrane, ranging from the extracellular ECL5 loop to the cytoplasmic NBD2 over a distance of nearly 75 Å, leaving the protein in a non-functional state, and provide the putative allosteric pathway. The identified allosteric structural pathway is not only in agreement with experimental data but enhances our mechanical understanding of MRP1, thereby facilitating the rational design of chemosensitizers toward the success of chemotherapy treatments.

Neurobiology of Disease

General Biochemistry

Biophysics

ATP-Binding Cassette (ABC)

MRP1

NBD2

cryo-EM

Multidrug resistance-associated protein 1 (MRP1/ABCC1) is a member of the superfamily of ATP-binding cassette (ABC) transporters. In healthy tissues, MRP1 fulfills important roles such as efflux of xenobiotics and endogenous metabolites and transport of inflammatory mediators ¹. However, in cancer cells it is overexpressed, resulting in resistance to anti-cancer drugs. It is especially prevalent in neuroblastoma and cancer cells found in the lung, breast, and prostate ². The 1,531 amino acids human transporter is composed of three transmembrane domains (TMD0, TMD1, and TMD2), two nuclear binding domains (NBD1 and NBD2), and an interfacial domain known as the Lasso motif ³ (Fig. 1A). TMD1 and TMD2 are composed of six helices each and together form the transport path. NBD1 and NBD2 are structurally and functionally nonequivalent, with only one competent ATPase site ⁴. It has been demonstrated that while both ATPase sites bind ATP, hydrolysis occurs only at NBD2. Each of the ATPase sites consists of Walker A and B motifs from one NBD and the signature motif from the other (Fig. 1B). Walker A and the signature motif are responsible for ATP binding while Walker B has a role in ATP hydrolysis. In addition, the Q-loop is known for mediating TMD-NBD interactions ^5,6. Another noticeable element is the s5/h2 loop that is part of the TMD-NBD2 interface, where a cytoplasmic helix from the TMD settles into a socket on the NBD surface (Fig. 1B) ³.

In this study we used the inward-facing state of the apo-protein, in which the protein is open to the cytoplasm, waiting for ligand recruitment (PDB: 5uj9), for investigating how a single mutation, F583A, affects the dynamics of the protein. In the initial resting state, the NBDs are widely separated, and the translocation pathway is continuous with the cytoplasm, ready for ligand binding ³. As previously described in the Cryo-EM studies, MRP1 shifts between multiple inward-facing conformations in the absence of ligand ⁷. Ligand binding induces movement of the TMDs towards each other, to form a high-affinity ligand-binding site ³. The TMDs movement brings the NBDs closer and aligns them for proper future dimerization. NBD1 binds one ATP molecule and then NBD2 binds the second one, forming an NBD-dimer. The dimerization initiates a further conformational change of the protein to an outward-bound-facing state. Consequently, the affinity of the ligand reduces, and it is released to the extracellular side of the membrane. In the final step of this process, ATP hydrolysis occurs in NBD2 followed by dissociation of ATP (from NBD1) and ADP and Pi (from NBD2) and the protein returns to its initial inward-facing state ⁴.

While the schematic formalism of the translocation is clear, the details of the reaction mechanism are yet to be elucidated. Our previous experimental data showed that the point mutation F583A essentially eliminates organic anion transport and impairs the release of ADP, thus effectively locks the transporter in a low-affinity binding state. F583 is located in the extracellular side of the protein and is at the C-terminus end of ECL5 as it turns into helix 11 (Fig. 1C). It is the first outward-facing amino acid in MRP1 found to be critical for the protein’s transport. The F583A mutation disrupts MRP1’s function through a long-range allosteric effect (over 75 Å) ranging from the ECL5 loop to the cytoplasmic NBDs ⁸. It is still unclear how the mutation affects the intracellular NBDs by propagating across the membrane.

Allostery, a well-known notion in biology, largely involves conformational and functional alterations in individual proteins. Because of the paramount role allostery takes place in the function and dynamics of proteins, it is not surprising it is associated with pathologies and diseases. Although an allosteric effect is mostly associated with effector molecules, it can also be impaired through mutation, in the regulatory site, and even in residues spread along the protein. One of the most used methods to study the allosteric-communication pathways is molecular dynamics simulations. Examining the dynamics of the system enables mapping the potential allosteric site and its communication with the functional site ⁹.

Herein, we tackled previously unanswered fundamental questions regarding the protein’s mechanism. We explain how a single mutation can substantially decrease the transport activity of MRP1, on an atomic resolution scale. We suggest that the exofacial F583A mutation induces a long-range allosteric effect ranging from the ECL5 loop to the cytoplasmic NBD2, through helix 10. The mutation unravels the secondary structure of the TMD/NBD2 interface, thus disrupting the signaling between the two functional parts of the protein. Moreover, we show unrevealing of the well-known Walker A, Walker B, and Q-loop elements, that was previously shown to be critical for ATP binding and hydrolysis as well as in the s5/h2 loop element. As a consequence, the appropriate conditions for hydrolysis may not occur and the protein will be stuck in an inactive state. Our findings exemplify how an allosteric effect affects the protein dynamics, consequently hampering its translocation activity.

Homology Modeling. Homology modeling of membrane proteins, such as MRP1, is similar to that of water-soluble proteins, and modeling membrane proteins results are as adequate as those obtained when modeling their water-soluble counterparts, for example ^10,11. Given that the human and bovine MRP1 share a high sequence identity (91 %) and similarity (99 %) (E-value lower than approximately 5∙10^− 324), the latter was selected as a template for the homology modeling process. This sequence resemblance between the proteins ensures a high success of the modeling process.

Model Quality Assessment. Before the MD simulations, the human MRP1 model structure was evaluated for its quality. The Ramachandran plot showed that 99 % of the residues lie in allowed regions, from which 87 % in the most favored regions. No major stereochemical clashes or bad contacts for main-chain or side-chain parameters were detected. All bond lengths, backbone, and rotamer angles were in good agreement with standard values. Thus, the resulted model of the homology modeling procedure was a relatively reliable structure.

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 ¹², 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 ¹³.

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 ¹⁴. 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 ¹⁵. 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 nm²) in the first eigenvector compared with the WT protein (116 ± 7 nm²), 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 ¹⁶. 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 ³. 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 ¹⁹.

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 ²⁰. 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.

ABC transporters are one of the largest and oldest protein superfamilies found in all phyla taxa, coupling ATP binding and hydrolysis with translocation of a wide array of substrates across the membrane. Although a vast amount of structural and biochemical data was accumulated about their mechanism of action (e.g., ^3,4,8,21), a few questions are yet to be answered, for example, how does a specific point mutation exert its influence on these transporters? To tackle that question, we extensively studied the effect of the nonconservative substitution F583A, located on the exofacial segment of MRP1, on the activity and structural dynamics of this transporter. Previously, it was experimentally shown that the F583A mutation causes a decrease in MRP1 transport activity by impairing the release of ADP and locking the transporter in a low-affinity solute binding state ⁸. Thus, rather than merely affecting the substrate selectivity of MRP1, the NBD’s interaction of the protein was substantially altered as well. These observations were highly surprising because the exofacial location of F583 at ECL5 of MRP1 is ~ 75 Å far-off from the NBDs. However, the mechanism of the F583A mutation remained vague and elusive, and it was not known how the loss of an aromatic side chain at position 583 disrupts the ATPase activity of MRP1, even in the absence of substrate.

Since neither the current structures and models of MRP1 nor the biochemical and mutagenesis studies offer substantive mechanistic insights with respect to how the loss of the exofacial F583 alters the nucleotide interactions of MRP1, further in-depth investigation was needed to understand precisely how the mutation effect is transmitted across the membrane. It was still to be determined how long-range interactions ranging from a single amino acid span the two bundles of helices (TMD1 and TMD2) of MRP1 and affect the two NBDs. To that end, herein we used modeling and MD simulations, combined with thorough macroscopic and microscopic-scale analyses, and successfully tracked the mechanistic progression of the mutation from the ECL5 loop, via the helices-bundles, up to the NBDs. To the best of our knowledge, this is the first MD study that is based on the newly determined bovine MRP1 cryo-EM structure ³.

We offer an atomic-scale resolution explanation for how the loss of an aromatic side chain at position 583 predisposes MRP1 to favor a conformation with lower substrate affinity. This conformation also disrupts the functional cooperativity between the NBDs, such that post hydrolysis release of ADP is impaired. Thus, F583 appears to be critical for the transition of MRP1 from its low-affinity substrate-binding state to a high-affinity binding state during its transport cycle.

Allostery is an intrinsic characteristic of dynamic macromolecules, representing a phenomenon in which a functional or structural change occurs at a given location while it is generated by a perturbation produced by a distant site. Based on our findings we were able to follow the step-by-step structural pathway of the mutation effect, discovering a hidden allosteric pathway embedded in MRP1, and revealing key cryptic residues. This pathway (Fig. 9) is evoked by the F583A mutation that triggers small conformational changes, that in-turn transmit the allosteric effect in a wave-like manner along with structural elements and specific amino acids.

The F583A mutation (Fig. 9, circle 1) induces a curvature change in helix 10’s angularity. Since helix 10 is composed of two rod-like elements (speared by W553 and V554) and the angle between them is not rigid, the mutation affects the opening of the angle, making it more obtuse and consequently, the helix becomes straighter (Fig. 9, circle 2). This generates a change to the helix’s orientation in relation to NBD2 (Fig. 9, circle 3). As a result, a disruption in the TMD-NBD2 interface communication is evoked by influencing the orientations of residues E521 (helix 10), H1364, and R1367 (NBD2). Lastly, important secondary structure elements in NBD2 which are responsible for ATP binding and hydrolysis are unraveled (4). The TMD-NBD interaction interface is influenced by the mutation, and proper communication between the TMDs and the NBDs is broken. Proper binding of ATP is hindered, and therefore the hydrolysis mechanism is possibly damaged. This exemplifies how the allosteric signal travels from the ECL5 loop to the NBDs and transmits its effect across the membrane. Having portrayed the entire putative embedded structural allosteric pathway that propagates across the membrane to a distance of 75 Å enabled us to identify the chronological sequence of the molecular events. Notably, this long-range coupling chain of effects occurs at the nanosecond timescale, and MD simulations appeared to be the appropriate tool for its discovery.

Several key residues that play a significant role by mediating the allosteric propagation were identified. These residues can be classified into two groups based on their positions: E521, W553, and V554 located at helix 10, and H1364 and R1367 located at NBD2. Not surprisingly, most of these residues are known to be involved in MRP1’s function and expression: E521A mutant is characterized by a significant decrease at MRP1 level, altered kinetic interactions with ATP, and shows 50% deficiency of transport activity ²²; W553 has a role in controlling the affinity of LTC4 to MRP1²³; R1367 is necessary for stable MRP1 expression, H1364A and R1367A mutations resulted in reduced MRP1 levels ²⁴.

In this study, we harnessed the abilities of computational modeling and MD simulations to investigate a well-studied protein. We were able to explain the molecular structural mechanism of the F583A mutation effect by discovering a long-range conformational transmission across the membrane and an unknown putative allosteric pathway in MRP1.

Our work does not only provide a mechanistic structure-function relationship framework for MRP1 but may contribute valuable insights for the design of specific chemosensitizers and chemotherapeutic agents for MRP1. From a broader perspective, we showed that the effect of point mutations at membrane proteins, and specifically ABC transporters, can be tracked and portrayed over a long-range of tenths of angstroms by computational means, and identified the allosteric hotspots along the pathway. We hope that our study will promote further works that couple computational mechanistic molecular explanations with biochemical data.

Homology Modeling. The template used to model the human MRP1 (UniProt: P33527) was the cryo-EM bovine MRP1 (PDB: 5uj9). The high sequence identity (91%) ensures a high success of the modeling process. The human and the bovine amino acid sequences were aligned using MAFFT ²⁵. Homology models were generated with Modeller 9.4 ²⁶. The model with the best score was assessed for its quality concerning its energy and stereochemical geometry using Procheck 3.5 ²⁷, ERRAT ²⁸, and WHATCHECK ²⁹.

Systems Setup. The model with the best score was inserted into a preequilibrated 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine (POPC) bilayer (adopted and modified from ³⁰). The insertion was done so that the protein’s rough axis was perpendicular to the membrane surface plane, and all colliding lipids and water molecules, within 2 Å of the protein, were removed. The lipids membrane was scaled to 448 lipids by GROMACS tools to fit the protein size ³¹. The simulation box was filled with water using the SPC216 model and the system’s total charge was neutralized by adding Na⁺ and Cl⁻ ions to a final concentration of 0.1 M. Then, the apo human MRP1 model was embedded in a membrane by the g_membed package in GROMACS ³², removing colliding lipids and water molecules. To construct the F583A mutant the Wizard > Mutagenesis option in PyMOL was used ³³.

MD Details. All MD simulations were performed using GROMACS 5.1.4 package (http://www.gromacs.org/) ³¹ with AMBER99sb-ildn force field ³⁴. A total of six systems were created, three for WT apo-protein and three for mutant F583A apo-protein. The systems were subjected to rigorous energy minimization using the steepest descent algorithm with a sequential decreasing tolerance, from 1000 to 200 kJ‧mol^− 1‧nm^− 1. Then, the systems were equilibrated under NPT for 500 ps for position restraint, allowing the solvent to equilibrate around the protein without disturbing the protein structure. Each system was simulated further for 10 ns of unconstrained equilibrium simulation, using Berendsen’s coupling algorithm ³⁵ for keeping the pressure constant. Afterward, the systems were submitted for unbiased MD production runs of 100 ns each. All simulations were conducted using the LINCS algorithm ³⁶ to constrain bond lengths and angles of hydrogen atoms, allowing a time step of 2 fs. The pressure was kept constant using the Parrinello-Rahman coupling algorithm at 1 bar by applying a semi-isotropic coupling constant of τ = 1 ps. Each simulation had a different initial velocity for every atom that was randomly generated from the Maxwell-Boltzmann distribution at 310 K employing a coupling constant of τ = 0.1 ps. A cutoff of 1.2 nm was used for van der Waals interactions, and long-range electrostatic interactions were computed using the PME method ³⁷.

Visualization and Analyses. The analyses were conducted using in-house purpose-written python scripts and the GROMACS analysis package utilities ³⁸. Three WT trajectories and three F583A mutant trajectories were used for our analyses, to ensure statistical significance. All snapshots were prepared using PyMOL ³³.

Cluster Analysis. A geometric clustering algorithm was performed to identify similar structures, based on the respective RMSD between all conformations sampled during the MD simulation ³⁹. A 3 Å predefined RMSD cut-off was chosen to determine the cluster group, in order to generate a sufficient number of clusters while maintaining statistical significance, e.g., ^40,41. The cluster number for each frame was calculated, summarized to a table consisting of the number of frames for each cluster. The average of the three WT simulations and the F583A mutant’s simulations were calculated separately. Then, the corresponding conformation was determined.

Evolutionary Conservation Score. Conservation analysis was performed using the Consurf server ⁴² and multiple sequence alignment by CLUSTALW ⁴³.

Principal Component Analysis (PCA). The covariance matrix was calculated, diagonalized, and analyzed by GROMACS utility tools. The corresponding eigenvalues are plotted against the first five corresponding eigenvector indices obtained from the Cα covariance matrix constructed from the MD trajectory. Visualization projection of the obtained results is shown by arrows representing the principal components of the protein’s motion throughout the simulations using a PyMOL-designed script.

RMSD Matrix. RMSD matrix was produced in two-dimensional format to compare values of the WT structure in the trajectory with respect to the F583A mutant structure.

Interhelical Opening Angle Between Helices 10 and 11. Interhelical angles were calculated according to the following procedure: One residue was selected at the C-terminal of helix 10 (I571), one residue at the N-terminal of helix 11 (I577), and one residue at the ECL5 loop (N574). The Cα atoms’ center of mass of the residues was determined. A vector-defining helix was drawn between these centers of mass along each helix. The interhelical angle was calculated between every two successive vectors.

DSSP (Define Secondary Structure of Proteins) Algorithm. Secondary structure content was determined with the STRIDE algorithm ⁴⁴.

Energy Calculations. To discover key residues in the communication between helix 10 and NBD2 we filtered residues in the upper side of NBD2 that are located within an average distance shorter than 3 Å from the lower part of helix 10, both in the WT and the F583A mutant. Energy calculations (vdW + Electrostatic) were done using the MM-PBSA method for GROMACS for each of the filtered residues ⁴⁵. Only two residues presented energy differences between the WT protein and the F583A mutant, H1364 and R1367. We chose to use a dielectric constant of ε = 20 because the residues we analyzed are located on the protein’s surface and electrically charged ⁴⁶.

Dihedral angles. Atoms defined for dihedral calculations were: Cα-Cβ-Cγ-Cδ2 for H1364, Cδ-Nε-Cζ-Nη2 for R1367 and Cβ-Cγ-Cδ-Oε1 for E521.

Acknowledgments

All authors acknowledge Prof. Suzan Cole (Queen's University, Canada) for seminal and inspiring studies of MRP1 over many years and fruitful collaboration. Y.B.K. acknowledges the financial support of the Teva Pharmaceuticals bio-innovation fellowship and Prof. Dan Peer (Tel-Aviv University, Israel).

Author Contributions

Y.B.K. performed the simulations and the analyses. Y.B.K. and A.G. processed the data, conceptualized and conceived the analytical methods, and drafted the manuscript. Y.T. supervised the study and was in charge of overall professional scientific direction and planning. All authors discussed the results, provided critical feedback, helped shape the research and analysis, and finalized the manuscript.

Competing interests

The authors declare no competing interests.

Li, Q., Kato, Y., Sai, Y., Imai, T. & Tsuji, A. Multidrug Resistance—Associated Protein 1 Functions as an Efflux Pump of Xenobiotics in the Skin. Pharm. Res. 22, 842–846 (2005).
Munoz, M., Henderson, M., Haber, M. & Norris, M. Role of the MRP1/ABCC1 Multidrug Transporter Protein in Cancer. IUBMB Life. 59, 752–757 (2007).
Johnson, Z. L. & Chen, J. Structural Basis of Substrate Recognition by the Multidrug Resistance Protein MRP1. Cell. 168, 1075–1085 (2017).
Johnson, Z. L. & Chen, J. A. T. P. Binding Enables Substrate Release from Multidrug Resistance Protein 1. Cell. 172, 81–89 (2018).
Wilkens, S. Structure and mechanism of ABC transporters. F1000prime reports. 7, 14 (2015).
Szöllősi, D., Rose-Sperling, D., Hellmich, U. A. & Stockner, T. Comparison of mechanistic transport cycle models of ABC exporters. Biochimica et Biophysica Acta (BBA) - Biomembranes. 1860, 818–832 (2018).
Wang, L. et al. Characterization of the kinetic cycle of an ABC transporter by single-molecule and cryo-EM analyses.eLife9, (2020).
Weigl, K. E. et al. An Outward-Facing Aromatic Amino Acid Is Crucial for Signaling between the Membrane-Spanning and Nucleotide-Binding Domains of Multidrug Resistance Protein 1 (MRP1; ABCC1). Mol. Pharmacol. 94, 1069–1078 (2018).
Protein Allostery in Drug Discovery. vol. 1163 (Springer Singapore, 2019).
Ganoth, A., Alhadeff, R. & Arkin, I. T. Computational study of the Na+/H + antiporter from Vibrio parahaemolyticus. Journal of Molecular Modeling. 17, 1877–1890 (2011).
Ganoth, A., Alhadeff, R., Kohen, D. & Arkin, I. T. Characterization of the Na+/H + Antiporter from Yersinia pestis. PLoS ONE. 6, e26115 (2011).
Capra, J. A. & Singh, M. Predicting functionally important residues from sequence conservation. Bioinformatics. 23, 1875–1882 (2007).
Poteete, A. R., Rennell, D. & Bouvier, S. E. Functional significance of conserved amino acid residues. Proteins: Structure, Function, and Genetics. 13, 38–40 (1992).
Franco, R. & Zavala-Flores, L. in Encyclopedia of Cancer 1–5 (Springer, Berlin Heidelberg, 2012). ABCC Transporters doi:10.1007/978-3-642-27841-9_7076-2.
Papaleo, E., Mereghetti, P., Fantucci, P., Grandori, R. & de Gioia, L. Free-energy landscape, principal component analysis, and structural clustering to identify representative conformations from molecular dynamics simulations: The myoglobin case. Journal of Molecular Graphics and Modelling. 27, 889–899 (2009).
Yang, R., Scavetta, R. & Chang, X. Interaction between the Bound Mg·ATP and the Walker A Serine Residue in NBD2 of Multidrug Resistance-Associated Protein MRP1 Plays a Crucial Role for the ATP-Dependent Leukotriene C4 Transport. Biochemistry. 47, 8456–8464 (2008).
Dhasmana, D., Singh, A., Shukla, R., Tripathi, T. & Garg, N. Targeting Nucleotide Binding Domain of Multidrug Resistance-associated Protein-1 (MRP1) for the Reversal of Multi Drug Resistance in Cancer. Sci. Rep. 8, 11973 (2018).
Szöllősi, D., Rose-Sperling, D., Hellmich, U. A. & Stockner, T. Comparison of mechanistic transport cycle models of ABC exporters. Biochimica et Biophysica Acta (BBA) - Biomembranes. 1860, 818–832 (2018).
Kumar, S. & Nussinov, R. Close-Range Electrostatic Interactions in Proteins. ChemBioChem. 3, 604–617 (2002).
He, J., Han, Z., Farooq, Q., ul, A. & Li, C. Study on functional sites in human multidrug resistance protein 1 (hMRP1). Proteins: Structure, Function, and Bioinformatics prot.26049(2021) doi:10.1002/prot.26049.
Amram, S. et al. Structural Characterization of the Drug Translocation Path of MRP1/ABCC1. Israel Journal of Chemistry. 54, 1382–1393 (2014).
Iram, S. H. & Cole, S. P. C. Mutation of Glu521 or Glu535 in Cytoplasmic Loop 5 Causes Differential Misfolding in Multiple Domains of Multidrug and Organic Anion Transporter MRP1 (ABCC1). Journal of Biological Chemistry. 287, 7543–7555 (2012).
Conseil, G., Arama-Chayoth, M., Tsfadia, Y. & Cole, S. P. C. Structure-guided probing of the leukotriene C 4 binding site in human multidrug resistance protein 1 (MRP1; ABCC1). The FASEB Journal. 33, 10692–10704 (2019).
Iram, S. H. & Cole, S. P. C. Expression and Function of Human MRP1 (ABCC1) Is Dependent on Amino Acids in Cytoplasmic Loop 5 and Its Interface with Nucleotide Binding Domain 2. Journal of Biological Chemistry. 286, 7202–7213 (2011).
Katoh, K., Rozewicki, J. & Yamada, K. D. MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization. Briefings in Bioinformatics. 20, 1160–1166 (2019).
Webb, B. & Sali, A. Comparative Protein Structure Modeling Using MODELLER.Current Protocols in Bioinformatics54, (2016).
Laskowski, R. A., MacArthur, M. W., Moss, D. S. & Thornton, J. M. PROCHECK: a program to check the stereochemical quality of protein structures. Journal of Applied Crystallography. 26, 283–291 (1993).
Colovos, C. & Yeates, T. O. Verification of protein structures: Patterns of nonbonded atomic interactions.Protein Science2, (1993).
Dunbrack, R. & Whatcheck in Dictionary of Bioinformatics and Computational Biology (John Wiley & Sons, Ltd, 2004). doi:10.1002/9780471650126.dob0791.pub2.
Janosi, L. & Gorfe, A. A. Simulating POPC and POPC/POPG Bilayers: Conserved Packing and Altered Surface Reactivity.Journal of Chemical Theory and Computation6, (2010).
van der Spoel, D. et al. GROMACS: Fast, flexible, and free.Journal of Computational Chemistry26, (2005).
Abraham, M. J. et al. High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 1–2. GROMACS, 19–25 (2015).
The PyMOL Molecular Graphics System, Version 2.0 Schrödinger, LLC.
Lindorff-Larsen, K. et al. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins: Structure, Function, and Bioinformatics. 78, 1950–1958 (2010).
Berendsen, H. J. C., Postma, J. P. M., van Gunsteren, W. F., DiNola, A. & Haak, J. R. Molecular dynamics with coupling to an external bath. The Journal of Chemical Physics. 81, 3684–3690 (1984).
Hess, B. P. L. I. N. C. S. A Parallel Linear Constraint Solver for Molecular Simulation. Journal of Chemical Theory and Computation. 4, 116–122 (2008).
Essmann, U. et al. A smooth particle mesh Ewald method. The Journal of Chemical Physics. 103, 8577–8593 (1995).
van der Spoel, D. et al. Fast, flexible, and free. Journal of Computational Chemistry 26. GROMACS, 1701–1718 (2005).
Keller, B., Daura, X. & van Gunsteren, W. F. Comparing geometric and kinetic cluster algorithms for molecular simulation data. The Journal of Chemical Physics. 132, 074110 (2010).
Amber-Vitos, O., Chaturvedi, N., Nachliel, E., Gutman, M. & Tsfadia, Y. The effect of regulating molecules on the structure of the PPAR-RXR complex. Biochimica et Biophysica Acta (BBA) - Molecular and Cell Biology of Lipids. 1861, 1852–1863 (2016).
Amber-Vitos, O., Kucherenko, N., Nachliel, E., Gutman, M. & Tsfadia, Y. The Interaction of FABP with Kapα. PLOS ONE. 10, e0132138 (2015).
Ashkenazy, H. et al. ConSurf 2016: an improved methodology to estimate and visualize evolutionary conservation in macromolecules. Nucleic Acids Res. 44, W344–W350 (2016).
Thompson, J. D., Gibson, T. J. & Higgins, D. G. Multiple Sequence Alignment Using ClustalW and ClustalX.Current Protocols in Bioinformatics00, (2003).
Heinig, M. & Frishman, D. STRIDE: a web server for secondary structure assignment from known atomic coordinates of proteins. Nucleic Acids Res. 32, W500–W502 (2004).
Homeyer, N. & Gohlke, H. Free Energy Calculations by the Molecular Mechanics Poisson – Boltzmann Surface Area Method. Mol. Inf. 31, 114–122 (2012).
Li, L., Li, C., Zhang, Z. & Alexov, E. On the Dielectric “Constant” of Proteins: Smooth Dielectric Function for Macromolecular Modeling and Its Implementation in DelPhi. Journal of Chemical Theory and Computation. 9, 2126–2136 (2013).

No competing interests reported.

PDF

Editorial decision: Major revision
05 Apr, 2021
Reviews received at journal
16 Mar, 2021
Reviewers agreed at journal
05 Mar, 2021
Reviewers invited by journal
04 Mar, 2021
Editor assigned by journal
03 Mar, 2021
Editor invited by journal
03 Mar, 2021
Submission checks completed at journal
03 Mar, 2021
First submitted to journal
01 Mar, 2021

You are reading this latest preprint version

Extracellular Mutation Induces an Allosteric Effect Across the Membrane and Hampers the Activity of MRP1 (ABCC1)

Status:

Version 1

Abstract

Figures

Introduction

Results

Macroscopic View: Structural Characterization of the Protein’s Global Dynamics

Microscopic view: Structural characterization of helices 10–11 and NBD2

Discussion

Materials And Methods

Declarations

References

Additional Declarations

Status:

Version 1