Pirkhezranian et al. BMC Genomics (2020) 21:60 https://doi.org/10.1186/s12864-020-6458-7 RESEARCH ARTICLE Open Access Interaction of camel Lactoferrin derived peptides with DNA: a molecular dynamics study Zana Pirkhezranian1, Mojtaba Tahmoorespur1*, Xavier Daura2,3, Hassan Monhemi4 and Mohammad Hadi Sekhavati1* Abstract Background: Lactoferrampin (LFampin), Lactoferricin (LFcin), and LFchimera are three well-known antimicrobial peptides derived from Lactoferrin and proposed as alternatives for antibiotics. Although the intracellular activity of these peptides has been previously demonstrated, their mode of action is not yet fully understood. Here, we performed a molecular dynamics simulation study to understand the molecular interactions between camel Lactoferrin derived peptides, including CLFampin, CLFcin, and CLFchimera, and DNA as an important intracellular target. Results: Our results indicate that all three peptides bind to DNA, albeit with different propensities, with CLFchimera showing the highest binding affinity. The secondary structures of the peptides, modeled on Lactoferrin, did not undergo significant changes during simulation, supporting their functional relevance. Main residues involved in the peptide-DNA interaction were identified based on binding free energy estimates calculated over 200 ns, which, as expected, confirmed strong electrostatic interactions between DNA phosphate groups and positively charged peptide side chains. Interaction between the different concentrations of CLFchimera and DNA revealed that after binding of four copies of CLFchimera to DNA, hydrogen bonds between the two strands of DNA start to break from one of the termini. Conclusions: Importantly, our results revealed that there is no DNA-sequence preference for peptide binding, in line with a broad antimicrobial activity. Moreover, the results showed that the strength of the interaction between DNA and CLFchimera is concentration dependent. The insight provided by these results can be used for the rational redesign of natural antimicrobial peptides targeting the bacterial DNA. Keywords: Antimicrobial peptide, DNA binding, Lactoferrin, Molecular dynamics simulation, CLFchimera Background net positive charge and an amphipathic structure [2–4]. Antibiotic resistance is becoming a serious global health One subgroup of AMPs includes peptides derived from problem, as infections by multidrug-resistant pathogens large proteins. Lactoferrampin (LFampin) and Lactoferri- are increasing at an alarming pace. There is thus an ur- cin (Lfcin) are two well-known antimicrobial peptides gent need to introduce new and safe antimicrobial derived from the Lactoferrin protein (LF) [5, 6]. These agents, including antimicrobial peptides (AMPs), as al- two cationic antimicrobial peptides have activity against ternatives to current antibiotics [1]. AMPs have evolved a broad spectrum of microorganisms including bacteria, as a natural defense mechanism for fighting microbial fungi and viruses [5, 6]. infections [1]. They are a diverse group of innate im- We have recently reported that a camel Lactoferrin mune system molecules that exist in all organisms [1]. chimera (CLFchimera) resulting from the fusion of the AMPs usually contain 12–50 amino acid residues, have a C-terminal ends of camel Lactoferricin 17–30 (CLFcin) and camel Lactoferrampin 265–284 (CLFampin) using * Correspondence: Tahmoores@um.ac.ir; sekhavati@um.ac.ir the side chain of lysine as linker to the second peptide, 1 Department of Animal Science, Faculty of Agriculture, Ferdowsi University of Mashhad, Mashhad, Iran has a broad-spectrum activity against both Gram- Full list of author information is available at the end of the article positive and Gram-negative bacteria [7–9]. Furthermore, © The Author(s). 2020 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Pirkhezranian et al. BMC Genomics (2020) 21:60 Page 2 of 14 Reyes-Cortes et al. (2016) showed that this chimeric experimentally, which relates directly to the thermody- peptide mediated its antibacterial activity by entering the namics and kinetics of the system. Herein, homology- cytoplasm through translocation across the bacterial based models were constructed for camel LFcin, LFam- membrane and possibly interacting with internal organ- pin and LFchimera and their interaction with DNA was elles [10]. To date, there has been no precise explanation analyzed using MD simulation, as a means to under- for the mechanisms underlying the antimicrobial peptide stand a reported intracellular mechanism of action of function, but it is known that DNA is one of the most these peptides. The findings of this study provide basic important intracellular targets for AMPs [11]. Thus, nu- directions for future studies regarding the function of cleic acids have been proven as intracellular targets for AMPs with intracellular activity and their potential re- some antimicrobial peptides such as MDpep9 [11], design with therapeutic purposes. Buforin I [12], Indolicidin [13], Cecropin PR39 [14], and NK18 [15]. Previous computational studies also showed Methods that Buforin II (from the stomach tissue of the Asian Molecular structure models toad bufo garagrizans) and Lasioglossin II (derived from An arbitrary 12-bp DNA sequence adopting a canonical bee venom) had considerable affinity for DNA [16, 17]. BDNA structure (entry 1BNA from the Protein Data Considering these reports, Uyterhoeven et al. (2008) Bank) was initially chosen for this study (Table 1). Since showed using MD simulation that Arg 2, Arg 14 and the simulations using this sequence indicated that the Arg 20 of Buforin II were mainly responsible for the interaction with the peptides is DNA-sequence inde- interaction with DNA and using Fluorescent Intercalator pendent (through the backbone phosphates), no add- Displacement (FID) assay indicated that disrupting itional sequences were used in the study. The length of Buforin-DNA interactions generally decreased the anti- the DNA allows the interaction with more than one pep- bacterial activity of the peptide [16]. In another study, tide, as shown in the Results and Discussion. The native Tang et al. (2009) demonstrated that MDpep9, a recently structure of camel Lactoferrin was also retrieved from discovered antimicrobial peptide derived from larvae of the PDB (entry 1DTZ). Camel Lactoferrin was used as a housefly (Musca domestica), a traditional food source in template structure for peptide modeling. Lactoferricin, China, is able to form bonds with DNA phosphate Lactoferrampin and CLFchimera were modeled with groups and insert between the base pairs of the DNA Modeller 9.2 [19] and PEP-fold server (http://bioserv. helix [11]. rpbs.univ-paris-diderot.fr/services/PEP-FOLD/) [20] Although the intracellular activity of Lactoferrin de- using default parameters. The quality of the models was rived peptides has been previously demonstrated [10], examined with PROCHECK (http://servicesn.mbi.ucla. the exact mechanism of action has not been yet estab- edu/PROCHECK/) [21]. lished. As recognition of specific DNA sequences by proteins is highly complex, involving structural, ener- Molecular dynamics simulations getic and dynamic aspects, the interaction cannot be eas- The complexes BDNA-CLFampin, BDNA-CLFcin, and ily characterized at the atomic level by experimental BDNA-CLFchimera were studied by molecular dynamics approaches alone [18]. The use of computational tech- simulation with the GROMACS 2016.1 package [22–24] niques such as molecular dynamics (MD) simulations and CHARMM27 force field [25]. Peptides and DNA provides complementary information, inaccessible were solvated in a cubic box using the Simple Point Table 1 Details of sequences, simulation lengths and replicates System Composition Simulation Length (ns) Replicates Box size (nm3) BDNA CGCGAATTCGCG 200 3 4.28 CLFampin DLIWKLLVKAQEKFGRGKPS 200 3 4.98 CLFcin KKCAQWQRRMKKVR 200 3 4.89 CLFchimera DLIWKLLVKAQEKFGRGKPS 200 3 5.19 KRVKKMRRQWQACKKS 1-CLFampin / BDNA 200 3 6.65 1-CLFcin / BDNA 200 3 6.38 1-CLFchimera / BDNA 200 3 6.90 2-CLFchimera / BDNA 200 3 7.81 3-CLFchimera / BDNA 200 3 8.36 4-CLFchimera / BDNA 200 3 8.96 Pirkhezranian et al. BMC Genomics (2020) 21:60 Page 3 of 14 Charge (SPC) water model [26]. To neutralize the overall each peptide and BDNA, including root mean square de- charge of the systems, Na and Cl ions were added as ap- viation (RMSD), root-mean-square-fluctuations (RMSF), propriate. Periodic boundary conditions were applied solvent accessible surface area (SASA), contacting sur- from this step onward. The system was first energy min- face area (CSA), hydrogen bonds, salt bridges, and center imized using the steepest descent algorithm to relax of mass distance were analyzed during the simulation high-energy contacts. After energy minimization, the using GROMACS built-in tools. An RMSD-based con- system was simulated under the NPT ensemble for 500 formational clustering algorithm, using the gmx-cluster ps, with initial velocities taken from a Maxwell- module of GROMACS, was applied to extract represen- Boltzmann distribution corresponding to 100 K. During tative structures. The clusters were obtained using a cut- this initial simulation time, the peptide and DNA atoms off of 1.5 Å for the RMSD to the centroid. were positionally restrained while the temperature was gradually increased from 100 K to 300 K at 1 atm. Bond Binding free energy estimates lengths were constrained for all atoms using the LINCS Binding free energies were estimated for BDNA- algorithm (SETTLE for water), allowing a time step in CLFampin, BDNA-CLFcin, and BDNA-CLFchimera the leap-frog integrator of 2 fs. Temperature and pres- complexes using molecular mechanics energies in com- sure were couple to the reference values using the Nosé- bination with Poisson-Boltzmann and surface area con- Hoover and Parrinello-Rahman algorithms, respectively tinuum solvation (MM/PBSA). The calculations were [27–29]. Additional 100 ps at 300 K and 1 atm, without performed with the g_mmpbsa program (https://rashmi- position-restraints, were subsequently run. In the pro- kumari.github.io/g_mmpbsa/) [30], using the single tra- duction phase, the equilibrated systems were run in the jectory approach. The solute dielectric constant was set NPT ensemble at 1 atm and 300 K for 200 ns. Long- to 8 [31] and the ionic strength was chosen to corres- range electrostatics were evaluated using the Particle pond to a NaCl concentration of 150 mM. The calcula- Mesh Ewald (PME) algorithm [28]. The real space com- tion of the Gpolar solvation term was performed with the ponent of PME and the van der Waals interactions were linearized Poisson-Boltzmann (PB) equation using a grid calculated with a cutoff of 1.0 nm. Three replicates of resolution of 0.05 nm and the bondi set of atomic radii. 200 ns were run per system, with different initial config- The Gnonpolar term was calculated with the SASA model urations generated by insertion of the peptides at ran- using default parameters [30]. The entropic component dom positions. The simulations performed and their of the binding free energy was disregarded. The average lengths are detailed in Table 1. Dynamics and stability of binding energy and its standard deviation were Fig. 1 Structural fluctuation analysis. a RMSD as a function of time; b RMSF per residue; c Cartoon structure of CLFchimera (C1), CLFampin (C2) and CLFcin (C3) at 0, 100 and 200 ns (red, blue and green, respectively); d Sequence alignment of the three peptides. RMSD and RMSF quantities were computed for structures at 0.1-ns intervals from the 200-ns simulations after least square fitting to the initial structure using the backbone atoms Pirkhezranian et al. BMC Genomics (2020) 21:60 Page 4 of 14 calculated with the MmPbSaStat.py python script were relaxed along independent 200-ns simulations, per- (http://rashmikumari.github.io/g_mmpbsa/) using the formed in triplicate (Table 1). To that end, the homology second half of the simulations production phase (100 to models obtained for the peptide structures were first ex- 200 ns), by taking 1000 snapshots at 100-ps intervals. To amined for overall quality. The Ramachandran plot for estimate the contribution of each residue to the total CLFampin, CLFcin and CLFchimera revealed that 93.3, binding free energy, the MmPbSaDecomp.py python 100.0 and 93.5% of the residues were situated within the script was used [30, 32]. It should be noted that this ap- most favored region, respectively, whereas the remaining proach represents a crude estimate of the binding free residues were found within the additional allowed energy that, most certainly, severely overestimates the region. real value, as noted by several authors [33]. However, the limitations of the approach are likely to affect the re- Structural fluctuation analysis lated systems studied here in similar ways and are there- Root-mean-square deviations from the initial structure fore expected to allow for a qualitative comparison. of the peptide as a function of simulation time and root- mean-square fluctuations of peptide residues are pre- Results and discussion sented, for one of the 200-ns replicates, in Fig. 1. The Molecular dynamics simulation in aqueous solution of the behavior of these quantities in the remaining replicates isolated peptides and BDNA is consistent with the observations made here (see Add- Before simulating the interaction between the different itional file 1: Figure S1 and Additional file 2: Figure S2). peptides and BDNA, the individual model structures The RMSD values are stable after the initial 100 ns, the Fig. 2 COM distance analysis. a COM distances between CLFcin, CLFampin and CLFchimera and DNA along 200 ns. b Structures at times t = 0 (cyan) and t = 200 ns (purple): (B1) CLFcin-DNA, (B2) CLampin-DNA, and (B3) CLFchimera-DNA Pirkhezranian et al. BMC Genomics (2020) 21:60 Page 5 of 14 larger peptide CLFchimera showing higher RMSD and Center of mass distances RMSF values. CLFchimera was obtained from the C- The center of mass distance between peptide and DNA term-C-term fusion of CLFampin 265–284 and was calculated as a function of time (Fig. 2a). Side view CLFcin17–30, using a lysine (Lys21) as linker [7, 8]. Fig- of snapshots of the first and last configurations are ure 1b shows that the global fluctuations of the corre- shown in Fig. 2b. In all three replicates, COM distances sponding sequences in the shorter peptides are lower in were initially around 3, 3 and 4 nm for CLFcin, CLFam- general than in the fusion peptide, as expected in light pin and CLFchimera, respectively. The peptides instantly of the structures shown in Fig. 1c. It is worth noting that moved toward the DNA grooves and COM distances de- the shorter CLFcin adopts a more stable helical structure creased rapidly. The three replicates show some differen- than CLFampin when isolated in solution, to become tial behavior in terms of final distance and convergence more flexible in the fusion peptide. Structures from the (Additional file 3: Figure S3A and Additional file 4: Fig- stable part of the 200 ns simulations with all residues in ure S4A), as well as in terms of position and orientation the most favored regions of the Ramachandran plot were (Additional file 3: Figure S3B and Additional file 4: Fig- used as initial structures for the corresponding simula- ure S4B), suggesting that the binding is not specific, as tion of peptide-DNA systems. demonstrated further below. Molecular dynamics simulation of the peptide-DNA Number of hydrogen bonds and salt bridges systems The number of hydrogen bonds between peptide and Simulations between CLFchimera, CLFcin and CLFam- DNA showed significant variation during simulation pin and BDNA were performed for 200 ns in triplicate. (Fig. 3; Additional file 5: Figure S5 and Additional file 6: To construct the system, the peptide was introduced in Figure S6). The average number of hydrogen bonds in the BDNA box at a random position and orientation. the second half of the three simulation replicates (100– Center of mass distance (COM), hydrogen bonds, salt 200 ns, 300 ns in total) was 5.66 ± 0.23, 4.61 ± 0.55 and bridges and contacting surface area between peptide and 2.63 ± 0.27 for CLFchimera, CLFcin and CLFampin, re- DNA were analyzed. spectively (see also Additional file 7: Table S1 for Fig. 3 Number of hydrogen bonds with DNA as a function of simulation time (200 ns). a CLFampin, b CLFcin, c CLFchimera. d Snapshot at t = 135 ns of the CLFchimera-DNA system, indicating hydrogen bonds (red lines) and salt bridges (yellow dashed) Pirkhezranian et al. BMC Genomics (2020) 21:60 Page 6 of 14 Fig. 4 Contacting surface area between peptide and DNA along a 200 ns MD simulation details), suggesting that CLFchimera establishes more and 4.76 ± 0.36 nm2 for the CLFchimera, CLFcin and stable interactions with DNA. CLFampin systems, respectively (see Additional file 7: A representative snapshot of the CLFchimera-DNA Table S1 for details). The CSA is higher for CLFchimera interaction is illustrated in Fig. 3d. In this frame, it can than for the other two peptides, in line with the observed be seen that hydrogen-bonding interactions are mainly interactions. established between positively charged residues of the peptide and the DNA-backbone phosphate groups, MM/PBSA binding free energy estimate which constitute also salt bridges. The binding free energy was estimated using the MM/ Salt bridges also play a fundamental role in protein- PBSA method. The results for the period 100–200 ns in ligand interactions [34, 35]. In several studies, a cutoff of one of the replicates are presented in Table 2. As indi- 4 Å between N-O atom pairs has been used to define salt cated in the Methods section, particularly for this type bridge formation [36, 37]. Here, we calculated salt bridges of systems (high charge density), the single-trajectory between P atoms from the nucleic-acid backbone and N MM/PBSA approach represents a very crude estimate of atoms from Lysine and Arginine residues, and thus used the binding free energy that, most certainly, severely 5 Å as cutoff. The average number of salt bridges in the overestimates the real value. Nevertheless, the calcula- second half of the three simulation replicates between tions will be used here to qualitatively compare and rank DNA and CLFchimera, CLFcin, CLFampin were 4.09 ± the different systems, which should be relatively safe 016, 3.17 ± 0.28 and 1.71 ± 0.44 (see Additional file 7: given that the nature of the interactions is the same in Table S1 for details). Again, CLFchimera establishes more all cases. The results indicate that CLFchimera has the salt bridges with DNA than the other two peptides. lowest DNA-binding energy. The plot of the binding free energy along the period 100–200 ns in one of the repli- Contacting surface area cates is shown in Fig. 5 (see Additional file 9: Figure S8 The solvent-accessible surface area was calculated with for the other two replicates). No significant differences the Gromacs library [38]. The contacting surface area can in the obtained binding free energy values were observed be then calculated using the following formula: CSA = ( among replicates (− 786 ± 2.545, − 731 ± 3.521 and − SASA Peptide(s) + SASA DNA – SASA Peptide(s)-DNA)/2 712 ± 7.801 kJ/mol for CLFchimera; − 340 ± 4.437, − [39]. Initially, the CSA was close to zero due to the dis- 352 ± 4.437 and − 316 ± 7.215 kJ/mol for CLFcin; − 71 ± tance between peptides and DNA. The evolution of the 3.063, − 78 ± 5.103 and − 62 ± 2.202 kJ/mol for CSA is shown in Fig. 4 for one of the simulation replicates CLFampin). (see Additional file 8: Figure S7 for the other two). In all The free energy values for the CLFchimera-DNA sys- three replicates, the CSA is stable after the initial 100 ns, tem were decomposed into residue contributions using indicating a stable interaction has been reached. The aver- the MmPbSaDecomp.py python script. The results, pre- age CSA in the period 100–200 ns is 5.92 ± 0.41, 4.9 ± 0.1, sented in Fig. 6 for one of the simulation replicates, Table 2 Binding free energy for the three peptide-DNA systems calculated by the MM/PBSA method (one simulation replicate) Peptides van der Waal (kJ/mol) Electrostatic (kJ/mol) Polar solvation (kJ/mol) Non-Polar solvation (kJ/mol) Binding energy (kJ/mol) CLFcin − 141 ± 1 − 1885 ± 2 1707 ± 7 −21 ± 0.1 −340 ± 4 CLFampin − 120 ± 1 − 825 ± 1 891 ± 3 −18.1 ± 0.1 −71 ± 3 CLFchimera − 152 ± 1 − 2396 ± 2 1781 ± 3 −20.75 ± 0.1 −786 ± 3 Pirkhezranian et al. BMC Genomics (2020) 21:60 Page 7 of 14 Fig. 5 Estimated binding free energy for the peptide-DNA systems. Calculated with the MM/PBSA method on the 100–200 ns period of one of the simulation replicates indicate that residues LYS5, LYS9, LYS13, ARG16, bovine Lactoferricin indicated that the core hexapeptide LYS18, ARG27, LYS34 and LYS35 are more relevant for “RRWQWR” in this peptide has a significant role in binding. On the other hand, GLU12 and SER36 have a antimicrobial activity [42]. The first two amino acids detrimental effect. The contributions in the other two from this central core in CLFchimera (ARG27 and simulation replicates follow the same trends (Add- ARG28) made a considerable contribution to the inter- itional file 10: Figure S9). action with DNA in our simulations (Fig. 6); however, Previous experimental studies revealed that substitu- they were not as effective as other positively charged res- tion of positively charged residues such as LYS269, idues. Investigation of minimum distances (averaged LYS277 and LYS282 with alanine in bovine Lactoferram- over the three replicates) showed that LYS5 and LYS35 pin (LYS9, LYS13 and LYS18 in CLFchimera) resulted in were closest to DNA, 0.13 ± 0.03 nm and 0.12 ± 0.02 nm, a dramatic decrease in antimicrobial activity [40, 41], a respectively (see Additional file 11: Figure S10). finding consistent with our in silico results (Fig. 6). As shown in Fig. 6, GLU12 and SER36 play a major However, Karn et al. (2006) showed that substitution of inhibiting role in the interaction with DNA. Additional GLU276 (GLU12 in CLFchimera) with glycine in bovine file 10: Figure S9 shows that they displayed also the lar- Lactoferrampin had no effect on increasing antimicrobial gest minimum distance to DNA, with 0.64 ± 0.13 nm activity [40]. Several experimental studies regarding and 0.57 ± 0.09 nm, respectively. Fig. 6 Contribution to DNA binding free energies of amino-acid residues in CLFchimera Pirkhezranian et al. BMC Genomics (2020) 21:60 Page 8 of 14 Fig. 7 Number of hydrogen bonds with DNA at different concentrations of CLFchimera. a 1-CLFchimera, b 2-CLFchimera, c 3-CLFchimera, d 4-CLFchimera Effect of CLFchimera concentration on DNA binding chains and DNA phosphate groups (see Additional file 14: Based on the previous sections, CLFchimera showed a Table S2). As these groups are equal in all nucleic-acid considerably higher affinity than CLFcin and CLFampin bases, we may conclude that the binding of CLFchimera for DNA. Therefore, we chose the former peptide as to DNA is DNA-sequence unspecific. For this reason, a candidate and evaluated the effect of its concentration potential DNA-related antimicrobial activity of this pep- on DNA binding by performing simulations with 1, 2, 3 tide could consist in the disruption of replication. Uyter- and 4 CLFchimera molecules and one DNA helix. hoeven et al., 2008 reported that Buforin II, another antimicrobial peptide, would probably target nucleic Number of hydrogen bonds and salt bridges acids in a non-sequence-specific manner [16]. Sim et al., As expected, when peptide concentrations rose, the 2017 also examined the interaction of Buforin II and number of hydrogen bonds between DNA and peptides increased but showed saturation (Fig. 7 for one of the replicates of each system; Additional file 12: Figure S11 and Additional file 13: Figure S12 for the other two rep- licates). The average number of hydrogen bonds in the 100–200 ns period of the three replicates was 5.66 ± 0.23, 9.66 ± 0.56, 11.75 ± 0.45, and 12.4 ± 0.37 for 1, 2, 3 and 4 peptides, respectively (see Additional file 7: Table S1 for details). The corresponding average number of salt bridges was 4.09 ± 0.16, 5.79 ± 0.72, 6.50 ± 0.24, and 7.03 ± 0.06, respectively (see Additional file 7: Table S1 for details). The largest change in the number of hydro- gen bonds and salt bridges occurs in the transition be- tween one and two peptides, and shows saturation after three peptides. The average percentage of hydrogen bonds involving Fig. 8 Contacting surface area at different concentrations the DNA phosphate group revealed that nearly all of CLFchimera hydrogen bonding was formed between peptide side Pirkhezranian et al. BMC Genomics (2020) 21:60 Page 9 of 14 Table 3 Binding free energy estimated by the MM/PBSA method, for one replicate at each CLFchimera concentration Peptides van der Waal (kJ/mol) Electrostatic (kJ/mol) Polar solvation (kJ/mol) Non-Polar solvation (kJ/mol) Binding energy (kJ/mol) 1-CLFchimera −152 ± 1 −2396 ± 2 1782 ± 3 −20.75 ± 0.1 −786 ± 3 2-CLFchimera − 176 ± 1 − 2398 ± 3 1650 ± 7 −21.6 ± 0.1 −946 ± 4 3-CLFchimera − 137 ± 2 − 2171 ± 6 1323 ± 10 −17.5 ± 0.2 −1004 ± 4 4-CLFchimera −137 ± 4 − 2001 ± 11 1137 ± 24 −36 ± 0.5 −1027 ± 16 DesHDAP1 with DNA and suggested that a large per- in CSA occurs for the transition between one and two centage of hydrogen bonds were formed between pep- peptides. However, the difference between 3 and 4 pep- tide side chains and phosphate groups in the nucleic tides is slightly larger than could be expected from the acid backbone (96.3 and 81.7%, respectively). Then, they previous analysis of interactions, indicating that the in- demonstrated experimentally that Buforin II and DesH- crease in contact area is not proportionally translated in DAP1 did not show signs of sequence-specific DNA additional specific interactions. binding [43]. In another study on a particular AMP de- rived from Chinese traditional edible housefly larvae by MM/PBSA binding free energy estimate Tang et al., 2009, it was revealed that the phosphate As shown in Table 3 (one of the replicates), the anion of the DNA double helix is one of the binding CLFchimera-DNA binding energy decreases with in- sites in DNA-peptide interaction [11]. The results of the creasing peptide concentration. The decrease in binding present study are in accordance with all these previous energy per peptide is again largest in the transition from findings. On the other hand, an experimental study on one to two peptides. The results for the different repli- the interaction between bovine Indolicidin and DNA cates are in this regard consistent: − 946 ± 4.504, − 927 ± demonstrated that Indolicidin bound tightly to the ds 9.31 and − 894 ± 5.056 kJ/mol for 2-CLFchimera; − [AT], ds [GC] and ds [AG] sequences, but formed loose 1004 ± 4.007, − 1020 ± 9.23 and − 975 ± 8.09 kJ/mol for bonds with ds [GT] [44]. 3-CLFchimera; − 1027 ± 15.718, − 1071 ± 12.21 and − 1022 ± 10.12 kJ/mol for 4-CLFchimera. The saturation Contacting surface area effect observed on going from three to four peptides Results from the CSA analysis at different concentra- could be partially due to repulsion between the highly tions of CLFchimera revealed an increase in CSA with charged peptides, in addition to crowding at the sites fa- peptide concentration (Fig. 8 for one of the replicates of vorable for the interaction with DNA phosphate groups. each system and Additional file 15; Figure S13 for the The binding energy contribution per amino-acid resi- other two replicates). The average CSA (three replicates) due (Fig. 9 for one replicate and Additional file 16: Fig- for 1, 2, 3 and 4 peptides was 5.92 ± 0.41, 10.30 ± 0.83, ure S14 for the other two) indicates that the residues 11.79 ± 1.14 and 14.55 ± 1.39 nm2, respectively (see Add- with important contributions are the same at the differ- itional file 7: Table S1 for details). The highest increase ent concentrations. It also shows that they are primarily Fig. 9 Contribution of amino-acid residues to DNA-binding energy at different concentrations of CLFchimera Pirkhezranian et al. BMC Genomics (2020) 21:60 Page 10 of 14 Fig. 10 Trends in the number of hydrogen bonds and salt bridges, and electrostatic, van der Waals and binding energy in simulations of DNA with CLFampin, CLFcin and different concentrations of CLFchimera Fig. 11 Top seven, most-populated conformations of CLFchimera in the simulations Pirkhezranian et al. BMC Genomics (2020) 21:60 Page 11 of 14 Fig. 12 Monitoring of BDNA conformation during MD simulation. a changes in BDNA conformation in the absence of the peptide. b change of BDNA conformation with presence of four copies of the peptide. c RMSD analysis of the DNA without and with four copies of the peptide responsible for the decrease in binding energy with in- Effect of CLFchimera concentration on DNA conformation creasing concentration. Analysis of RMSD and hydrogen bonds between the two Tables 2 and 3 demonstrate that complex formation strands of the DNA revealed that the presence of 1, 2 and stability was highly correlated with electrostatic and 3 molecules of CLFchimera had no significant effect interaction, as expected. Our MM/PBSA results conform on DNA conformation during MD simulations. How- to the results obtained by Pandey et al., 2018 on similar ever, four copies of the peptide caused DNA partial de- systems [18]. naturation (Fig. 12a and b), significantly increasing the A computational study by Khabiri et al., 2017 reported RMSD values (Fig. 12c) and decreasing the number of that changes in binding free energy do not correlate backbone hydrogen bonds (Fig. 13) after about 50 ns of strongly with salt bridges or hydrogen bonding in simulation. This result was observed for two (out of protein-DNA interactions [45]. Contrary to their results, three) replicates of four copies of concentrations. we show that for our peptide-DNA systems the electro- static energy and binding free energy are correlated with Conclusions the number of hydrogen bonds and salt bridges (Fig. 10), To improve our understanding of the possible intracellu- with an R2 between binding free energy and the number lar mechanisms of antimicrobial peptides (AMPs) of hydrogen bonds and salt bridges of − 0.92 and − 0.95, respectively, and an R2 of − 0.95 and − 0.97 between electrostatic energy and hydrogen bonds and salt brid- ges, respectively. Structural characterization of CLFchimera To assess the conformational heterogeneity of the pep- tide, clustering analysis was performed for all simulation trajectories. CLFchimera was thus distributed into 22 clusters. The conformation of the top seven clusters, covering 97% of all structures, is shown in Fig. 11. Clus- ters one and two comprised 65% of all structures. The analysis shows that the helical segments of CLFchimera Fig. 13 Analysis of the number of DNA interstrand hydrogen bonds are stable along the simulations, with major conform- in the presence of zero and four copies of CLFchimera ational changes affecting the turn between them. Pirkhezranian et al. BMC Genomics (2020) 21:60 Page 12 of 14 derived from camel lactoferrin, we planned an in silico Additional file 6: Figure S6. Third replicate: Number of hydrogen study based on molecular dynamics simulation. In this bonds with DNA as function of simulation time (200 ns). (A) CLFampin, study, DNA was selected as a well-known intracellular (B) CLFcin, (C) CLFchimera. target for AMPs. Overall, the simulation results indi- Additional file 7: Table S1. The value of hydrogen bond, salt bridge and contacting surface area for three replicates. cated that the chimeric peptide CLFchimera has a higher Additional file 8: Figure S7. Second and Third replicates: Contacting affinity for DNA than its component peptides CLFcin surface area between peptide and DNA along a 200 ns MD simulation. and CLFampin. It is shown that the interaction between Additional file 9: Figure S8. Second and Third replicates: Estimated DNA and CLFchimera follows a saturation curve with binding free energy for the peptide-DNA systems. Calculated with the increasing peptide concentration, the largest gain in MM/PBSA method on the 100–200 ns period of one of the simulation replicates. binding free energy corresponding, for a 12-bp DNA, to Additional file 10: Figure S9. Second and Third Replicates: the transition between one and two peptide copies. Sat- Contribution to DNA binding free energies of amino-acid residues in uration is reached after three peptide copies, with four CLFchimera. copies inducing the denaturation of the DNA helix. Additional file 11: Figure S10. Minimum distances during 200 ns Binding free energy analysis revealed that the number simulation. of lysine and arginine residues (12 in CLFchimera, 4 in Additional file 12: Figure S11. Second replicate: Number of hydrogen bonds with DNA at different concentrations of CLFchimera, (A) CLFampin, 7 in CLFcin) is largely responsible for the CLFchimera, (B) 2-CLFchimera, (C) 3-CLFchimera, (D) 4-CLFchimera. higher affinity of the CLFchimera-DNA complex. It is in Additional file 13: Figure S12. Third replicate: Number of hydrogen fact concluded that electrostatic interactions and par- bonds with DNA at different concentrations of CLFchimera, (A) ticularly salt bridges play a pivotal role in the interaction CLFchimera, (B) 2-CLFchimera, (C) 3-CLFchimera, (D) 4-CLFchimera. between CLFchimera and DNA. Moreover, nearly all Additional file 14: Table S2. Interaction results from MD simulations of CLF-chimera concentration with DNA. hydrogen bonds between CLFchimera and DNA in- Additional file 15: Figure S13. Second and Third replicates: Contacting volved the backbone phosphate groups of DNA, indicat- surface area at different concentrations of CLFchimera. ing that CLFchimera targets the nucleic acids in a Additional file 16: Figure S14. Second and Third replicates: sequence-unspecific manner. Contribution of amino-acid residues to DNA-binding energy at different A conformational clustering analysis suggests that concentrations of CLFchimera. CLFchimera mostly performs its activity through the interaction of its two α-helices with DNA. These sec- Abbreviations ondary structure elements, which are stable throughout AMP: Antimicrobial peptides; CLFampin: Camel lactoferrampin; CLFcin: Camel Lactoferricin; COM: Center of mass; CSA: Contacting surface area; the simulations, facilitate an optimal presentation of DA: Deoxyadenosine; DC: Deoxycytidine; DG: Deoxyguanosine; charge density and orientation for the interaction with DNA: Deoxyribonucleic acid; DT: Deoxythymidine; LFampin: Lactoferrampin; the phosphate groups. The outcomes of this study pro- LFcin: Lactoferricin; LINCS: Linear constraint solver; MM/PBSA: Poisson- Boltzmann and surface area continuum solvation; PB: Poisson-Boltzmann; vide insight into the structural and dynamic aspects of PDB: Protein data bank; PME: Particle mesh ewald; RMSD: Root mean square the interaction between Lactoferrin-derived AMPs and deviation; RMSF: Root-mean-square-fluctuations; SASA: Solvent accessible DNA as well as new directions for the design of novel surface area; SPC: Simple point charge AMPs with intracellular activity. Acknowledgments We are grateful to Faezeh Fatemizadeh and Oscar Conchillo Solé for their contribution to manuscript preparation. Supplementary information Authors’ contributions Supplementary information accompanies this paper at https://doi.org/10. ZP carried out the project, performed bioinformatics analyses and wrote the 1186/s12864-020-6458-7. manuscript. MT was the corresponding author and contributed to manuscript preparation. XD designed the experiments and revised the final draft. HM conceived the study and commented on the manuscript. MHS was Additional file 1: Figure S1. Second replicate: Structural fluctuation the corresponding author and contributed to manuscript preparation. All analysis. (A) RMSD as function of time; (B) RMSF per residue. authors have read and approved the manuscript. Additional file 2: Figure S2. Third replicate: Structural fluctuation analysis. (A) RMSD as function of time; (B) RMSF per residue. Funding Additional file 3: Figure S3. Second replicate: COM distance analysis. No funding was obtained for this study. (A) COM distances between CLFcin, CLFampin and CLFchimera and DNA along 200 ns. (B) Structures at times t = 0 (cyan) and t = 200 ns (purple): Availability of data and materials (B1) CLFcin-DNA, (B2) CLampin-DNA, and (B3) CLFchimera-DNA. Data sharing is not applicable to this article as the datasets that support the Additional file 4: Figure S4. Third replicate: COM distance analysis. (A) findings of this study are included within the article and its Additional files. COM distances between CLFcin, CLFampin and CLFchimera and DNA along 200 ns. (B) Structures at times t = 0 (cyan) and t = 200 ns (purple): Ethics approval and consent to participate (B1) CLFcin-DNA, (B2) CLampin-DNA, and (B3) CLFchimera-DNA. Not applicable. Additional file 5: Figure S5. second replicate: Number of hydrogen bonds with DNA as function of simulation time (200 ns). (A) CLFampin, (B) CLFcin, (C) CLFchimera. Consent for publication Not applicable. Pirkhezranian et al. BMC Genomics (2020) 21:60 Page 13 of 14 Competing interests 17. Bandyopadhyay S, Lee M, Sivaraman J, Chatterjee C. Model membrane The authors declare that they have no competing interests. interaction and DNA-binding of antimicrobial peptide Lasioglossin II derived from bee venom. Biochem Biophys Res Commun. 2013;430(1):1–6. Author details 18. Pandey B, Grover A, Sharma P. Molecular dynamics simulations revealed 1 Department of Animal Science, Faculty of Agriculture, Ferdowsi University of structural differences among WRKY domain-DNA interaction in barley Mashhad, Mashhad, Iran. 2Institut de Biotecnologia i de Biomedicina, (Hordeum vulgare). BMC Genomics. 2018;19(1):132. Universitat Autònoma de Barcelona, Bellaterra, Spain. 3Catalan Institution for 19. Eswar N, Webb B, Marti-Renom MA, Madhusudhan M, Eramian D, My S, Research and Advanced Studies (ICREA), Barcelona, Spain. 4Department of Pieper U, Sali A. Comparative protein structure modeling using Modeller. Chemistry, Faculty of Science, University of Neyshabur, Neyshabur, Iran. Curr Protoc Bioinformatics. 2006;15(1):5.6. 1–5.6. 30. 20. Maupetit J, Derreumaux P, Tuffery P. PEP-FOLD: an online resource for de Received: 8 June 2019 Accepted: 7 January 2020 novo peptide structure prediction. Nucleic Acids Res. 2009;37(suppl_2): W498–503. 21. Laskowski RA, MacArthur MW, Moss DS, Thornton JM. PROCHECK: a program to check the stereochemical quality of protein structures. J Appl References Crystallogr. 1993;26(2):283–91. 1. Wang S, Zeng X, Yang Q, Qiao S. Antimicrobial peptides as potential 22. Berendsen HJ, van der Spoel D, van Drunen R. GROMACS: a message- alternatives to antibiotics in food animal industry. Int J Mol Sci. 2016;17(5): passing parallel molecular dynamics implementation. Comput Phys 603. Commun. 1995;91(1–3):43–56. 2. Linde A, Ross C, Davis E, Dib L, Blecha F, Melgarejo T. Innate immunity and 23. Van Der Spoel D, Lindahl E, Hess B, Groenhof G, Mark AE, Berendsen HJ. host defense peptides in veterinary medicine. J Vet Intern Med. 2008;22(2): GROMACS: fast, flexible, and free. J Comput Chem. 2005;26(16):1701–18. 247–65. 24. Abraham MJ, Murtola T, Schulz R, Páll S, Smith JC, Hess B, Lindahl E. 3. Skalickova S, Heger Z, Krejcova L, Pekarik V, Bastl K, Janda J, Kostolansky F, GROMACS: high performance molecular simulations through multi-level Vareckova E, Zitka O, Adam V. Perspective of use of antiviral peptides parallelism from laptops to supercomputers. SoftwareX. 2015;1:19–25. against influenza virus. Viruses. 2015;7(10):5428–42. 25. MacKerell AD Jr, Banavali N, Foloppe N. Development and current status of 4. Pirkhezranian Z, Tanhaeian A, Mirzaii M, Sekhavati MH. Expression of the CHARMM force field for nucleic acids. Biopolymers. 2000;56(4):257–65. Enterocin-P in HEK platform: evaluation of its cytotoxic effects on cancer 26. Berendsen HJ, Postma JP, van Gunsteren WF, Hermans J. Interaction models cell lines and its potency to interact with cell-surface glycosaminoglycan by for water in relation to protein hydration. In: Intermolecular forces. molecular modeling. Int J Pept Res Ther. 2019;19(4):1–10. Dordrecht: Springer; 1981. p. 331–42. 5. Sang Y, Blecha F. Porcine host defense peptides: expanding repertoire and 27. Nosé S, Klein M. Constant pressure molecular dynamics for molecular functions. Dev Comp Immunol. 2009;33(3):334–43. systems. Mol Phys. 1983;50(5):1055–76. 6. Baloglu S, Boyle SM, Vemulapalli R, Sriranganathan N, Schurig GG, Toth TE. 28. Essmann U, Perera L, Berkowitz ML, Darden T, Lee H, Pedersen LG. A Immune responses of mice to vaccinia virus recombinants expressing either smooth particle mesh Ewald method. J Chem Phys. 1995;103(19):8577–93. listeria monocytogenes partial listeriolysin or Brucella abortus ribosomal L7/ 29. Evans DJ, BLJTJocp H. The nose–hoover thermostat. J Chem Phys. 1985; L12 protein. Vet Microbiol. 2005;109(1–2):11–7. 83(8):4069–74. 7. Tanhaiean A, Azghandi M, Razmyar J, Mohammadi E, Sekhavati MH. 30. Kumari R, Kumar R, Consortium OSDD, Lynn A. g_mmpbsa—A GROMACS Tool for Recombinant production of a chimeric antimicrobial peptide in E. coli and High-Throughput MM-PBSA Calculations. J Chem Inf Model. 2014;54(7):1951–62. assessment of its activity against some avian clinically isolated pathogens. 31. Cuervo A, Dans PD, Carrascosa JL, Orozco M, Gomila G, Fumagalli L. Direct Microb Pathog. 2018;122(5):73–8. measurement of the dielectric polarization properties of DNA. Proc Natl 8. Tanhaeian A, Ahmadi FS, Sekhavati MH, Mamarabadi M. Expression and Acad Sci. 2014;111(35):E3624–30. purification of the Main component contained in camel Milk and its 32. Kollman PA, Massova I, Reyes C, Kuhn B, Huo S, Chong L, Lee M, Lee T, antimicrobial activities against bacterial plant pathogens. Probiotics Duan Y, Wang W. Calculating structures and free energies of complex Antimicrob Proteins. 2018;10(4):1–7. molecules: combining molecular mechanics and continuum models. Acc 9. Daneshmand A, Kermanshahi H, Sekhavati MH, Javadmanesh A, Ahmadian Chem Res. 2000;33(12):889–97. M. Antimicrobial peptide, cLF36, affects performance and intestinal 33. Homeyer N, Gohlke H. Free energy calculations by the molecular mechanics morphology, microflora, junctional proteins, and immune cells in broilers Poisson− Boltzmann surface area method. Mol Inform. 2012;31(2):114–22. challenged with E. coli. Sci Rep. 2019;9(1):1–9. 34. Xie N-Z, Du Q-S, Li J-X, Huang R-B. Exploring strong interactions in proteins 10. Reyes-Cortes R, Acosta-Smith E, Mondragón-Flores R, Nazmi K, Bolscher JG, with quantum chemistry and examples of their applications in drug design. Canizalez-Roman A, Leon-Sicairos N. Antibacterial and cell penetrating PLoS One. 2015;10(9):e0137113. effects of LFcin17–30, LFampin265–284, and LF chimera on 35. Bosshard HR, Marti DN, Jelesarov I. Protein stabilization by salt bridges: enteroaggregative Escherichia coli. Biochem Cell Biol. 2016;95(1):76–81. concepts, experimental approaches and clarification of some 11. Tang Y-L, Shi Y-H, Zhao W, Hao G, Le G-W. Interaction of MDpep9, a novel misunderstandings. J Mol Recognit. 2004;17(1):1–16. antimicrobial peptide from Chinese traditional edible larvae of housefly, 36. Kumar S, Nussinov R. Relationship between ion pair geometries and with Escherichia coli genomic DNA. Food Chem. 2009;115(3):867–72. electrostatic strengths in proteins. Biophys J. 2002;83(3):1595–612. 12. Park CB, Kim HS, Kim SC. Mechanism of action of the antimicrobial peptide 37. Barlow DJ, Thornton J. Ion-pairs in proteins. J Mol Biol. 1983;168(4):867–85. buforin II: buforin II kills microorganisms by penetrating the cell membrane 38. Gromiha MM. Protein bioinformatics: from sequence to function: academic and inhibiting cellular functions. Biochem Biophys Res Commun. 1998; press; 2010. 244(1):253–7. 39. Zeng S, Zhou G, Guo J, Zhou F, Chen J. Molecular simulations of 13. Hartmann M, Berditsch M, Hawecker J, Ardakani MF, Gerthsen D, Ulrich AS. conformation change and aggregation of HIV-1 Vpr13-33 on graphene Damage of the bacterial cell envelope by antimicrobial peptides gramicidin oxide. Sci Rep. 2016;6:24906. S and PGLa as revealed by transmission and scanning electron microscopy. 40. van der Kraan MI, van der Made C, Nazmi K, van’t Hof W, Groenink J, Antimicrob Agents Chemother. 2010;54(8):3132–42. Veerman EC, Bolscher JG, Amerongen AVN: Effect of amino acid 14. Ho Y-H, Sung T-C, Chen C-S. Lactoferricin B inhibits the phosphorylation of substitutions on the candidacidal activity of LFampin 265–284. Peptides the two-component system response regulators BasR and CreB. Mol Cell 2005, 26(11):2093–2097. Proteomics. 2012;11(4):M111 014720. 41. van der Kraan MI, Nazmi K, Teeken A, Groenink J, van't Hof W, Veerman EC, 15. Yan J, Wang K, Dang W, Chen R, Xie J, Zhang B, Song J, Wang R. Two hits Bolscher JG, Amerongen AVN. Lactoferrampin, an antimicrobial peptide of are better than one: membrane-active and DNA binding-related double- bovine lactoferrin, exerts its candidacidal activity by a cluster of positively action mechanism of NK-18, a novel antimicrobial peptide derived from charged residues at the C-terminus in combination with a helix-facilitating mammalian NK-lysin. Antimicrob Agents Chemother. 2013;57(1):220–8. N-terminal part. Biol Chem. 2005;386(2):137–42. 16. Uyterhoeven ET, Butler CH, Ko D, Elmore DE. Investigating the nucleic acid 42. Chan DI, Prenner EJ, Vogel HJ. Tryptophan-and arginine-rich antimicrobial interactions and antimicrobial mechanism of buforin II. FEBS Lett. 2008; peptides: structures and mechanisms of action. Biochim Biophys Acta. 2006; 582(12):1715–8. 1758(9):1184–202. Pirkhezranian et al. BMC Genomics (2020) 21:60 Page 14 of 14 43. Sim S, Wang P, Beyer BN, Cutrona KJ, Radhakrishnan ML, Elmore DE. Investigating the nucleic acid interactions of histone-derived antimicrobial peptides. FEBS Lett. 2017;591(5):706–17. 44. Hsu C-H, Chen C, Jou M-L, Lee AY-L, Lin Y-C, Yu Y-P, Huang W-T, Wu S-H. Structural and DNA-binding studies on the bovine antimicrobial peptide, indolicidin: evidence for multiple conformations involved in binding to membranes and DNA. Nucleic Acids Res. 2005;33(13):4053–64. 45. Khabiri M, Freddolino PL. Deficiencies in molecular dynamics simulation- based prediction of protein–DNA binding free energy landscapes. J Phys Chem B. 2017;121(20):5151–61. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.