Introduction

Methicillin-resistant Staphylococcus aureus (MRSA) is a significant global healthcare challenge due to its resistance to methicillin and other beta-lactam antibiotics1. MRSA infections range from mild skin conditions to severe, life-threatening diseases such as septicemia, pneumonia, and endocarditis, contributing to prolonged hospital stays, increased healthcare costs, and limited treatment options2. The high morbidity and mortality rates associated with MRSA infections underscore the urgent need for new therapeutic strategies and antimicrobial agents. Given the continual evolution of MRSA strains, identifying new drug targets and mechanisms to inhibit the pathogen is essential for reducing its public health impact3,4.

Natural products (NPs) have re-emerged as pivotal agents in the search for new antibacterial strategies, despite their inherent limitations in drug development. Quercetin, a well-known flavonoid with broad-spectrum antibacterial properties5, exemplifies this potential. Although its poor oral bioavailability presents challenges for direct clinical use, understanding its mechanisms of action at the molecular level could lead to the development of optimized derivatives with enhanced therapeutic efficacy. Despite its recognized efficacy, the precise mechanisms by which quercetin exerts its antibacterial effects remain under-explored, particularly against MRSA. Recent studies propose multiple possible mechanisms of action for quercetin, highlighting its complex interaction with bacterial cells. Some reports suggest that quercetin inhibits β-lactamase, an enzyme used by resistant bacteria to inactivate β-lactam antibiotics6. This mode of action suggests a potential for quercetin to restore the efficacy of these antibiotics against resistant strains. Other research indicates that quercetin interacts with penicillin-binding proteins (PBPs), which are critical for bacterial cell wall synthesis and maintenance7,8. Moreover, the role of quercetin in modulating bacterial virulence factors has also been highlighted. A particularly intriguing study demonstrates that quercetin targets the ClpP protease, a key virulence factor in MRSA, significantly reducing the virulence of MRSA. This finding is critical as it suggests that quercetin not only inhibits growth but also diminishes the pathogenicity of MRSA9. However, despite these promising findings, the absence of crystal structures or structural data showing quercetin bound to MRSA-specific target limits the ability to visualize interactions at the molecular level. This limitation is crucial, as structural insights are essential for the rational design of quercetin-based therapies or their derivatives.

Given the absence of experimentally determined structures of quercetin-bound MRSA proteins, we leveraged a “targeted subtractive proteomics” approach to identify bacterial proteins that could interact with quercetin in a manner analogous to their plant counterparts10. The term “subtractive proteomics” refers to the strategic screening of an organism’s proteome (MRSA in the current study) to identify and prioritize essential proteins as viable drug targets, narrowing the potential targets to those most impactful on bacterial survival when inhibited11. Within this framework, molecular docking predicts the most probable quercetin binding sites in the identified protein, providing initial insights into binding affinity and geometric compatibility. Subsequent molecular dynamics (MD) simulations offer a dynamic representation of this interaction over time, simulating the effect of quercetin cellular-like conditions. These simulations assess the stability of the complex and monitor the conformational changes of the enzyme, crucial for understanding the functional implications of this interaction10. This approach not only advances our understanding of quercetin’s antibacterial mechanism but also refines hypotheses before experimental validation, ensuring that computational predictions focus on strategically relevant targets.

In the current study, we explored potential bacterial targets of quercetin within MRSA using an integrated computational approach, aiming to elucidate its possible mode of action. We employed targeted subtractive proteomics to identify MRSA proteins structurally similar to known quercetin-binding enzymes in other organisms. The identified proteins were further refined against the Database of Essential Genes (DEG)12 and DrugBank13 to identify essential, druggable targets. MurG, a crucial enzyme in the peptidoglycan biosynthesis pathway, was found to be a significant target for quercetin binding and interactions between quercetin and MurG were investigated using structural modelling and MD simulations to understand the stability and dynamics of the quercetin-MurG interaction. These results suggest that quercetin can inhibit S. aureus by potentially binding to MurG, leading to defects in peptidoglycan biosynthesis. While quercetin itself may not be an ideal therapeutic, this study provides critical insights that will guide the development of more effective, modified derivatives that target MRSA.

Materials and methods

MRSA proteome retrieval

In this study, the complete proteome of S. aureus strain MRSA252 (Proteome ID: UP000000596) was retrieved from UniProt14. This strain was selected due to its clinical relevance as a methicillin-resistant S. aureus strain associated with nosocomial outbreaks, as well as the availability of a complete and annotated proteome dataset for computational analyses.

Filtering proteome to paralog sequences

To ensure a manageable and non-redundant dataset for downstream analyses, the MRSA252 proteome was refined using the CD-HIT tool15 (with a threshold of 90%). This step clusters similar protein sequences and retains only unique representatives, effectively reducing redundancy in the dataset. By eliminating paralogous sequences, this process ensures a streamlined and accurate dataset, facilitating the identification of essential and druggable proteins during subsequent analyses.

Filtering essential proteins

Essential proteins, critical for bacterial survival, were identified by querying the non-redundant proteome against the DEG database12. This process involved stringent criteria, including an e-value threshold of < 0.0001 and bit score > 100, to ensure the selection of high-confidence essential proteins. By focusing on these essential proteins, we aimed to prioritize targets whose inhibition could significantly disrupt bacterial viability.

Druggability analysis

To further refine the list of essential proteins, druggability analysis was conducted by querying the filtered essential proteins against the DrugBank database. This process utilized stringent criteria, including bit scores > 100 and e-values < 0.005, to identify high-confidence druggable proteins based on their alignment with known drug targets. These proteins were prioritized for their potential to interact with quercetin, a compound with an established role as a glycosyltransferase ligand.. This database encompasses a comprehensive collection of drug targets and includes 4323 non-redundant protein sequences related to drug targets, enzymes, transporters, and carriers associated with 6712 drug entries. These entries comprise 131 FDA-approved biotech (protein/peptide) drugs, 1,448 FDA-approved small-molecule drugs, 5080 experimental drugs, and 85 nutraceuticals13.

Identification of quercetin co-crystallized proteins in the PDB

The initial phase of the study involved identification of proteins co-crystallized with quercetin in the Protein Data Bank (PDB)16. This process began by searching for the PDB for entries where quercetin is recorded as a ligand. Search queries that incorporated keywords related to “quercetin” and its chemical identifiers or ligand IDs refined the search to isolate only those protein entries that feature quercetin in their crystallized structures. Relevant data such as protein names, PDB IDs, structural resolutions, and details concerning the binding sites were extracted using the PDB advanced search function (https://www.rcsb.org/). Structural analysis and visualization were performed using UCSF ChimeraX17 to assess factors such as quercetin’s orientation within the binding site and key amino acids involved in the interaction. Subsequently, proteins confirmed to be co-crystallized with quercetin were categorized and organized according to protein family, biological function, and the nature of their interaction with quercetin.

Selection of analogs of quercetin co-crystallized proteins from a filtered MRSA proteome

Sequence comparison techniques were used to identify structural analogs of proteins known to interact with quercetin and druggable proteins from the MRSA proteome. This process began by transforming the identified set of druggable proteins into a custom database using the "makeblastdb" command from BLAST+18. The database was then queried using "blastp" against a curated list of protein sequences known to bind quercetin, which included structural PDB entries such as 6IJD, 1H1I, 8SFW, 3BPT, 2C9Z, and 8WRJ. The goal was to identify proteins within the MRSA252 proteome that not only exhibit significant sequence similarity to known quercetin-binding proteins but also possess potential binding sites for quercetin. This step aimed to narrow down potential targets to those with a higher likelihood of interaction with quercetin, based on their structural and functional similarities to previously characterized proteins.

Molecular docking

Prior to docking, MurG and potential ligands, quercetin and UDP, underwent preparatory steps using Schrödinger’s Maestro software suite V. 202019 to optimize their structures and ensure precise interaction modeling. The preparation of the protein structure was executed using the Protein Preparation Wizard module in Schrödinger Maestro19, which included adding hydrogen atoms, removing water molecules, and optimizing hydrogen bond assignments. The structure was then minimized using the OPLS3e force field to refine the geometry20. The quercetin SDF file, retrieved from PubChem, was prepared with the LigPrep module in Schrödinger Maestro.

Active sites on MurG were identified by superimposing its structure, retrieved from AlphaFold-221, onto the structure of 8SFW—a glycosyltransferase known to interact with quercetin. Subsequent double docking experiments involved separately docking quercetin and UDP into their predicted active sites on MurG using the Glide module of Schrödinger Maestro in the extra precision (XP) docking mode. The Docked complex was exported in ‘.pdb’ format and 2D representations of protein–ligand interactions were generated using Discovery Studio Visualizer V. 2022.

MD simulations

The stability and dynamic behavior of the MurG protein in complex with quercetin and UDP, as well as its apo form, were explored through MD simulations spanning 250 ns. These simulations were conducted using the Desmond module of Schrödinger Maestro V. 202019, employing the OPLS3e force field20. MD simulations are crucial for elucidating the conformational flexibility of protein–ligand complexes, identifying critical interactions that contribute to binding, and potentially uncovering allosteric effects or induced-fit phenomena22. For the simulations, the protein–ligand complex and apo form of protein were solvated in an orthorhombic box filled with TIP3P water molecules, and counter ions were introduced to achieve electrical neutrality of the system. The preparation of the simulation system involved energy minimization followed by a heating phase and equilibration under NPT conditions, which maintained a constant pressure (1 atm) and temperature (300 K) using the Nose–Hoover thermostat23. The production phase of the MD simulations was executed over 250 ns, yielding trajectories consisting of 1000 frames. These comprehensive simulations provide valuable insights into the molecular interactions and structural dynamics critical for understanding the mechanistic basis of MurG inhibition by quercetin.

Post MD simulation assessment

The trajectories from the MD simulations were meticulously analyzed using the Simulation Interaction Analysis module in Schrödinger Maestro V. 202019. Initially, the trajectories were aligned to the original protein structure to remove any translational and rotational motions of the system. This step is crucial for accurate assessment of structural changes. The structural stability of each protein–ligand complex was evaluated by calculating the root-mean-square deviation (RMSD) of both the protein backbone and the ligand heavy atoms across simulation time. Additionally, the structural flexibility was assessed by computing the root-mean-square fluctuation (RMSF) of the protein residues, with RMSF values highlighting regions of variable mobility24. Analysis of hydrogen bonding, hydrophobic contacts, and pi-stacking interactions were conducted to identify key residues involved in ligand binding. The frequency and duration of these interactions were meticulously monitored throughout the simulation using the Protein–Ligand Interaction module. These interaction profiles were subsequently compared with those obtained from initial docking studies to evaluate the consistency and reliability of the predicted binding modes25. Post-simulation analyses further focused on evaluating several critical parameters to elucidate the structural and dynamic nuances of the MurG-quercetin complex. These parameters included the radius of gyration (rGyr) for evaluating protein compactness and the molecular surface area (MolSA) alongside the polar surface area (PSA) to analyze the interaction surface. Additionally, the solvent-accessible surface area (SASA) was examined to gain insights into hydration shell dynamics and ligand accessibility.

Following the MD simulation, Molecular Mechanics Generalized Born Surface Area (MM-GBSA) analysis was performed on the trajectory to evaluate interaction of quercetin with MurG, utilizing Python script ‘thermal_mmgbsa.py’ from Schrödinger. This script is designed to calculate the binding free energies of ligand–protein complexes, providing insights into the interaction dynamics. The general structure of the command used for running the script included various parameters to specify the input file, frame range, and computational resources, formatted as follows: “run thermal_mmgbsa.py [input_file] -j energy-calc -start_frame [start_frame] -end_frame [end_frame] -step_size [step] -HOST localhost:[cores] -NJOBS [jobs] > [output_file]”. The resultant data encompassed various energy components, including Binding energy (ΔGBind), Coulomb energy (ΔGCoulomb), Covalent bond energy (ΔGCovalent), Hydrogen-bonding correction (ΔGHbond), Lipophilic energy (ΔGLipo), Pi-Pi packing correction (ΔGPacking), Solvation energy (ΔGSolvation), and Van der Waals energy (ΔGvdW). To effectively visualize and interpret this data, violin plots were created using Plotly Studio (https://chart-studio.plotly.com/).

Lastly, Principal Component Analysis (PCA) was employed to discern collective movements within the protein–ligand complex, focusing specifically on the Cα backbone atoms. This analysis was facilitated by the Bio3D package, which computed the covariance matrix of positional changes to capture the essence of protein dynamics during the simulations26. The eigenvectors and eigenvalues derived from PCA elucidated the direction and magnitude of atom movements, offering a granular view of structural variations over the course of the simulation, thus enhancing our understanding of the dynamic behavior of the protein–ligand complex in comparison to apo form of the protein under study.

Analysis of MurG conservation across S. aureus strains

Protein sequences of MurG from multiple S. aureus strains were retrieved from the UniProt database and aligned with Clustal Omega to identify conserved regions and any potential variability among the strains. In addition to the sequence alignment, 3D structures of MurG proteins were obtained from the AlphaFold-2 database for the various strains. These structures were superimposed using ChimeraX17 to assess the structural conservation of the binding site residues identified in the docking studies. This approach allowed for a detailed visualization of the degree of conservation at the molecular level across different strains.

Results

Identification of MurG as a potential quercetin target in S. aureus MRSA252 by subtractive proteomics

S. aureus MRSA252 is a representative of the epidemic methicillin-resistant S. aureus EMRSA-16 clone, which is known for its role in nosocomial infections globally. This strain’s unique genetic profile, with about 6% of its genome containing novel elements not present in other sequenced strains14, makes it an ideal model for studying the evolution of virulence and antibiotic resistance. As such, subtractive proteomics was employed to identify potential quercetin binding targets within the S. aureus MRSA252 proteome. Initially, non-redundant sequences from among the 2640 proteins of the S. aureus MRSA252 proteome were selected using CD-HIT15, which reduced the protein count to 2,625. Subsequent screening involved the use of the DEG database12, which contains a list of essential genes for bacterial survival and are thus potential targets for antibacterial strategies. This strategy reduced the number of proteins to 1637 (S1 Table), which were subsequently screened against the DrugBank database13, which includes information on proteins that are known drug targets or have been implicated in drug interactions. This stepwise approach produced a list of 816 essential and druggable proteins (S2 Table) from the S. aureus MRSA252 proteome that were used for further analyses (Fig. 1).

Fig. 1
figure 1

Outline of targeted subtractive proteomics methodology employed for filtering proteins from the proteome of S. aureus MRSA252.

Full size image

We hypothesized that quercetin-binding proteins in S. aureus could share significant structural similarities with known quercetin-binding proteins from other organisms. To investigate this, the PDB was searched for proteins co-crystallized with quercetin, resulting in six identified proteins (Table 1). Notably, four of these proteins—PDB entries 6IJD, 8SFW, 2C9Z, and 8WRJ—belong to the glycosyltransferase family and all showed significant homology in their active sites, where quercetin and UDP bind (Fig. 2). This alignment highlights their structural and potential functional similarities, underlining their importance in the glycosylation of natural products. Of the remaining two, one was identified as quercetin 2,3-dioxygenase (PDB entry 1H1I), while the other was human 3-hydroxyisobutyryl-CoA hydrolase (PDB entry 3BPT).

Table 1 Protein structures co-crystallized with quercetin from the PDB.
Full size table
Fig. 2
figure 2

Superimposition of glycosyltransferases from plants co-crystallized with quercetin retrieved from protein data bank (PDB).

Full size image

Using the identified quercetin-binding glycosyltransferases as a reference, we performed BLASTP analysis to identify homologous proteins within the 816 essential and druggable candidates from the S. aureus MRSA252 proteome. The analysis revealed MurG (NCBI accession: Q6GGZ0) as the top match (Table 2). MurG, an essential glycosyltransferase involved in bacterial cell wall biosynthesis, displayed notable structural similarities to quercetin-binding glycosyltransferases, particularly within its active site. These similarities were confirmed through structural alignment of the AlphaFold-2 predicted structure of MurG with the glycosyltransferase 8SFW, highlighting the quercetin-binding region (Fig. 3a)27. Despite only having 26% identity and 46% similarity by alignment, comparison of the AlphaFold-generated 3D structure of MurG and 8SFW showed a high degree of similarity in their active sites, in particular in the quercetin binding region (Fig. 3b). In addition to MurG, MenB (sp|Q6GI37|MENB_STAAR) was identified as an essential protein via DEG screening (Table 2). However, MenB was deprioritized for further analysis due to (i) the absence of a structurally resolved active site, (ii) lack of experimental evidence linking it to quercetin binding, and (iii) insufficient functional data clarifying its specific role in bacterial survival despite its essentiality annotation in DEG.

Table 2 BLASTP search results highlighting protein homology between quercetin-Interacting proteins and S. aureus proteins.
Full size table
Fig. 3
figure 3

(a) Role of MurG in Biosynthetic pathway of Peptidoglycan in S. aureus (b) Structural Alignment of MurG (AF-Q6GGZ0 from AlphaFold-2) with 8SFW (from PDB), suggestive of structural similarity amongst both proteins. This image was made with BioRender (https://biorender.com/).

Full size image

Molecular docking reveals interaction between quercetin and MurG

Given the identification of MurG as the most likely quercetin binding target by structural alignment, we sought to further explore MurG-quercetin interactions using molecular docking. MurG acts as a glycosyltransferase that catalyzes the transfer of N-acetylglucosamine (GlcNAc) from UDP-GlcNAc to a growing lipid-linked oligosaccharide precursor, specifically converting Lipid I to Lipid II. Lipid I consists of a C55-undecaprenyl-phosphate linked to N-acetylmuramic acid (MurNAc) that is attached to a pentapeptide. MurG adds a GlcNAc residue to MurNAc, resulting in the formation of Lipid II. This addition is a critical step because Lipid II serves as the primary building block for the synthesis of the bacterial cell wall peptidoglycan layer27.

To investigate the quercetin and MurG interaction further, molecular docking was performed using the high-resolution MurG AlphaFold-2 model (AF-Q6GGZ0) with quercetin and UDP as ligands. Both ligands were optimized using Maestro software suite of Schrödinger and docking was targeted to the enzyme active site. Molecular docking showed interaction between quercetin, UDP and MurG, with quercetin shown to fit snugly into the binding pocket (Fig. 4). An XP Glide score of − 6.4 kcal/mol indicates a strong binding affinity for quercetin with MurG and multiple residues in the active site were seen to hydrogen bond with quercetin, including His14, Asp125, Gln285 and Arg287 (Fig. 4). In addition, a pi-anion interaction with Asp289, along with carbon-hydrogen bonds with Ser263, contributes to the robustness of the binding. To validate this interaction, we compared it with the binding affinity of MurG’s natural substrate, UDP-GlcNAc. Docking of only the GlcNAc moiety at the same site as quercetin yielded a comparable XP Glide score of − 6.16 kcal/mol, confirming that quercetin’s binding affinity is within the relevant range. These data suggest a strong interaction and binding affinity between MurG and quercetin.

Fig. 4
figure 4

Interaction of quercetin and UDP with MurG (AF-Q6GGZ0 from AlphaFold-2) predicted from computationally.

Full size image

MD simulations of MurG in complex with quercetin and UDP

Due to the strong binding affinity between MurG and quercetin, we sought to further investigate the stability and conformational changes of these complexes through MD simulations. MD simulations of the MurG protein with quercetin and UDP were conducted for 250 ns using Schrödinger Desmond (Fig. 5). The RMSD analysis provides insights into the structural stability of the protein and ligand complexes over the course of the simulation. For MurG, RMSD values indicated that after initial equilibration, the protein backbone stabilized at an average value of 2.5 Å after 25 ns, suggesting that the system reached a stable conformation relatively early in the simulation (Fig. 5a). For quercetin, the RMSD values relative to the protein ('Lig fit Prot’ in Fig. 5a) initially fluctuated but eventually stabilized around 5 Å, while the internal fluctuations of quercetin ('Lig fit Lig’ in Fig. 5a) were around 2 Å, indicating a stable interaction within the binding pocket. UDP exhibited a similar pattern with stable ligand RMSD values at approximately 2 Å; however, a notable peak between 75 to 90 ns in the 'Lig fit Prot’ RMSD suggests a transient conformational change in the ligand or its interaction with MurG before restabilizing (Fig. 5b).

Fig. 5
figure 5

Quantitative Assessment of Dynamic Interactions and Structural Stability in MurG Complexed with Ligands During MD Simulations. Subfigures (a, c) detail the RMSD and RMSF for MurG in complex with quercetin, respectively, while subfigures (b, d) display the RMSD and RMSF for MurG in complex with UDP. These metrics elucidate the conformational stability and residue flexibility within MurG throughout the ligand binding process.

Full size image

RMSF was utilized to characterize local dynamics along the MurG protein chain. Peaks on the RMSF plot reveal regions with the highest fluctuations during the simulation, typically observed in the N- and C-terminal regions (Fig. 5c). Secondary structure elements like alpha-helices and beta-strands, highlighted against red and blue backgrounds respectively, demonstrated less fluctuation, indicating higher structural rigidity compared to loop regions. Persistent secondary structure elements, such as alpha-helices and beta-strands, were maintained through over 70% of the simulation time, underscoring the structural integrity of MurG under dynamic conditions (Fig. 5c). Furthermore, residues interacting with the ligands are marked with green vertical bars, denoting crucial contact points that could influence ligand binding and efficacy. This comprehensive analysis suggests a complex interplay between MurG and its ligands, quercetin and UDP. It highlights the potential of quercetin to stably bind within the active site of MurG, potentially influencing its function and paving the way for further experimental validation of its inhibitory effects on the peptidoglycan biosynthesis pathway in S. aureus.

Interaction dynamics of MurG with quercetin and UDP during a 250 ns MD simulation, highlighting their stability and engagement at their respective binding sites. A comprehensive timeline of these interactions is presented, demonstrating that quercetin maintains an average of approximately 8 contacts, while UDP maintains around 16 contacts, indicating both ligands remain securely bound within their respective pockets without dissociation (Fig. 6a,b). However, it is important to note that the interactions made by quercetin and UDP cannot be directly compared as they occur at different sites of MurG. When looking at which MurG residues interact with the ligands across every simulation frame, interactions involving quercetin are stronger with residues such as His14, Asp284, Asp289, and Gln90 from around 35 ns onward, suggesting that the ligand achieves a better conformation for interaction as the simulation progresses (Fig. 6c,d). When interaction type was investigated (including hydrogen bonds, hydrophobic interactions, ionic bonds, and water bridges), hydrogen bonds appeared to dominate at the strongest interacting residues, with values exceeding 1.0 indicating that multiple contacts of the same type are occurring simultaneously at different interaction sites (Fig. 6c,d).

Fig. 6
figure 6

MD Simulation Interaction Analysis. Time-resolved interaction profiles of MurG with (a) quercetin and (b) UDP, showing the fluctuation in binding stability and key residues involved in ligand binding over the course of the simulation; interaction fraction profiles of MurG with (c) quercetin and (d) UDP.

Full size image

Investigation of ligand-atom interactions with protein residues further shows that certain residues engage in multiple interactions of a single type with the same ligand atom (Fig. 7a) and emphasizes the complex and dynamic nature of protein–ligand interactions within the enzyme. In addition, quercetin and UDP were shown to achieve a stable conformation within the binding pocket and maintain consistent compactness, according to measurements of average RMSD values and rGyr data (Fig. 7b). Further evidence for the maintenance of structural integrity of quercetin and UDP is found in their intraHB values, which measure internal hydrogen bonds within each ligand (Fig. 7). Taken together, Fig. 7 provides a detailed visualization of the interactions between the MurG protein and its ligands, quercetin and UDP, during a 250 ns MD simulation. The top panel displays a comprehensive schematic of ligand-atom interactions with protein residues, showing interactions that occur more than 20% of the time throughout the simulation. Notably, some interactions exceed 100%, indicating that certain residues engage in multiple interactions of a single type with the same ligand atom, highlighting the complex and dynamic interaction landscape within the enzyme’s active site. The bottom panels quantify various aspects of the ligands’ behavior throughout the simulation, illustrating the stability and interaction dynamics of Quercetin and UDP with the MurG protein. The RMSD measures the average deviations of the ligands from their initial positions, indicating that after initial adjustments, both ligands achieve a stable conformation within the binding pocket. The radius of gyration data suggests that both ligands maintain consistent compactness, indicative of stable intramolecular interactions. The count of internal hydrogen bonds within each ligand, reflected in the intraHB values, shows that both Quercetin and UDP preserve their structural integrity throughout the simulation, maintaining consistent internal hydrogen bonding. MolSA and SASA provide insights into how each ligand interacts with its environment, reflecting their respective hydrodynamic properties and interaction potentials. These values suggest how the ligands expose themselves to the surrounding solvent, influencing their interaction dynamics. PSA focuses on areas contributed by oxygen and nitrogen atoms, offers insights into the potential for hydrogen bonding and other polar interactions, which are crucial for the specificity and stability of ligand binding. Figure 7 elucidates the significant, stable interactions and the complex interplay of forces that potentially influence MurG’s enzymatic activity. It underscores the dynamic interactions of quercetin and UDP within the active site of MurG, supporting their roles as effective modulators of the enzyme’s function. This comprehensive analysis highlights the intricate and active engagement between the ligands and MurG, showing multifaceted interactions within the enzyme’s active site. Lastly, S1 Video provides a visual representation of the stable and consistent interaction between quercetin, UDP, and MurG.

Fig. 7
figure 7

Percentage of Protein–Ligand contacts and Ligand properties of quercetin and UDP while interacting with MurG.

Full size image

MM-GBSA assessment of quercetin interacting with MurG during MD simulation

To further interrogate the MurG-quercetin interaction and investigate the distribution of binding energies, we performed MM-GBSA analysis assessing the interaction between quercetin and MurG over a 250 ns MD simulation (Fig. 8). The violin plots depicted in the figure show the distribution of various energy components that are crucial for understanding the binding affinity between quercetin and MurG. These plots are particularly useful for visualizing the density and distribution of data points, with wider sections indicating a higher frequency of occurrences.

Fig. 8
figure 8

MM-GBSA Energy Component Analysis of MurG-Ligand Interactions. Violin plots illustrate the distribution of various energy components. (ΔGBind = Binding energy, ΔGCoulomb = Coulomb energy, ΔGCovalent = Covalent bond energy, ΔGHbond = Hydrogen-bonding correction, ΔGLipo = Lipophilic energy, ΔGPacking = Pi–Pi packing correction, ΔGSolvation = Solvation energy and ΔGvdW = Van der Waals energy).

Full size image

The analysis revealed that the total binding energy, ΔGBind, which is a primary indicator of binding affinity, becomes more negative throughout the simulation. This suggests that the interaction between quercetin and MurG strengthens over time, pointing to an increased binding affinity as the MD simulation progresses. Dividing the simulation into two halves, from 0 to 125 ns and from 125 to 250 ns, provides a granular view of how this interaction evolves. In the first half, ΔGBind values are somewhat dispersed, ranging from − 20 to − 60 kcal/mol, but they begin to cluster more densely towards the lower end of this range in the second half, indicating a growing stabilization of the ligand-enzyme complex. Electrostatic contributions, represented by ΔGCoulomb, exhibit considerable variability, ranging broadly from about − 20 to − 40 kcal/mol throughout the simulation. This variability underscores that electrostatic interactions play a significant but fluctuating role in the binding process. Conversely, Van der Waals interactions, denoted by ΔGvdW, maintain relatively stable values between − 20 to − 40 kcal/mol, consistently contributing to the interaction’s stability. Solvation energy, ΔGSolvation, shows minor fluctuations but remains largely in the negative range, suggesting that solvation dynamics generally favor the binding of quercetin to MurG. Other energy components such as hydrogen bonds, lipophilic interactions, and packing energies exhibit less dramatic variations but align with the overall trend of increasing binding affinity across the simulation. Overall, the data from the MM-GBSA analysis indicates that the interaction between quercetin and MurG not only strengthens over time, as evidenced by the increasingly negative ΔGBind, but also achieves greater stability due to contributions from various interaction forces including electrostatic, Van der Waals, and solvation dynamics. The increasingly negative values of ΔGBind throughout the duration of the MD simulation significantly reinforce the inference that the binding affinity of quercetin to MurG improves and stabilizes over time.

Comparative analysis of MurG conformational dynamics in apo and ligand-bound states

Having shown that MurG and quercetin exhibit a strong and increasingly stable interaction over time, we then sought to compare the behaviour of MurG in its bound and unbound (apo) state. MD simulations comparing the behavior of holo and apo forms of MurG (AF-Q6GGZ0) when complexed with quercetin and UDP are shown in Fig. 9. The evaluation highlights key structural and dynamic properties of the protein across different states. RMSD over time, shown in Fig. 9a, tracks the average deviations of the protein’s atomic positions from a reference conformation throughout the simulation. The results indicate that the protein in complex with quercetin and UDP displays lower RMSD values (2.26 ± 0.19 Å) than in its apo form (2.82 ± 0.33 Å), suggesting greater conformational stability when bound with the ligands. Figure 9b presents the RMSF, which measures the flexibility of each residue within the protein structure. It reveals that RMSF values are generally lower in the protein–ligand complex, indicating that the interaction with quercetin and UDP reduces residue fluctuations and stabilizes the protein structure. MolSA, depicted in Fig. 9c, measures the total molecular surface area of the protein. The analysis shows that MolSA remains largely consistent between the apo (14,702.17 ± 230.82 Å2) and ligand-bound (14,568.09 ± 201.44 Å2) states, implying that ligand binding does not significantly alter the protein’s overall molecular footprint. Figure 9d assesses SASA, which quantifies the surface area of the protein accessible to solvent molecules. This measurement is crucial for understanding protein-solvent interactions. The results indicate that the apo-protein exhibits higher SASA values (17,635.33 ± 248.97 Å2) compared to the ligand-bound state (16,406.43 ± 215.73 Å2). This reduction in SASA upon ligand binding suggests a decrease in solvent exposure, correlating with increased structural compactness and stability. PSA, evaluated in Fig. 9e, measures the area of the protein that can interact with polar molecules. Findings indicate that the apo-protein exhibits a higher PSA (8689.02 ± 165.84 Å2) compared to the ligand-bound state (8289.28 ± 161.32 Å2). This reduction in PSA in the ligand-bound complex suggests fewer exposed polar interactions, which may contribute to the protein’s stability upon ligand binding. Finally, the rGyr in Fig. 9f evaluates the compactness of the protein’s structure, reflecting how tightly the mass is distributed around the center of gravity. The ligand-bound protein shows a lower rGyr value (21.04 ± 0.09 Å) than the apo-protein (21.89 ± 0.23 Å), indicating a more compact and stable structure when complexed with quercetin and UDP. Figure 9 demonstrates that MurG in complex with quercetin and UDP is structurally more stable than in its apo form. This enhanced stability is evidenced by lower RMSD and RMSF values, decreased solvent exposure as shown by lower SASA, a more compact structure indicated by reduced rGyr, and reduced polar surface area, collectively suggesting that ligand binding significantly stabilizes the protein.

Fig. 9
figure 9

Comprehensive Post-MD Simulation Assessment of MurG in Apo and Ligand-Bound States. This figure evaluates the structural and dynamic properties of the MurG enzyme in its unbound (apo) form compared to its complexed state with ligands. The panels provide insights into various aspects of protein behavior: (a) RMSD over time showcasing the conformational stability, (b) RMSF illustrating residue flexibility, (c) Molecular Surface Area (MolSA) reflecting the overall molecular footprint, (d) Solvent Accessible Surface Area (SASA) indicating exposure to the solvent, (e) Polar Surface Area (PSA) measuring interactions with polar molecules, and (f) Radius of Gyration (rGyr) depicting the compactness of the protein structure.

Full size image

PCA analysis of apo form versus ligand-complexed states of MurG

The PCA of the MurG enzyme (AF-Q6GGZ0) from a post-MD simulation highlights differences between its behavior in the unbound (apo) form and when complexed with ligands (Fig. 10). PCA is used here to emphasize variation and extract prominent patterns from the dataset.

Fig. 10
figure 10

Post MD Simulation PCA of MurG Conformational Variability. PCA results show the variance in protein conformations of MurG in both Apo- and ligand-bound states over the course of a MD simulation. This analysis underscores the impact of ligand binding on the conformational behavior of MurG.

Full size image

The upper section of Fig. 10 examines MurG in its apo form. The PCA scatter plots display the first three principal components (PC1, PC2, PC3), with data points transitioning from blue to red, likely indicating progression through the MD simulation or different conformational states within the simulation. The scree plot shows the proportion of variance explained by each principal component. In this case, PC1 and PC2 show a broad spread in conformations, highlighting substantial conformational flexibility or diversity in MurG’s apo form. PC3 also displays noticeable variation but to a lesser extent. The scree plot emphasizes that the first few components account for a significant portion of the total variance, with nearly 54% captured by PC3, suggesting these components capture major movements within the protein structure.

The lower section of Fig. 10 focuses on MurG in complex with the ligands. The layout and color progression are similar to the upper section but likely denote interaction dynamics over time or various binding states of the ligands. This section reveals a more clustered distribution of data points compared to the apo form, indicating reduced conformational flexibility due to ligand binding. The narrowed range of motion in PC1 and PC2, along with the trend in PC3, suggests a decrease in the diversity of protein conformations compared to the apo form. The corresponding scree plot reveals a shift in variance distribution, with less variance explained by the first three components (38.4% by PC3), highlighting that ligand binding influences the dynamics and potentially stabilizes specific conformations over others.

In summary, Fig. 10 demonstrates that the apo form of MurG exhibits greater conformational diversity and flexibility, as evidenced by the broader distribution of PCA scores. In contrast, the complexed form with ligands displays a more restricted conformational space, suggesting that ligand binding stabilizes specific protein structures, potentially offering a more predictable and consistent interaction interface for therapeutic targeting. This difference could be critical in understanding how ligand binding impacts the structural and functional dynamics of MurG.

Conservation of MurG binding site Across S. aureus strains

Having shown that MurG and quercetin are predicted to have a strong binding interaction, similar to the native substrate UDP, and having investigated key binding residues in MurG from one S. aureus strain, we sought to investigate the presence of these residues in MurG from other key S. aureus clones. Multiple sequence alignment of MurG protein sequences from various S. aureus strains revealed a significant degree of conservation, particularly within the binding site residues (Fig. 11).

Fig. 11
figure 11

(a) MSA of MurG protein sequences from various S. aureus strains. The alignment reveals a significant degree of conservation within the binding site residues, which are highlighted in red boxes. Only minor variations are observed in other parts of the protein, indicated by white backgrounds, (b) Superimposition of 3D structures of MurG proteins from various S. aureus strains using ChimeraX. The structural overlay demonstrates a remarkable overlap of the binding site residues across different strains, confirming the structural conservation of the binding pocket. This high level of conservation suggests that the MurG binding site is preserved across the global S. aureus population.

Full size image

To further substantiate these findings, the 3D structures of MurG from various strains were superimposed using ChimeraX. The resulting structural overlay (Fig. 11) demonstrated a remarkable overlap of the binding site residues across the different S. aureus strains, confirming the structural conservation of the binding pocket. This high level of conservation suggests that the MurG binding site is preserved across the global S. aureus population, reinforcing the potential of MurG as a universal therapeutic target. These findings underscore the relevance of MurG as a broadly applicable target for antibacterial strategies, supporting the validity of targeting this enzyme across diverse S. aureus strains, including those beyond the MRSA252 strain studied initially. This broad applicability highlights the potential of MurG inhibitors, such as quercetin, in the development of new treatments for a wide range of S. aureus infections.

Discussion

The rise of multidrug-resistant S. aureus (MRSA) poses a critical challenge to global health, necessitating innovative strategies to identify novel antibacterial targets1,2,3,4. While natural products like quercetin have re-emerged as promising candidates due to their broad-spectrum activity6, their precise molecular targets in MRSA remain poorly characterized. Prior studies attributed quercetin’s antibacterial effects to β-lactamase inhibition6, allosteric modulation of PBPs7,8, or suppression of virulence via ClpP protease targeting 9. However, these mechanisms lack structural validation, and co-crystallization data for quercetin-bound S. aureus proteins remain absent9.

Quercetin exhibits well-documented antibacterial activity against S. aureus6,9, yet its precise molecular target(s) in this pathogen remain unresolved. To address this gap, we employed a targeted subtractive proteomics strategy to systematically identify quercetin-binding proteins within the S. aureus proteome. Our approach was guided by an evolutionary hypothesis: given that quercetin serves as a natural substrate for plant glycosyltransferases (e.g., UGT89C1 [PDB:6UD] and UGT202A2 [PDB:8SFW]), which catalyze its glycosylation as part of flavonoid metabolism6,27,28, we postulated that S. aureus might harbor a structurally analogous glycosyltransferase capable of binding quercetin. This rationale led us to prioritize MurG, a bacterial glycosyltransferase essential for peptidoglycan biosynthesis, due to its unexpected structural homology with plant quercetin-binding enzymes. By integrating subtractive proteomics—which filtered the MRSA252 proteome for essential, druggable targets—with molecular docking and dynamics simulations, we identified MurG as a high-affinity quercetin target. This strategy diverges from conventional subtractive proteomics workflows10,11,29,30,31 by focusing on cross-kingdom enzyme mimicry, a novel paradigm for antibiotic discovery that leverages evolutionary conservation to repurpose plant-natural product interactions against bacterial pathogens.

The choice of S. aureus MRSA252 (ST36/CC30) as the model strain was driven by its global epidemiological significance and unique genomic features. MRSA252, a representative of the EMRSA-16 clone, is a major cause of hospital-acquired infections and carries a genome enriched with novel virulence and resistance genes absent in other strains14,32,33. Its proteome (UniProt ID: UP000000596) is one of the most comprehensively annotated among MRSA isolates, facilitating large-scale computational analyses14,32,33. Furthermore, MRSA252’s resistance profile—including mecA-mediated methicillin resistance and blaZ-encoded β-lactamase activity34,35 exemplifies the therapeutic challenges posed by multidrug-resistant pathogens. By focusing on MRSA (rather than methicillin-susceptible S. aureus), this study addresses the urgent need for therapies that circumvent conventional resistance mechanisms.

Subtractive proteomics has emerged as a powerful tool for prioritizing therapeutic targets in pathogenic bacteria10,11,29,30,31. In this study, a subtractive proteomics approach was employed to identify potential drug targets within the MRSA252 proteome, focusing on their interaction with quercetin. This method effectively screened the MRSA proteome to pinpoint essential and druggable proteins, leading to the identification of MurG, a key enzyme in the peptidoglycan biosynthesis pathway, as a prime target. MurG is essential because it catalyzes a critical step in cell wall formation, converting Lipid I to Lipid II, a necessary component for bacterial cell wall integrity and survival35. The significance of glycosyltransferases, like MurG, lies in their role in catalyzing the transfer of sugar moieties to various substrates, a process crucial in both bacterial and plant systems35,36. To achieve this, a multi-tiered filtering strategy was implemented to identify MRSA252 proteins that are both essential and structurally homologous to plant quercetin-binding enzymes. First, the proteome was clustered using CD-HIT15, removing paralogs and generating a non-redundant dataset. Next, cross-referencing with the DEG database12 isolated proteins critical for survival. Finally, screening against DrugBank13 further refined the list to those with recognized or potential druggable pockets. The subsequent structural alignment with plant glycosyltransferases co-crystallized with quercetin revealed MurG (UniProt: Q6GGZ0) as the sole MRSA252 protein exhibiting significant homology to PDB-8SFW. MurG’s role in Lipid II biosynthesis27, combined with its structural resemblance to plant homologs, positioned it as a high-priority target. This approach contrasts with prior subtractive proteomics studies in Staphylococcus saprophyticus11 or Acinetobacter baumannii30, which focused on generic essential pathways without leveraging evolutionary mimicry.

The evolutionary parallels between plant glycosyltransferases and bacterial MurG, while compelling at the structural level, demanded rigorous validation under dynamic, near-physiological conditions. MD simulations served as a critical bridge between static docking poses and biologically relevant interactions10,37, offering insights into how quercetin’s binding mimics its natural role in plant systems while destabilizing MurG’s function in S. aureus. By extending beyond the static snapshots provided by docking, MD accounts for conformational flexibility and intermolecular forces, thereby highlighting whether quercetin’s plant-derived functional groups can stably occupy MurG’s UDP-GlcNAc-binding site in S. aureus. This approach parallels investigations into synthetic MurG inhibitors such as murgocil, which similarly depend on sustained interactions within the UDP binding pocket; however, in this study, quercetin occupies the GlcNAc binding pocket38. In analyzing the quercetin–MurG complex via MD, attention focuses on the persistence of hydrogen bonds and the spatial alignment of critical residues, such as His14 and Asp125, with comparable binding motifs in plant glycosyltransferases27,28. These interactions can suggest that quercetin effectively mimics UDP-sugar moieties, thus impeding Lipid II biosynthesis. Moreover, assessing contributions from van der Waals and electrostatic forces reveals how quercetin’s hydroxyl substituents might replicate key features of native UDP-sugars, potentially causing MurG to adopt a less catalytically active conformation. When PCA methods indicate a more rigid enzyme-ligand complex, it points to a “lock” mechanism reminiscent of plant glycosyltransferases in flavonoid glycosylation28. This cross-kingdom resemblance helps validate the original hypothesis that quercetin’s antibacterial activity leverages structural features evolved for flavonoid processing in plants. Such mechanistic evidence contrasts with prior claims that quercetin’s anti-MRSA effects stem solely from β-lactamase or ClpP inhibition6,9, thereby refining our understanding of its mode of action.

Lastly, the clinical relevance of this conservation is underscored by sequence alignment of MurG across S. aureus strains. In the last section of result section we found that residues critical for quercetin binding (His14, Asp125, Arg287) are 100% conserved, even in community-associated MRSA clones (e.g., USA300). This universality positions MurG as a viable target for quercetin derivatives against diverse MRSA isolates, circumventing strain-specific resistance mechanisms35,39.

Conclusion and future directions

Our findings demonstrate that quercetin acts as a promising inhibitor of MurG in S. aureus, supported by both subtractive proteomics and comprehensive molecular dynamics simulations. By exploiting structural similarities between plant glycosyltransferases and bacterial MurG, quercetin effectively competes with UDP-GlcNAc, disrupting a critical step in peptidoglycan biosynthesis. This mechanism diverges from previously proposed targets for quercetin, such as β-lactamases or ClpP proteases, and underscores the potential advantages of repurposing naturally occurring molecules against bacterial pathogens. The observed stability of the quercetin–MurG complex, alongside reduced protein flexibility, suggests that bacterial cell wall assembly could be impeded by quercetin in a manner akin to synthetic MurG inhibitors, but with the added advantage of utilizing a natural scaffold. This novel insight into quercetin’s antibacterial mode of action thus paves the way for developing enhanced therapeutic strategies against MRSA and potentially other resistant pathogens.

Future work should prioritize experimental validation of quercetin’s inhibitory effects on MurG through rigorous in vitro and in vivo assays. High-resolution structural studies, such as X-ray crystallography or site-directed mutagenesis of key residues, would confirm the critical interactions identified in silico. Such data would also guide rational design of quercetin-based analogs with improved pharmacokinetic properties, possibly by leveraging structural parallels between plant and bacterial glycosyltransferases. Assessing the synergistic potential of quercetin in combination with existing antibiotics may further illuminate its role in mitigating resistance. Overall, this research highlights the utility of integrating computational screening with targeted enzymatic assays, offering a blueprint for discovering new antibacterial agents that harness evolutionarily conserved architectures across distinct biological kingdoms.