S. aureus strains used were the quality control strain ATCC29213 and the strains preserved in our laboratory (234 S. aureus strains, including strains NM98, XF10, XF12, Wld19, JY4, JY39, and JY41, were isolated from mastitis milk samples from dairy farms in different Ningxia areas in 2017–2018)). AP (Sigma-Aldrich, United States) was dissolved in DMSO (Solarbio Technology Co., Ltd., China) to prepare a stock solution at 40,960 μg/ml.

The 234 S. aureus strains isolated and the strain ATCC29213 were removed from the −80°C refrigerator and screened using CHROMagar^™ Staph aureus (CHRO Magar, France). Single purple-red colonies were selected for purification and culture. DNA was extracted using the TIANamp Bacteria DNA Kit (Tiangen Biotech Co., Ltd., China). Bacterial genomic DNA was used as the template, and hla primers (F: 5′-GGT​TTA​GCC​TGG​CCT​TC-3′; R: 5′-CAT​CAC​GAA​CTC​GTT​CG-3′), 2× Taq PCR Master Mix (Vazyme Biotech Co., Ltd., China), and deionized water (Tiangen Biotech Co., Ltd., China) to prepare a 20-μL reaction system. Gene fragment amplification was performed at 58°C, used as the optimum annealing temperature. The sequencing results were uploaded to the NCBI website for sequence comparison by using the BLAST tool.

The hla sequences sequencing results were submitted to NCBI, and the gene sequences of known bacteria with high similarity to hla sequences were selected through BLAST homology search. The hla sequences detected in our region were used to construct phylogenetic trees by using CLASTALX 1.83 and MEGA X 10.2.6 and the neighbor-joining method. The gene trees were embellished through the ITOL online website (https://itol.embl.de/).

hla-containing S. aureus cells were cultured to the log phase. Then, 5 μL of the bacterial solution was removed, added dropwise vertically on dried rabbit blood agar plates, and incubated overnight at 37°C to observe the size of the hemolytic rings. Meanwhile, the hla-containing S. aureus cells were cultured to the log phase, and the supernatant was collected after centrifugation (10,000 g, 10 min). Then, the phosphate buffered saline (PBS), bacterial culture supernatant, and rabbit defibrinated erythrocytes were added to a 1.5-ml centrifuge tube, gently inverted and mixed, placed in a biochemical incubator at 37°C, and incubated overnight. The morphology and number of erythrocytes were observed under a biological microscope (Motic China Group Co., Ltd., China).

Several S. aureus strains with high hemolytic capacity were screened and cultured to the log phase. The bacterial solution, culture medium, and drug were mixed in 10-ml centrifuge tubes to achieve final drug concentrations of 0, 16, 32, 64, 128, and 256 μg/ml, respectively. After overnight incubation, the culture supernatant was collected through centrifugation (10,000 g, 10 min). Then, PBS, bacterial culture supernatant, and rabbit defibrinated erythrocytes were added to a 1.5-ml centrifuge tube, mixed gently and completely, and incubated for 15 min at 37°C. After incubation, the supernatant was centrifuged (5,500 g, 4°C, 1 min) and collected. The hemolytic activity was determined by measuring the absorbance of the supernatant solution at OD₅₄₃ nm. The hemolytic activity of the control group supernatant without AP was regarded as 100% and was used as a reference to calculate the hemolysis ratio. Changes in the hemolytic capacity of S. aureus before and after AP treatment were analyzed through the hemolysis assay.

Gene ontology (GO) of Hla protein was performed using online databases such as DAVID (https://david.ncifcrf.gov/) and STRING (https://cn.string-db.org/), with a view to enriching and annotating Hla protein in terms of biological processes (BP), molecular functions (MF), and cellular components (CC). STRING was used to search for related proteins that interact with Hla protein and contains experimental data, text mining results from PubMed abstracts, and comprehensive data from other databases, as well as results predicted using bioinformatics methods. Protein–protein interaction (PPI) analysis was performed using this database. Cytoscape plugin Cytohubba was used to select the top 10 highly connected genes (Hub) genes by association.

Two strains with more significant hemolysis and ATCC29213 were selected, and different AP concentrations were added to make the final drug concentration for each strain reach 0, 64, and 256 μg/ml. Their total RNA was isolated using TriZol reagent (Takara, Japan). The isolated total RNA was then treated with DEPC water (Tiangen Biotech Co., Ltd., China) and reverse transcribed using the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, United States). PCRs were performed using the SYBR^®Green Real-Time Fluorescence PCR Premix (Qiagen, Germany) in a 20-μL reaction system. The reactions were run using an ABI 7500 fluorescence quantitative PCR instrument (Applied Biosystems, United States), and the results were processed using the ΔΔCT method, in which gyrB was used as an internal reference gene. The primers used were gyrB forward, 5′-GCC​GAT​TGC​TCT​AGT​AAA​AGT​CC-3′, gyrB reverse, 5′-GAT​TCC​TGT​ACC​AAA​TGC​TGT​G-3′; hla forward, 5′-GGT​TTA​GCC​TGG​CCT​TCA​GC-3′, hla reverse, 5′-ACC​AGT​AAC​ATT​ACC​GTT​GAA​TCC​A-3′; clfA forward, 5′-CGA​TTC​TGA​CTC​CGA​CAG​TGA​TTC​C-3′, clfA reverse, 5′-CGC​TGT​CTG​AAT​CTG​AGT​CGC​TAT​C-3′; fnbA forward, 5′-AAG​ATC​AGC​AGA​TGT​AGC​GGA​AGC-3′, fnbA reverse, 5′-CTC​GTT​GTC​CTG​CAT​GAG​GTT​CTA​C-3′.

The secondary structure was predicted and analyzed using the online software Predict Protein (http://www.predictprotein.org/) as well as Jpred4 (http://www.compbio.dundee.ac.uk/jpred4/). Homology modeling and conservativeness analysis were performed using AlphaFold (Q2G1X0, Identity: 99.7%), SWISS-MODEL (A0A0D3Q6S7), and the online server Consurf Web Serve (http://consurf.tau.ac.il/). The conserved amino acids and solvent accessibility were mapped to the three-dimensional structure of Hla protein using the software PyMOL to analyze the linkage between amino acid conservativeness and protein structure and function. Using the programs ProtScale (http://us.expasy.org/tools/protparam.htmL) and Compute PI/MW (http://www.expasy.org/cgi-bin/protscale.pl) in the online software ExPASy, the physicochemical property analysis was performed.

The hla sequences were converted into amino acid sequences, and the tertiary structures of Hla protein were obtained through homology modeling by using AlphaFold (Q2G1X0, Identity: 99.7%) and SWISS-MODEL (A0A0D3Q6S7) (Tanaka et al., 2011), with the help of the molecular dynamics program GROMACS 2020.03 based on the AMBER99SB force field for 50-ns molecular dynamics simulations. The processed dominant conformation was used as the initial conformation for Hla protein docking. Small molecule structure files of AP were obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). Energy minimization based on the MMFF94 force field was performed, and structure files were generated and used for docking. Subsequently, the AutoDockTools 1.5.6 program was used to hydrogenate the receptor protein Hla. In the meantime, the number of rotatable bonds of ligand molecules and docking boxes (Grid Box) were set. Semi-flexible molecular docking to receptor proteins was performed sequentially using AutoDock Vina and Smina programs. The obtained results were verified on the basis of the consistency of the aforementioned three procedures. The final results, which were selected with a balance of scoring and conformational fit, were visualized using PyMOL (https://pymol.org/2/) and Discovery Studio 2020.

To further determine the stable binding of AP to Hla protein. After molecular docking, the Hla–AP complex was used as the initial conformation for 70-ns all-atom dynamics simulations performed using the classical molecular dynamics simulation software GROMACS 2020.03. The stability of the system was also measured on the basis of the root mean square deviation (RMSD), which indicates the degree of structural changes in protein molecules. Root mean square fluctuation (RMSF) was used to reflect the degree of freedom of motion of each atom in the molecule. The radius of gyration (Rg) reflects the degree of tightness of the protein before and after binding. The free energy landscape is used to characterize changes in free energy and the stability of conformation experienced by the protein during the simulation.

The g_mmpbsa calculation script (Kumari et al., 2014) for the Hla–Ap complex was applied using the Molecular Mechanics Poisson-Boltzmann surface area (MM-PBSA) method to calculate the binding free energy. This approach defines the binding free energy as follows: Δ G bind = G complex − G free - protein − G free - ligand

By definition, in solution, we can write the free energy of a molecule as follows: G = E gas − TS gas + G solvation where the solvation free energy can be further decomposed into polar and non-polar fractions as follows: G solvation = G polar + G nonpolar

In the MM-PBSA method, energy, and entropy contributions of the gas phase are calculated according to the MM method: Egas = E MM = E bond + E angle + E dihedral + E vdw + E coulomb S gas = S MM where E_bond, E_angle, and E_dihedral correspond to bond, bond angle, and dihedral interactions, respectively, and E_vdw and E_coulomb represent van der Waals and Coulomb electrostatic interactions, respectively.

The solvation energy in the MM-PBSA method contains two components: polar and non-polar solvation energies. The polar solvation energy is derived from the electrostatic interaction between the solute and solvent molecules. It is calculated using an implicit solvent model, where the solvent is considered a continuous medium, and the corresponding Poisson Boltzmann equation is linearized and solved numerically. G polar = G PB

The non-polar solvation free energy can be calculated using the empirical surface area method and is therefore also referred to as the surface solvation energy. The calculation requires knowledge of the solvent-accessible surface area A of the molecule with two empirical parameters γ and b: G nonpolar = G surface = γ A + b

Combining these aforementioned terms gives the free energy equation for MM-PBSA: G = E MM − TS MM + G PB + G surface

All assays were analyzed in triplicate. Statistical significance was analyzed using one-way ANOVA with IBM SPSS Statistics 25.0 statistical software, and a p value of <0.05 was considered statistically significant.

DNA of the 234 S. aureus strains and the quality control strain ATCC29213, isolated from mastitis milk samples from dairy farms in various Ningxia regions and kept in our laboratory, were used as template DNA for hla testing. The corresponding fragments of hla amplification are presented in Figure 1. The hla detection rate in the 234 strains was 98.29% (230/234).

The hla sequences detected in this region were uploaded to the NCBI database for blast comparison, and the results showed that the sequences had high homology with S. aureus SA1428, S. aureus 6538P, and S. aureus 13. The phylogenetic analysis showed that the detected sequences were the same as those carried by S. aureus BMS-2 α-hemolysin KX951417.1, S. aureus SA1428 CP048431.1, and S. aureus 2030RH1 CP039848.1, and the bootstrap value was 100 (Figure 2), indicating that the evolutionary tree constructed was highly reliable and reflects the high homology between the detected hla genes and those carried by the aforementioned strains.

To detect the hemolytic capacity of S. aureus, a hemolysis test was conducted using rabbit erythrocytes as they are more sensitive to the Hla hemolytic activity. Seven S. aureus strains with high hemolytic capacity were screened using rabbit blood agar plates, and their hemolytic capacity was further determined by microscopically observing the status of the rabbit blood erythrocytes treated with the bacterial solution. As shown in Figure 3A, a larger hemolysis ring indicates the stronger hemolysis ability of the strain. Further tests confirmed that strains with different sizes of hemolytic rings have different degrees of effects on the status of the rabbit blood erythrocytes. Figure 3B presents the blank control. In this case, rabbit blood erythrocytes with a diameter of approximately 7 μm were morphologically intact and more numerous in the field of view. Figure 3C shows rabbit blood erythrocytes treated with a strain with a medium-sized hemolysis loop; some erythrocytes in the field of view were morphologically incomplete and so appeared mutilated. Moreover, the number of erythrocytes in the field of view was considerably reduced compared with that in the blank control group. Figure 3D shows erythrocytes treated with a strain with a larger hemolysis loop; only few erythrocytes were scattered in the field of view, and the edges of the erythrocytes were rough and mutilated. Based on the aforementioned criteria, the strains NM98, XF10, XF12, Wld19, JY4, JY39, and JY41 were confirmed to be the most prominent hemolytic S. aureus strains in the region. The standard strain ATCC29213 was also found to have some hemolytic ability.

After qualitative detection of the S. aureus hemolytic ability, the two strains NM98 and XF10, which were more significant in hemolysis among the 7 strains, and the standard strain ATCC29213 were selected for quantitative detection. The effect of different AP concentrations on the hemolytic activity of the aforementioned three strains was further investigated through quantitative assays. The results showed (Figure 4) that the inhibitory effect of AP on the hemolytic activity of these three strains was concentration-dependent, with a highly significant (p < 0.01) decrease in hemolysis after AP treatment compared with the control without AP. Among them, the decrease in hemolysis was more in the NM98 strain than in the other two strains. Hemolysis was almost undetectable when the AP concentration reached above 64 μg/ml.

GO revealed (Figure 5A) that Hla protein in the BP was highly associated with S. aureus pathogenesis, hemolysis, and colonization adhesion. In terms of CCs, Hla protein are mainly localized in the extracellular region and cellular anatomical entities. In terms of MF, Hla is mainly involved in toxin activity. The STRING database was used to synthetically mine the results from PubMed abstracts, synthesize data from other databases, and predict results using bioinformatics methods to synthetically mine other proteins that have interactions with Hla protein (Figure 5B) and filter the top 10 genes by using Cytoscape plugin Cytohubba Hub genes. The proteins encoded by sarA, clfA, hlb, clfB, sigB, fnbA, saeS, saeR, and eno are strongly associated with Hla protein. Hla protein is significantly associated and interact more closely with several global regulators, virulence proteins, and biofilm formation proteins in S. aureus.

Fluorescence quantitative PCR was used to detect the effect of AP at concentrations of 0, 64, and 256 μg/ml on the transcriptional levels of hla, clfA, and fnbA genes of S. aureus, aiming to investigate the transcriptional and translational activities of different AP concentrations on these genes. There are also some relevant differences. The results showed (Figure 6) that the transcript abundance of hla genes in NM98, XF10, and ATCC29213 was highly significantly downregulated after AP treatment at 64 μg/ml compared with the control (p < 0.01), with the NM98 strain showing the largest fold downregulation. Moreover, the transcript abundance of the three strains continued to be slightly downregulated when the AP dose was increased to 256 μg/ml. The aforementioned results are in good agreement with those of the hemolysis test analysis. The transcript abundance of clfA and fnbA genes in the aforementioned three strains were highly significantly downregulated (p < 0.01) after AP treatment at 64 and 256 μg/ml compared with the control. Among them, the NM98 strain exhibited a greater downregulation of the transcript abundance of clfA and fnbA genes after AP treatment at 256 μg/ml. In summary, AP can highly significantly downregulate the hla, clfA, and fnbA transcript levels.

As shown in Figure 7A, Hla protein can bind to form heptameric aqueous proteins, which are secreted by bacteria in the monomeric form (Figure 7B). The Hla protein contains 1.36% α-helix, 48.81% β-fold, and 49.83% irregular coiled. The association of conserved amino acids with structural stability and functional performance was analyzed through homology modeling and conservativeness analysis of Hla protein and projection of conserved assignments onto tertiary structures according to color (Figures 7B, C). The regions SER41-ASP45, ARG55-THR60, ASP152-LEU157, and TYR191-MET197 were found to be highly conserved, as revealed through the analysis. The prediction of amino acid solvent accessibility by the convolutional neural network algorithm reveals that multiple residues used to stabilize the highly conserved nature of the specific protein conformation are wrapped inside: ILE14, ILE43, THR57, SER82, VAL168, TYR182, ASN193, LEU195-MET197, LEU219 PHE224, TYR256, and ILE284. Amino acids located on the surface and highly conserved play crucial roles in biological functions, mainly ASP44, ASP45, LYS50, PRO103, ASN105, TYR148, GLN150, LYS154, THR155, LEU157, ASN176, ASP185, SER186, GLN194 ASN209, SER225, ASP254, ASN269, LYS271, and ASP285.

The data of three dimensions of RMSD, Rg, and Gibbs free energy were used to plot the free energy morphology (Figure 8A), obtain the time point of the low-energy structure, and finally determine the stable conformation of Hla–AP at 3650 ps (Figure 8B). AP could bind to the “triangular” region that connects the “stem” and “cap” regions in the Hla monomer structure. The triangular region of the Hla monomer plays a major role in the conversion of the Hla monomer to a heptamer (Walker and Bayley, 1995; Song et al., 1996). Therefore, we can infer that the binding of the triangular region to AP leads to a conformational change in the Hla monomer, probably hindering the formation of the Hla protein heptamer. Conformational visualization analysis revealed that the AP molecule can form three hydrogen bonds with residues ASN105, SER106, and THR155 of Hla protein; bind with PRO103 through alkyl intermolecular forces; and form carbon hydrogen bonds with LYS154. Among them, PRO103 and LYS154 are highly conserved residues that play pivotal roles in biological functions, and when AP binds to these proteins, the Hla biological functions may be affected. In this binding mode, AP binds to the Hla protein with lower energy, which fully reflects the reliability of AP binding to Hla. Therefore, a strong interaction between the aforementioned residues and AP is likely. These results are also confirmed through binding free energy analysis.

Molecular dynamics simulations of the Hla protein and its Hla–AP complex were performed for 70 ns by using GROMACS 2020.03. The stability of the complex conformation was measured on the basis of the RMSD of the complex, and as shown in Figure 9A, the complex reached equilibrium after 30 ns. The mean RMSD value was 0.328 nm during the molecular simulation of 30–70 ns, and the RMSD value fluctuated less in this time period, reflecting the stability of the system and the reliability of binding of AP and Hla. The RMSF was used to reflect the degrees of freedom of motion of the atoms in the molecule and to assess the flexibility of the residues (Figure 9B). The RMSF values of the Hla protein residues bound to AP were lower than those of the residues in the free Hla protein. The RMSF values of most residues were <3 Å, indicating that AP has a conformational stabilization effect on the Hla protein. Meanwhile, the Rg analysis also confirmed the aforementioned results; when Rg is smaller, the protein structure is more dense and more stable (Figure 9C). The Rg value of the Hla–AP complex was in the range of 2.20–2.29 nm during the simulation from the beginning to 30 ns, with a wave at 30 ns, and then gradually tended to be stable during the subsequent simulation (mean = 2.25 nm), further confirming the stability of AP and Hla protein. The role of hydrogen bonds in complex formation is also critical; the more the number of hydrogen bonds formed, the more stable the complex is. The number of hydrogen bonds formed between Hla and AP fluctuates between 1 and 2 during the 70-ns simulation (Figure 9D). Up to four hydrogen bonds can be formed, and a certain number of hydrogen bonds will further stabilize this binding conformation.

The total binding free energy of the Hla–AP complex during the simulation was calculated using the MM-PBSA method (Table 1), where the sum of van der Waals free energy, electrostatic free energy, polar solvation free energy, and solvent-accessible surface area free energy was the total binding free energy. The total binding free energy of the complex was −81.479 ± 13.502 kJ/mol (Table 1), where the binding of AP to Hla protein was influenced by the van der Waals free energy and polar solvation free energy. Meanwhile, the binding free energy was decomposed to each residue (Figures 10A, B), and the key residues in the region of AP to Hla protein binding were identified, in which the two residues contributing most to the binding were LYS154 and PRO103 (binding energy contributions: −9.07 ± 0.40 kJ/mol and −5.62 ± 0.20 kJ/mol, respectively), followed by PRO151 (−2.47 ± 0.18 kJ/mol), LYS168 (−1.83 ± 0.10 kJ/mol), ILE156 (−1.46 ± 0.07 kJ/mol), and PHE228 (−1.35 ± 0.11 kJ/mol).