Potential Inhibitors of Lumpy Skin Disease’s Viral Protein (DNA Polymerase): A Combination of Bioinformatics Approaches

1. Introduction

The Lumpy skin disease (LSD) is a contagious viral disease that affects ruminant animals such as cattle, water buffaloes, and giraffes, posing a significant threat across borders [1,2,3,4]. This disease is attributed to the Lumpy skin disease virus (LSDV), a member of the Capripoxvirus genus within the Poxviridae family that exhibits close kinship with the sheep pox virus (SPPV) and goat pox viruses (GTPV) [5]. It has a genome spanning 151 kilobase pairs (kb) and 156 putative genes [6]. Similar to other poxviruses, LSDV possesses several conserved genes involved in its fundamental replicative mechanisms. Seven homologues of genes found in chordopoxviruses are associated with DNA replication. These genes, namely, LSDV039, LSDV077, LSDV082, LSDV083, LSDV112, LSDV133, and LSDV139, are either known or potentially involved in DNA replication processes [6]. However, the virus is most likely spread mechanically by blood-sucking arthropods such as flies, mosquitoes, and ticks. Direct contact between animals can also transmit the virus to a lesser extent [7,8]. Additionally, infected animals can spread the virus through milk, blood, nasal secretions, and saliva, which serve as alternative transmission routes through feeding or drinking [9]. Affected animals primarily exhibit symptoms such as fever, nodular skin lesions, dramatic decrease in milk production, and weight loss [10]. The introduction of LSDV into a herd can result in a high incidence rate, ranging from 5 to 45%, and a case death rate from 0.5% to 7.0% [11,12]. Consequently, LSD poses a significant economic threat to the livestock industry worldwide because it causes vast economic losses, including abortions in females and sterility in males [13].

It was first reported in Africa (Zambia) in 1929 and was identified as a communicable disease in the 1940s [14]. It has also spread to countries in Central and Eastern Africa, the Middle East, Asia, and Eastern Europe [15]. LSD outbreak data were accessed from the World Organization for Animal Health (WOAH) World Animal Health Information System (WAHIS) database in Southeast Asian countries during the study period between October 2020 and October 2021. During the epidemic period, 866 LSD outbreaks were reported in six Southeast Asian countries, including 1,758,923 susceptible cattle, 93,465 cases, 5,936 deaths, and 1,117 culled cattle [16]. In Bangladesh, the outbreak first emerged in the Chattogram region in July 2019 and then quickly spread throughout the entire country [17].

To date, the scientific community has preferred some homologous and heterologous vaccines to develop immunity owing to their cross-protection against lumpy skin diseases, such as Kenyan sheep and goat pox (KSGP) O-180 strain vaccines and Gorgan goat pox (GTP) vaccines. While the KSGPO-180 vaccine failed to protect cattle against LSDV, the Gorgan GTP vaccine successfully prevented the clinical symptoms of LSD in all vaccinated calves [18]. Certain SPPV vaccines provide only partial or incomplete protection against LSDV. Despite the use of goat pox vaccine virus KSGPO-240 and Romanian SPPV strains, their attenuation levels were insufficient to ensure the safety of cattle. Consequently, vaccinated animals remain ill [19]. The Neethling vaccine was developed in 2016 and tested in six Balkan countries. Its effectiveness was found to be average at 79.8% (95% CI: 73.2–84.7) in the Balkan region, with varying rates in different countries. Albania had a 62.5% effectiveness rate, whereas Bulgaria and Serbia had effectiveness rates of 97% [20].

Over the past three decades, many diseases have been treated using repurposed drugs. Zidovudine (AZT), the first successful repurposed drug used to treat Human Immunodeficiency Virus (HIV) in 1987, was originally developed to treat cancer [21,22]. Recently, a drug named Remdezivir (RDV), a repurposed drug approved by Food and Drug Administration (FDA) in 2020, was used to treat Covid-19 patients [23]. Additionally, a recent in vitro study showed that ivermectin, an antiparasitic drug, is effective against capripoxviruses [24]. This study found that ivermectin demonstrated significant inhibitory effects on viral replication and the attachment and penetration stages of the LSDV virus. Specifically, the results showed that ivermectin reduced viral replication by 99.82 and 99.87% at the replication stage, respectively. It also exhibited inhibitory effects of 68.38 and 25.01% at the attachment stage and 57.83 and 0.0% at the penetration stage [24]. Nevertheless, plant-based phytocompounds with antiviral activity are suitable natural solutions and show high efficacy in inhibiting several viral diseases, such as HIV and Covid-19 [25]. For the alternative treatment approach, propolis-alginate nanoparticles (Propolis-ALg NPs) have a potential therapeutic approach through different routes, including eye drops, oral routes, and topical spray. Transmission electron microscopy was used to characterize the propolis-ALg NPs. Propolis-ALg NPs effectively treat infected animals by reducing fever and boosting overall health [13]. Furthermore, colchicine, a potent natural alkaloid, has the potential to effectively treat a range of skin diseases, either as a standalone therapy or in combination with other drugs [26]. In this study, ivermectin was used as a control drug to serve as a benchmark for comparison with repurposed drugs and plant-based compounds investigated for their antiviral properties against LSDV. Although these have demonstrated some level of effectiveness, the efficacy and safety of these treatments for LSDV remain unclear. Considering these constraints, it is crucial to identify a specific drug that can effectively inhibit or reduce LSDV infection [27,28].

Currently, there is no cure for LSD, and available treatments mainly focus on relieving the symptoms of the disease [13]. Researchers are exploring new drug targets for LSD to develop more effective treatments. Computer-aided drug design is a valuable tool in this process because it allows for the design of new drugs based on the structure of protein targets within the virus [29]. Computational methods have been employed to analyze the protein targets of the LSD virus and design novel lead molecules or repurposed drugs for drug development. Target selection is a critical step in drug discovery, particularly for viruses with evasion mechanisms [29]. The identification of highly conserved and essential targets for viral survival and replication is crucial for successful drug development. Therefore, we targeted the DNA replication associated with DNA replication, including the DNA polymerase LSDV039, which is potentially involved in DNA replication.

Considering the morbidity and mortality, a truly effective antiviral drug against LSDV has emerged. Repurposed drugs [27,28,30], which can be rapidly used, are readily available, and show significance as antivirals, are nothing but a great choice to eradicate LSDV. Hence, this study aimed to identify LSDV inhibitors that could effectively inhibit viral replication. Virtual screening and molecular modeling techniques enable the identification of compounds (plant-based compounds and repurposed drugs) with high binding affinities for viral proteins. Further experimental validation will provide valuable insights into the effectiveness and safety of potential inhibitors and pave the way for the development of new therapeutic strategies against LSDV.

2. Materials and Methods

2.1. Target Selection and Validation

In this study, target proteins were screened after an extensive literature review (Supplementary Table S1). These proteins are essential for LSDV replication and survival of the LSDV [6]. Interestingly, most of these proteins share sequence similarities with sheeppox virus (SPPV) and goatpox (GTPV) [31]. The amino acid sequences of these proteins were obtained from UniProt KB (https://www.uniprot.org/), a widely used protein database, and their structures were predicted using AlphaFold2, an advanced protein structure prediction tool [32]. To validate the predicted structures, we used the PROCHECK [33] and ERRAT plot [34] programs available from the Structural Analysis and Verification server (SAVES) (http://nihserver.mbi.ucla.edu/SAVES), which is commonly used to assess the quality of protein structures [35]. The conserved domains and superfamilies of the selected proteins were analyzed using Chimera X software [36], InterProScan [37], and Pfam (Proteins Families Database) [38]. InterProScan provides structure-based categorization, domain and homologous superfamily predictions, and other high-level results. The Pfam database stores a curated set of sequence-aligned and curated data on protein families. Additionally, the active sites of the selected proteins were predicted using the SiteMap module of Schrödinger suite 2021-2 (Schrödinger, New York, United States), a suite of computational tools widely used in drug discovery and molecular modeling [39].

2.2. Ligand Selection

In the present study, we investigated the potential use of antiviral compounds as LSDV inhibitors. To create a dataset of active compounds, we conducted an extensive literature search of various databases and compiled a dataset comprising active phytocompounds and repurposed drug compounds. We created a library of active antiviral compounds by conducting a comprehensive search of related literature in various databases, including PubMed, Google Scholar, Web of Science, and Scopus [40]. We retrieved the 3D structures of active phytocompounds and drug compounds in SDF (Structure-Data File) format from the PubChem [41] and DrugBank [42] databases. The SDF format is a file format commonly used to represent chemical structures and associated data. Open Babel software was used to convert the SDF into the PDB format [43]. After obtaining the 3D structures of the phytocompounds, their molecular interactions, binding modes, and potential as LSDV inhibitors were investigated (Figure 1).

2.3. Ligand Preparation

LigPrep, a tool in the Schrödinger suite (Maestro 12.8) [44], was utilized for ligand preparation, which accounted for the metal-binding states, desalting, and generation of tautomers and stereoisomers. It also facilitated the generation of various possible ionization states at a targeted pH of 7.0 ± 2.0 using Epik. An OPLS4 force field was used during the preparation.

2.4. Protein Preparation

The first step in protein processing was to extract the desired protein structure from the AlphaFold2. Subsequently, the Protein Preparation Wizard (PrepWizard), part of the Schrödinger suite (Maestro 12.8) [45], was used to perform several preprocessing steps on the protein structure. These steps involve the addition of missing hydrogen atoms, side-chain optimization, correction of the ionization state of the protein, and assignment of bond orders and formal charges. In addition, advanced algorithms have been used to refine protein structures, such as the removal of water molecules and energy minimization using the OPLS4 force field.

2.5. Site map Analysis

The selected protein was deposited on a Site Map, a module of the Schrödinger suite (Maestro 12.8), for binding-site analysis [46]. This software employs advanced algorithms to identify binding sites, evaluate the location of binding sites with a high degree of confidence, and predict the druggability of these sites.

2.6. Receptor Grid Generation

The interaction region between the protein and ligand was determined by generating a receptor grid using the receptor grid generation tool in Maestro 12.8 [47,48]. This grid defined the specific area surrounding the active site of the protein, and its dimensions were established based on the x-, y-, and z-coordinates. The receptor grid box had a resolution centered at coordinates 4.2, 15.18, and 1.51 along the x-, y-, and z-axes, respectively.

2.7. Molecular Docking

2.7.1. Virtual Screening and Ligand Docking

During drug development, a computational approach known as virtual screening or virtual ligand screening is used to search small-molecule libraries and identify drug-like compounds that can bind to therapeutic targets [49]. It has been used to discover novel chemical entities in structure-based drug design. We assigned 718 compounds, including 381 phytocompounds and 337 drug compounds, to the virtual screening workflow (VSW) in Maestro 12.8 [50,51,52].

The high-throughput virtual screening (HTVS) mode of Glide is executed during the initial step. The top 50% of the resulting ligands were retained for further analysis in the subsequent stage, which involved Glide Standard Precision (SP). Once again, the top-scoring 50% of the ligands were selected to proceed with the Glide Extra Precision (XP) modes [53].

The screened compounds were again examined in the ligand-docking module in the extra precision (XP) mode, as it provides greater accuracy. The docking process used a flexible docking mode that automatically generated conformations for each ligand input. The G Score of the glide was examined as follows:

G Score = a × vdW+b × Coul + Lipo + H-bond + Metal + BuryP + RotB + site

Here, vdW is the Van der Waals energy, Coul is the Coulomb energy, Lipo is the lipophilic contact, HBond denotes hydrogen bonding, Metal means metal binding, BuryP indicates the penalty for buried polar groups, RotB means the penalty for freezing rotatable bonds, Site denotes polar interactions in the active site, and a=0.065 and b=0.130 are the coefficients of vdW and Coul, respectively.

Furthermore, additional settings were configured, including the incorporation of Epik state penalties into the docking score, execution of post-docking minimization, and computation of the root-mean-square deviation (RMSD) with respect to the input ligand geometries. The protocol provides the top-scoring ligands in the XP description based on the glide score and glide energy. Glide calculates the energies of a wide variety of interactions between ligands and proteins, such as hydrophobic interactions, hydrogen bonds, internal energies, pi-stacking interactions, salt bridges, desolvation, and RMSDs.

2.7.2. Free Energy Calculation by MM-GBSA

Following docking in Glide, both the receptor and ligands were considered Prime MM-GBSA. The Prime MM-GBSA module of the Schrodinger suite 2021-2 was employed to calculate the binding free energies of the protein-ligand complexes with high accuracy [54,55]. After selecting the receptor and ligand, Prime MM-GBSA uses the pose viewer file (pv.maegz) to show the appropriate descriptions [56]. The ligand-protein binding energy (G bind) was estimated using the following equation:

ΔG bind = G complex ₋ (G protein + G ligand)

In the given context, the terms G complex, G protein, and G ligand represent the minimum free energies associated with the protein-ligand complex, free protein, and free ligand, respectively.

2.7.3. ADMET Profiling of Novel Antiviral Compounds against LSDV

Ligands with the highest glide scores were carefully selected after virtual screening and ligand docking. To ensure that these selected ligands have the potential to be effective drug candidates, we performed a thorough evaluation using the Lipinski rule and ADMET analysis. We conducted this analysis using the QikProp module of the Schrödinger suite, specifically Maestro 12.8 [57].

2.7.4. Molecular Dynamics (MD) Simulation

MD simulations were performed using the Schrödinger LLC Desmond software for a duration of 100 nanoseconds [58]. Prior to MD simulations, a crucial docking step was performed to predict the static binding position of the ligand at the active site [59]. The MD simulations employed Newton's classical equation of motion to simulate atomic movements over time and predict the ligand-binding status in a physiologically relevant environment [60]. The ligand-receptor complex was prepared using Maestro's Protein Preparation Wizard, which involved the optimization, minimization, and filling of missing residues, if necessary. The system was created using a built-in tool.

The MD simulations were performed with the TIP3P (Intermolecular Interaction Potential 3 Points Transferable) solvent model in an orthorhombic box, which allowed for a 10 Å buffer region between protein atoms and box sides while maintaining a temperature of 300 K and pressure of 1 atm. Optimized Potentials for Liquid Simulation (OPLS 2005) force field were employed to describe the interactions within the system [61].

To mimic physiological conditions, counterions and 0.15 M sodium chloride were added to neutralize the total charge of the system. The models were loosened prior to the actual simulation. Trajectories were stored and inspected at regular intervals of 100 ps for further analysis. Furthermore, the trajectories obtained from the simulations were utilized in the Veusz software v. 3.6.2 (https://veusz.github.io/) to analyze structural insights and stability through root-mean-square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), and solvent-accessible surface area (SASA).

3. Results

3.1. LSDV DNA Polymerase Protein Has a Groove-Like Active Site

Identifying the key proteins associated with diseases is a fundamental step in understanding the underlying mechanisms and developing targeted therapeutic interventions [62]. In this study, we focused on the ability of the protein to replicate DNA during the LSDV life cycle.

Based on their DNA replication ability, we found 12 out of 156 putative genes that could replicate the DNA of the LSDV (Table S1) [6,63]. From these genes, we selected only LSDV039 genes that code DNA polymerase enzymes, and its principal function is to faithfully duplicate the genome, ensuring that genetic information is preserved and passed on from generation to generation [6,64]. Furthermore, DNA polymerase is important for the growth and survival of all living organisms, from tiny viruses and bacteria to more complex organisms such as humans [65,66]. Unregulated DNA replication is often linked to various medical conditions such as cancer, autoimmune diseases, and viral and viral/bacterial infections [67].

Analysis using Chimera X software revealed that the DNA polymerase belongs to the DNA-directed DNA polymerase family B. It has four functional domains: DNA polymerase B exonuclease N-terminal (residues 1-22), DNA-directed DNA polymerase family B exonuclease (residues 64-349), DNA-directed DNA polymerase family B viral insert (residues 350-481) and DNA-directed DNA polymerase family B multifunctional (residues 493-989) (Figure 2). Additionally, the findings from the InterProScan and Pfam analyses supported the previously mentioned information. According to InterProScan, the conserved site within the 2–515 residues and 520-1003 residues of amino acid residues indicates that it belongs to the Ribonuclease H-like superfamily and the DNA/RNA polymerase superfamily.

DNA polymerase protein structure prediction was carried out using AlphaFold 2, and analysis of its binding site through the site map revealed that the active site structure closely resembles a groove. Additionally, a site map with a high DScore of 1.048 was selected as the active site. Proteins with Dscores between 0.7 and 0.8 are moderately druggable, while those with Dscores greater than 1.0 are considered very druggable [68]. Active-site residues were used for receptor grid generation, followed by molecular docking (Table 1, Figure 3).

3.2. Structure Validation

To assess the quality of the 3D model, a Ramachandran plot was constructed using the PROCHECK software (Figure 4). Ramachandran plot analysis of the predicted model revealed that 91% of the residues were located in the most favorable regions, whereas 8.0% were in the allowed regions. This confirmed that the predicted model was of excellent quality. ERRAT is a measure of the overall quality of non-bonded atomic interactions, where a higher score denotes superior quality. A high-quality model is generally considered to have a range of >50. According to the ERRAT server, the current 3D model's overall quality factor was 88.777 (Figure 5). The overall structural validation results are presented in Table 2.

3.3. Binding Profile Analysis of Bonded Interactions

Four potential antiviral compounds were screened from 1098 compounds using a Virtual Screening Workflow (VSW), followed by ligand docking in Glide Extra Precision (XP) modes. Among them, 380 were phytocompounds and 718 were drugs. Virtual Screening Workflow (VSW) separated 160 conformers out of 1098, and further interaction analyses of the protein-ligand complex yielded the 10 best compounds from the result of ligand docking in Glide Extra Precision (XP) modes. Compounds with a higher Glide Gscore (Gscore is a specific scoring function used to quantify the strength of the binding interaction between the ligand and protein) and binding interactions were evaluated as the best conformers.

Ligand docking performed by Maestro (Schrodinger) showed that canagliflozin had the best Glide Gscore of −9.858 kcal/mol, with a binding affinity of −45.68 kcal/mol. Compound tepotinib had the Glide Gscore −8.856 kcal/mol and its binding energy was −47.99 kcal/mol. Rhein and taxifolin got the Glide Gscore −8.965 kcal/mol and −7.195 kcal/mol, respectively, showing binding affinity −44.72 kcal/mol and −44.48 kcal/mol (Table 3). The positive control was also depicted with a Glide Gscore of −5.630 kcal/mol, which has a binding affinity of −26.28 kcal/mol.

The interaction of amino acid residues with selected ligands showed that canagliflozin establishes two non-covalent bonds with the amino acid residues Lys 383 and Asn 659 (Table 4, Figure 7 a). Canagliflozin formed two hydrogen bonds with Asn 659 with a bond distance of 1.80 and 1.73 Å, respectively. In addition, a salt bridge was formed between Lys 383 and canagliflozin, with a bond distance of 2.61 Å. Tepotinib formed two non-covalent interactions involving hydrogen bonds with Glu 339 residues with a bond distance of 1.93 Å and a salt bridge with Glu 399 with a bond distance of 2.94 Å. Moreover, it also establishes four Pi interactions in which Tyr 477, Phe 499, and Tyr 663 amino acid residues are involved in pi-pi stacking with a bond distance of 5.49, 4.89, and 5.38 Å, respectively, whereas Pi cation is only formed with amino acid residue Lys 483 having a bond distance of 4.01 Å (Figure 7 b). Rhein formed two non-covalent interactions involving two hydrogen bonds with Val 498 and Thr 664 residues having a bonding distance of 2.12 and 2.09 Å, respectively (Figure 7 c). It also established a salt bridge with Lys 337 with a bond distance of 2.71 Å. Taxifolin is involved in forming five non-covalent hydrogen bonds with amino acid residues Ser 338, Glu 339, Ser 656, Lys 483, and Thr 664 having a bond distance of 2.10, 1.78, 2.02, 2.00, and 2.05 Å, respectively, and a Pi interaction specifically Pi-cation with Lys 483 having a bond distance of 3.90 Å (Figure 7 d). In case of positive control, it formed three non-covalent hydrogen bonding with Asp 337, Ser 338, and Glu 339 residues where the bond distance between them was 1.55, 2.16, and 1.83 Å, respectively (Figure 6).

3.4. Binding Free Energy Calculations for Non-Bonded Interactions

Non-bonded interactions are more abundant than bonded interactions during protein-ligand docking. One of the significant non-bonded interactions is the Van der Waals energy. Coulomb energy is also one of the major forms of non-bonded energy. In addition, lipophilic energy, which is an important non-bonded interaction, is a ligand–receptor complex. Binding free energy calculations by Prime MMGBSA module reveals that canagliflozin formed a Van der Waals interactions having energy −42.95 kcal/mol. Tepotinib formed Van der Waals energy having energy −58.20 kcal/mol. Rhein and taxifolin are involved in interaction with Van der Waals energy having −34.55 kcal/mol and −24 kcal/mol, respectively. Positive control formed that energy with −63.90 kcal/mol. Other forms of non-bonded interactions and their energies are listed in Table 5.

To determine the non-bonded interaction energy in the protein ligand complex, the following equation was employed: MMGBSA dG Bind (NS) = MMGBSA dG Bind - Rec Strain - Lig Strain.

3.5. ADMET Profiling of Novel Antiviral Compounds against LSDV

We identified the top ten potential antiviral compounds based on their high Glide Gscore, indicating their promising binding affinity to viral targets. These compounds were subjected thorough ADMET (Absorption, Distribution, Metabolism, Excretion, and Toxicity (ADMET) analyses using the QikProp module, which assesses crucial drug properties. The selected compounds exhibited favorable physicochemical and ADMET properties, making them excellent candidates for antiviral drug studies (Table 6).

3.6. Molecular Dynamics (MD) Simulation

MD simulations were performed to study the stability and behavior of the protein-ligand complexes over time. In this study, we examined the control complex and the four newly chosen compounds to understand their characteristics through a range of analyses, including RMSD, RMSF, Rg, and SASA. As shown in Figure 8a, the selected compound Rhein showed the lowest RMSD values compared to the positive control and the other selected compounds. A lower RMSD value indicates a higher level of system compactness [69]. The average RMSD values of ivermectin B1a, canagliflozin, tepotinib, rhein and taxifolin were 0.88 Å, 0.88 Å, 1.44 Å, 0.33 Å and 0.51 Å, respectively (Table 7). Canagliflozin exhibited a stable conformation from 15 to 63 ns. Tepotinib showed a stable conformation from 1 to 77 ns, after which fluctuations were observed from 78 to 100 ns. Rhein exhibited structural stability from 30 to 63 ns, whereas the remainder fluctuated. Taxifolin showed an initial fluctuation until 10 ns, after which a specific stabilization could be observed post-10 ns up to 89 ns; again, fluctuation occurred from 90 to 91 ns, and the rest showed structural stability. The control drug Ivermectin B1a showed stability and compactness at 16–36 ns and 44–100 ns, respectively.

Additionally, RMSF analysis was employed to examine the average atomic displacement from the mean positions in the residues during molecular dynamics simulation (Figure 8(b)). The canagliflozin showed the highest average RMSF value (1.69 Å), while the tepotinib showed the lowest RMSF value (1.26) among all systems. Higher RMSF values indicate that the system has greater flexibility during the molecular dynamics simulations [69]. The average RMSF value of Rhein (1.62 Å) was quite similar to the RMSF value of taxifolin (1.63 Å). However, Ivermectin B1a (control) showed average RMSF values of 1.58 Å (Table 7).

The radius of gyration (Rg) was used to measure the compactness of the system throughout the simulation (Figure 8(c)). Herein, the positive control ivermectin B1a showed the lowest Rg value of 32.72 Å among all the systems. In contrast, canagliflozin showed an average effect.

The Rg value (32.97 Å) that was quite similar to tepotinib (32.99 Å) and taxifolin (32.89 Å). However, rhein showed a slightly higher Rg value of 33.18 Å (Table 7). Figure 8d illustrates the SASA analysis for all systems. The SASA value provides insights into the compactness of the systems, with lower values indicating higher compactness, and higher values suggesting openness [69]. The selected complexes had average SASA values of 47186.72, 47677.70, 47724.08, 48127.33, and 47439.90 Å for ivermectin B1a, canagliflozin, tepotinib, rhein, and taxifolin, respectively, as indicated in Figure 8. Among all systems, the positive control, ivermectin B1a, displayed the lowest average SASA value of 9829.01 Å. The average SASA values of canagliflozin (47677.70 Å), tepotinib (47724.08 Å), and taxifolin (47439.90 Å) were quite similar, while rhein showed slightly higher SASA values of 48127.33 Å among all the systems.

4. Discussion

LSD is a viral disease that affects livestock and is considered a transboundary animal disease because of its ability to cross national borders, affecting multiple countries and regions. LSDV has significantly expansion of LSDV beyond its original endemic regions in Africa. LSDV have been reported in various geographical locations worldwide, including Asian, European, and Middle Eastern countries [70]. This disease was first reported in Bangladesh in July 2019 [71], China in August 2019 [72], India in November 2019 [73], Nepal in June 2020 [74], and other Asian countries [75]. In recent years, multiple cluster outbreaks of LSD have occurred on the Asian subcontinent, leading to substantial cattle mortality and posing a serious concern for the agricultural and livestock sectors [76]. For instance, in India alone, over 155,000 cattle deaths will be reported by 2022, underscoring the urgent need for disease control and management [77].

According to a report by the Department of Livestock Services (DLS) in 2019, 553,528 cattle were affected by LSD in all eight divisions of Bangladesh. The highest incidences were observed in Chattogram (8.26%) and Khulna (6.52%), whereas the lowest was observed in Sylhet (0.01%) (DLS, 2019) [78]. In one study, the LSDV collected from two Bangladeshi strains, BD-V392.1 and BD-V395.1, were distinct from contemporary field strains found in other Asian countries such as Hong Kong, China, Taiwan, and Vietnam. This suggests a unique genetic lineage or origin for Bangladeshi LSDV strains, which differs from strains circulating in neighboring Asian countries [79]. The Livestock Research Institute (LRI) in Bangladesh provides vaccines that are used as facultative measures, meaning that it is not mandatory but optional for farmers to vaccinate their livestock against LSDV. However, vaccination coverage is limited, resulting in many small holdings remaining unvaccinated [79]. Furthermore, the genetic variation of LSDV is one of the reasons why vaccines are not effective in a specific geographical region. Currently, no subunit vaccinations or chemotherapeutic medications are available for the treatment of LSD. It is currently feasible to prevent or treat this disease using other therapies such as phytochemical-based medicines or repurposed medications.

In this study, we used computer-aided methods to identify potential DNA polymerase inhibitors. The top four compounds, two drug compounds (canagliflozin and tepotinib) and two phytocompounds (rhein and taxifolin), were identified as highly promising LSDV DNA polymerase protein inhibitors. Molecular docking results reveal that canagliflozin poses the highest Glide Gscore (−9.858 kcal/mol) and the binding affinity (−45.68 kcal/mol) among the derived compounds. tepotinib, rhein, and taxifolin exhibited higher Glide Gscores and binding affinities than those of the control, Ivermectin B1a (−5.630 kcal/mol) (Table 3). These compounds form hydrogen bonds, pi-pi interactions, salt bridges, and cation-pi bonds with the LSDV DNA polymerase protein.

The formation of hydrogen bonds between the ligand and protein is highly selective and specific [80]. This depends on the spatial arrangement of the atoms in both the binding site of the protein and the ligand. The complementarity of these structures allows precise interactions to occur. The presence of hydrogen bonds stabilizes the ligand-protein complex, enhancing binding specificity [81]. Our results showed that each of the selected compounds formed hydrogen bonds with specific active-site residues of the LSDV DNA polymerase protein (Table 4), maintaining a specific distance. Taxifolin formed the highest number of hydrogen bonds. The positive control, ivermectin B1a, formed three hydrogen bonds. The hydrogen bond distance between Asn 659 and canagliflozin is 1.80 Å and 1.73 Å, respectively. The difference in the bond distance was due to the differences in the molecules within the same amino acid residue to which they bonded. The minimum bond distance between the ligand and specific amino acid residues indicates the maximum hydrogen bond energy.

Pi-Pi interactions are frequent in protein crystal structures, which help in the interaction of proteins and small molecules. The geometry of aromatic compounds and electrostatic interactions play a role in Pi-Pi interactions [82]. Besides, Pi-Pi stacking is also involved in the binding energy of receptor ligands [83]. From our findings, only Tepotinib formed three Pi-Pi interactions with the Tyr 477, Phe 499, and Tyr 663 amino acid residues, maintaining bond distances of 5.49, 4.89, and 5.38 Å, respectively.

Salt bridges are one of the strongest non-covalent interactions in nature. It plays crucial roles in protein folding, mediates protein-protein interactions, and facilitates molecular recognition [84]. Our results suggest that canagliflozin, tepotinib, and rhein, but not taxifolin, form salt bridges. Canagliflozin forms a salt bridge with Lys 483 with a bond distance of 2.61 Å. Tepotinib forms a salt bridge with Glu 339 with a bond distance of 1.93 Å. Moreover, rhein formed a salt bridge with Lys 337 by maintaining a bond distance of 2.71 Å. The Positive control ivermectin B1a did not exhibit such bonding interactions.

Cation-pi bond interactions play a crucial role in protein stability and structure. It involves attractive forces between the aromatic ring of a small molecule and a positively charged cation in a protein [85,86,87,88]. In terms of energy, this interaction is comparable to or sometimes stronger than a hydrogen bond [89]. Research suggests that cation-pi interactions are essential for protein-ligand recognition and have valuable implications for predicting drug-receptor interactions [90]. Our study reveals that tepotinib and taxifolin from our selected compounds only formed cation-pi bonded interactions with Lys 483 having bond distances of 4.01 and 3.90 Å, respectively. In contrast, positive control did not form these bonds.

A similar binding energy pattern was observed for the non-bonded interactions (Table 5). Coulomb energy is a form of electrostatic energy and a significant non-bonded energy found in protein-ligand complexes. Rhein has the highest coulomb energy (−77.35 kcal/mol), and ivermectin B1a possesses the lowest coulomb energy without any protein-ligand strains. Lipophilic energy is important for the evaluation of drug uptake and metabolism. Canagliflozin possesses the highest lipophilic energy (−25.21 kcal/mol) without strain. The Van der Waals interaction energies showed the opposite results. The Van der Waals interaction energy without receptor and ligand strains for the positive control was the highest (−63.90 kcal/mol) compared to the selected compounds.

The drug-like properties of the compounds satisfied Lipinski's Rule of Five, indicating their potential as drug candidates. Pharmacokinetics and toxicology offer important insights into the interactions between drug molecules in the human body [91]. The properties of these compounds were calculated using ADMET analysis tools such as QikProp. These tools predict the results of various chemical and physical properties of drug candidates, including molecular weight (MW), surface area, hydrogen bond (HB) interactions, lipophilicity, and oral absorption rates in humans [92]. Similarly, Qplog Po/w, QPlogS, Qplog HERG, and Human oral absorption were the major parameters used for the assessment of chemical compounds and for identifying the pharmacokinetic features of drugs (Table 4).

In addition, most drugs on the market have a molecular mass between 200 and 600 Da, with the majority being less than 500 Daltons [93]. To determine whether Canagliflozin, Tepotinib, rhein, and taxifolin are viable options, Lipinski's Rule of Five was considered. Lipinski's Rule of Five is a set of criteria used to evaluate a compound's drug-like behavior, including assessments such as the AMES test, Veber rule, and bioavailability [94]. Canagliflozin, tepotinib, rhein, and taxifolin satisfied most of the requirements, specifically Lipinski's Rule of Five, ensuring their drug-likeliness behavior (i.e., AMES test, Veber rule with no violations, better bioavailability value, and other parameters), which could be used as drug candidates (Table 1 and Supplementary Table S1).

Molecular simulations are a powerful approach for understanding the stability and dynamics of protein-ligand complexes [95,96]. Where the Mean Square Deviation indicates the average atomic displacement, with higher RMSD values indicating structural deviation or flexibility and a lower value indicating structural stability. The Square Fluctuation quantifies the average atomic positional variance from their mean positions within a molecule, with higher RMSF values suggesting greater atomic fluctuation or flexibility, and lower RMSF values indicating less atomic fluctuation and greater stability. The radius of gyration in molecular dynamics is a measure of the compactness or spread of a molecule's structure, where a higher Rg indicates a more extended or less compact structure, and a lower Rg signifies a more compact or folded structure. The solvent-accessible surface area quantifies the surface area of a molecule that is accessible to solvent molecules, where a higher SASA value suggests a greater exposed surface area, and a lower SASA value indicates a reduced exposed surface area. These parameters collectively provide essential insights into the molecular behavior in simulations [63,97,98].

The stability and dynamics of the compounds in protein-ligand complexes were studied. From the average mean value of the molecular dynamic simulation trajectory and analysis of the whole simulation data within 100 ns, it can be determined that the phytocompound rhein showed the highest stability in RMSD and RMSF. This fluctuated slightly in comparison with that of the positive control. The Rg and SASA values indicated less compactness and less exposure to the surface area compared to the positive control. It exhibited structural stability from 30 to 63 ns, whereas the rest of the simulation fluctuated. Taxifolin maintained its compactness in the RMSD. It showed stability at 10–89 ns; after that, it fluctuated from 1 to 10 ns and 91 to 100 ns, which is comparatively less than that of the positive control and gives a slightly higher RMSF value. The Rg and SASA values provided more compactness and a greater exposed surface area, respectively, than the other three compounds, but were marginally higher than those of the positive control. Tepotinib had the highest RMSD value among all compounds, indicating a lower level of compactness and even much higher compactness than the positive control, although this drug gave a minimum value in RMSF compared to the positive control and the rest of the compounds. It exhibited compactness in the range of 1–77 ns, and instability in the range of 78–100 ns. Canagliflozin showed the same RMSD score as the positive control. The RMSF of this compound was the highest among all the other compounds, including the positive control. Rg and SASA values followed the same trends. Canagliflozin exhibited compactness in 15 to 63 ns.

Canagliflozin and tepotinib are repurposed drugs, which means they have already been approved for use in treating specific diseases, and currently, there is interest in exploring their potential use for different diseases, including LSD. Canagliflozin (Invokana) is primarily used to the treatment of type 2 diabetes. It belongs to a class of drugs known as sodium-glucose co-transporter 2 (SGLT2) inhibitors. Canagliflozin is typically prescribed as a third-line treatment option after metformin (first-line medication) and other second-line alternatives in cases where these options are insufficient to effectively manage blood sugar levels [99,100,101]. Tepotinib belongs to a class of medications known as kinase inhibitors. It works by blocking the action of an abnormal protein that signals cancer cells to multiply [102,103,104].

Similarly, in this study, two specific phytochemicals, rhein and taxifolin, were identified as top-screened compounds with drug-like properties in inhibiting the macro domain of the chikungunya virus [105], with a focus on finding phytochemical inhibitors (drug-like properties of the top-screened phytochemicals rhein and taxifolin) against a unique macro domain found in the conserved N-terminal region of the non-structural protein nsP3 of the chikungunya virus. Medicinal plants have long been used as a source of potential therapeutic compounds and may offer a promising approach to combat LSD. This study investigated the binding affinities of two phytocompounds (rhein and taxifolin) and two repurposed drugs (canagliflozin and tepotinib) to the active sites of the DNA polymerase protein LSDV039, which are important targets for antiviral drug development. Post-MM/GBSA analysis, which calculated the binding free energy, indicated that the complexes exhibited a better binding free energy at 100 ns than the pre-MM/GBSA binding free energy. This indicates that the binding of phytocompounds to target proteins became more favorable over time. This suggests that they have the potential to interact effectively with these proteins, possibly by inhibiting their activity and interfering with viral replication.

However, to develop these selected compounds as LSD-specific antivirals, it is crucial to understand their behavior in their natural hosts, confirm their ability to inhibit replication, and provide antiviral prophylaxis. Notably, the identified compounds (canagliflozin, tepotinib, rhein, and taxifolin) are publicly available, which facilitates their rapid utilization and further research. Further experimental studies are necessary to validate their efficacy as novel compounds against LSD.

Overall, LSD continues to pose significant challenges to the livestock industry, and more comprehensive control measures, including improved vaccination coverage and alternative treatment options, are required to mitigate its impact. This study emphasized the importance of pharmacokinetic properties, efficacy, and safety levels in determining therapeutic efficacy and releasing new drugs. Based on in-silico screening of natural compounds, that is, rhein and taxifolin, and repurposed drugs, that is, canagliflozin and tepotinib, combinatorial docking, molecular dynamic simulation, drug-likeliness analysis, and all other experimental data were used. Taxifolin and canagliflozin have been suggested as drug candidates for the treatment of LSD.

5. Conclusions

This study provides valuable insights into the potential of specific phytocompounds and repurposed drugs against LSD. These findings suggest that canagliflozin, tepotinib, rhein and taxifolin have strong binding affinities for important viral proteins, exhibit favorable dynamics and hydrogen-bonding patterns, and possess favorable drug-like properties. This study serves as a foundation for future research and for the design of specific drugs targeting LSD. Experimental studies are required to confirm the efficacy of these phytocompounds as potential therapeutic agents against viruses.