Integrated LC-MS/MS and network pharmacology approach for predictingactive ingredients and pharmacological mechanisms of Tribulus terrestris L. against cardiac diseases

Abstract Tribulus terrestris L. (Gokshura) is a medicinal herb used for treating cardiac diseases and several other diseases. However, the active ingredients and the possible mechanism of action for treating cardiac diseases remain unclear. Hence, the study was designed to identify the active ingredients and to explore the potential mechanism of action of Tribulus terrestris L. for treating cardiac diseases by an integrated approach of metabolomics and network pharmacology. We performed HPLC-QTOF-MS/MS analysis to identify putative compounds and network pharmacology approach for predictive key targets and pathways. Using molecular docking and molecular dynamics simulation, we identified the active ingredients in Tribulus terrestris L. that can act as putative lead compounds to treat cardiac diseases. A total of 55 putative compounds were identified using methanolic extract of Tribulus terrestris L. using HPLC-QTOF-MS/MS analysis. Network pharmacology analysis predicted 32 human protein targets from 25 secondary metabolites, which have shown direct interaction with cardiac diseases. Based on the degrees of interaction, the hub targets such as TACR1, F2, F2R, ADRA1B, CHRM5, ADRA1A, ADRA1D, HTR2B, and AVPR1A were identified. In silico molecular docking and simulation resulted in the identification of active ingredients such as Kaempferol 3-rutinoside 7-glucuronide, Keioside, rutin, moupinamide, aurantiamide, quercetin-3-o-α-rhamnoside, tribuloside, and 3ˈˈ,6ˈˈ- Di-O-p-coumaroyltrifolin against hub protein targets. Hence, these compounds could be potential lead compounds for treating cardiac diseases. A further assessment of its efficacy can be made based on in vivo and in vitro studies for better understanding and strong assertion. Communicated by Ramaswamy H. Sarma


Introduction
Herbal medicines are being used in ayurveda to treat several diseases since time immemorial.In recent years, there has been a rise in the search for new drugs from medicinal plants to combat the most common metabolic disorder such as cardiac diseases which leads to heart attack, myocardial infarction, hypertension, hypotension, cardiac failure, venous thrombosis, angina pectoris, coronary artery disease, cardiac arrhythmias, and many more (Reiner et al., 2019).Herbal medicine comprises several bioactive compounds that are being used to treat cardiac diseases with multiple target strategies instead of single target strategies.Several high throughput approaches such as liquid chromatography-mass spectrometry are being used for the identification of phytoconstituents in crude extracts and the reverse network pharmacology approach in understanding the mechanism of action in terms of identification of active ingredients against multiple targets related to cardiac diseases (Noor et al., 2022).
In ayurveda, Tribulus terrestris L. is a popular medicinal herb commonly known as gokshura or puncture vine.It belongs to Zygophyllaceae and the presence of several bioactive compounds attributes to various pharmacological activities such as antimicrobial, antitumor, antihypertension, and anti-cardiovascular diseases (Chhatre et al., 2014).Many studies have shown that Tribulus terrestris L. saponins are the main active compounds present in fruits that showed protective activity against myocardial infarction, ischemic disease, and coronary disease (Guo et al., 2007;Mohd et al., 2012).However, the identification of active ingredients from Tribulus terrestris L. other than steroidal saponins and the mechanism of action of Tribulus terrestris L. against cardiac diseases remain elusive.
The study aims to identify putative compounds using the HPLC-QTOF-MS/MS technique and network pharmacology approach coupled with in silico molecular docking and molecular dynamics (MD) simulation to identify active ingredients as potential leads, hub targets, and key pathways against cardiac diseases to understand the molecular mechanisms of action of Tribulus terrestris L. with multi-target strategy.

Sample collection and identification
According to the ayurvedic pharmacopeia of India, Tribulus terrestris L. (whole plant) was collected from Southern India and assigned a unique accession number.The herbarium was prepared and taxonomically identified to preserve the specimen for future reference.Germplasm of plant species was maintained at the greenhouse at Manipal School of Life Sciences, MAHE, Manipal, India.

DNA barcoding
The genomic DNA was isolated for Tribulus terrestris L. using a modified CTAB protocol (Tiwari et al., 2012).The isolated DNA was visualized using 0.8% agarose gel and molecular authentication was performed with a universal nuclear internal transcribed spacer region marker (nrITS) followed by Sanger DNA sequencing with Genetic analyzer 3500.

Metabolite fingerprinting
The whole plant of Tribulus terrestris L. sample was dried and ground with liquid nitrogen with a mortar and pestle.The ice-cold methanolic extraction protocol was followed to extract the metabolites present in Tribulus terrestris L. (De Vos et al., 2007).The extract was dissolved with 100% Methanol and stored at À 80 � C before use.Metabolite fingerprinting was performed using Agilent LC-QTOF-6530 with HPLC-QTOF-MS system (Agilent Technologies, USA).The separation was carried out on the Xbridge C18 column (4.6 � 250 mm, 1.8 mm) with 0.1% formic acid in deionized water and 0.1% formic acid in acetonitrile as mobile phases.The LC-MS sample injection volume was 5 ml.The gradient elution program was performed for 60 min run (5-35%B; 0-45 min, 35-75%B; 45-47 min, 75%B; 47-52 min, 75-5%B; 52-54 min, 5%B; 54-60 min) with a flow rate of 0.5 ml/min.The LC-MS analysis was performed in duplicates to obtain the accuracy of the chromatogram peaks and MS/MS fragmentation analysis was operated in the positive ionization mode with the mass range from m/z 50-1000.The mass spectrometry conditions such as spray voltage at 3 kV, nitrogen gas temperature at 350 � C, nebulizer pressure at 35 psi, drying gas flow rate at 8 L/min, and collision energy used in MS/MS experiment was ramped from 10-60 eV at the interval of 5 eV.

Data analysis
Agilent MassHunter Qualitative Analysis (version B.06.01) for processing of LC-MS data and molecular feature extraction algorithm was used to find the metabolite features.The raw data of HPLC-QTOF-MS/MS was exported to MS-DIAL software 3.50 (Tsugawa et al., 2015) for peak analysis with a tolerance of 0.01 Da and 0.40 Da for MS1 and MS2 respectively in centroid mode.All the other parameters from MS-DIAL were kept at default for data processing and exported the ions obtained to MS-Finder version 3.50 (Tsugawa et al., 2016) to calculate the molecular formula and structural compound identification.The putative compounds were identified with the top score assigned and further screened using the dictionary of natural products (http://dnp.chemnetbase.com) and a literature search to identify putative plant compounds.The 2D structures for each putative compounds were built using 2D sketcher tool (Schrodinger 2021-2, LLC, New York) and PubChem compound IDs were retrieved using the PubChem database (Kim et al., 2016).

Target prediction
The secondary metabolites of Tribulus terrestris L. were screened by the molecular structure's ontology and the literature search.The 3D structures were retrieved from the PubChem database (Kim et al., 2016) and searched against the human targets based on a similarity score (>0.7) in the BindingDB database (Gilson et al., 2016).The gene symbol for the human protein targets was retrieved from the UniProtKB/Swiss-Prot database (Bateman et al., 2021).

Compound-target-disease network construction
The protein targets related to cardiac diseases were analyzed from the Therapeutic target database (TTD) (Wang et al., 2020).In TTD, the successful disease targets which are enrolled in clinical trials above phase 1 were included.The tripartite network was constructed by linking the secondary metabolites of Tribulus terrestris L., successful targets of cardiac diseases and several types of cardiac diseases.The network visualization was performed using Cytoscape v3.7.2 (Shannon et al., 2003) to predict the pharmacological actions of compounds and their targets related to cardiac diseases.

Protein-protein interaction network and module analysis
Protein-protein interaction (PPI) network was constructed for human cardiac disease targets with StringDB (version 10.5) database by setting the confidence score greater than 0.9 (Szklarczyk et al., 2019).The network was retrieved and imported into the MCODE Cytoscape plugin version 1.6.1 (Bader & Hogue, 2003) to identify the highly interconnected protein network for module analysis.Finally, the core target proteins were selected based on the degree level within the PPI and module analysis network using Cytoscape plugin cyto-Hubba (version 0.1) in Cytoscape.

Functional enrichment
The KEGG pathway and Gene Ontology (GO) enrichment analysis for predicted targets from module analysis and the overall target proteins related to cardiac diseases was performed using the following packages clusterProfiler (Yu et al., 2012), org.Hs.eg.db, andDOSE (Yu et al., 2015) in R version 3.6.1 (Team, 2013).

Ligand and protein preparation
The 3D structure minimization for 25 putative compounds that showed interaction with cardiac diseases was performed using the LigPrep tool (Schrodinger 2021-2, LLC, New York).The Protein Data Bank (Burley et al., 2021) and Alpha Fold databases (Jumper et al., 2021) were used to retrieve the crystal structures of hub targets.Maestro 13.2's protein preparation wizard was used to minimize the crystal structures of hub targets.

QikProp ADME analysis
The ADMET properties of minimized ligand structures were predicted using the QikProp module in the Schrodinger suite (Schrodinger 2021-2, LLC, New York).The ADMET properties such as the dipole moment, rule of five, molecular weight, octanol-water partition coefficient, solubility, blockage of HERG K þ channels, caco2 cell permeability, blood-brain partition coefficients, and human oral absorption were calculated.

Molecular docking and post docking energy minimization studies
The SiteMap module was used to identify the probable binding sites of ligands to hub targets and the Glide-Receptor Grid generation module was used to generate the grid.The glide module of Schrodinger suite 2021-2 in extra precision mode was used to perform molecular docking with the minimized ligands and the selected active site of the protein.Molecular mechanics with generalized Born and surface area solvation (MM-GBSA) of Schrodinger 2021-2 was used to calculate the free energy of the ligand in a complex with a protein receptor.,The post docking energy minimization calculations were performed for the minimized XP docked poses of ligand-receptor complexes of hub targets using variable surface generalized born as a solvent model and OPLS3e as the force field.

MD simulation
MD was performed using the Desmond module of the Schrodinger suite 2021-2.In MD, an orthorhombic box was generated with SPC solvent model system.The salt concentration was kept constant at 0.15 M, and the model was optimized.The constant temperature and pressure (NPT) were maintained and simulation was performed for 100 ns.A total of 1000 frames were generated during the MD simulation.These simulation results were analyzed by calculating the root mean square deviation, root means square fluctuation, and ligand interaction profile obtained using the simulation interaction diagram.

Sample collection and molecular authentication
The fresh sample of Tribulus terrestris L. was collected from Chitradurga, Karnataka, India (Latitude: 14.2251 � N, Longitude: 76.3980 � E).The herbarium was stored at the Department of Biotechnology, Manipal School of Life Science, MAHE, Manipal, India.The germplasm was maintained in the greenhouse at Manipal School of Life Sciences, MAHE, Manipal, India.The DNA barcode sequence of Tribulus terrestris L. with nrITS marker was deposited to the GenBank database under the accession number MW591544 for future reference.

Putative compound identification using HPLC-QTOF-MS and MS/MS analysis
The total ion chromatogram (TIC) was obtained from HPLC-QTOF-MS/MS for Tribulus terrestris L. extract run in a positive mode of ionization (Supplementary Figure 1).HPLC-QTOF-MS and MS/MS fragmentation data resulted in the identification of 55 putative metabolites.The mass spectrometry data of the identified putative compounds including retention time, experimental m/z, error in parts per million (<15 ppm), molecular formula, score obtained (max ¼ 10), and the ontology of the compounds were shown in Supplementary Table 1.The 2D structures for each putative compound and PubChem CID were shown in Supplementary Table 2.The experimental mass spectra matching the in-silico mass spectra for each putative compound obtained using MS-DIAL integrated MS-Finder software were shown (Supplementary Figure 2).

Prediction of cardiac disease-related targets
Based on the compound class, 14 primary and 41 secondary metabolites were identified.A total of 210 human protein targets were predicted from 34 secondary metabolites.Seven metabolites that do not show any interaction with human targets were removed from the network construction.Out of 210 targets, 32 targets related to cardiac diseases were screened using the Therapeutic Target Database (TTD).The compound-target-disease network related to cardiac diseases was constructed using Cytoscape 3.7.2(Figure 1).A total of 25 secondary metabolites showed interaction with 32 potential targets of cardiac diseases for protein-protein interaction network analysis.

Identification of hub targets
The protein-protein interaction (PPI) network consists of 32 nodes and 114 undirected edges with a network density of 0.115 and an average number of targets of 3.56.The module analysis identified the two clustering modules: Module1 resulted in an MCODE score of 8.0 containing 9 nodes with 36 edges and Module 2 resulted in an MCODE score of 4.0 containing 5 nodes and 10 edges.Module 1 showed the higher average connectivity, and it was further analyzed using Cytohubba plug-ins to identify the hub targets.The top nine nodes (TACR1, F2, ADRA1B, F2R, CHRM5, AVPR1A, ADRA1A, ADRA1D, and HTR2B) with degree values greater than 8 were selected as hub targets from module 1.These nodes change color from red to yellow with decreasing degree values (Figure 2).

Identification of key pathways and gene ontology for hub targets and overall targets against cardiac diseases
Gene ontology enrichment analysis for hub targets resulted in 39 biological processes (BP), 5 cellular components (CC), and 2 molecular functions (MF) for a p-value threshold of 0.01.Gene ontology related to cardiac diseases was screened and obtained a total of 9 BP, 4 CC, and 2 MF.In BP, targets are involved in calcium ion homeostasis, vasoconstriction, activation of phospholipase C activity, calcium ion transport into the cytosol, cytosolic calcium ion transport, G proteincoupled serotonin receptor, excretion, negative regulation of  autophagy, and cell growth (Figure 3).In CC, targets are involved in caveola, membrane raft, membrane microdomain, and postsynaptic membrane.In MF, the targets are mainly involved in G protein-coupled amine receptor activity and Gprotein alpha-subunit binding.Similarly, gene ontology enrichment for 32 target proteins resulted in 132 biological processes (BP), 20 cellular components (CC), and 10 molecular functions (MF).After ranking based on gene count, a total of 35 relevant biological functions were screened and includes mainly calcium ion homeostasis (GO: 0055074), Gprotein coupled receptor signaling pathway (GO: 0007187), activation of adenylate cyclase activity (GO: 0007190), and many more (Supplementary Table 3).

ADME analysis of ligands
The assessment of 25 minimized ligand structures was performed using various descriptor calculations such as dipole moment, rule of five, QPlog Po/w, QPPCaco, QPlogHERG, QPlogS, PSA and, human oral absorption (Supplementary Table 5).All the ligands satisfied the dipole moment, the logarithm of the partition coefficient of the compound between n-octanol and water (QPlog Po/w), and predicted aqueous solubility (QPlogS).Out of 25 ligands, 16 ligands were within the acceptable range for % of oral absorption and 17 ligands were predicted Caco-2 cell permeability within the acceptable range.The QlogHERG values for 14 ligands were less than À 5 and hence, these ligands can block the hERG potassium channel.

Molecular docking analysis
Based on the least resolution, the 3D structures of hub targets were selected.(Mahmood et al., 2005) were retrieved from RCSB PDB (Burley et al., 2021) and ADRA1A (AlphaFold ID: AF-P35348-F1), ADRA1B (AlphaFold ID: AF-P35368-F1), ADRA1D (AlphaFold ID: AF-P25100-F1) were retrieved from AlphaFold database (Jumper et al., 2021).Using site map analysis, the top five docking sites for the hub target protein crystal structures were determined.The binding site with top site score and D-score were further considered for extra precision molecular docking analysis (Supplementary Table 6).Extra precision molecular docking for cocrystallized ligands with hub targets were shown in Supplementary Table 7.However, the cocrystallized ligands bound with hub targets (F2 and HTR2B) didn't show any poses in extra precison mode.Extra precision molecular docking showed a good binding affinity for all the compounds with hub targets except for ADRA1D protein (Table 1).In which, several flavonoids such as rutin, Keioside, Kaempferol 3-rutinoside 7-glucuronide, 3'',6''-Di-O-p-coumaroyltrifolin, tribuloside, and quercetin 3-rhamnoside showed higher glide score against hub targets.The 3D surface view of the top hit ligands against the hub target was represented (Figure 4).
The free ligand binding energy was analyzed using MM-GBSA for 25 ligands against hub targets except for ADRA1D protein.The free energy for the top three hit compounds from each hub target and molecular interactions were represented (Table 2).Based on glide score, free energy, and   molecular interactions, the top three hit compounds from each hub target were subjected to molecular dynamics simulation to identify the stable interactions.

MD simulations
MD simulation helps to monitor the stability and compatibility of protein-ligand complexes.The interaction of TACR1 with Kaempferol 3-rutinoside 7-glucuronide (Complex 1 A) was stable throughout the simulation of 100 ns (Figure 5).
The MD simulation was performed for the top hit compounds against CHRM5, F2R, and F2 protein respectively.The RMSD plot, RMSF plot, and protein-ligand interactions for different complexes against CHRM5, F2R, and F2 protein were represented (Figure 5; Supplementary Figures 3A and 4A).After comparing the RMSD and RMSF plots, rutin has the potential to inhibit F2R and CHRM5 proteins and tribuloside can inhibit F2 protein respectively.
The MD simulation was performed for the top hit compounds against ADRA1B protein (Complex 5 A, 5B, and 5 C).The RMSD plot, RMSF plot, and protein-ligand interactions for different complexes against ADRA1B protein were represented (Figure 6; Supplementary Figures 3B and 4B).After comparing the RMSD and RMSF plots, rutin has the potential to inhibit ADRA1B protein more effectively.Similarly, MD simulation was performed for other proteins (HTR2B, ADRA1A, and AVPR1A).Based on RMSD and RMSF plots along with protein-ligand interactions, we identified that Keioside can inhibit ADRA1A protein effectively whereas rutin has more potential to inhibit HTR2B and AVPR1A proteins respectively (Figure 6; Supplementary Figures 3B and 4B).

Discussion
Tribulus terrestris (L.) has multi-therapeutic potential and different parts of the plant are being used for treating several diseases including cardiac diseases.Several bioactive compounds are present in T. terrestris L extract, and the identification of these compounds from the plant extract can help us in developing new drugs and pharmacological agents (Katiyar et al., 2012).Hence, in the present study, the methanolic extraction of T. terrestris was analyzed using LC-MS and MS/MS analysis in a positive mode of ionization, resulting in an identification of putative compounds by matching the experimental spectra with in silico fragmentation databases such as MassBank, HMDB, FooDB, UNPD, KNaPSacK, NANPDB present in MS-DIAL integrated with MS-Finder software.
To identify new drug candidates for treating cardiac diseases, a network pharmacology approach was utilized to study the interactions between the active compounds and the disease targets.We identified the key targets include tachykinin receptor (TACR1), adrenergic receptors (ADRA1A, ADRA1B, ADRA1D), proteinase-activated receptors (F2, F2R), acetylcholine receptor (CHRM5), 5-hydroxytryptamine receptor (HTR2B) and Vasopressin and oxytocin receptors (AVPR1A).Adrenergic receptors (ADRA1B) regulate the systemic arterial blood pressure, and positive regulation of heart rate (Kim et al., 2021).The proteinase-activated receptor (F2R) is involved in vascular, hematopoietic development, and various physiological responses including inflammation, platelet activation, and intimal hyperplasia (Leger et al., 2006).The prothrombin (F2) receptor maintains vascular integrity, angiogenesis, and tissue repair.Tachykinin receptor (TACR1) plays a key role in the regulation of blood pressure, heart frequency, and stretching of blood vessels and as well as angiogenesis, ischemia, inflammation, and cardiac response to stress (Mistrova et al., 2016).
Structure-based drug discovery is a promising technique for the optimization of lead compounds rapidly and costeffectively.Molecular docking and dynamics studies were incorporated for the hub targets based on network pharmacology analysis.The top hit ligands such as Rutin, Tribuloside, Quercetin-3-O-a-rhamnoside, moupinamide, Kaempferol 3-rutinoside 7-glucuronide, Aurantiamide, 3",6"-Di-O-p-coumaroyltrifolin and Keioside were obtained that can act as potential lead compounds against cardiac disease- related targets.Previously, these compounds have reported antioxidant, anti-inflammatory, and cardioprotective activities (Ganeshpurkar & Saluja, 2017;Li et al., 2016).Although these putative lead compounds violate the Lipinski rule of Five for predicting drug-likeness properties, they showed better binding interaction and stability with the corresponding targets.The approved cardiac drugs such as digitoxin derived from the plant Digitalis lanata and reserpine from the roots of Rauwolfia serpentina also violate the drug-likeness properties (Shaito et al., 2020).Hence, an alternate strategy to the Lipinski rule (Lipinski et al., 2001) needs to be developed for natural products to identify the lead compounds.Further, the identified putative lead compounds from Tribulus terrestris L. need to be validated through in vivo or in vitro experiments to confirm their activity against cardiac diseases.

Conclusions
In the current study, HPLC-QTOF-MS/MS was used to identify putative compounds from Tribulus terrestris L. extract.Network pharmacology analysis combined with in silico molecular docking simulations to investigate the active ingredients and pathways of Tribulus terrestris L. in treating cardiac diseases.Based on molecular docking and MD simulation, we identified the potential lead compounds such as rutin, keioside, tribuloside, kaempferol 3-rutinoside 7-glucuronide, and moupinamide can act as promising drug candidates to treat cardiac diseases.However, further validation of these drug compounds should be assured based on in vivo and in vitro experiments.

Figure 1 .
Figure 1.Compound-target-disease network of Tribulus terrestris L. plant extract related to cardiac diseases.

Figure 2 .
Figure 2. (A) Protein-protein interaction network of targets related to cardiac diseases using the StringDB app installed in Cytoscape v3.7.2.(B) Module identification using MCODE plugin installed in Cytoscape.(C) Hub targets with the color codes representing from red to yellow (in terms of degree).

Figure 3 .
Figure 3. Heat plot representation of KEGG pathways and Gene Ontology for hub targets (A) Gene ontology analysis enriched in biological processes (B) cellular processes (C) molecular functions.(D) KEGG pathway analysis of hub targets.
Top three hit ligands for each hub targets represented in red color.

Figure 4 .
Figure 4. 3D surface view of the top hit ligand against hub targets for treating cardiac diseases.

Table 1 .
Results of extra precision molecular docking with docking score for 25 ligands from Tribulus terrestris L. plant extract.