In this study, we started from 27 peptide molecules that were previously identified in the venom of the wild bee Hylaeus signatus and which were screened for antimicrobial activity (Nešuta et al., 2016). Three peptides were excluded as they contained D-amino acids. Additionally, one further peptide was excluded as its amino acid sequence was incompletely characterized. The PEP-FOLD 3.5 webserver was used to predict the peptide structures from the amino acid sequences of the peptides (Lamiable et al., 2016). This webserver uses a Hidden Markov Model suboptimal sampling algorithm to predict the peptide structures. The resulting peptide structures were used as the starting point for 20 ns molecular dynamics simulations, and the root mean square deviations of the alpha carbon atoms were calculated. The final frames of these molecular dynamics simulations were used for further studies.

The three-dimensional structure of the spike receptor-binding domain (RBD) of SARS-CoV-2 at 2.43 Å resolution (PDB ID: 6M0J) was retrieved from the Protein Data Bank. The protein structure was prepared by removing heteroatoms and water molecules using Discovery Studio (Discovery Studio, 2009). Additionally, energy minimization of the protein structure was performed using the AMBER14 (Case et al., 2014) force field in YASARA software (Krieger et al., 2013). Molecular docking and dynamics studies used this energy minimized protein structure.

The peptides and RBD protein were uploaded as ligand and receptor molecules respectively to the HDOCK, HPEPDOCK, and ClusPro web servers. HDOCK uses a combined template-based and template-free algorithm in an automatic manner (Yan et al., 2017), while HPEPDOCK uses a hierarchical algorithm (Zhou et al., 2018). After docking the peptides to the protein using HDOCK and HPEPDOCK, ClusPro was used to calculate binding energies (Comeau et al., 2004). The top ten peptides with the highest docking scores were selected for further evaluation. These ten peptides were docked to the RBD using the HawkDock web server. For further analysis, the four highest-scoring peptides were chosen based on their binding free energy and docking scores from HawkDock. On the HawkDock server, the HawkRank scoring system, the ATTRACT docking algorithm improved in groups, and MM/GBSA free energy decomposition analysis are implemented on a multipurpose platform (Weng et al., 2019). PyMOL and Discovery Studio (Discovery Studio, 2009) were used for structural analysis of the top four protein-peptide complexes.

Molecular dynamics simulations was performed in YASARA dynamics (Land and Humble, 2018) using the AMBER14 force field (Wang et al., 2004). The docked peptide-protein complexes were initially cleaned, optimized and the hydrogen bond network was oriented. A cubic simulation cell was created with periodic boundary conditions and the TIP3P water model was used (Harrach and Drossel, 2014). The simulation cell was extended by 20 Å in each direction beyond the protein-peptide complexes. The physical conditions of the simulation cells were set at 298 K, pH 7.4, and 0.9% NaCl (Krieger and Vriend, 2015). The initial energy minimization of the simulation cells were conducted by steepest gradient approaches with simulated annealing methods (5,000 cycles). The time step of the simulation was 2.0 fs. Long-range electrostatic interactions were calculated by the Particle Mesh Ewald (PME) method with a cut-off radius of 8.0 Å (Essmann et al., 1995; Krieger et al., 2006; Harvey and De Fabritiis, 2009). The simulation trajectories were saved every 100 ps. The simulations were run for 100 ns at constant pressure and temperature, using a Berendsen thermostat. The simulation trajectories were used to calculate root mean square deviations, root mean square fluctuations, solvent accessible surface areas, radii of gyration, and hydrogen bonds (Adji et al., 2021; Dutta et al., 2021; Obaidullah et al., 2021).

The per residue energy contribution of the peptide-protein structures was calculated using the pyDockEneRes webserver. The first, last, and average structure was extracted from simulations trajectories and utilized as input entry, and the average higher energy from the hotspot residues were tabulated. This tool can be utilized for the identification of the hotspot residues (Romero-Durana et al., 2020). E r e s i p y D o c k S C e l e _ v d w = ∑ s i d e c h a i n i s c E i s c e l e + E i s c v d w = − Δ Δ G ( Romero-Durana et al ., 2020

AllerTOP (Dimitrov et al., 2014a) and AllergenFP (Dimitrov et al., 2014b) webservers were used to predict the allergenicity of the peptides. AllergenFP implements five E-descriptor-based fingerprinting, whilst AllerTOP uses both k-nearest neighbor (kNN) and amino acid E-descriptors to predict the allergenicity of peptides. The ToxinPred (Gupta et al., 2013) webserver was utilized to predict the toxicity of the peptides.

Venom extracted from hymenopteran insects, including the solitary bee Hylaeus signatus, is a prominent source of antimicrobial peptides (AMPs). Several α-helical amphipathic AMPs, referred to as HYL, have been identified in the solitary bee, with antimicrobial activity against distinct strains of pathogenic bacteria and fungi as well as the ability to lyse cancer cells (Slaninová et al., 2011; Slaninová et al., 2012; Nešuta et al., 2016). Additionally, HYL had low hemolytic activity, suggesting that it may be safe for use in humans (Nešuta et al., 2016).

Notably, the earlier-discovered antimicrobial peptides, including cecropin-A from Hyalophora cecropia, exhibited antiviral activity (Both, 1980). Cecropin-A was shown to have antiviral activity against HSV-1, HIV-1, and JV, whereas melittin from A. mellifera had antiviral activity against influenza A viruses, HSV-1, HIV-1, JV, TMV, RSV, VSV, enterovirus, and coxsackievirus (Albiol Matanic and Castilla, 2004). In addition, alloferon 1 and alloferon 2 from C. vicina were effective against influenza viruses (Chernysh et al., 2002), TnGlv1 and TnGlv2 from Trichoplusia ni were effective against AcMNPV (Moreno-Habel et al., 2012), attC and dptB from Drosophila melanogaster were effective against SINV (Huang et al., 2013), and a myristoylated peptide from H. virescens was effective against HIV-1 and HSV-1 (Ourth, 2004; Feng et al., 2020). Therefore, we hypothesized that the HYL antimicrobial peptides from the solitary bee could have antiviral activity and consequently devised an in silico study to test this hypothesis. We found that several HYL peptides bind the RBD of SARS-CoV-2 in silico with favorable binding energy and stable conformations through 100 ns of molecular dynamics simulation.

Molecular docking can be utilized to rapidly assess the binding affinities and modes between a target substrate, such as a protein, and diverse ligands, including small peptides, to assist antiviral drug design (Willett and Glen, 1995). Virtual screening can identify the preferred binding orientation, optimal conformation, and binding sites of protein and peptide molecules (Alonso et al., 2006; Liu et al., 2018). Our computational approach confers improved structural precision and rapid prediction of both the existence and strength of binding through the use of multiple diverse docking algorithms aimed at achieving accuracy (Bartuzi et al., 2017).

The predicted peptide structures from PEP-FOLD 3.5 are shown in Supplementary Figure 1. The peptide structures were optimized by molecular dynamics simulations in which the root mean square deviations of the peptides were found to be small (Supplementary Figure 2). After predicting the peptide structures, all 23 peptides were docked to the RBD of SARS-CoV-2 by using the HDOCK, ClusPro, and HPEPDOCK webservers (Supplementary Table 1). Binding energies and/or docking scores were recorded for the ten highest-scoring peptides (Supplementary Table 2). Moreover, the binding free energies and docking scores from HawkDock (Supplementary Table 3) were used to select the four highest-scoring peptides (Figure 1). The highest docking score was −2,848.16 kcal/mol for peptide P20, followed by −2,815.13, −2,735.49, and −2,660.34 kcal/mol for P18, P21, and P2, respectively. The maximum binding free energy was −40.38 kcal/mol for the P20 peptide, followed by −37.2, −31.17, and −28.15 for P18, P21, and P2, respectively (Table 1 ).

The interactions between the RBD protein and the four highest-scoring peptides are shown in Figure 2. The P2:RBD complex had five hydrogen bonds at the RBD residues ALA363, VAL367, SER371, SER373, and ASN343, and six hydrophobic interactions at ASP364, LEU368, PHE374, TRP436, PHE342, and LEU441. The P18:RBD complex formed six hydrophobic interactions at TRP436, PHE374, LEU441, LEU368, VAL367, and LEU335, and three hydrogen bonds at SER371, ASN343, and GLU340. The P20:RBD complex had five hydrophobic interactions at TRP436, PHE374, PRO337, VAL367, and LEU335, and four hydrogen bonds at ASN440, ASN343, PHE342, and SER371. The P21:RBD complex had seven hydrophobic interactions at VAL503, TYR508, LYS378, ARG408, ALA411, TYR380, PRO412 position, one unfavorable bond at SER373 and three hydrogen bonds at SER375, VAL407, ILE410 position ( Table 2; Figure 2 ).

The interactions between the RBD protein and the remaining six of the top ten peptides are shown in supplementary figure 3. The P3:RBD complex formed two hydrogen bonds at ALA344 and ASN343, five hydrophobic interactions at LEU441, LEU335, PHE374, VAL362, and VAL367, and one unfavorable bond at ARG509. The P6:RBD complex had five hydrophobic interactions at LEU441, PHE374, LEU368, VAL367, and LEU335, and six hydrogen bonds at ASN440, ARG509, TRP436, ASN343, PHE342, and GLY339 position. The P16:RBD complex had four hydrophobic interactions at CYS391, CYS525, ALA522, and VAL362, two hydrogen bonds at ASP389 and THR523, one salt bridge at ASP389 and one unfavorable bond at CYS361. The P17:RBD complex had five hydrophobic interactions at LEU335, VAL367, LEU368, PHE342, and TRP436, and three hydrogen bonds at SER371, ASN343, and ASN440. The P22:RBD complex formed five hydrogen bonds at ASN440, ASN343, SER371, GLU340, and ASP364, three hydrophobic interactions at TRP436, VAL367, and LEU335, and one salt bridge at GLU340. The P23:RBD complex formed six hydrophobic interactions at ALA522, THR333, CYS391, CYS525, VAL362, and LEU390 and two hydrogen bonds at GLY526 and ASP389 (Supplementary table 4, Supplementary figure 3).

Molecular dynamics simulations were conducted to explore the structural stability of the docked peptide-protein complexes. The root mean square deviations (RMSD) from the peptide-protein systems were calculated from the simulation trajectories. Figure 3A indicates that the peptide-protein complexes involving the P2, P18, P20, and P21 peptides had an initial upwards RMSD trend, suggesting flexibility of the complexes. The peptide-protein systems subsequently stabilized after 30 ns and maintained their integrity for the remaining 70 ns of the simulations. The P21 complex had a comparatively higher RMSD than the other three peptide complexes, which indicates the more flexible nature of this complex. However, all four peptide-protein complexes had RMSD lower than 2.5 Å, which indicates the stable nature of the complexes over the whole simulations.

The solvent-accessible surface area (SASA) of the complexes was also explored. Figure 3B indicates that the P20 and P21 complexes increased in SASA upon binding with the protein target, and after 30 ns the P21 complex decreased in SASA, maintaining a similar SASA trend to the rest of the complexes. The P20 complex stabilized after 50 ns and maintained its integrity. The radius of gyration (Rg) of the complexes was also explored to understand the compactness of the complexes. High Rg indicates an extended nature whereas a low Rg indicates a more compact structure. Figure 3C indicates that P20 had slightly higher Rg than other complexes which indicates that this complex is less compact. The other complexes had smaller deviations, which indicate the rigid nature of these complexes. Hydrogen bonding can play a crucial role in determining the stability of protein complexes. Figure 3D indicates that all complexes had a stable hydrogen bonding profile over the 100 ns of simulation.

The root mean square fluctuations (RMSF) of the complexes were explored to understand the flexibility on a residue-by-residue basis. Figure 4 indicates that almost every residue of the complexes had an RMSF lower than 1.4 Å, indicating the stable nature of the protein-peptide systems.

After 100 ns of simulation time, the P2:RBD complex was stabilized by five hydrogen bonds at ASN343, SER371, ASN370, ASP364, and ALA363, as well as five hydrophobic interactions at LEU441, TRP436, LEU368, VAL367, and LEU335. The P18:RBD complex had two hydrogen bonds at SER371 and GLU340, as well as six hydrophobic interactions at TRP436, PHE374, PHE342, LEU368, VAL367 and LEU335. The P20:RBD complex formed two hydrogen bonds at ASN440 and SER371, in addition to three hydrophobic interactions at TRP436, LEU368 and LEU335. The P21:RBD complex had four hydrogen bonds at ILE410, LYS378, TYR508 and PHE374, in addition to three hydrophobic interactions at PRO412, ALA411 and VAL407 (Supplementary table 5, Supplementary figure 4).

In the P2:RBD simulation, the interactions at LEU441, TRP436, ASN343, SER371, LEU368, VAL367, ASP364, and ALA363 were occupied throughout the whole simulation period. Similarly for TRP436, PHE374, SER371, GLU340, LEU368, VAL367, and LEU335 in the P18:RBD simulation, ASN440, TRP436, SER371, and LEU335 in the P20:RBD simulation, and PRO412, ALA411, ILE410, LYS378, VAL407, and TYR508 in the P21:RBD simulation ( Table 2; Figure 2 , Supplementary table 5, Supplementary figure 4).

The active sites of the RBD protein include amino acid residues 340, 374, 375, 378, 403, 420, 477, 478, 499, 543, and 546 (Lan et al., 2020). At the start of the simulation, the P2 peptide bound to the active site of RBD via a hydrophobic interaction with the conserved residue PHE374 (bond distance 3.31360), whereas the P18 peptide interacted with the RBD active site at both PHE374 (Bond distance 3.41730) and GLU340 (bond distance 2.70389) with a hydrophobic interaction and a hydrogen bond respectively. The P20 peptide bound to the active site of RBD at PHE374 (bond distance 3.44294) via a hydrophobic interaction, and the P21 peptide interacted with the RBD active site at SER375 (bond distance 2.49462) and LYS378 (bond distance 3.83247) position with a hydrogen bond and a hydrophobic interaction respectively (Table 2; Figure 2).

After 100 ns of simulation time, the P18 peptide was bound to the active site of RBD at two conserved residues, PHE374 (bond distance 3.44205) and GLU340 (bond distance 2.73681), via a hydrophobic interaction and a hydrogen bond respectively. The P21 peptide interacted with the active site of RBD via LYS378 (bond distance 2.81702) and PHE374 (bond distance 2.87497) via two hydrogen bonds (Supplementary table 5, Supplementary figure 4). These interactions of the peptides with the active site of the RBD indicate that these peptides may inhibit the RBD with high binding affinity.

Previously suggested multifunctional peptides, including ALPEEVIQHTFNLKSQ and DIENLIKSQ from Bacillus-fermented soybean, failed to bind to the active site of the RBD. However, our top four peptides bound to the active site. In addition, P18 and P20 had higher binding energy than DIENLIKSQ when employing the same HawkDock server (Padhi et al., 2021). The computationally designed peptide APASMFLGKGDHEILM made no interactions with the active site when docked using the same HPEPDOCK server that we used for our top four peptides (Priya et al., 2021). Another peptide modeling and screening study suggested that AVP0671 can bind to the RBD, albeit not at the active site, and the HDOCK scores were not reported. Meanwhile, our top four suggested peptides not only exhibited active site binding but also had a more favorable HDOCK docking score (Rathod et al., 2020). Antiviral peptides S2P25 and S2P26 were predicted to form one bond at the active site, but had lower docking interaction energies than our P21 and P18 peptides when comparing the Cluspro scores (Chowdhury et al., 2020). Additionally, AVP1795, identified from the computational screening of 645 antiviral peptides, exhibited only one interaction (ARG403) at the active site, which was occupied for ≥90% of the time in MD simulations, whereas our peptides were stably bound over the whole simulation period (Hossain et al., 2021).

Multiple methodology developed to understand the hotspot residues in protein-protein interactions which are based on snapshot from MD simulation sampling (Gohlke et al., 2003). In addition to pyDockEneRes, mm_pbsa.pl tools successfully implemented for hotspot resdiues identifications from cytokines and receptor interface (Du et al., 2020) as well as RBD from SARS-CoV-2 and antibodies (Yang et al., 2021). Moreover, the pyDockEneRes tools enables the hotspot residues identifications from the protein-protein interactions which is key to understand the biological process at molecular level (Romero-Durana et al., 2020). In our study, the P2 peptide and RBD complex, four important residues contribute to overall energy contributions: Leu335, Glu340, Val367, and Phe374; in the P18-RBD complex, hotspot residues include Asp164, Ser373, Ser375, and Trp380 (Table 3). Leu335, Asn343, Asp364, Val367, Ser373, Phe374, and Glu484 were major energy contributing residues in the P20-RBD complex, whereas Leu335, Phe342, Asn343, Asp364, Val367, Ser373, Phe374, and Glu484 were key energy contributing residues in the P22-RBD complex.

An allergenic antigen can activate Th2 cells resulting in stimulation of B cells to generate immunoglobulin E (IgE) that binds to FcεRI and activates eosinophils leading to inflammation and tissue shrinkage (Dimitrov et al., 2013). The online AllerTOP tool, which assesses allergens using E-descriptors affined with amino acid attributes, was utilized to predict the allergenicity of our peptides (Dimitrov et al., 2014a). Three distinct web-based approaches were utilized to assess toxicity (Table 4). According to the AllergenFP v1.0 webserver, P2 is a probable allergen, whereas P18, P20, and P21 are probably non-allergenic. The same result was obtained using the AllerTOP v2.0 webserver. In addition, all four peptides (P2, P18, P20, and P21) were categorized as non-toxic by the ToxinPred web server.