Molecular dynamics simulationsThe closed state x-ray structure (PDBID: 1EK9)12, and the open state cryo-electron microscopy (cryo-EM) structure (PDBID: 5NG5)14 were used as the initial TolC coordinates for MD simulations. The D371G mutant which we hypothesized destroyed key salt bridges on the periplasmic entrance side of TolC was also investigated via MD simulations; the initial coordinates of this mutant were created using the Pymol34. These coordinates were subjected to MD simulations to understand the time dependent behavior of the protein under physiological conditions. Since TolC is a Gram-negative bacterium protein spanning the outer membrane and the periplasmic space, the orientation of TolC was arranged according to the Orientations of Proteins in Membranes (OPM) database35,36. The overall system was prepared on the CHARMM-GUI platform37,38,39,40. The protein was placed in a bilayer that contains 256 POPE lipids and solvated with TIP3P water molecules. Na+ and Cl- ions were added to fulfill the isotonic condition of physiological environment (0.15 M).MD simulations were performed using NAMD41,42 and CHARMM36 force field parameters with periodic boundary conditions38. The total number of atoms for all our simulations are around 150,000 with box dimensions of 94 × 94 × 173 Å. The particle-mesh Ewald method with a grid spacing of maximal 1 Å per grid in each dimension was used43. 12 Å cutoff radius with a switching function turned on at 10 Å was applied to calculate van der Waals (vdW) interactions. The RATTLE algorithm44 was utilized to set the step size of 2 fs during numerical integration with the Verlet algorithm. Temperature regulation was managed with Langevin using a damping coefficient of 5 ps−1. Pressure regulation was achieved through a Nose-Hoover Langevin piston, and volume variations were set to remain uniform in all directions. Systems were minimized for 50,000 steps to remove undesired van der Waals contacts. The energy-minimized systems were subjected to 5 ns MD run with constant temperature control (NVT). The production parts of the MD simulations of length 500 ns were run in the isothermal-isobaric (NPT) ensemble at 310 K and 1 atm. The data from MD simulations were analyzed by the open-source Python package ProDy45,46,47,48 as well as custom Python and tcl scripts.We carried out nine ‘free MD simulations’ by placing carbenicillin molecule at different locations along the channel. In three of these simulations, carbenicillin molecule located at the periplasmic entry site which corresponds 0–25 Å region of the protein (see Fig. 1B), three are in the central region where the transfer of the drug molecules to the exit region occurs (25-100 Å), and the other three are in the exit region (100–140 Å). We label them entry/transfer/exit and we report results from the averages over 100 ns simulations in each region in Supplementary Fig. S6. The simulation parameters mentioned earlier for closed and open structures were also used for these simulations. The RMSD of Cα atoms for these simulations, and the number of hydrogen bonds in each region during the simulations are shown in the Supplementary Fig. S6. We found that TolC-carbenicillin hydrogen bond densities are the largest in the entry region, very few occur in the transfer region and approximately the one hydrogen bond is maintained by carbenicillin throughout the simulation in the exit region. These occurrences are similar to what we had observed in SMD simulations. Residues display larger than 5% hydrogen bond occupancy at least one of the simulations are displayed in Supplementary Table S2. The only one with significant overlap which was not observed in our SMD simulations (Fig. 5D) is R143 which is at the periplasmic tip of TolC and is below our docking position for SMD simulations. Therefore, the drugs pulled in SMD simulations are never in contact with this residue.Principal Component Analysis (PCA)Our MD simulations cover the sub-μs time scale dynamics of the protein. Normal mode analysis (NMA) models protein as a harmonic oscillating system and determines the protein dynamics around an energetically stable state of the protein. Low-frequency vibrations, also known as low-energy modes, are associated with collective motions, whereas higher-frequency modes are linked to local deformations49,50,51,52,53,54,55. Principal component analysis (PCA) is an approach for comprehending protein dynamics that involves analyzing an ensemble of conformations obtained from MD simulations21. The mathematical principles underlying PCA and NMA share similarities: both require solving linear equations. In PCA, the variance-covariance matrix (C) is diagonalized; C approximates the inverse of the Hessian matrix employed in NMA49. The covariance matrix is constructed using the analysis of multiple snapshots from a MD trajectory, rather than through the calculation of the second derivative near a selected energy minimum. In both instances, the outcomes provide insights into the potential motions a particular structure can undergo, while in a stable conformation. Reducing the dimensions obtained from MD simulations can be advantageous for an unbiased global analysis and to isolate specific configurations. We used principal components to discern the dominant motions of TolC during the MD simulations. PCA was carried out using ProDy45,46,47,48 and the Cα atoms were subjected to analysis. In our previous work we have shown that 40 ns chunks of equilibrated trajectories can sufficiently sample the elastic motions around a selected energy minimum of a protein56,57. We have therefore utilized 240–280 and 460–500 ns fragments of the MD trajectories. Extracellular loops and a part of the equatorial domain which displays large motions obscuring the motions of the overall TolC were excluded in the construction of the C matrix.Perturbation Response Scanning (PRS)The PRS method is a powerful tool we previously developed to identify key residues that are functional in the interconversion of one conformer of a protein into another. PRS uses linear response theory to shed light on the possible conformational changes of protein might undergo by applying external forces on residues in silico23. In previous studies58,59, we and others used PRS to study water-soluble proteins. However, TolC is a transmembrane protein with much larger size and requires embedding the TolC homotrimer in a lipid bilayer. Studying TolC dynamics via PRS also provided us the opportunity to evaluate the utility of this powerful method on membrane proteins and solve possible problems arising from this relatively complex environment.PRS requires experimentally obtained initial and target coordinates of the Cα atoms of a protein, Sinitial and Starget, respectively. In our case, these structures were the closed and open forms of TolC. R0 represents the unperturbed state coordinates (initial), R1 represents the perturbed state coordinates (target) of the Cα atoms. The external forces vector (ΔF) in random directions are sequentially applied on each residue of the initial structure. Each predicts displacements in the initial structure, ΔR1 by$${\Delta {\bf{R}}}_{1}=\langle {{\bf{R}}}_{1}\rangle -\langle {{\bf{R}}}_{0}\rangle \approx \frac{1}{{k}_{B}T}{\bf{C}}\Delta {\bf{F}}$$
(1)
where kB is the Boltzmann constant, and T is temperature in Kelvin units.A total of 250 random forces are applied to each Cα atom in the protein and results in 250 N predicted displacement vectors, where N is the number of residues. The agreement between experimentally observed and predicted displacements of each residue i in response to the applied force is calculated by Pearson correlation coefficient (Ci) between predicted and experimental displacements. Thus, each ΔR1 for applied random forces is compared with the experimentally determined difference between Sinitial and Starget coordinates (ΔS), averaged over all affected residues k:$${C}_{i}=\mathop{\sum }\limits_{k=1}^{N}\frac{\left[{\left(\Delta {R}_{k}\right)}^{i}-{\left(\overline{\Delta R}\right)}^{i}\right]\left(\Delta {S}_{k}-\overline{\Delta S}\right)}{\left(N-1\right){\sigma }_{R}{\sigma }_{S}}$$
(2)
∆Sk is the experimental displacement between initial MD frame and the target conformation where \(\overline{\Delta S}\) the overbar represents average displacement, and σR and σS denote the standard deviations of experimental and predicted structures, respectively. In this equation, the magnitudes of the displacements, ∆Rk, are compared with ∆Sk. Ci takes on values in the range [0,1] where a correlation greater than 0.7 refers to strong correlation with the experimentally determined conformational changes. PRS was calculated from 240-280 and 460-500 ns chunks of these 500 ns simulations. Extracellular loops and a part of the equatorial domain were excluded in PRS calculations. Since TolC is a homotrimeric outer membrane protein, we performed PRS analysis using each chain as a replica. In the end, obtained results were averaged over 6 sets of data.Steered Molecular Dynamics (SMD)We performed SMD calculations using the closed-form TolC structure (PDB ID: 1EK9) as the initial point12 in separate SMD runs, carbenicillin, piperacillin, or oxacillin was positioned at the periplasmic tip of the protein by docking, CB-Dock2 tool60. The charges of carbenicillin, piperacillin, and oxacillin in water are −2, −1, and −1, respectively. While these charges might shift depending on the position of the drugs in the protein environment, we assumed the pKa of all positions and the charges remained constant in our simulations. The simulation boxes for the membrane, protein, and three different antibiotic molecules were optimized using the CHARMM-GUI web server as described in the “MD simulations” section. Force fields for the antibiotic molecules used in our study were generated by CHARMM general force field module of CHARMM-GUI39,61,62. The systems obtained with the three different antibiotics were minimized for 50,000 steps and equilibrated for 0.5 ns. A 50 kcal mol−1 Å−2 harmonic spring constant was applied to a carbon atom on the hexane structure of antibiotic molecules (depicted as purple spheres in Fig. 1C) for forcing the drug molecule pass through the TolC channel, a value that satisfies the stiff spring approximation63. Antibiotics were pulled along the z-axis until they reached the extracellular loops of the protein at a constant velocity of 5 Å ns−1. The molecules traveled approximately 140 Å; each SMD simulation thus lasts nearly 30 ns. No restraints were applied to the substrates in the x and y directions. The pulling simulations of each antibiotic molecule were repeated 45 times to increase the signal-to-noise ratio. The average work done to reach the extracellular side was then calculated for each antibiotic molecule. We note that the general trends and the energy levels in the work curves did not change when the pulling atom was modified to the nitrogen atom in the middle of the carbenicillin molecule from the carbon atom on the hexane ring. Also, modification of the pulling velocity to 0.5 Å ns−1 from 5 Å ns−1 did not change the profiles in the first 30 Å of the pump. We found that residues that are in contact with the carbenicillin above the threshold we have set in Fig. 5D overlap with our 0.5 Å ns−1 pulling simulations. We thus concluded that doing 45 pulling simulations at 5 Å ns−1 provides the information necessary to compare our simulation findings with our mutational scanning experiments.Bacterial Strains and Growth ConditionsEscherichia coli (E. coli) cells were grown at 37 °C in a M9 minimal medium supplemented with 0.4% glucose and 0.2% amicase. The BW25113 E. coli strain (CGSC No.: 7636) and the ΔtolC732::kan E. coli strain (CGSC No.: 11430) were obtained from Coli Genetic Stock Center31. Using the method from reference64, we removed the kanamycin resistance marker in the ΔtolC732::kan E. coli strain. In the article, we refer to the E. coli strain tolC gene deletion as ΔtolC strain.The tolC gene from BW25113 (wild type, WT) strain was PCR amplified using 5′-ATTCAAAGGAGGTACCCACCATGAAGAAATTGCTCCCCATTC-3′ (forward), and 5′-AGAAATCGATTGTATCAGTCTCAGTTACGGAAAGGGTTATGAC-3′ (reverse) primers. PCR amplified tolc gene was then cloned into the pSF-Oxb14 plasmid which was obtained from Oxford Genetics (OGS557, Sigma) using the NEBuilder HiFi DNA Assembly kit (New England Biolabs); by strictly following the protocols recommended by the manufacturers. The resulting plasmid (pSF-Oxb14-tolC) contains kanamycin resistance cassette and an Oxb14 constitutively open promoter region. The introduced E. coli tolc gene sequence was confirmed by Sanger sequencing.Whole gene saturation mutagenesis library (SML) for each codon in the tolC gene used in our experiments was carried out by two separate PCR reactions for each codon in the tolC gene. First, a portion of tolC gene in the pSF-Oxb14-tolC plasmid was amplified and the target codon was randomized using a primer containing NNS nucleotide sequence (where N represents A, C, G, or T, and S represents G or C nucleotides). The first PCR products were referred as insert. Second PCR reaction amplifies the remainder of pSF-Oxb14-tolC plasmid which is referred as backbone. The custom software to design mutagenesis primers is available (See https://github.com/ytalhatamer/DMS_PrimerDesignTool). The PCR products (insert and backbone) were assembled using NEBuilder HiFi DNA Assembly Kit (E5520, New England Biolabs), and the resulting plasmids were introduced into NEB-5-alpha cells (C2987, New England Biolabs). Plasmids were extracted from the transformed cells using the Nucleospin Plasmid kit. This assay resulted in libraries for each residue. Equimolar amounts of the library were pooled into 5 sublibraries and these pooled sublibraries were transformed into ΔtolC strain for selection experiments31. We refer to this library as TolC-SML. We refer the ΔtolC strain supplemented with the pSF-Oxb14-tolC plasmid as ΔtolC+ptolC. All selection experiments were conducted with minimal M9 media containing 50 µg ml−1 of kanamycin to prevent loss of the pSF-Oxb14-tolC plasmid.Determination of minimum inhibitory concentrations (MIC)Bacterial cells were overnight grown and optical densities (OD600) of cultures were measured in a 1 cm pathlength cuvette. OD600 was then adjusted to 0.001 for which corresponds ~5 × 105 colony forming unit (CFU) per milliliters. Bacterial cells were then grown in 96 well plates ( ~ 200 µl per well) in which drug gradients were created using serial dilutions (Fig. 2A). Highest starting antibiotic concentrations on these plates were 100 µg/ml for carbenicillin, 20 µg/ml for piperacillin, and 1875 µg/ml for oxacillin. The dilution factor for carbenicillin and piperacillin was 1/2 while it was 1/5 for oxacillin. These 96-wells were incubated for 18 hours in a humidity-controlled shaker operated at 37 °C and 400 RPM. After the overnight incubation, cell densities were measured using a plate reader (VictorX, PerkinElmer). The experiments were performed in triplicates for each antibiotic molecule. For each strain and antibiotic compounds, the MIC value was defined as the lowest antibiotic concentration at which the final OD600 was below ~0.04 after background correction.Selection AssayWe used three antibiotics: carbenicillin, piperacillin, and oxacillin as selection agents. We used the saturation mutagenesis library (SML) for tolC31. In this library, all residues except start codon (471 residues in the mature TolC protein, and the 22 residue-long signal sequence) were mutated to all other possible residues, and a pool of mutant library was generated (Fig. 8). We measured the fitness effects of mutations under selection using a liquid-based sequencing assay. Our assay does not allow us to evaluate the expression levels or stability of thousands of TolC mutants. However, by completing the screen with three different antibiotics, we can confidently deduce that a TolC mutant is expressed and stable enough if the relative fitness of the mutant under the selection with at least one antibiotic, with respect to the wild-type TolC, is non-deleterious, suggesting that the efflux for at least one antibiotic is still occurring.Fig. 8: Selection assay protocol.We conducted deep mutational scanning on TolC, systematically mutating all 493 residues, which includes the 21-amino acid long signal sequence preceding the start codon (excluding it). Residues were mutated to all the 19 other amino acids, along with the stop codon. Obtained mutants were pooled together and grown overnight. Overnight grown mutant pool and ΔtolC strains were split into cultures: untreated (UT), CBC, PIP, OXA; the experiment was replicated once. Cultures other than the untreated were treated with the relevant selection factors. Cultures were exposed to drugs for 3 h for selection. This was followed by 6 hours long incubation in fresh growth media without any selection agent for the recovery of cells (gray abbreviations on the tubes represent previously applied antibiotics). We isolated tolC mutant carrying plasmids and PCR amplified the tolC gene for amplicon sequencing. Mutation frequencies were calculated by deep sequencing and fitness values of each mutation were calculated using the formula in the figure. (Created with BioRender.com).We grew the mutant library in M9 minimal medium supplemented with 0.4% glucose and 0.2% amicase overnight, diluted to a final optical density of 0.001, and then exposed these cultures to one of the selection factors: carbenicillin, piperacillin, and oxacillin for 3 h at 37 °C, 400 rotation per minute (RPM). In parallel, we grew cultures without any selection (untreated: UT) for separating drug-induced fitness effects from other potential fitness effects that are not related with the bactericidal effects of antibiotics. The concentrations for these antibiotics in this selection assay were 20 µg/ml and 10 µg/ml for carbenicillin, 10 µg/ml for piperacillin, and 460 µg/ml for oxacillin. These concentrations were close to the MIC of the related compound (Fig. 2B–D). We note that for carbenicillin we used two carbenicillin concentrations because the dynamic range of the fitness values was not large enough when the selection was done with the lower concentration. The dose-response curve for carbenicillin is very steep (higher Hill coefficient), making it difficult to adjust the selection strength.At the end of the selection period, these cultures were spun down at 5000 × g for 5 minutes. The pellets obtained from spinning were resuspended in fresh M9 minimal medium supplemented with 0.4% glucose, 0.2% amicase and 50 µg/ml kanamycin and incubated at 37°C, 400 RPM for 6 hours to increase cell density and ensure harvesting of enough tolC carrying plasmid for sequencing. We recorded cell densities during this recovery period as depicted in Supplementary Fig. 1. Cells were harvested by spinning cultures at 5000 × g for 5 minutes. Pellets were collected for plasmid purification. tolC gene variants on these plasmids were amplified with polymerase chain reaction (PCR) using 5’-TACCCGGCAGATCTTTGTCGATCCTA-3’ (forward), and 5’-GTGAGCTGAAGGTACGCTGTATCTCA-3’ (reverse) primers.These products were sequenced on Illumina NovaSeq 6000, producing 151 bp reads. Sequence reads were compared with the wild type tolC sequence and mutations were listed. Sequence reads that had mutations in more than one residue were excluded from the analysis. Synonymous mutations yielding the same amino acid replacement were grouped together and used as a reference for relative fitness calculations. The frequency of each mutation was calculated by dividing the number of counts for that mutation with the number of all reads, including alleles with multiple mutations31. We then used the formula displayed in Fig. 8 to calculate relative fitness effects of each mutation.Statistics and ReproducibilityMD simulations were performed for 500 ns. The average RMSF graph was plotted for the homotrimeric TolC protein. The data points for chain RMSFs were plotted as pale dots. The principal components (PC) and perturbation response scanning (PRS) analyses were performed on two chunks from the related trajectory. The p-value for PRS data were calculated by one-sample two-tailed t-test, results were reproducible (p < 0.01). SMD simulations have 45 replicates for each antibiotic molecule, antibiotic-work curves for each run are available. Work data were collected for each 200 fs. Error bars for work graphs were determined by using block averaging for the data in each bin, and by applying the bootstrapping method (See Fig. 5A and Fig. 7A). Significance threshold for fitness measurements is 2.5σ (99.4% confidence interval).Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
Structural shifts in TolC facilitate Efflux-Mediated β-lactam resistance
