Identification of antidiabetic constituents in Polygonatum odoratum (Mill.) Druce by UPLC-Orbitrap-MS, network pharmacology and molecular docking

2.1 Sample preparation

Dried P. odoratum rhizomes were purchased from Beijing Tong Ren Tang Co., Ltd (Lot Number: 20210802). 2.0 g of P. odoratum rhizome powder was extracted by ultrasonic with methanol (1:25, w/v) for 30 min. The extraction solution was forced to pass through 0.22 μm filter membrane (Millipore, Merck). The filtered methanol extraction solution was ready for injection.

2.2 UPLC-Orbitrap-MS analysis of methanol extract of P. Odoratum

UPLC-Orbitrap-MS analysis was conducted on the methanol extract of P. odoratum via the Thermo Fisher Scientific UPLC system (Dionex UltiMate 3000) coupled with an LTQ Orbitrap Velos mass spectrometer. XCalibur software (version 2.0.7) was used for instrument control, data acquisition, and data analysis. The chromatographic separation was performed on an Agilent ZORBAX SB C18 column (2.1 mm × 150 mm, 3.5 μm). 0.1 % formic acid in water (solvent A) and methanol (solvent B) as mobile phase were used at a flow rate of 0.3 mL/min. The gradient program was set as follows: 0–5 min, 5–20% B; 5–15 min, 20–45% B; 15–25 min, 45–80% B; 20–30 min, 80–100% B; equilibration time of 5 min at 5% B. A total of running time was 35 min. The column temperature was set at 35 °C, and the injection volume was 10 μL. The MS data were acquired using electrospray ionization (ESI) in the negative ionization modes, spray voltage, 4 kV (-4 kV in ESI-); sheath gas (N2, >95%), 35 bar; auxiliary gas (N2, >95%), 15 bar; heater temperature, 275 °C; capillary temperature, 275 °C. MS Scanning mode: full MS scan ranging from 50 m/z to 1500 m/z and the resolution was 35,000.

2.3 Disease targets of the compounds in methanol extract of P. odoratum

Traditional Chinese Medicine Systems Pharmacology Database (TCMSP, https://tcmspw.com/tcmsp.php) was initially searched by using chemical names, InChIKeys, and CAS numbers of all the compounds identified in the methanol extract of P. odoratum via UPLC-Orbitrap-MS (Ru et al., 2014). A total of 183 components were identified; which corresponded to 313 targets. For the compound that cannot be found in TCMSP, we queried the Simplified Molecular Input Line Entry System (SMILES) numbers of the corresponding components in the PubChem database https://pubchem.ncbi.nlm.nih.gov/ and import them into Swiss Target Prediction https://swisstargetprediction.ch/ platform with screening threshold > 0.1. A total of 94 components were obtained, which corresponded to 673 targets. Taken together, we obtained 277 chemical components and 897 drug targets. The names of drug targets were then converted into gene symbols through the UniProt database https://www.uniprot.org/ (Consortium, 2020).

2.4 Disease targets of T2DM

To obtain the target genes related to diabetes among all identified targets, we used “type 2 diabetes mellitus” as the keyword to obtain 15,776 gene targets related to T2DM from DrugBank https://go.drugbank.com/, Online Mendelian Inheritance in Man (OMIM) https://www.omim.org/ and Gene Cards https://www.genecards.org databases (Hamosh et al., 2005, Safran et al., 2021, Wishart et al., 2018). All obtained targets were validated and calibrated against the UniProt database (Consortium, 2020). In order to obtain T2DM-related drug targets in terms of P. odoratum methanol extract, we used the online analysis tool Venny (https://bioinfogp.cnb.csic.es/tools/venny/) to get the intersection of identified targets of P. odoratum and T2DM. A total of 703 common targets were collected in this study.

2.5 Construction of the compounds-targets network

703 targets for T2DM and the corresponding 277 components in the methanol extract of P. odoratum were inputted into CytoScape (version 3.7.1) for the construction of the “component-target” network (Killcoyne et al., 2009). In order to identify core bioactive components, we used the “Analyze Network” function in CytoScape to perform a topology analysis of the network and calculate the degree value of each node in the network (Chin et al., 2014). Based on the degree numbers, the top five core active ingredients were selected. To further demonstrate the key “ingredient-target” network, only T2DM targets associated with the top five active ingredients were kept drawing a sub-network of the interaction between the core active ingredients and disease targets.

2.6 Construction and analysis of protein–protein interaction network

In order to understand the relationships among T2DM targets, we imported the targets corresponding to the active components of methanol extract of P. odoratum into the STRING database https://cn.string-db.org/ in order to obtain the protein–protein interaction (PPI) network. “Multiple proteins” were selected for the specific input method and “Homo sapiens” for species to remove the independent nodes in the network. TSV files were downloaded in the analysis results, which were then imported into CytoScape (version 3.7.1) for data visualization and network analysis. The Maximal clique centrality (MCC) algorithm in the Cytohubba plugin was used as the screening tool for the core target (Chin et al., 2014). All parameter settings were in default, and the nodes with the top five scores were selected as the hub nodes of the network. The first neighbors of each hub node in the network were grouped into one category. In the PPI network, nodes represented proteins and edges represented interactions between proteins.

2.7 Functional classification of target proteins in T2DM

In addition to studying interactions among targets in T2DM, we also used the online classification system PANTHER https://pantherdb.org/to classify T2DM-related proteins for their functions (Mi et al., 2021). First of all, we imported the proteins into PANTHER, selected “ID List” for List Type, “Homo sapiens” for Organism, and “Functional classification in gene list” for analysis mode to display protein classification in pie chart. Based on the number of proteins (genes) covered, the top five protein categories were selected and the list information was downloaded for each protein category. In order to visualize the relationship between the active components of P. odoratum and protein classification, we used the online drawing platform HIGHCHARTS (https://www.highcharts.com.cn/) to draw a Sankey diagram so as to show the relationship between the five core active components of P. odoratum and the top five protein classes (PC).

2.8 Enrichment analysis

Disease Ontology (DO) function in the DOSE package (version 3.22.0) was used to enrich disease-related target proteins into pathways (Yu et al., 2014). In order to screen the enriched pathways with statistical difference, the specific parameters of DO were set to ont = “DO”, P-value Cutoff = 1 and Q-value Cutoff = 1. According to the size of Gene Ratio, a bubble chart was drawn for the result. For Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses, UniProtKB IDs obtained from the UniProt database for all the target proteins were imported into the database DAVID https://david.ncifcrf.gov/summary.jsp with the identifier selected as “UNIPROT_ACCESSION” to perform GO and KEGG enrichment analysis, the results of which were downloaded for further analysis (Huang da et al., 2009; Sherman et al., 2022). For the GO enrichment analysis based on the number of differentially expressed genes, the top 10 GO terms based on P-values in Biological Process (BP), Molecular Function (MF), and Cellular Component (CC) were selected to draw a bubble chart. For KEGG pathway enrichment analysis, the top 30 terms with based on P-values were selected to draw a bubble chart. The smaller the P-value, the more significant the enrichment. Redder nodes and larger node sizes in bubble chart indicated that the pathways were more important in the entire network.

2.9 Molecular docking

Molecular docking is a convenient method for screening the genuine interactions between active ingredients and target proteins, which also facilitated the understanding of the binding mechanism (Morris and Lim-Wilby, 2008). In this study, we conducted molecular docking to test whether the core active components in the ethanol extract of P. odoratum identified by network pharmacology were genuinely combined with the core target proteins. First of all, we downloaded the 3D Conformer small molecule files in sdf format corresponding to the core components from the PubChem (https://pubchem.ncbi.nlm.nih.gov) database and used Open Babel (version 3.1.1) software to convert sdf files into mol2 file. For the 3D structures of top 5 target proteins, their pdb files were sourced from the RCSB https://www.rcsb.org database, that is, actin Cytoplasmic 1 (ACTB) (PDB ID: 2OAN), transcription factor activator protein AP-1 (JUN) (PDB ID: 5VPC), signal transducer and activator of transcription 3 (STAT3) (PDB ID: 6NJS), hypoxia-inducible factor 1 alpha (HIF1A) (PDB ID: 1LQB), and matrix metalloproteinase 9 (MMP9) (PDB ID: 5UE3). PyMol (version 2.4.0) software was used to remove water molecules and redundant ligands for pre-docking of small molecule components to target proteins. AutoDock4 (version 1.5.7) was then used to hydrogenate proteins that were set as acceptors, the results of which were generated as pdbqt files. Small molecules were also conducted for hydrogenation operation, and torsion trees were automatically configured in the software, saving as a ligand file in pdbqt format. The pdbqt structure files of receptor and ligand were imported into AutoDock4 to define the molecular docking range. The target protein was taken as the center of the grid, while the parameters of the center coordinates (center X/Y/Z) and box size (size X/Y/Z) were adjusted to ensure that the protein macromolecules were completely covered by the docking box before molecular docking was conducted.

In AutoDock4, molecular docking was accomplished by detecting protein macromolecules and inserting tiny drug molecules, as well as configuring operating methods and docking parameters. The pdbqt format was used to compute the minimal binding energy. We converted the pdbqt file to PDB format to demonstrate hydrogen bonds between ligands and proteins, calculated the distance between the protein residue stick structure and the hydrogen bond. The final generated receptor and ligand complex PDB files were imported into LigPlot++ (version 2.2.5) software to visualize and analyze the two-dimensional interaction diagram of ligands and target proteins. In addition, we used the “pheatmap” package to draw a heat map of docking binding energy to visualize the binding energy between different receptors and ligands, and used the binding energy of these active components and target proteins as the evaluation criteria for the possibility of molecular docking. The smaller the binding energy, the greater the possibility of docking. In this study, the binding energy threshold was set to −4 kcal/mol. If the docking binding energy was less than −4 kcal/mol, it was considered that the binding ability between the ligand and the receptor was in good condition.

3.1 Chemical composition

The total ion chromatogram (TIC) of UPLC-Orbitrap-MS in the negative ion mode is shown in Fig. 1. A total of 334 compounds belonging to five categories were identified, which were flavonoids, alkaloids, terpenoids, saponins, and organic acids. Supplementary Table S1 contains the retention time, molecular formula, molecular ion peak, and adduct ion of all the identified compounds. The results showed a relatively good mass spectrometry response in the negative ion mode.

Close

Fig. 1

UPLC-Orbitrap-MS total ion current chromatogram of the methanol extract of P. odoratum in the negative ion mode. Representative components were labelled by numbers in the chromatogram. 1: 5,7,3′-Trihydroxy-6,4′,5′-trimethoxyflavone; 2: 9-Methoxycamptothecin; 3: 9-Aminocamptothecin; 4: Diacerein; 5: 5-Hydroxy-1-tetralone.

Export to PPT

3.2 Analysis of compound-target network

Due to the multi-component and multi-target characteristics of TCM functions, we constructed a component-target network for the methanol extract of P. odoratum in Fig. 2. The degree value of the component was calculated in the network according to the number of targets linked to each component. According to the degree values, the top 5 components were highlighted in the center as the five hub components with unique protein targets clustered in a circle and connected individually (see Table 1), while the shared protein targets formed the outermost circles; in addition, the remaining active components formed the purple circles, and the green circles represented the remaining targets (Fig. 2A). To highlight the hub components of P. odoratum, we also plotted a sub-network of hub component-target protein interaction (Fig. 2B), in which the shared components were displayed in the central area and separated the shared nodes of adjacent hub components to show their association.

Close

Fig. 2

Network visualization of multi-components and multi-targets relationships for the methanol extract of P. odoratum. (A) Interactional network of active components in methanol extract of P. odoratum and T2DM disease target proteins. (B) Network of hub components and corresponding T2DM target proteins. Node sizes in the network indicated the importance of the nodes in the network.

Export to PPT

Table 1 Chemical structures of the five representative components with antidiabetic effects identified from the methanol extract of Polygonatum odoratum (Mill.) Druce.

No.	Compound Name	Chemical Structure
1	5,7,3′-Trihydroxy-6,4′,5′-trimethoxyflavone
2	9-Methoxycamptothecin
3	9-Aminocamptothecin
4	Diacerein
5	5-Hydroxy-1-tetralone

3.3 PPI network construction and analysis

A PPI network was created from the obtained 703 common targets via STRING database, which was then imported into CytoScape for further visualization. The sizes represent the node degree in the network (Fig. 3A). In specificity, the PPI network contains 703 nodes and 16,494 edges, the average node degree is 46.9, and the expected number of edges is 7580, which means that these proteins have more functional and physical interactions with each other and can be biologically linked as a group. We used the size and color of the node to represent the interaction between proteins, where the larger the size and the redder the color, the more the node is associated with other nodes. In addition, Five core targets including Actin Cytoplasmic 1 (ACTB), Transcription factor AP-1 (JUN), Signal transducer and activator of transcription 3 (STAT3), Hypoxia-inducible factor 1 alpha (HIF1A), and Matrix metalloproteinase 9 (MMP9) and their directly connected targets in Fig. 3B. The common targets of 5 core targets were clustered in the center, and the associated targets between two of the core targets were shown in the adjacent area, and the targets associated with a single core target were shown in the outermost layer.

Close

Fig. 3

PPI network of active components in methanol extract of P. odoratum and T2DM disease targets. (A) PPI network of intersection targets of active components of P. odoratum. (B) Core target network. The size of a node in the network indicates how important the node is in the network.

Export to PPT

3.4 Classification of core target proteins

Functional classification of all the target proteins were performed and the results were shown in the pie chart in Fig. 4A, according to which the top five subclasses were Protein Modifying Enzyme (PC00260), Metabolite Interconversion Enzyme (PC00262), Transmembrane Signal Receptor (PC00197), Transporter (PC00227) and Gene-Specific Transcriptional Regulator (PC00264) (Mi et al., 2021). In order to check which protein subclasses of the potential key active components of P. odoratum played a role, we plotted the top five protein classes Sankey diagram in Fig. 4B. The associations between components and proteins were shown according to the connections among hub components-core targets-protein classes, and it could be seen in Fig. 4B that Protein Modifying Enzyme (PC00260) was associated with 80 core targets and Metabolite Interconversion Enzyme (PC00262) was associated with 56 core targets. The results indicated that these two protein subclasses were enriched with more genes (proteins) related to T2DM and might play key roles in the treatment with P. odoratum active ingredients.

Close

Fig. 4

Pie chart of protein classification and Sankey diagram for top five protein classes. (A) Pie chart of PANTHER protein classes. (B) Sankey diagram of top five protein classes containing the most amount of core targets. The five hub components were shown on the left, while the top five protein classes with gene hits were displayed on the right, which were linked by all core targets of the five hub components.

Export to PPT

3.5 Enrichment analysis

As for the enrichment analysis, a total of 703 target proteins were first analyzed via DO enrichment, according to which 706 DO terms were significantly enriched (P-value < 0.05). In the DO enrichment (Fig. 5A), the most significant disease pathways included urinary system disease, over-nutrition, kidney disease, and obesity. Furthermore, DO enrichment analysis also revealed multiple cancer-related disease signals. Since DO enrichment only suggested diseases associated with the input target proteins, we also performed GO functional enrichment and KEGG pathway enrichment analyses to further investigate the mechanisms of these core target proteins. For BP, CC, and MF, the core targets for the treatment of T2DM were analyzed using the DAVID database. The top 10 items of GO with statistical significance (P-value < 0.05) for each category were shown in Fig. 5B. The plasma membrane, an integral component of the plasma membrane, response to drugs, protein phosphorylation, protein kinase activity and enzyme binding and other metabolic aspects were most relevant to GO functional enrichment. Similarly, for the KEGG enrichment analysis, P-value < 0.05 was used as the threshold to obtain the signaling pathways in Fig. 5C. The top 30 items of KEGG results were mainly composed of signaling pathways such as pathways in cancer, neuroactive ligand-receptor interaction, lipid and atherosclerosis and age-rage signaling pathway in diabetic complications.

Close

Fig. 5

Enrichment analysis of target proteins. (A) DO disease enrichment analysis. (B) GO functional enrichment analysis. (C) KEGG pathway enrichment analysis. The abscissa ratio in the figure represented the percentage of the number of proteins enriched in each term or pathway to the total core target proteins, and the ordinate was the enriched term or pathway.

Export to PPT

3.6 Molecular docking study

In order to validate the core components that have potential roles in the treatment of T2DM, we calculated the binding affinities of the top five bioactive components in P. odoratum and target proteins of T2DM. According to the heat map of the binding energy (Fig. 6), it could be found that the binding energies between active components and target proteins were less than 0 kcal/mol, indicating that the top 5 active components had the potential to bind to the top 5 target proteins. In particular, the binding energy between JUN and CID 5496475 was the highest at −0.44 kcal/mol, which suggested that the binding between the protein and the compound was very weak. As for other interactions, STAT3 and CID 1236175 had the binding energy of −4.12 kcal/mol, HIF1A and CID 1236175 at −4.02 kcal/mol, and MMP9 and CID 72402 at −4.56 kcal/mol. In particular, MMP9 and CID 1236175 had the smallest value of the binding energy value at −4.69 kcal/mol, indicating that a strong binding occurred between the two molecules.

Close

Fig. 6

Heatmap of binding energies between potential active ingredients and core disease targets. The left side of the figure is the ID number of the active ingredient in PubChem, and the top of the figure is the core target name. The darker the color (red), the lower the binding energy and the easier the binding between the receptor and the ligand. CID 72402: 9-Aminocamptothecin, CID 123617: 9-Methoxycamptothecin, CID 14057196: 5-Hydroxy-1-tetralone, CID 5496475: 5,7,3′-Trihydroxy-6,4′,5′-trimethoxyflavone, CID 26248: Diacerein.

Export to PPT

The optimal molecular docking diagrams of receptors and ligands were exhibited at Fig. 7. It could be seen that 9-Methoxycamptothecin formed two hydrogen bonds with GLN in STAT3 with bond lengths of 2.1 and 3.5, respectively (Fig. 7A). In addition, 9-Methoxycamptothecin formed a hydrogen bond with both PRO-496 and ILE-498 in HIF1A with bond lengths of 2.8 and 2.2, respectively (Fig. 7B). 9-Aminocamptothecin formed two hydrogen bonds with LYS-92 in MMP9 with bond lengths of 2.1 and 2.5 and formed one hydrogen bond with both ARG-95 and ASP-185 with bond lengths of 3.1 and 2.4, respectively (Fig. 7C). 9-Methoxycamptothecin formed one hydrogen bond with GLY-215 in MMP9 with a bond length of 2.0 and interacted with PHE-250 via two hydrogen bonds with bond lengths of 1.9 and 3.6 (Fig. 7D). For each receptor and ligand combination, we also mapped the protein–ligand interaction. For molecular docking diagrams between the remaining target proteins and active components, please see Supplementary Figure S1.

Close

Fig. 7

Molecular docking between active components in P. odoratum and core target proteins of T2DM disease. (A) 9-Methoxycamptothecin and STAT3. (B) 9-Methoxycamptothecin and HIF1A. (C) 9-Aminocamptothecin and MMP9. (D) 9-Methoxycamptothecin and MMP9. Yellow dotted lines indicated hydrogen bonds, while the numbers indicated bond lengths.

Export to PPT