Molecular rationale for the impairment of the MexAB-OprM efflux pump by a single mutation in MexA.

Efflux pumps of the Resistance-Nodulation-cell Division (RND) superfamily contribute to intrinsic and acquired resistance in Gram-negative pathogens by expelling chemically unrelated antibiotics with high efficiency. They are tripartite systems constituted by an inner-membrane-anchored transporter, an outer membrane factor protein, and a membrane fusion protein. Multimerization of the membrane fusion protein is an essential prerequisite for full functionality of these efflux pumps. In this work, we employed complementary computational techniques to investigate the stability of a dimeric unit of MexA (the membrane fusion protein of the MexAB-OprM RND efflux pump of Pseudomonas aeruginosa), and to provide a molecular rationale for the effect of the G72S substitution, which affects MexAB-OprM functionality by impairing the assembly of MexA. Our findings indicate that: i) dimers of this protein are stable in multiple µs-long molecular dynamics simulations; ii) the mutation drastically alters the conformational equilibrium of MexA, favouring a collapsed conformation that is unlikely to form dimers or higher order assemblies. Unveiling the mechanistic aspects underlying large conformational distortions induced by minor sequence changes is informative to efforts at interfering with the activity of this elusive bacterial weapon. In this respect, our work further confirms how molecular simulations can give important contribution and useful insights to characterize the mechanism of highly complex biological systems.

Introduction

Bacterial Multi-Drug Resistance (MDR), i.e., the tolerance of bacteria to several antibiotics belonging to different classes, has become one of the most serious threats to human health [1], [2], [3], [4], [5]. The multidrug efflux pumps of the Resistance-Nodulation-cell Division (RND) superfamily [6], [7], [8] confer intrinsic and acquired MDR in Gram-negative pathogens [9], [10], [11], [12], [13]. These pumps are tripartite machineries spanning the whole periplasmic space [14], [15], [16], [17], [18], [19], [20], [21], [22], and are composed of: i) a secondary RND antiporter embedded in the inner membrane (IM) [23], [24], [25], [26], [27], [28], which is the engine of the pump and is responsible of drug recognition and selectivity [29], [30], [31], [32], [33], [34]; ii) a channel embedded in the outer membrane (OM) [35], [36], called outer membrane factor (OMF), through which toxic molecules are extruded out of the bacterium; and iii) six membrane fusion proteins (MFP), linking the IM component to the OM one and building up a trans-envelope efflux system [37], [15], [16], [38], [39], [19], [40], [41], [42], [43], [20], [44].

Biochemical and structural studies demonstrated that functioning efflux pumps such as AcrAB-TolC and MexAB-OprM (constitutively expressed in Escherichia coli and Pseudomonas aeruginosa, respectively) assemble with stoichiometry 3(RND transporter):6(MFP):3(OMF) [15], [16], [38], [19], [43], [45], [20], [44]. Nevertheless, it is still unclear if the minimal functional oligomerisation unit of the MFP component is constituted by a monomer or by a dimer. In vivo and in vitro studies suggested a dimeric structure as the minimal unit assembling into a trimer of dimers in the whole efflux pump [44], [22].

However, some doubts have been cast that such an oligomerization might be an artifact of experimental conditions, since hydrodynamic and size exclusion chromatography studies found that only monomers were present in solution [46], [28]. Another study showed that the degree of oligomerization of MFPs and the stoichiometry of the MFP-transporter-OM channel complex might depend on the experimental condition e.g., the pH value [47].

In this work we focused on the MFP MexA, part of the major RND pump MexAB-OprM of P. aeruginosa [37], [48], [40], [51], [52], [53], [49], [50], one of the critical and high-priority Gram-negative pathogens. The monomeric structure of MexA as determined from X-ray crystallography [52] features an elongated shape in which four domains are connected by semi-flexible linkers (Fig. 1). Several cryo-EM structures confirmed an elongated morphology of MexA in the whole tripartite pump, where a hexameric arrangement of these proteins bridges the gap between MexB and OprM as a part of the duct formed by the MexB-MexA-OprM complex [15], [17], [54], [38], [43], [20].

Fig. 1

X-ray crystal structure of the MexA protein (PDB ID: 2V4D [52], chain L, residues 36 to 362 according the full Uniprot sequence numbering employed by Nehme et al. [51], or 13 to 339 according to the sequence numbering employed in the X-ray crystal structure numbering. The front view of the protein is shown on the left. The four domains of the protein are highlighted with different colors: α-hairpin (red; residues 96–158), lipoyl (green; residues 61–95, 159–194), β-barrel (blue; residues 50–60, 195–285) and Membrane Proximal (MP, mauve; residues 36–49, 286–362). The black sphere in the lipoyl domain identifies the mutation G72S considered in this work. A Zooming on the lipoyl domain is repoerted on the right side of the picture.

The flexibility of MexA, has been suggested to be crucial for the assembly and overall dynamics of the pump, particularly at the hinge between the α-helical hairpin and lipoyl domain [40], [55], [22]. Among the various mutations in this protein that compromise the functionality of MexAB-OprM, Poole and co-workers [51] identified only one impairing MexA-MexA interactions but not interfereing with MexA-MexB association. Such a mutation, G72S, occurs in the lipoyl domain (at the interface with the α-hairpin domain, see Fig. 1), leading to the hypothesis that the native structure of the protein, probably altered by the substitution, could be necessary for MexA-MexA association. Note that, in this work, we followed the residue numbering used by Nehme et al. [51], who employed the full protein sequence deposited in the Uniprot databank (). This sequence includes 23 residues from the periplasmic signal sequence at the N terminus, which were removed by Symmons et al. [52] to obtain the X-ray structure 2V4D (thus, the mutated residue becomes G49 in this structure). Importantly, both the presence of a glycine at position 72 and the key involvement of the lipoyl domain in MexA multimerization are conserved features across the MFP family [16], [38], [19], [49], [51], [43], [44], [22]; for instance, the lipoyl domains of the MFP AcrA of E. coli are key to its oligomerization [16], [44]. Moreover, the disruption of the very same domain in DevB (MFP associated to a bacterial ABC efflux pump) hampers MFP hexamerization [56].

In this study, by employing complementary computational methods, we assessed the propensity of the wild-type (WT) and of the G72S variant of MexA to dimerize in solution in the absence of MexB and OprM. We found that dimers of WT MexA can form with relatively high stability, supporting the hypothesis that at least the first stage of their oligomerization can occur independently of interactions with partner proteins. In particular, the resulting dimer retains an elongated shape and a high flexibility in the peripheral α-hairpin and the MP domains, which are both crucial to form stable and functional tripartite assemblies together with the MexA cognate components. On the contrary, the G72S substitution was found to induce dramatic changes in MexA structure and dynamics, significantly reducing the propensity of the protein to dimerize and ultimately hampering the formation of a stable tripartite complex. Our findings support the hypothesis that is the dimer is the minimal multimerization unit of MexA and highlight the pivotal role of the lipoyl domain in the assembly process. Furthermore we provide a molecular rationale for the available experimental data [49] and important insights into the possible drug targeting of MFPs for more effective inhibition of the major efflux system of P. aeruginosa [57], [58], [59].

Systems and methods

Choice of the monomeric structure of MexA

The re-refined crystal structure of the WT MexA protein (PDB ID: 2V4D [52], containing residues 13 to 339 (36 to 362 according to the numbering in [51] out of 360 aminoacids reported in the sequence), was the most recent and reliable structure available at the time this investigation began. In this crystallographic structure, twelve monomers of MexA are arranged into a bi-hexameric structure featuring two main protein conformations that differ in the orientation of the MP domain relatively to the rest of the system. The so-called unrotated state (e.g., chain L in 2V4D) features an ordered β-ribbon linker between the β-barrel and the MP domains, which gets distorted in the so-called rotated state (e.g., chain M in 2V4D). As a result, the MP domain twists by ∼85° relatively to the rest of the protein (Figure S1A). The unrotated state has been suggested to represent the preferred resting state of the protein in the absence of stabilizing contacts from other copies of MexA, thus it was used as reference structure to build homology models of the G72S variant of MexA (hereafter MexAG72S, see Fig. 1). The recently published cryo-EM structure of the full MexAB-OprM assembly [43] further supported our choice of the L chain as template. Indeed, chain L has an overall root-mean-square displacement (RMSD) lower than chain M from the monomeric MexA protein in the structure of the full assembly (PDB ID: 6IOK, see Figures S1B, C, E), although the overall RMSD difference between the conformations of chains L and M reduces from 7.1 Å in the X-ray structure (2V4D) to 3.5 Å in the cryo-EM one (6IOK). Only the β-barrel domain displays a significantly different arrangement between the structures of the L chain from PDB IDs 2V4D and 6IOK (Figures S1D, E).

Homology modelling

The monomeric unit of MexAG72S was built by structure-based homology modelling using the software MODELLER 9.16 [60] and assuming as template the chain L of the re-refined crystal structure of the WT MexA protein (PDB ID: 2V4D [52]. As for the WT protein, 326 residues were modelled in the G72S variant. Fifty models of the mutant protein were generated, from which two structures were chosen for further studies. The first structure (hereafter MexAG72S1) had the highest score among all the generated models [61] and a C α-RMSD of 2.2 Å with respect to the chain L in the X-ray crystal structure 2V4D. The score was evaluated using the Discrete Optimized Protein Energy (DOPE), a statistical potential frequently employed to assess homology models in protein structure prediction. The second (hereafter MexAG72S2) featured the lowest C α-RMSD (0.8 Å) from the experimental X-ray structure. All models were further validated for their overall quality using the ERRAT [62] and VERIFY3D [63] programs available in the SAVES server (; see Table S1).

MD simulations

Two independent all-atom molecular dynamics (MD) simulations were performed for the WT and G72S-substituted MexA monomers, respectively (Table S2). As stated above, the chain L from the re-refined crystal structure of MexA [52] was chosen as starting structure for the WT system (hereafter MexAWT), while the starting structures for MexAG72S were the two homology models previously described. Concerning the dimers of MexA, a first MD simulation was performed on a dimer model of the WT protein (chains L and M from the X-ray structure 2V4D, hereafter MexA2X-ray), in order to address its stability in the solvent and in the absence of partner proteins [44]. In addition, since as discussed above the unrotated state should represent the conformation of the isolated monomer, twelve MD simulations of the WT dimer were performed starting from the corresponding twelve top ranked structures obtained from ensemble-docking calculations (see Molecular Docking section). Regarding the MexAG72S dimers, no refinement of the docking poses was performed through MD simulations, as the dimeric structures strongly pointed to unlikeliness of arrangements compatible with the proper functional oligomerization within the MexAB-OprM assembly (see Results and Discussion). Hydrogen atoms were added as needed to all the systems using the tleap module of the AMBER18 package [64]. Solvated systems were prepared by inserting the structures of the proteins into truncated octahedron boxes filled with a 0.1 M NaCl water solution. The resulting systems contained ∼200 K and ∼220 K atoms for the monomeric and for the dimeric units, respectively (Table S2). MD simulations were performed using AMBER18. The AMBER force fields ff14SB [65], TIP3P [66], and the one reported in Ref. [67] were used to represent respectively the proteins, the water molecules and the monovalent ions. The structure of each system was first relaxed for 2000 steps (500 using the steepest descent algorithm, followed by 1500 steps of relaxation using the conjugate gradient method) in presence of soft restraints (force constant of 1.0 kcal·mol−1·Å−2) on all non-hydrogen atoms of the protein. The restraints were then removed, and the system was further optimized for 10,000 steps. Next, heating up to 310 K was performed in 1 ns under NVT conditions using the Langevin thermostat, and pressure equilibration of 2 ns was finally performed under the NPT ensemble (using the Berendsen barostat). Subsequently, production runs of 4 μs in length each were performed for the WT (two replicas) and each of the G72S variants of MexA. Concerning the dimers of MexA, the productive runs were of 1 μs in length each, except for MexA2X-ray, which was simulated for 2 μs. Conformational snapshots were saved every 100 ps, resulting in 10,000 frames per μs. The analysis of the simulations was performed using VMD1.9.3 [68], the cpptraj module of AMBER18, and utilities of the GROMACS 5.1.4 package [69].

Post-processing of MD simulations

Different analyses were performed with the aim to highlight structural and dynamical implications deriving from the G72S substitution. Cluster analysis of MD trajectories was performed using the average-linkage algorithm as implemented in cpptraj, using the RMSD of the backbone as metric and a cut-off of 4 Å. Secondary structure was calculated using the secstruct command of cpptraj on the L chain of the X-ray structure 2V4D [52]) and on the equilibrium trajectories. H-bonds were calculated using the hbond command of cpptraj with cut-offs of 3.5 Å and 145° for the distance between donor and acceptor atoms and the acceptor–donor-hydrogen angle, respectively. The number of (pseudo) hydrophobic contacts between each domain of MexA was calculated using the hbond command of the GROMACS2019 package [69], using the flag -contact to record the number of interatomic distances lower than 4.5 Å between atoms of two group of residues containing hydrophobic substituents (namely amino acids A, L, V, I, P, F, M, G). The fraction of native contacts between the two monomers in the dimeric structures of MexA was calculated using the nativecontacts command of cpptraj. As reference structure we chose the WT dimer formed by the chains L and M in the X-ray structure with PDB ID 2V4D [52],a contact was recorded for every two C α atoms on different monomers and distant less than the default cut-off (7 Å) from each other.

(Pseudo) binding free energies of MexAWT dimerization were estimated with the Molecular Mechanics – Generalized Born Surface Area (MM–GBSA) method [70] using the MMPBSA.py tool of AMBER18 [64]. To this aim, an additional cluster analysis was performed on the MexA2Xray trajectory, using the RMSD of the backbone as metric and a cut-off of 2 Å.; the MM-GBSA calculations were performed on 100 frames extracted from the most populated conformational cluster. According to the MM–GBSA theory, the free energy of binding is evaluated through the following formula:where , , and are the absolute free energies of complex, receptor, and ligand, respectively, averaged over the equilibrium trajectory of the complex (so-called single trajectory approach). can be decomposed as:where is the difference in the molecular mechanics energy, is the solvation free energy, and is the solute conformational entropy (not evaluated here). The first two terms were calculated with the following equations:

is the molecular mechanics energy change, contributed by the bonded (, , and ) and by the non-bonded ( and ) terms of the force field. can be modeled as the sum of an electrostatic contribution (, evaluated using the MM-GBSA approach) and of a non-polar term (, proportional to the difference in solvent-exposed surface area ). was calculated using the implicit solvent model in [71] (igb = 8 option) in combination with mbondi3 and intrinsic radii. Partial charges were taken from the AMBER force field, and relative dielectric constants of 1 for the solute and 78.4 for the solvent (0.1 M KCl water solution) were used. was approximated by the LCPO6 method implemented within the sander module of AMBER. The MM-GBSA method offers a computationally cheap platform to evaluate pairwise residue contributions to .

To compare the large-amplitude motions of WT and G72S-substituted MexA monomers during the MD simulations, a Principal Component Analysis (PCA) [72] was performed with the covar and anaeig tools of the GROMACS 2019 package [69]. The analysis of the latter system was restricted to the portion of the MD simulation prior to the structural collapse. PCA is based on the covariance matrix C with entries , where and are the vector coordinates of the ith and jth Cα atom and and represent the average coordinates calculated over the equilibrium trajectory after least-squares fitting of the lipoyl and β-barrel domains to a reference (average) structure. The eigenvectors of C represent the directions of the principal (collective) motions, whose amplitudes are determined by the corresponding eigenvalues. We restricted our analysis to the first five modes, which account for more than the 90% of the total variance (Table S3). Gibbs free energy surfaces (FES) on the different pairs between the first three PCs as reaction coordinates were estimated using the GROMACS 2019 tool sham with the “-ls” flag [73], [74]. These PCs encompass the motions featuring the maximum overlap (Table S4; overlaps were quantified by means of the anaeig tool of GROMACS 2019 – flag “-over”) between the WT1/2 and G72S1/2 trajectories. Plotting the FES along maximally overlapping PCs allows estimating the “thermodynamic cost” associated to similar motions in the two systems, mirroring changes in their intrinsic flexibility. A similar estimation was recently performed on the homologous MFP AcrA from E. coli [75].

Most analyses were performed on the equilibrium trajectories, defined by the time range greater than 0.5 μs for all systems but MexAG72S (monomeric G72S variant of MexA), where it was set to greater than 2 μs (see Fig. 2C). Exceptions were the cluster analysis and the root-mean-square fluctuation (RMSF) profiles, for which the MD trajectories sampled after 0.5 μs were used.

Fig. 2

Conformational dynamics of the WT and G72S variant of MexA. A) Conformations of the top five structural clusters extracted from the equilibrium trajectories of MexAWT1 (left) and MexAWT2 (right) compared to the X-ray structure of the protein (chain L from PDB ID 2V4D [52], coloured as in Fig. 1). The clusters are shown as grey ribbons, the top one being solid and the 2nd to 5th increasingly transparent. B) Same as in A for MexAG72S1 (left) and MexAG72S2 (right). C) Profiles of the RMSD (upper panel), radius of gyration (RoG, middle), and RMSF (lower) extracted from the simulations of WT and G72S variant of MexA. RMSD and RMSF were calculated after alignment of the lipoyl and β-domains of the protein. G72S1/21 and G72S1/22 refer to the RMSF profiles calculated respectively before and after protein collapse.

Molecular docking

Model structures of WT and mutant dimers of MexA were generated by means of protein–protein docking simulations using the software ATTRACT [76], [77]. As demonstrated by experimental data [52], the high flexibility of MexA is important in setting the possible structural arrangements assumed by dimers and oligomers of this protein. In our docking calculations we accounted for protein flexibility by considering an ensemble of monomer conformations of the protein [78]. For both MexAWT and MexAG72S the conformational ensemble was generated as follows: first, we concatenated the two productive trajectories of each system; next we performed a cluster analysis as described in the section above and imposed that the selected clusters covered up to 90% of the conformational space sampled by each protein. This resulted in five and twenty-five cluster representatives for MexAWT and MexAG72S, respectively. Each docking calculation produced 10,000 structures, leading to a total of 250,000 (5·5·10,000) and 6,250,000 (25·25·10,000) complex structures generated for MexAWT and MexAG72S, respectively. To limit the analysis to putative functional structures of the dimer we added soft upper restraints to the distance between selected C α atoms on each chain. Namely, we applied a semi-harmonic potential (k = 2.0 kcal·mol−1·Å−2) centered at 25 Å along the distance between the atoms nearest to the center of mass of the C α-hairpin domains on both chains (residue L133), as well as to the distance between atoms nearest to the center of mass of the MP domains (residue A331; Figure S2). These restraints filtered out, for instance, structures featuring the two MexA monomers not aligned along the same direction. Structures resulting from docking simulations were ranked according to the overall RMSD with respect to the dimer formed by the chains L and M of the X-ray structure [52].

Results and discussion

In this section we discuss the main results from the analysis of the homology models and of the MD simulations performed on the WT and the G72S variant of MexA. We first discuss the results for MexA monomers and subsequently the outcome of docking and MD simulations of the dimers.

MexA monomers

Structure-based homology modelling was used to introduce the point mutation G72S in MexA, As starting structures for subsequent MD simulations (vide infra), we selected the model featuring the highest score among all those generated by MODELLER (hereafter MexAG72S1). In addition, we also selected the model featuring the lowest Cα-RMSDs from the experimental structure (hereafter MexAG72S2 see Systems and Methods). The overall C α-RMSDs of MexAG72S1 and MexAG72S2 from the L chain in the X-ray structure amounted to 2.2 Å and 0.8 Å respectively (Figure S3). The same comparison performed at the level of each of the four protein domains of MexA (after structural alignment of each domain), resulted in values lower than 0.4 Å and 0.3 Å, respectively, for MexAG72S1 and MexAG72S2. Thus, the small differences between the X-ray structure of the WT protein and those generated in silico for the G72S mutant are due to slightly altered inter-domain connections. Before undergoing MD simulations, both models were validated as detailed in Systems and Methods and found to be characterized by values of the key parameters well within the typical ranges assumed for accurate structures (Table S1).

. MexAWT maintained the elongated shape assumed in the X-ray crystal structures and in the functional assembly [38], [52], [43] in two independent MD simulations (hereafter MexAWT1 and MexAWT2, respectively), each of 4 μs in length (Fig. 2). The values of the C α-RMSD with respect to the L chain of the X-ray crystal structure 2V4D [52] (after fitting the whole protein) were 6.8 ± 2.4 Å for MexAWT1 and 5.3 ± 1.6 Å for MexAWT2. Moreover, the gyration radius (RoG) remained fairly constant (around 37 Å) over both trajectories. Nonetheless, the protein displayed a large degree of flexibility, especially in the peripheral α-hairpin and MP domains (similarly to the homologous protein AcrA in E. coli [75], [79]), which together with the β-barrel domain, interact with partner proteins in the MexAB-OprM complex [43].

In detail, MexAWT1 showed nearly regular oscillations in the global RMSD, while MexAWT2 displayed an overall smoother profile (Fig. 2C). Note that, in all cases, these oscillations are due to reversible changes in the mutual arrangement of the four domains, which remain otherwise very close to the conformations assumed in the X-ray structure 2V4D (Figure S3). Indeed, upon structural alignment of the two central domains (lipoyl and β-barrel) of MexA, it becomes evident that the peaks in the overall RMSD profile are due to changes in the orientation of the α-hairpin, the MP, or both domains with respect to their neighbour ones. The MP domain displayed a slightly higher intradomain flexibility, which could be explained by its peripheral location coupled to a much larger number of unstructured amino acids compared to the α-hairpin domain. Indeed, in the L chain of the X-ray structure 2V4D this domain features 21 residues within unstructured regions (a number comparable to those of the lipoyl and β-barrel domains – 24 and 19 respectively), while only 2 are present in the α-hairpin domain. This peculiarity of the MP domain also affects its secondary structure, which in both MD simulations features the largest number of residues changing their preference towards a specific secondary structure element (mostly turn-to-bend, bend-to-turn, coil-to-β-sheet, or β-sheet-to-turn conversion; see Figure S4).

The oscillations seen in MexAWT1 are coupled to the alternation of two possible H-bonds patterns formed between residues of the β-barrel and of the lipoyl domains (Figure S5). The first one, involving residues R167 on the lipoyl domain and N270 on the β-barrel domain, is associated with a conformation of the lipoyl and β-barrel domains very close to that found in the X-ray structure. The second one, involving residues K86 and S165 on the lipoyl domain and again N270 on the β-barrel domain, induces a 35° rotation of the former with respect to the latter (high RMSD values). No significant variation was recorded in the number of hydrophobic contacts between these two domains (data not shown).

Our findings are consistent with a functional role for the elongated conformation of the protein, demonstrated both by structural data on the full MexAB-OprM complex [43], and by biochemical experiments performed on the homologous protein AcrA of E. coli [80].

In addition, our results agree with previous data from ns-long MD simulations performed on a truncated model of MexA (lacking the MP domain not available at the time [55]) and on AcrA [79], as well as with μs-long simulations on the latter protein [75].

. Unlike the WT protein, MexAG72S displayed large structural distortions after about 1.5 μs in two independent simulations, reaching a stable conformation only after ∼ 2 μs (Fig. 2C). In both cases, these large structural rearrangements were correlated to an increased compactness of MexAG72S as compared to the structure of MexAWT (Fig. 2C), mainly due to a collapse of the α-hairpin and MP domains onto the lipoyl and β-barrel domains (Fig. 2A-B and Fig. 3). Importantly, the α-hairpin domain and the β-barrel and MP domains are responsible for the functional interactions of MexA respectively with the porter domain of the cognate transporter MexB and with the coiled-coils helices of the outer membrane factor OprM [75], [42], [43], [22]. Thus, their collapse onto the central domains of the protein will drastically hamper the ability of the latter to establish flexible but functional interactions that are mandatory for the biological activity of MexAB-OprM. Note that the high RMSF values of the α-hairpin and MP domains displayed in Fig. 2C also derive from their collapse onto the central domains mentioned above, resulting in large distortions of the protein from the average (reference) structure both before and after the collapse. Indeed, the profiles calculated separately on the trajectories representing the states before and after this event showed that the functional movements of the protein initially resemble those of MexAWT, except for a slightly larger mobility of the MP domain. As expected, an overall flatter profile is associated to the dynamics of MexA after the folding of the peripheral domains onto the central core.

Fig. 3

Main rearrangements occurring in the G72S variant of MexA after ∼ 1.5 μs in two independent MD simulations. The structure evolved from an elongated shape (left) into a more compact conformation (right), in which the α-hairpin and the MP domains bend towards the β-barrel domain. In the right image the two conformations were aligned by matching their lipoyl and β-barrel domains in order to highlight the rotations of the peripheral domains (indicated by black arrows).

To further investigate if the collapse of the mutant is related to changes in the intrinsic flexibility of the protein, a PCA was performed on MexAWT1/2 and MexAG72S1/2 (restricted in the latter case to the portion of the MD trajectory before the collapse of the protein). We analysed the first five PCs, virtually covering all possible global motions in both systems. The peripheral domains display wider-amplitude motions in MexAG72S (Figure S6), as seen from the sizeably larger eigenvalues associated to this system (Table S3). The different flexibilities of the α-hairpin and MP domains in the two models are also mirrored in the overall larger RMSF values associated to MexAG72S1/2 when comparing PCs showing the largest overlap (that is, the more collinear dynamics) with MexAWT1/2 (Table S3, Figure S7). Consistently with these findings, the free energy surfaces calculated using the pairs PC1-2, PC1-3, and PC2-3 as reaction coordinates all display significantly shallower and wider basins in MexAG72S1/2 than in MexAWT1/2 (Figure S8). Overall, these results suggest that the mutation increases protein flexibility along pre-existing modes, namely anticorrelated motions between the α-hairpin and MP domains, which are those involved in the main structural changes due to the G72S substitution.

After the collapse of the protein, the interaction of the α-hairpin domain with the lipoyl/ β-barrel domains is enhanced by the onset of a few hydrophobic contacts established by residues nearby the hinge of the former domain and up to 6 (lipoyl), or 1 (β-barrel), residues of the latter (Figure S9). Moreover, about 1 μs after the collapse, a persistent H-bond network is formed between residues Q132, D136 (near the tip of the α-hairpin), E246, and R262 (in the middle of the β-barrel).

Concerning the lipoyl/ β-barrel interaction, upon movement of the MP domain, the number of interdomain hydrophobic contacts drops to about half the initial value, while no significant changes in the H-bond pattern were detected (data not shown). Two H-bonds between N51 and E53 on the β-barrel and R326 on the MP domain were formed a few hundred ns after the collapse of the α-hairpin domain (Figure S9). Interestingly, the number of hydrophobic contacts between these domains did not show any relevant change during the trajectory (data not shown).

The G72S substitution also significantly affects the secondary structure of the protein compared to WT, the largest changes being localized in the β-barrel domain (and of the same kind of those discussed for the MP domain in MexAWT; see Figure S4). This is perhaps not surprising as this domain is involved in the formation/rupture of several hydrophobic and H-bond interactions with all the other domains of MexA (Figure S9). On average, the mutation increases by about dozen units the total number of amino acids undergoing significant secondary structure changes with respect to the experimental structure.

Importantly, all the rearrangements occurring in MexAG72S do not involve significant intra-domain conformational changes, as reflected in the low RMSD values obtained for each domain after self- aligning the MD trajectory (Figure S10). The only exception is represented by the MP domain, which is rich not only in β-sheets, but also in turns and coils regions.

In fact, the value of the RMSD calculated only for the β-sheet domain is almost flat and features values comparable to those found for the other domains. Therefore, the rearrangements occurring in this domain do not alter the structure of its main core.

Clearly, our work does not rule out the possibility that the WT monomer could also sample compact conformations like those detected here for MexAG72S. Indeed, while our findings further highlight the need for extensive MD simulations to investigate the effect of point mutations on protein structure and dynamics, we are aware that our simulations could be statistically insufficient in length and number [81]. Nonetheless, the data point sharply to a drastic effect by this single mutation, which is perhaps not too surprising in view of previous work reporting changes in fold and functionality on proteins smaller than MexA [82].

MexA dimers

In the functional tripartite pump MexAB-OprM, MexA proteins assemble into hexamers formed by trimers of dimers [14], [15], [48], [38], [51], [43], [44], [22]. However, the role of the RND transporter and of the OMF proteins in stabilizing long-lived oligomeric assemblies of MFP proteins (such as MexA) is still subject to investigation [38], [42], [43], [22]. To assess the possibility to recover stable dimers of MexA in the absence of interacting partners and in aqueous solvent, we first performed a 2 μs-long MD simulation of the dimer formed by chains L and M from the X-ray crystal structure in Ref. [52] (hereafter MexA2Xray). As expected, X-ray structure, is stable along the trajectory (RMSD oscillations of about 4.5 Å and 2.5 Å in size around the experimental and average structures, respectively; see Figure S11). Chain L was significantly more stable than chain M with respect to the geometry assumed in the experimental structure. Moreover, as for the WT monomer, all domains but the MP were pretty rigid, with RMSD values below 4.5 Å after single-domain alignment; in contrast, the MP domains featured values around 10 Å (data not shown).

Given the different conformations assumed by chain L and M in hexameric assemblies of MexAWT [37], [15], [38], [43], [22], we investigated the possibility to obtain additional putative structures of MexA dimers, possibly encountered during the early stages of the oligomerization process. To this aim, we performed ensemble-docking calculations [83] as described in Systems and Methods. Namely, we selected the representative structures of the top ten conformational clusters extracted from the cumulative equilibrium MD trajectories of MexAWT1 and MexAWT2. In both systems, these clusters cover more than 2/3 of the conformations sampled by the protein during the whole simulation time (Figure S12). Several thousands of docking poses were generated (see Systems and Methods); the twelve top ranked docking models, shown in Fig. 4A, reproduced fairly well the overall structure of the experimental dimer, with values of the C α RMSD ranging from 5.4 to 7.4 Å for the whole assembly and lower than 4.2 Å for all regions but the MP domain (Table S5), as expected due to the different arrangement of this region between chains L and M in the experimental structure. Consistently with the good superposition between this structure and the docking models, the latter reproduced a significant fraction of the native contacts (Fig. 4C). We further investigated the stability of these docking models by performing, for each of them, a MD simulation of 1 μs in length (hereafter MexA2dock,i, with i = 1,…,12; see Table S2). It turns out that these models are also all relatively stable, with C α RMSD values at the end of the simulation ranging from ∼4.5 to ∼10 Å for the whole assembly. Moreover, the fraction of native contacts remained virtually constant during all MD simulations (Figure S13).

Fig. 4

A) Conformations of the twelve top ranked docking poses of MexAWT-MexAWT (left). The reference X-ray structure (2V4D) is shown as black ribbons, while the poses are shown in semi-transparent dark (chain L) and light (chain M) ribbons coloured from red to white to blue according to their docking score. The arrow on the right points to the conformations assumed by these poses at the end of MD simulations of 1 μs in length. B) Conformations of the twelve top ranked docking poses of MexAG72S-MexAG72S. C) Fraction of native contacts fNC detected in the twelve top ranked docking poses for and . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Overall, these results support the hypothesis that WT MexA proteins can assemble into (meta)stable dimeric units with lifetimes beyond the μs timescale, without the participation of other components of the MexAB-OprM pump. The latter are likely required in order to stabilize the individual components over longer timescales and/or to induce structural rearrangements matching the steric requirements needed for the assembly of the tripartite pump [15], [38], [46], [84], [43], [22].

Regarding MexAG72S, the large structural changes due to the G72S substitution impair the propensity of MexA towards dimerization. Indeed, the association time between MexA monomers is most likely longer than a few μs, implying an unfavourable interaction between collapsed geometries of MexAG72S. Nonetheless, to further investigate the propensity of MexAG72S towards dimerization, we performed protein–protein docking calculations between elongated conformations of MexAG72S, extracted from the part of the trajectories preceding the structural collapse of the protein (Figure S12). The top ranked twelve dimeric structures are shown in Fig. 4B. Comparing them to Fig. 4A, it is evident that they feature a larger displacement from the experimental structure compared to dimers of MexAWT. Indeed, the values of the C α RMSD range from 6.9 to 13.1 Å, significantly larger than those obtained for the WT dimers. Furthermore, consistently larger values are also obtained for virtually all clusters and all MexA regions (Table S5). In agreement with these data, a relatively low fraction of native contacts was found in all poses (Fig. 4C). Intriguingly, a minor change in the relative arrangement of protein domains in the elongated conformations of MexAG72S had a drastic impact on the propensity to form native-like dimeric structures. Summarizing, our findings indicate that dimeric structures of MexAG72S are most likely unstable and thus not functional.

To identify protein residues possibly important for dimerization, we estimated pairwise residue contributions to the (pseudo) binding free energy via MM-GBSA calculations performed on the equilibrium trajectory of the WT dimer. Four pairs were found to contribute by more than 10 kcal/mol to the stabilization of the dimer, all involving interactions between charged or polar residues (Table S6 and Figure S14). The arrangement of these pairs in both the elongated and collapsed conformations of the MexA monomer are shown in Figure S15. Two of them, contributing to lipoyl – lipoyl and β-barrel – lipoyl interactions, are also present in recent structures of the full MexAB-OprM pump (Table S6).

For each residue listed in Table S6 we recorded the intramolecular H-bonds occurring along the trajectories of MexAWT and MexAG72S (Table S7). The latter were split in two separated trajectories, corresponding to the conformations assumed before and after the collapse of the protein. About half of the residues featured markedly different interactions in MexAWT vs. MexAG72S after the collapse. Several new and stable interactions between the α-hairpin and β-barrel domains appeared in the latter, involving that are residues that are key for dimerization and would be thus unavailable for molecular recognition between two MexA proteins.

Conclusions and perspectives

In this work we performed homology modelling, ensemble-docking calculations and multiple μs long all-atom MD simulations on the WT and the G72S variant of the membrane fusion protein MexA from P. aeruginosa. The main goals of the study were to assess the stability of WT MexA dimers in water solution and in the absence of protein partners, and moreover to shed light on the molecular mechanism behind the impaired functionality of the MexAB-OprM pump experimentally ascribed to the G72S substitution.

First, we showed that the WT protein keeps its overall elongated structure during μs-long MD simulations and retains however the functional flexibility of the peripheral α-hairpin and MP domains.

Second, we showed that the G72S substitution impairs the stability of the linker between the lipoyl and the α-hairpin domains, resulting in enhanced oscillations around the hinge. These oscillations induce large inter-domain rearrangements in the structure of the protein, leading to an overall collapse of its peripheral domains on the β-barrel domain in a μs timescale. This collapse is accompanied by several changes in the secondary structure of the domains involved in both MexA-MexB and MexA-OprM binding, and the resulting compact structure is stabilized by the formation/rupture of several interdomain interactions. Our findings suggest that the mutation alters the balance between structural flexibility and stability, which makes possible to assemble a stable complex in the WT protein. In MexAG72S there is a shift toward increased flexibility, leading in turn to the sampling of conformational states that are incompatible with dimerization.

Third, we showed that dimers of WT MexA are stable in aqueous solvent on the μs timescale, even in the absence of partner proteins. Importantly, this stability coexists with a retained plasticity of the peripheral domains of the protein, responsible for the recognition of cognate sites on MexB and OprM. Clearly, the interaction with other dimers as well as with partner proteins can: i) affect stability on timescales longer than those simulated in this work; ii) induce further conformational changes (particularly in the MP domains) favouring the interaction with MexB and/or OprM. Nonetheless, the study of the dimerization propensity of MexA (and homologous proteins) in the absence of other pump components is relevant to understand the dynamics underlying the full assembly of MexAB-OprM, and homologous tripartite systems. Indeed, it is quite likely that in these conditions MexA will have a low propensity towards oligomerization to more than two units, as the formation e.g., of a hexamer would imply a quite substantial thermodynamic gain, that is a structure perhaps too stable to engage the pump.

Fourth, we showed by extensive ensemble-docking calculations followed by MD simulations that the mutation has a drastic effect on the capability of MexA to assemble into functional dimers.

Overall, our results provide a possible explanation for the experimental findings by Poole and co-workers [49], who discovered that the G72S substitution significantly affects the functionality of the pump by interfering with MexA oligomerization. In addition, the present work provides precious insights for structure-based rational design of drugs targeting the assembly of the major efflux pump MexAB-OprM of P. aeruginosa [85]. In particular, the different dynamical behaviour of the four protein domains in the WT, mutated, monomeric and dimeric forms of the protein could be exploited to support studies aimed at inhibiting the assembly of efflux systems [57], [86], [8].

Disclaimer

The Authors received no support, financial or otherwise, for the research, authorship, and/or publication of this article from Zimmer Biomet nor was the research conducted on behalf or in the interests of Zimmer Biomet.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.