Molecular Networking, Network Pharmacology, and Molecular Docking Approaches Employed to Investigate the Changes in Ephedrae Herba before and after Honey-Processing

Raw Ephedrae herba (REH) and honey-processed Ephedrae herba (HEH) were the different decoction pieces of Ephedrae herba (EH). Honey-processing that changes REH into HEH has been shown to relieve cough and asthma to a synergistic extent. However, the chemical markers and the synergistic mechanism of HEH need to be further studied. In this study, the ultra-high performance liquid chromatography coupled with hybrid quadrupole time of flight mass spectrometry (UPLC-Q-TOF-MS) and molecular networking (MN) were used to investigate the chemical composition of REH and HEH, which led to the identification of 92 compounds. A total of 38 differential chemical markers for REH and HEH were identified using principal component analysis (PCA) and orthogonal partial least squares discriminant analysis (OPLS-DA). Network pharmacology suggests that the synergistic effect of HEH in relieving cough and asthma may be due to 31 differential chemical markers acting through 111 biological targets. Among them, four compounds and two targets probably played an important role based on the results of molecular docking. This study enriched our knowledge about the chemical composition of REH and HEH, as well as the synergistic mechanism of HEH.


Introduction
According to Leigong Paozhi Lun (a Chinese ancient book), "honey-processing allows drugs become sweet and sluggish with effects of benefiting Qi, moistening lung, relieving cough, and stopping pain and dysentery" [1]. Modern research has shown that honeyprocessing playing a synergistic role to enhance the therapeutic effect, but also moderating the medicinal properties to reduce side effects, such as honey-processed astragalus [2], honey-processed licorice [3], etc.
Ephedrae herba (EH), deriving from the dried stems of Ephedra sinica Stapf, is an herbal plant that commonly grows in northern China [4], which has been widely used in traditional Chinese medicine for the treatment of common cold, coughs, asthma, and edema for thousands of years [5]. The main components of EH include alkaloids, flavonoids, tannins, polysaccharides, organic acids, volatile oils, and many other active compounds [6]. Figure 1A shows the total ion flow diagram for fraction P1, and Table 1 shows the identified compounds. The MN was constructed based on MS/MS similar spectra for visualization as shown in Figure 2A. The constructed MN showed a total of 734 precursor ions visualized as nodes in the molecular map, which include 27 groups (node > 2) and 517 unique nodes. The proportion of the red and blue sectors at each node in MN represents the relative content of the compound in REH and HEH. All the compounds were determined from MS data combined with the MN chemical composition database and published literature. According to MS/MS fragment breakage patterns and MN derivation, 65 components were identified in eight major groups (groups a-h). All the structures are shown in Figure S1.     MN for the nodes annotated as "quinolinic acids" class (b). MN for the nodes annotated as "phenolic acids" class (c, d, and e). MN for the nodes annotated as "alkaloids" class (f, g, j and k). The pie chart for nodes was filled with different colors, red (raw Ephedrae herba) and blue (honey-processed Ephedrae herba), based on the proportion of the intensity of the ion peak corresponding to each metabolite in the two decoction pieces. Peak 6 loses a moiety of m/z 132 and then shows a molecular ion peak at m/z 152.0566, which is similar to adenosine. In MN, peak 6 as a node linking to cordycepin, adenosine, and adenine shows structural correlation and is 16 Da different from adenosine. Therefore, it is determined that peak 6 is isoguanosine.

Others
Three other compounds were also first identified from EH, including glucosyringic acid (

The Analysis of Fraction P2
The total ion flow diagram of the fraction P2 is shown in Figure 1B, and the identified chemical compositions are shown in Table 2. The constructed MN ( Figure 2B) was shown a total of 742 precursor ions visualized as nodes in the molecular map, which include 26 groups (node > 2) and 456 unique nodes. At each node in MN, the proportion of red and blue sectors represents the relative content of the compound in REH and HEH. The MS data along with the MN chemical composition database and published literature were used to determine these compounds. A total of 38 components including four major groups (group i-m) were identified based on the MS/MS fragment fracture pattern and MN derivation. The structures of all 38 compounds are shown in Figure S2.  Based on the MN, 92 components were identified by the MN method, and the structure and detailed mass spectral information of some of the markers are shown in Figure 2, Tables 1 and 2, respectively. Figure 2A shows that the red sectors of the alkaloid composition for group f are larger than the blue sectors. This indicates that the content of alkaloids, the main components of EH, decreases slightly after processing. In group a, vicenin-2, schaftoside, kaempferol-3-glucoside-7-rhamnoside, and quercetin-3-rhamnoside have smaller proportions in the red sectors than in the blue, indicating an increase in their levels after processing. At the same time, the red sectors of gallocatechin-(4 → 6"; 2 → O → 7")-(epi) gallocatechin is greater than the blue, indicating that the content of this component decreases after processing.

PCA and OPLS-DA Analyses for the Identified of Discriminatory Metabolites
To verify the accuracy of the MN data analysis, a comparative analysis of the REH and HEH data files was performed. The PCA score spots of REH and HEH were shown in Figure 3A,C. OPLS-DA analysis was shown in Figure 3B,D, indicating significant differences in the components of REH and HEH. Furthermore, the screening of differentially abundant metabolites in REH and HEH was performed based on the variable importance in projection value (VIP > 1) and p-value (p < 0.05); then, 38 differential chemical markers were screened. Heatmap analysis of these differential chemical markers ( Figure 3E) showed that 7b, 9-Dihydroxy-3-(hydroxymethyl)-1,1,6,8-tetramethyl-5-oxo-1,1a,1b,4,4a,5,7a,7b,8,9-decahydro-9aH-cyclopropa [3,4]benzo[1,2-e]azulen-9a-yl acetate, norephedrine, m-cumenyl methylcarbamate, L-norpseudoephedrine, L-ephedrine, dibutyl phthalate, pseudoephedrine, hexaxdecanoic acid and methylephedrine decreased after honey-processing, while other components increased, which is consistent with the results obtained previously by MN. This may be significantly correlated with the synergistic effect of HEH.

Prediction of the Target Proteins
In total, 1457 targets related to 31 components of the 38 differential chemical markers and 984 targets associated with diseases treated by HEH were identified. The 111 overlapping component and disease targets were obtained by Venn diagrams ( Figure 4A) that is intersecting the component targets and the disease targets.  the PCA score plot with 95% density ellipses (A,C); the OPLS-DA score plot with 95% density ellipses (B,D); the heatmap of differential compounds related to REH and HEH (E).

Prediction of the Target Proteins
In total, 1457 targets related to 31 components of the 38 differential chemical markers and 984 targets associated with diseases treated by HEH were identified. The 111 overlapping component and disease targets were obtained by Venn diagrams ( Figure 4A) that is intersecting the component targets and the disease targets.

Construction of the Component-Disease-Target Network
A component-disease-target network was developed to identify the relationships between the 31 active components and 111 overlapping component and disease targets, and Cytoscape V3.8.2 (La Jolla, CA, USA) was used for network visualization and topological analysis. The component-disease-target network contained 143 nodes and 394 edges ( Figure 4B).

Construction of the PPI Network
The 111 overlapping component and disease targets, which were considered the therapeutic targets of HEH in the treatment of cough and asthma, were imported into the STRING database, and the PPI network was constructed. The hub nodes of degree ≥ 2 D with Closeness Centrality ≥ CC and Betweenness Centrality ≥ BC were selected as the core PPI network ( Figure 4C). The sizes and colors of the nodes are proportional to the degree. The larger the node and the darker the color, the stronger is the interaction, indicating that the interaction plays a more central role in the PPI network. In total, 28 core targets and 376 edges representing protein-protein interactions were predicted. The top 10 core targets were TNF (degree = 85), IL6 (degree = 84), IL1B (degree = 79), ALB (degree = 76), AKT1 (degree = 74), IL10 (degree = 72), CXCL8 (degree = 71), CCL2 (degree = 67), INS (degree = 66), and IL4 (degree = 66) ( Table 3).

GO and KEGG Enrichment Analyses
The 111 overlapping components and disease targets were imported to OmicShare tools for GO functional enrichment and KEGG pathway enrichment analyses. In total, 833 GO terms with p < 0.05, including biological processes, cellular components, and molecular functions, were obtained. From these terms, the functional information and secondary classification histogram were drawn with OmicShare tools ( Figure 4D). The highest correlations in GO enrichment were associated with cell motility, movement of intracellular components, and inflammation, including positive regulation of cell migration, positive regulation of cell motility, positive regulation of cellular component movement, and positive regulation of locomotion. In total, 132 KEGG pathway items with p < 0.05 were obtained, and we selected the 20 top-ranking pathways based on their p-values and drew bubble plots with OmicShare tools ( Figure 4E).
The largest number of KEGG-enriched targets was Lipid and atherosclerosis, which included 15 targets. Cough and asthma treatment with EH was primarily associated with anti-inflammatory activity, and the enriched signaling pathways for TNF, MAPK, and PI3K-Akt were associated with inflammation. The PI3K-Akt pathway and asthma were associated with the main synergistic effects of HEH on cough and asthma [17,18], with the PI3K-Akt signaling pathway being enriched to a total of 10 targets. PI3K is a member of the lipid kinase family, and Akt is a downstream actor of PI3K. The PI3K-Akt pathway is involved in a variety of roles in vivo, such as antioxidant [19], anti-inflammatory [20], etc. MAPK signaling pathway is also a pathway of anti-inflammatory mechanism in vivo [21].
In summary, these differential chemical markers act on multiple targets and pathways to exert their pharmacological effects, and the variation in the levels of the differential chemical markers differentiates the efficacy of HEH from that of REH. GO and KEGG analyses suggest that HEH may have therapeutic and inflammatory effects as well as antiviral and antioxidant effects. This study provides more ideas for the pharmacological study of HEH.

Component-Target Molecular Docking
L-Ephedrine and pseudoephedrine, the main components of EH, quercetin and tricin, whose content increases after honey-processing, were selected to dock with the core target proteins TNF and IL6. The above compounds acting on these target proteins may exert anti-inflammatory and anti-cancer effects through biological pathways such as apoptosis and metabolic cell death. It is generally accepted that the more negative the docking affinity, the more likely it is to bind, and the result of the docking affinity is recorded in Figure 5.
The results show that all four compounds have the potential to bind to all two targets, and that the binding of components and amino acid residues is mainly through hydrogen bonds and van der Waals forces. The molecular docking results further validate that HEH may play an important role in the pathways screened out by network pharmacology and provide data to support further studies of HEH. The results show that all four compounds have the potential to bind to all two targets, and that the binding of components and amino acid residues is mainly through hydrogen bonds and van der Waals forces. The molecular docking results further validate that HEH may play an important role in the pathways screened out by network pharmacology and provide data to support further studies of HEH.

Plant Material and Chemicals
LC-MS grade methanol, acetonitrile, and formic acid were purchased from Fisher Scientific (Pittsburgh, PA, USA). Ultrapure water was prepared by a Synergy water purification system (Millipore, Billerica, MA, USA). The internal standard of hyperoside (≥98.0%, Lot no. Y20A9X59340) was purchased from Chengdu Push-Biotechnology Co. Ltd. (Chengdu, China). Other chemicals and reagents were of analytical grade. Ephedrae herba was acquired across the major pharmacies and identified as Ephedra sinica Stapf by Dr. Dan Zhang. The details of each sample are listed in Table 4. The specimens were stored at Hebei University of Chinese Medicine.

Sample Preparation for UPLC-Q-TOF-MS Analysis
An aliquot of 1.00 g of sample powder was immersed in 20 mL of 65% methanol (v/v), followed by ultrasonic extraction at room temperature for 30 min. The mixture was then centrifuged at 13,000 rpm for 10 min, and 1 mL of the extracts was separated into fraction P1 and P2 using a semi-preparative HPLC (Shimadzu Co., Kyoto, Japan). P1 was then diluted 60-fold in 65% methanol and P2 20-fold, 0.00125 mg/mL of hyperoside as an internal standard.

UPLC-Q-TOF-MS Analysis
The UPLC-MS analysis was performed on an Agilent 1290 Infinity II system coupled with an Agilent 6545 quadrupole time-of-flight mass spectrometer system (Q-TOF-MS) (Agilent Technologies, Santa Clara, CA, USA) equipped with an electrospray ionization interface.
A Waters Acquity UPLC HSS T3 column (2.1 × 100 mm, 1. The MS acquisition parameters were as follows: drying gas (N2) temperature, 320 • C; sheath gas temperature, 350 • C; drying gas (N2) flow rate, 10.0 L/min; sheath gas flow (N2) rate, 11 L/min; nebulizer gas pressure, 35 psi; capillary voltage, 3500 V; fragmentor voltage, 135 V; collision energy, 40 eV. The analysis was operated in a negative mode with the mass range of m/z 100-1000 Da. Data were analyzed by MassHunter Qualitative Analysis Software Version B.10.00 (Agilent Technologies, Palo Alto, CA, USA). Quality control (QC) samples were prepared by pooling the same amount of samples together, and every 5th was utilized as a QC sample.

Establishment of Molecular Networking
The MN was constructed using the UPLC-Q-TOF-MS/MS data from REH and HEH. All MS/MS data files were converted into 32-bit mzXML by using MSconvert software.

Mass Spectrometry Analysis and Metabolites Annotation
The LC-MS data acquisition was conducted on a MassHunter Workstation (Agilent Technologies). After normalization, the data set was introduced into SIMCA software (version 13.0, Umetrics, Umea, Sweden), for PCA and OPLS-DA. Heatmap analysis was generated by OriginPro 2019b software (Origin Lab Corporation, Northampton, MA, USA).

Network Pharmacology
To reveal the correlation between the differential compounds and the synergistic effect of relieving cough and asthma of HEH, and to predict the potential targets and pathways closely associated with the synergistic of HEH from a comprehensive perspective, a network pharmacology study was performed in this study. Secondly, with relieving asthma, asthma, cough, relieving cough, anti-inflammatory, inflammatory and antioxidant as the keywords, the related targets were collected using the following database: GeneCards (https://www.genecards.org, 10 April 2022), Online Mendelian Inheritance in Man (OMIM, https://www.omim.org, 10 April 2022), Therapeutic Target Database (TTD, http://db.idrblab.net/ttd, 10 April 2022), TCMSP, the Encyclopedia of Traditional Chinese Medicine (ETCM, http://www.tcmip.cn/ETCM/index.php/Home/ Index, 10 April 2022), and SymMap (http://www.symmap.org, 10 April 2022). After the collected targets were merged and removed duplication, the core targets were obtained by intersecting the targets of differential compounds and diseases based on the STRING database (https://cn.string-db.org, 10 April 2022).
Thirdly, each core target (enzymes, receptors, transporters, cytokines, proteins, and others) was classified to construct networks of differential compounds-core targets and the protein interaction (PPI). The network was generated by Cytoscape software (version 3.8.2) (http://www.cytoscape.org, 12 April 2022) with topological analysis. Through the topological analysis, the main differential compounds and the potential targets were screened out.

Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) Pathway Enrichment Analysis
The targets were input to the Database for Annotation, Visualization, and Integrated Discovery (DAVID) (https://david.ncifcrf.gov/summary.jsp, 13 April 2022) for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, and species were selected as "homo sapiens". OmicShare's platform tool (http://www.omicshare.com/tools, 13 April 2022) was used to draw the bubble plot. The color of the dot represents different p-values, and the size of the dot reflects the number of target genes expressed in the pathway. The rich factor represents the ratio of all target genes in a pathway to the number of all the annotated genes in the pathway. A higher rich factor represents a higher level of enrichment.

Component-Target Molecular Docking
AutoDock Tools (vision 4.2.6, La Jolla, CA, USA) is a receptor-ligand docking simulation program used for protein-ligand docking simulation and prediction of the docking affinity [22]. Molecular docking analysis was adopted to confirm the interactions between the core components and the core targets of REH and HEH in the treatment of asthma and to verify the accuracy of the network pharmacology prediction. The threedimensional (3D) structures of the target proteins were downloaded from the PDB database (https://www.rcsb.org/, 18 April 2022), and the MOL2 structures of the differential compounds were downloaded from the TCMSP database. We removed the water molecules, separated the proteins, added nonpolar hydrogens, calculated the Gasteiger charges for each structure, and saved the data as a PDBQT file. The target proteins were receptors, and the active components were ligands. AutoDock Vina 1.1.2 (La Jolla, CA, USA) was used to dock molecules with proteins. The Vina score is the core of the complex obtained by docking the receptor and ligand with the corresponding pocket parameters using the Vina procedure. The lower the Vina score is, the higher the affinity of the receptor and ligand. It is generally believed that, when the conformation of the ligand and receptor is stable, the lower the energy and the greater the possibility of interaction. When the docking affinity is less than 0, the ligand and the receptor can bind spontaneously. The conformation with the best affinity was selected as the final docking conformation and visualized with PyMOL 2.3 (New York, NY, USA).

Conclusions
In this study, UPC-Q-TOF-MS combined with the MN were used to identify the components of REH and HEH for the first time, and 92 compounds were tentatively identified, including 27 flavonoids, 9 alkaloids, 26 phenolic acids, 3 quinolinic acids, 6 fatty acids, and 21 others, and 38 differential chemical markers were screened out by MN and OPLS-DA. The synergistic mechanism of HEH was studied using network pharmacology and molecular docking. Thirty-one active components were probably acting through 111 biological targets, and the enriched signaling pathways for TNF, MAPK, and PI3K-Akt were associated with cough and asthma. Among them, four components and two targets probably played an important role. This study enriched our knowledge about the chemical composition and the synergistic mechanism of HEH.