Antifungal Activity and Computational Study of Constituents from Piper divaricatum Essential Oil against Fusarium Infection in Black Pepper

Fusarium disease causes considerable losses in the cultivation of Piper nigrum, the black pepper used in the culinary world. Brazil was the largest producer of black pepper, but in recent years has lost this hegemony, with a significant reduction in its production, due to the ravages produced by the Fusarium solani f. sp. piperis, the fungus which causes this disease. Scientific research seeks new alternatives for the control and the existence of other Piper species in the Brazilian Amazon, resistant to disease, are being considered in this context. The main constituents of the oil of Piper divaricatum are methyleugenol (75.0%) and eugenol (10.0%). The oil and these two main constituents were tested individually at concentrations of 0.25 to 2.5 mg/mL against F. solani f. sp. piperis, exhibiting strong antifungal index, from 18.0% to 100.0%. The 3D structure of the β-glucosidase from Fusarium solani f. sp. piperis, obtained by homology modeling, was used for molecular docking and molecular electrostatic potential calculations in order to determine the binding energy of the natural substrates glucose, methyleugenol and eugenol. The results showed that β-glucosidase (Asp45, Arg113, Lys146, Tyr193, Asp225, Trp226 and Leu99) residues play an important role in the interactions that occur between the protein-substrate and the engenol and methyleugenol inhibitors, justifying the antifungal action of these two phenylpropenes against Fusarium solani f. sp. piperis.


Introduction
Piper nigrum L. (Piperaceae), commonly known as black pepper, is one of the most famous and oldest spices in the world, with culinary and food preservative uses [1]. In the last century, the Japanese immigrants introduced the culture in Brazil, with excellent edaphic and climatic adaptation, driving the country to occupy the first place in world production. The Brazilian national harvest reached 44 thousand tons in 2011, with about 90% produced in the Pará State [2]. However, in last fifty years, the crop has been undermined by Fusarium solani f. sp. piperis, a biotrophic phytopathogen that causes the Fusarium infection with vessels obstruction and reduction of the culture, bringing serious social and economic losses for the regional producers [3].
Control measures for the Fusarium infection were initiated with research on the resistant black pepper specimens against the phytopathogen. However, the low genetic variability and the homogeneous plantations have also led to the spread of the fungus. Studies with Piper species from the Amazon, which were previously infected with F. solani f. sp. piperis, have identified plants with different response levels to the disease. Tolerant species have been used as rootstock, but, later, these plants exhibit incompatibility, and they die in the fourth year after cultivation [4]. Another alternative use for the native Piper species is the employment of their secondary metabolites as a natural fungicide. Piper essential oils and extracts have a high metabolic diversity, as the terpenes, phenylpropanoids, amides, lignans, neolignans, steroids, pyrones, piperolides, chalcones and flavones, which were previously reported as having antifungal properties [5,6].
The black pepper essential oil is associated with the numerous applications in the food and pharmacological industries and cosmetics and home remedies [7]. The aroma and flavor of black pepper essential oils are attributed to the presence of terpene compounds such as α-and β-pinene, sabinene, myrcene, limonene, β-caryophyllene and camphene [8]. Oil composition variability has been observed, according the major and popular cultivars. For instance, β-caryophyllene was identified in the leaf oil of all cultivars while α-and β-pinene, sabinene, myrcene, and limonene were not detected [9,10]. However, these terpene constituents are not defined as compounds having antifungal activity.
Previously, it was seen that the essential oil of Piper divaricatum occurring in the Marajó Island, Pará State, Brazil, is rich in eugenol and methyleugenol, two phenylpropenes with potent antifungal activity [11]. Seedlings of Piper divaricatum introduced in a nursery containing plants of black pepper infected with F. solani f. sp. piperis showed complete resistance to attack of the fungus, indefinitely. It was suspected that the presence of methyleugenol and eugenol in the oil of Piper divaricatum inhibits this fungus infection. Based on this finding, it was considered important to investigate if this other Piper species could be used in the future for controlling the Fusarium infection, either as a rootstock or in protection of the soil spreading the dried leaves in the culture areas of black pepper.
Many advances have been made in the knowledge of pathogenic plant fungi, at the cellular and molecular level. The β-glucosidase proteins from fungi are attractive targets for the design of new antifungal agents to control the Fusarium disease in plants [12,13]. Thus, the aim of this work was to provide information on possible roles of the main compounds of the P. divaricatum oil in the mechanisms of action against F. solani f. sp. piperis, evaluate the in vitro antifungal assays and perform the study of 3D structure of the β-glucosidase of F. solani (FsBglc). The 3D study was based on the structure of a homologous β-glucosidase protein by applying the method of comparative modeling by homology. Additionally, the molecular docking and calculation of the molecular electrostatic potential were applied to analyze the details of the protein-ligand interactions between the natural substrate (glucose/glc) and the two main phenylpropanoids of P. divaricatum, in the complex with FsBglc.

Antifungal Activity
The oil of Piper divaricatum showed high antifungal activity against F. solani f. sp. piperis, with inhibition above 90%. Methyleugenol and eugenol, the main constituents of P. divaricatum oil, were tested individually at concentrations between 0.25 and 2.5 mg/mL and they exhibited also strong antifungal activity against F. solani f. sp. piperis. The antifungal index for the oil of P. divaricatum ranged from 18.0% to 92.0%; the methyleugenol standard (95% pure) varied from 32.3% to 79.0%; and the eugenol standard (95% pure) inhibited totally the growth of the fungus, showing an antifungal index of 100%, at concentrations between 0.75 and 2.5 mg/mL. The MIC and IC50 values for the oil of P. divaricatum and the standards of methyleugenol and eugenol are shown in Table 2.  The antifungal activity of the oil of P. divaricatum has been previously reported. It showed germination inhibition for the fungal plant pathogens Marasmus perniciosa [18], Cladosporium cladosporioides and C. sphaerospermum [11], at the equivalent action of the miconazole, used as standard antifungal compound. Some studies have considered the inhibitory properties of methyleugenol against phytopathogens. It has strongly inhibited the growth of Alternaria humicola, Colletotrichum gloeosporioides, Rhizoctonia solani and Phytophthora cactorum, showing the best activity against this latter microorganism, after promoting its structural hyphae damage [19,20]. The vapor of eugenol has fungicidal activity against F. verticillioides, a phytopathogen that infects maize kernels [21]. Some fungi that promote the wood decay, as Lenzites betulina and Laetiporus sulphureus, have strong susceptibility to the eugenol, with antifungal index of 100% [22]. Nevertheless, eugenol showed weak fungicidal activity against Colletotrichum gloeosporioides, C. fragariae, C. acutatum and Botrytis cinerea, tested by direct bioautography [19,23].
Combined effects of methyleugenol (ME) and eugenol (E) were found to determine whether this combination has synergistic, additive or antagonistic effects against F. solani f. sp piperis, as previously noted by Benz (1971) [24]. This interaction was evaluated by comparing the effects of the isolated methyleugenol and eugenol at concentration of 0.75 mg/mL, and in combination at proportions of 4:1, 1:1 and 1:4. The results are shown in Figure 1. The antifungal index for methyleugenol was 76.3% and for eugenol was 100%. The antifungal index for the combination 4ME/1E decreased to 64.6%, indicating an antagonist effect. For the combinations 1ME/1E and 1ME/4E with values of 80% and 100%, respectively, an additive effect was observed, and it was highest than the isolated effect of methyleugenol. The significant antifungal activity of essential oil components arises from their hydrophobicity and could include attack on cell wall synthesis and retraction of cytoplasm in hyphae, which affects fungal growth and the morphogenesis, resulting in death of mycelium [24].

Homology Modeling and Others Structural Analysis
As explained in the previous section, the target was obtained using the Swiss-Model Workspace server [25] and the MODELLER-9v12 program [26]. The enzyme β-glucosidase from Kluyveromyces marxianus (KmBglI) (PDB entry 3AC0) [27], in complex with glucose (glc), was found as template for FsBglc, where the target primary sequence has 38% of identity with template. A sequential alignment between template and target is shown in Figure 2. The superimposition of 3D structures of the β-glucosidases is shown in Figure 3 and, as can be seen, the 3D structures of the target and template are highly similar (RMSD = 0.249 Å). The main difference is founded in a loop localized between the β14 and β5 β-sheets of the tertiary structures, where a residual gap can be found (see Figures 2 and 3). The results might be explained by the 3D model structure of FsBglc, only including the region corresponding to amino acids 1-340 in the amino-terminal, which is hypothesized for containing the catalytic domain.
In addition, the amino acids residues Asp45, Leu99, Arg113, Lys146, His147, Met190, Tyr193, Asp225 and Trp226, which are normally present in the catalytic region, were conserved in the target structure. Others important residues close to the catalytic region, such as Phe445, Phe508 and Glu590, were also conserved, but with different number positions, Phe436, Phe491 and Glu576, respectively.  The stereochemical quality of the homologue model and template were checked with the PROCHECK program [28] as implemented in the Swiss-Model Workspace [25]. The majority of the residues of the Swiss-Model and template were found occupying the most favored regions of the Ramachandran plot while the other residues occupied the additionally allowed regions. These plots represent the distribution of the Φ and Ψ angles in each residual of the amino acid. In homologue model, 89.4% of the residues were presented in the most favored regions, 8.0% in the additional regions, 1.3% in the generously allowed regions and 1.3% in the disallowed regions. For template, these values are 86.5%, 12.9%, 0.4% and 0.1%, respectively. Furthermore, the quality of the structure derived from homology modeling was validated by calculating the QMEAN Z-score, a combined scoring function consisting of three potential statistical terms and two additional terms, describing the agreement of predicted and observed secondary structure and the solvent accessibility, respectively [29]. The Z-score value for the homologue model was −0.839 while the Z-score value for the template was −0.150. These results suggest that the structural quality of homologue model is better than template, model in some aspects.
The electrostatic potential density maps (see Figure 4) might be indicators of electrophilic and nucleophilic centers, which govern the strength of bonds, the strength of non-bonded interactions, and molecular reactivity. They affect the strength of the interaction of the ligand with the receptor protein. It is possible to see that target and template models have different electrostatic potential densities in some regions outside of protein. However, the region corresponding to amino acids 1-340 in the amino-terminal has similar electrostatic potential density. Therefore, according the residual identity, where the principal residues were conserved in homologue model.

Molecular Docking and MEP Analysis
A wide method of evaluating a docking program is to gauge its performance in cognate self-docking and re-docking. In this process, a ligand is removed from the crystal structure with its target protein and the program is challenged to pose the ligand so closely as possible to the experimental coordinates. Here, this methodology was applied to validate the MVD software, which was used in the structures of the β-glucosidases homologue protein. Then, the template and homologue protein were overlapped, and the glucose (glc) was docked in 3D structure obtained by homology modeling. The results are shown in Figure 5A, where the molecular docking is completely superimposed in the glc experimental conformation. The MVD software can reproduce the top conformation of β-glucosidase ligands inside the binding pocket of the FsBglc catalytic site. Furthermore, previous computational studies applying molecular docking in β-glucosidase homologue models have been used successfully [30]. The docking study was employed to determine the binding energy of the natural substrates (glucose/glc), eugenol and methyleugenol. The scoring results of docking calculations using the MolDock Score function and the H-bond interactions observed for glc and other ligands are summarized in Table 3.
Eugenol and methyleugenol were docked against the FsBglc homologue protein (see Figure 5B,C). The same docking procedure used in crystal ligand was applied for these compounds. The search algorithm MolDock SE, in combination with the MolDock scoring function and the H-bond interaction, produced the best docking results. In Table 3, it is possible to see that glucose has the lowest MolDock score and the H-bond energy. It reflects the affinity of the FsBglc homologue by natural substrate of β-glucosidase enzymes. The H-bond interactions of eugenol and methyleugenol show a better agreement with the glucose results, when compared with the MolDock score. It can be explained by various possibilities of hydrogen interactions formed between each ligand and the catalytic site, where residues as Asp45, Arg113, Lys146, Tyr193 and Asp225 are induced in forming hydrogen interactions with the polar groups of ligands, according the experimental results ( Figure 4) [27]. Eugenol and methyleugenol can form 5 and 4 H-bond interactions, respectively ( Table 4). The H-bond interactions are in agreement with the IC50 values shown in the previous section of experimental results (see Table 4), where eugenol and methyleugenol have similar inhibitory activity values. Furthermore, an interesting new interaction found for these compounds was the hydrogen interaction with the residue His147. This interaction cannot be found in complex with glucose.  The nature of hydrogen bond interactions is still today the subject of many discussions and its importance in enzymatic environment has been shown in some drug design studies [31,32]. In this study, the hydrogen bond interactions were the main interactions found in protein-ligand complex as previously discussed. Thus, in order to understand why the eugenol and methyleugenol have showed similar interaction with the natural substrate, the molecular electrostatic potential (MEP) was built. MEP is a useful visual tool to understand the relative polarity of a molecule and serves to explain hydrogen bonding, reactivity, residual interaction, polarizability and the structure-activity relationship for biomolecules and drugs [32,33]. MEP was derived from DFT (B3LYP/6-31G*) calculations using Gaussian 09 package and visualized in Gaussview 3.07. These surfaces have an isodensity value of 0.002 a.u. In Figure 6, the most nucleophilic regions (negative electronic potential) are shown in red while the most electrophilic regions (positive electrostatic capability) have green color.   The natural substrate (glc) has an electronegative profile distributed around the whole molecule, this feature is important for the substrate be recognized by polar residues present in protein catalytic site, such as the residues Asp45, Arg113, Lys146, Tyr193 and Asp225. It can be observed the hydrogen bond interactions between OH groups of glc and the side chains of these amino acid residues (see Figure 6 and Table 3). Although eugenol and methyleugenol have fewer polar groups than glc, these compounds have shown a similar charge profile. This characteristic allows to seeing the same hydrogen bond interactions found in FsBglc-glc complex, where the phenyl group contributes to maintaining a similar electron density, as it was found in glc. A particular hydrophobic interaction can be found in all FsBglc-eugenol or FsBglc-methyleugenol complexes, where the aliphatic groups in these compounds interact with the side chain of residue Trp226 and the methoxy group (via methyl) interacts with the side chain of residue Leu99. It is in accordance with the hydrophobic properties found for these compounds in the experimental assays [24].

Plant Material
The leaves and thin stems (aerial parts) of Piper divaricatum G. Meyer were collected in the Marajó Island, municipality of Breves, Pará State, Brazil. The species was identified by Dr Elsie F. Guimarães, Jardim Botânico do Rio de Janeiro, Rio de Janeiro, Brazil, and a voucher (MG 165212) was deposited in the herbarium of Emílio Goeldi Museum, city of Belém, Pará State, Brazil.

Extraction and Analysis of Essential Oils
The aerial parts (100 g) were subjected to hydrodistillation using a Clevenger-type apparatus (3 h). The extracted oil was dried over anhydrous sodium sulfate, and its percentage content was calculated on the basis of the plant dry weight. The analysis of the oil was carried out on a THERMO DSQ II GC-MS instrument, under the following conditions: DB-5ms (30 m × 0.25 mm; 0.25 mm film thickness) fused-silica capillary column; programmed temperature: 60-240 °C (3 °C/min); injector temperature: 250 °C; carrier gas: helium, adjusted to a linear velocity of 32 cm/s (measured at 100 °C); injection type: splitless (2 mL of a 1:1000 hexane solution); split flow was adjusted to yield a 20:1 ratio; septum sweep was a constant 10 mL/min; EIMS: electron energy, 70 eV; temperature of ion source and connection parts: 200 °C. The quantitative data regarding the volatile constituents were obtained by peak-area normalization using a FOCUS GC/FID operated under conditions similar the GC-MS, except for the carrier gas, which was nitrogen. The retention index was calculated for all the volatiles constituents using an n-alkane (C8-C40, Sigma-Aldrich, St. Louis, MO, USA) homologous series. Individual components were identified by comparison of both mass spectrum and GC retention data with authentic compounds which were previously analyzed and stored in private library, as well as with the aid of commercial libraries containing retention indices and mass spectra of volatile compounds commonly found in essential oils [34]. The standards of eugenol and methyleugenol were purchased from Sigma-Aldrich, Brazil.

Fungal Strains and Culture Media
The phytopathogen Fusarium solani f. sp. piperis was obtained from an infected plantation of Piper nigrum in Pará State, Brazil. The fungus was cultured on potato dextrose agar (PDA) for ten days, at 27 °C, before the experiment to be carried out.

Effect of Essential oil on Mycelial Growth
The antifungal activity was performed by the agar dilution method [35] and expressed as percentage inhibition against the mycelia growth diameter. Aliquots of essential oil were incorporated in medium potato dextrose agar (PDA) sterilized, melting at concentration of 5 mg/mL. Mycelial discs (Ø, 0.9 cm) from Fusarium solani f. sp. piperis were transferred to the center of the plates and incubated in dark, at 27 °C. The mycelium growth was monitored during 7 days and the antifungal index was calculated according the Equation: Inhibition (%) = [1 − (∅growth of treatment)/(∅growth of control)] × 100 All treatments were tested in quadruplicate, and plates without the oil were used as control.

Determination of Minimum Inhibitory Concentration (MIC) and IC50
The oil of Piper divaricatum, eugenol and methyeugenol were incorporated in medium PDA at concentrations of 0.1, 0.25, 0.5, 0.75, 1.0 and 2.5 mg/mL. The minimum inhibitory concentration was defined as the lowest concentration able to cause 100% inhibition of mycelial growth. The IC50 values (concentration that inhibited 50% of the mycelium of fungi growth) were calculated by Probit analysis and nonlinear regression using the software GraphPad Prism 5.0.

Statistical Analysis
Antifungal experiments were performed in quadruplicate and data analyzed are mean ± SE subjected to one-way ANOVA. Means are separated by Tukey's multiple range tests when ANOVA was significant (p < 0.05) (GraphPad Prism 5.0).

Homology Modeling
Homology modeling is a well-known computational tool, which it is possible to build the secondary structure of a given protein based on knowledge of the primary structure [36]. This method is used to build protein structure, on the premise that specific functions are maintained during the evolutionary process between the conserved structurally homologous proteins [37]. As an initial step, the primary sequence of F. solani (NCBI code EEU48841.1) was obtained from the Genome Project, Nectria haematococca [38]. This sequence was subjected to BLAST search, in the Protein Data Bank (PDB), to identify suitable template structures for homology modeling. As homology modeling relies on sequence alignment, the target-template sequences were aligned using the Swiss-Model Workspace [25] and, after a careful examination of alignment errors, the automated comparative protein modeling program, MODELLER-9v12, was used to build the model [26]. The β-glucosidase from Fusarium solani (FsBglc) structure was analyzed and evaluated using PROCHECK program and QMEAN6 Z-score as implemented in Swiss-Model Workspace.

Molecular Docking
The 3D structure of FsBglc obtained by homology modeling was used as start point for molecular docking calculations. In order to validate the methodology used by Molegro Virtual Docker (MVD) software [39], a re-docking was performed using as reference the glucose (glc) coordinates, in complex with the crystal structure of Kluyveromyces marxianus (KmBglI) [27]. Then, the structures of eugenol and methyleugenol were tested in silico procedure. The MolDock SE (Simplex Evolution) search algorithm was used in the docking studies as implemented in MVD software. Furthermore, the MolDockScore, which is an adaptation of the Differential Evolution (DE) algorithm, and the H-Bond interactions were used as scoring function. The MolDock Score function was used because it yields a higher docking accuracy than other state-of-art docking algorithms (MVD: 87%, Glide: 82%, Surflex: 75%, FlexX: 58%) [39]. For the pose generation, 1500 maximum interactions were used to selecting a population size of 50. These poses can generate tests for different torsion angles, rotations and translations, which evaluate the affected part of the molecule and chooses the value that results in the lowest energy contributions. The MolDockScore function (Escore) is defined by the following energy terms: where Einter is the ligand-protein interaction energy and Eintra is the internal energy of the ligand. The Einter is determined by follow Equation: = + 332.0 4 The EPLP term is a piecewise linear potential using two different sets of parameters: one for approximating the steric term (van der Waals) between atoms and the other stronger potential for hydrogen bonds. The Eintra is calculated by follow Equation: The double summation calculates all the energy terms involving pairs of atoms of the ligand, except those connected by two bonds. The second summation calculates the torsional energy, where θ is the torsional angles of the bond. The average of the torsional energy term contributions is used if several torsions can be determined. The last term, Eclash, assigns a penalty of 1000, if the distance between two atoms (more than two bonds apart) is less than 2.0 Å. Thus, the Eclash term punishes infeasible ligand conformations [39].

Other Structural Analysis
Structural alignment of structures (template and target) and the electrostatic potential around the template and the modeled molecular systems were calculated by solving the nonlinear Poisson-Boltzmann equation, using finite difference method as implemented in Adaptive Poisson-Boltzmann Solver (APBS). The three-dimensional surfaces of the molecular electrostatic potential (MEP) were generated using the software Chimera and APBS [40], with an ionic strength of 12 mM. Mapping of the electrostatic potential onto the molecular surface of the protein was performed with a potential range from −10 eV to 10 eV. These calculations were carried using a PB2PQR Server [41]. After the docking calculation, the MEP was calculated using Gaussian-09 package, applying DFT (B3LYP/6-31G*) theory quantum level with CHELPG charges in vacuum, to verify the charge distribution profile for each molecule analyzed by molecular docking.

Conclusions
Methyleugenol and eugenol, the main constituents of the P. divaricatum oil, are responsible by the significant antifungal activity observed against Fusarium solanif. sp. piperis, which causes severe infection in black pepper. The interaction of FsBglc protein with the phenylpropanoid inhibitors is in agreement with the experimental data and, therefore, can provide valuable structural information for advances in rational drug design. The model also showed that the residues Leu99 and Trp226 can interact with eugenol and methyleugenol by hydrophobic interactions and, in this way justifying their antifungal action.