The Traditional Chinese Medicine Systems Pharmacology Database and analysis platform (TCMSP) (https://tcmsp-e.com/tcmsp.php) [7] is based on the traditional Chinese Medicine Systems pharmacology framework. The TCMSP includes 837 related disorders, 29,384 constituents, 3311 targets, and 499 varieties of Chinese therapeutic plants from the Chinese Pharmacopeia. The TCMSP also contains 12 crucial Absorption, Distribution, Metabolism, Excretion (ADME)-associated attributes of drugs, including drug similarity, human oral bioavailability (OB), Caco-2 permeability, half-life, blood–brain barrier integrity, and compliance with Lipinski’s rule-of-five for drug screening and assessment. Compositional information for Lycium barbarum was sourced from the TCMSP database.

ADME [8] refers to the Absorption, Distribution, Metabolism, Excretion, and Toxicity of drugs. Thus, the ADME content is applied to contemporary drug design and screening in pharmacokinetics research. This research is based on the characteristics of the ADME TCMSP databases and provides relevant data with a selected oral exploitation degree OB of $>$30%, and a drug-likeness (DL) value of $>$0.18, corresponding to a chemical composition as an efficient ingredient.

We used a two-step procedure to predict and screen the target proteins of effective Chinese medicine ingredients. First, the TCMSP database was utilized to predict protein targets for the active components of Lycium barbarum. Subsequently, the predicted targets in the TCMSP database were transformed into gene names using UniProtKB (http://www.uniprot.org) [9].

The GeneCards database (https://www.genecards.org/) offers extensive information on genes for individuals [10]. Therefore, we identified a total of 445 ferroptosis-related genes (FRGs) in the GeneCards database using the keyword “ferroptosis” and a relevance score of $>$1 as the screening criterion. Meanwhile, employing “ferroptosis” as the keyword in the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/) [11] resulted in obtaining an additional 64 FRGs from the reference gene set WP_FERROPTOSIS. Additionally, we downloaded “marker”, “suppressor”, and “driver” genes associated with ferroptosis from the FerrDb database (http://www.zhounan.org/ferrdb/current/) [12], from which 339 human protein-coding FRGs were screened; Supplementary Table 1 lists the FRGs from the three sources combined, yielding 626 FRGs for further analysis.

The search term “osteoradionecrosis of the jaw” was used to examine the GeneCards database. A total of 35 genes linked to ORNJ targets, with a relevance score $>$1, were found after screening (Supplementary Table 2).

The intersection of targeted proteins for screened Chinese herbal ingredients and screening genes associated with ORNJ (ORNJ-related genes, or ORNJRGs) yielded Lycium barbarum.

The screened genes relating to ORNJ and ferroptosis were intersected to identify any FRGs in Lycium barbarum therapy. These were designated as hub genes. Individual proteins that interact with one another form protein–protein interaction networks. The STRING database (https://cn.string-db.org) [13] can be employed to search for connections between anticipated and existing proteins. Using this database, the biological species was set to human with a confidence level of $\ge$0.400, and the hub genes were used to build a protein–protein interaction network that was visualized using Cytoscape 3.9.1 (Cytoscape Consortium, San Diego, CA, USA). Additionally, genes with similar functions were also predicted, and an interplay network was established through the GeneMANIA website (https://genemania.org).

Gene Ontology (GO) analysis is a general method for large-scale functional enrichment, including both molecular functions (MFs) and biological processes (BPs). The Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www.kegg.jp) is an extensive database that contains information on genomes, illnesses, biological pathways, and drugs [14]. The R package clusterProfiler [15] was employed to perform GO annotation analysis of hub genes, with entry screening standards of p.adjust 0.05 and a false discovery rate (FDR) value (q-value) of 0.05 using the Benjamini–Hochberg method for p-value correction.

The ENCORI database (https://rnasysu.com/encori/) [16] is version 3.0 of the starBase database. Interactions in ENCORI, including RNA-binding protein (RBP)–mRNA, RNA–RNA, ncRNA–RNA, RBP–ncRNA, miRNA–mRNA, and miRNA–ncRNA, are formulated based on CLIP-seq and degradome sequencing data specific to plants. The interplay between hub genes and miRNAs was predicted using the ENCORI database. Subsequently, the mRNA–miRNA data were screened. Three or more databases supported the establishment of this predicted mRNA–miRNA interplay network.

The CHIPBase database [17] (version 3.0) (https://rna.sysu.edu.cn/chipbase/) recognizes hundreds of thousands of binding motifs and binding sites from DNA-binding protein ChIP-seq information, while predicting transcriptional regulatory relationships between millions of genes and transcription factors (TFs). Meanwhile, TFs that bind to hub genes were identified using the CHIPBase database. Subsequently, to build an mRNA–TF interplay network, we filtered mRNA–TF data supported by four or more databases.

Additionally, the Comparative Toxicogenomics Database (CTD) [18] (http://ctdbase.org/) was used to predict underlying drugs or small-molecule compounds that interact with hub genes. Subsequently, mRNA–drug data supported by four or more databases were filtered to build an mRNA–drug interplay network.

The RCSB PDB database (https://www.rcsb.org) was used to extract crystal structures for the hub genes in pdb format [19]. Core small-molecule structures were obtained in mol2 format from the TCMSP database. For crystal structures downloaded from the RCSB PDB database, we used PyMOL (https://pymol.org/2/) to remove ligands, select chains, and extract them into PDB format. Protein molecule active docking sites were predicted using Discovery Studio 2.5.5 (Dassault Systèmes BIOVIA, San Diego, CA, USA). For the extracted protein structures mentioned above, the AutoDockTool-1.5.7 program (The Scripps Research Institute, La Jolla, CA, USA) [20] was utilized to eliminate water molecules, supplement hydrogen atoms, compute charges, and export files in both pdb and pdbqt formats. For the core small-molecule structures in mol2 format, we used the AutoDockTool-1.5.7 program to add hydrogen atoms, calculate charges, and determine rotatable bonds for the ligands. These formats were then exported as pdbqt format files.

Autodock Vina (1.2.0, The Scripps Research Institute, La Jolla, CA, USA) [21] was utilized for docking proteins and small molecules, generating nine complex conformations for each protein–molecule docking. The docking score Affinity in AutoDock Vina indicates the strength of binding: Affinity $>$–4 kcal/mol represents very poor or no binding, –7 kcal/mol Affinity –4 kcal/mol is defined as moderate binding, and Affinity –7 kcal/mol is considered strong binding. Conformations with a binding free energy of –5 kcal/mol were filtered through energy ranking, with the lowest binding free energy chosen as the target conformation for analysis. For all target conformations, PLIP 2021 (Freie Universität Berlin, Berlin, Germany) [22] was employed to analyze the intermolecular non-covalent interplay, such as hydrophobic interactions, hydrogen bonds, and the cation–$\pi$ interplay. All interplay was visualized using PyMOL (https://pymol.org/) and saved in png format.

Sprague-Dawley (SD) rats aged 6–8 weeks were randomly divided into the normal control (NC), ORNJ, and ORNJ + LbGP groups. To establish the ORNJ animal model, rats from the ORNJ and ORNJ + LbGP groups were irradiated once daily at doses of 7 Gy per fraction for a total of five fractions [23]. The NC rats did not receive radiation treatment. The ORNJ + LbGP rats were also administered oral LbGP at a dose of 400 mg/kg/d, whereas the NC and ORNJ rats were administered the same concentration of normal saline. After 3 weeks of radiotherapy, the first molar of the mandible was extracted under isoflurane local anesthesia. The rats were euthanized with excessive anesthesia via inhalation using isoflurane (2% concentration) delivered at 0.41 mL/min with a fresh gas flow rate of 4 L/min. The mandibular tissues were extracted, and jaw bone marrow mesenchymal stem cells (JBMMSCs) were isolated and cultured for 3–5 generations for the follow-up experiments. The identity of these cells was confirmed through morphology, surface marker expression, and differentiation ability. All the cultures were subjected to routine tests, and mycoplasma contamination was found to be negative. RT-qPCR was performed to evaluate the expression of the hub genes in the JBMMSCs. First, cDNA synthesis was performed under the following conditions: reverse transcription at 42 °C for 60 minutes, followed by enzyme inactivation at 70 °C for 5 minutes. The used primers are listed in Supplementary Table 3. The RT-qPCR reaction mixture consisted of 10 µL of 2$\times$ SYBR Green Master Mix, 0.5 µL each of the forward and reverse primers (10 µM), 2 µL of cDNA template, and nuclease-free water to a final volume of 20 µL. SYBR Green chemistry was employed to detect transcripts. The reactions were performed on a QuantStudio3 Real-Time PCR System (Applied Biosystems). The thermal cycling conditions used for qPCR were as follows: initial denaturation at 95 °C for 10 minutes, followed by 40 cycles of denaturation at 95 °C for 15 seconds, and annealing/extension at 60 °C for 60 seconds. GAPDH was used as the endogenous control gene against which gene expression levels were normalized. All reactions were performed in triplicate to ensure reproducibility and reliability of the data. The protein levels of TNF-$\alpha$, IL-1$\beta$, and IL-6 were determined by ELISA, using the following kits obtained from Elabscience (Wuhan, Hubei, China): rat TNF (E-EL-R2856c), rat IL-1$\beta$ (E-EL-R0012c), and rat IL-6 (E-EL-R0015c). The mandibular tissues of rats were fixed, paraffin-embedded, and sliced for preservation. The protein expression levels of hub genes in these tissues were compared by Ysubg immunohistochemical staining using the following Huabio antibodies: EGFR (1:100, 0407-21), TGFB1 (1:200, HA721143), IL-6 (1:100, R1412-2), IL-1$\beta$ (1:100, HA601036), and TNF (1:200, ER65189). Antibodies for the analysis of P53 (1:500, YT3528) and ESR1 (1:500, PT0632R) expression were purchased from Immunoway (Hangzhou, Zhejiang, China).

Statistical analysis was performed using R software (v4.3.0, R Foundation for Statistical Computing, Vienna, Austria) for bioinformatics and GraphPad Prism (v9.0.0.121, GraphPad Software (Dotmatics) Boston, MA, USA) for experimental data. Enrichment analyses (GO and KEGG) were conducted using the clusterProfiler package 4.2 (Guangzhou Medical University, Guangzhou, Guangdong, China), with significance defined as p.adjust 0.05 and FDR 0.05 (Benjamini–Hochberg correction).

Data from at least three independent experiments were expressed as the mean $\pm$ standard deviation. For comparisons of multiple groups (NC, ORNJ, and ORNJ + LbGP), one-way analysis of variance (ANOVA) was applied, followed by Tukey’s post-hoc test. Differences with a p-value 0.05 were considered statistically significant. This statistical approach was used to analyze data obtained by RT-qPCR, ELISA, and immunohistochemistry.

The flowchart for this study is shown in Fig. 1.

Fig. 1.

Technology roadmap. ORNJ, osteoradionecrosis of the jaw; FRGs, ferroptosis-related genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; TF, transcription factors.

A total of 188 chemical components from Lycium barbarum were extracted from the TCMSP database. Among these, 45 active compounds were selected through screening for values of OB $>$30% and DL $>$0.18. Detailed data for these compounds are shown in Table 1. Protein targets for the active compounds in Lycium barbarum were predicted using the TCMSP database. The targets were then translated into gene names using UniProtKB, yielding a total of 205 targets for the Lycium barbarum constituents. The component–target network diagram is shown in Fig. 2A, and specific target information is provided in Supplementary Table 4.

Fig. 2.

Target screening for Lycium barbarum and ORNJ. (A) Network diagram of predicted targets of active chemical constituents in Lycium barbarum. (B) Venn diagram of the predicted targets of active components in Lycium barbarum and ORNJ-related genes. (C) Active chemical components, ORNJ-related genes, and the KEGG pathway network map for Lycium barbarum. ORNJ, osteoradionecrosis of the jaw; ORNJRGs, osteoradionecrosis of the jaw related genes; KEGG, Kyoto Encyclopedia of Genes and Genomes. In panels (A) and (C), the green circles represent mRNAs, purple triangles are chemical components, blue hexagons are ORNJ, and orange quadrangles are the KEGG pathways.

The intersection of 205 potential target points for Lycium barbarum components and 35 disease target points was plotted in a Venn diagram (Fig. 2B). This yielded 11 potential Lycium barbarum, in traditional Chinese medicine formulation, treatment targets for ORNJ: RUNX2 (Runt-Related Transcription Factor 2), TGFB1 (Transforming Growth Factor Beta 1), TP53 (Tumor Protein p53), TNF (Tumor Necrosis Factor), IL-6 (Interleukin-6), VEGFA (Vascular Endothelial Growth Factor A), ESR1 (Estrogen Receptor 1), SPP1 (Secreted Phosphoprotein 1), IL-1$\beta$ (Interleukin-1 Beta), CXCL8 (C-X-C Motif Chemokine Ligand 8), and EGFR (Epidermal Growth Factor Receptor). A component-target-pathway network was constructed (Fig. 2C).

Seven intersection genes (TP53, EGFR, TGFB1, IL-6, IL1B, TNF, ESR1) were obtained from the intersection of FRGs and Lycium barbarum treatment ORNJRGs. These intersection genes were deemed as hub genes in the Venn diagram (Fig. 3A). Protein–protein interactions were examined in terms of the seven hub genes, resulting in the creation of a protein–protein interaction network (Fig. 3B), visualized through Cytoscape (Fig. 3C). Furthermore, for each of the seven hub genes, an interaction network of functionally related genes was predicted and built using the GeneMANIA website (Fig. 3D). This allowed the physical interactions, shared protein domains, gene interactions, and other relevant information for the hub genes to be observed.

Fig. 3.

Protein–protein interaction network. (A) Intersection Venn diagram of ferroptosis-related genes and 11 genes related to Lycium barbarum treatment of radiation-induced osteonecrosis of the jaw. (B) Protein interaction network diagram obtained from the STRING database. (C) Protein interaction network drawn by Cytoscape. (D) Interaction network of functionally similar genes of hub genes on the GeneMANIA website. In this network, black circles with white slash lines represent input hub genes, while other black circles without white slash lines represent predicted functionally similar genes.

The BP, MF, CC (Cellular Component), and biological pathways were analyzed for the seven hub genes (TP53, EGFR, TGFB1, IL-6, IL1B, TNF, ESR1). The results of the enrichment analysis are presented as bar graphs (Fig. 4A–C). Furthermore, the correlations between hub genes, GO results (Fig. 4D,E), and KEGG (Fig. 4F) enrichment analyses are presented as a circular network diagram.

Fig. 4.

Functional enrichment analysis (GO) and pathway enrichment (KEGG) analysis of hub genes. (A) Bar graph of the results of BP functional enrichment analysis of hub genes. (B) Ring bar chart of MF functional enrichment analysis results of hub genes. (C) Bar chart of KEGG pathway enrichment analysis results of hub genes. (D) Network diagram of BP functional enrichment analysis results of hub genes. (E) Ring network diagram of MF functional enrichment analysis results of hub genes. (F) Network diagram of KEGG pathway enrichment analysis results of hub genes. GO, Gene Ontology; BP, biological process; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes. (D–F): the blue circles are used to denote the BP entries, orange circles for MF entries, green circles for KEGG entries, and red circles for genes.

As shown in Fig. 4A–F, the hub genes were mainly enriched in regulating the generation of miRNAs involved in gene silencing through microRNAs, small RNA generation involved in gene silencing through RNA, miRNA-mediated gene silencing, gene silencing through RNA, post-transcriptional gene silencing, post-transcriptional gene silencing through RNA, miRNA generation involved in gene silencing through miRNA, primary sncRNA processing, and the negative and positive regulation of the generation of miRNAs involved in gene silencing through miRNA. The regulation of other BPs included cytokine receptor binding, signaling receptor activator activity, cytokine activity, receptor ligand activity, RNA polymerase II common transcription initiation factor binding, common transcription initiation factor binding, basal transcription machinery binding, basal RNA polymerase II transcription machinery binding, and ATPase binding, as well as MFs such as protease binding. The hub genes were mainly enriched in malaria, inflammatory bowel disease, cancer proteoglycans, cytomegalovirus infection, rheumatoid arthritis, the Advanced Glycation End-products-Receptor for Advanced Glycation End-products (AGE-RAGE) signaling pathway in diabetic complications, Chagas disease, Amoebiasis, the MAPK signaling pathway, and antifolate resistance pathway (Table 2).

The mRNA–miRNA data identified miRNAs that interact with the seven hub genes (TP53, EGFR, TGFB1, IL-6, IL1B, TNF, and ESR1). A visual representation of the mRNA–miRNA interaction network was created using the Cytoscape software (Fig. 5A). Six hub genes (TP53, EGFR, TGFB1, IL-6, IL1B, ESR1) and 73 miRNA molecules were included in this network, yielding 129 pairs of interactions between mRNA and miRNA. Detailed information on these interactions is shown in Supplementary Table 5.

Fig. 5.

Construction of mRNA–miRNA, mRNA–TF, and mRNA–drug interaction networks. (A–C) mRNA–miRNA (A), mRNA–TF (B), mRNA–drug (C) interaction networks of hub genes. Green squares are mRNAs; orange squares are miRNAs; blue squares are TFs; pink squares are small-molecule drugs. TFs, transcription factors.

We next searched the CHIPBase database for TFs that interact with the hub genes. Interaction information was extracted from the two databases and intersected with the seven hub genes, thereby demonstrating interactions between the seven hub genes and 81 TFs, as visualized using Cytoscape (Fig. 5B). Details of the mRNA–TF interactions are shown in Supplementary Table 6.

Possible medications or chemical compounds for the seven hub genes were identified using the CTD database. The mRNA–drug interaction network (Fig. 5C) revealed 73 possible pharmaceuticals or chemical substances that interact with six hub genes (TP53, TGFB1, IL-6, IL1B, TNF, ESR1) (Supplementary Table 7).

Based on the disease target–Lycium active ingredient–pathway network shown in Fig. 2C, three active ingredients related to the seven hub genes were identified: $\beta$-sitosterol (MOL000358), glycitein (MOL008400), and quercetin (MOL000098). By using Autodock Vina for the molecular docking of hub genes with core active ingredients and then selecting optimal conformations, eight protein–active ingredient docking structures were obtained, all with binding free energies –5 kcal/mol (Table 3). Molecular interactions for the eight binding conformations were subsequently analyzed using PLIP and visualized with PyMOL (Fig. 6A–H). From the docking results and a molecular mechanics perspective, TGFB1 was capable of forming hydrogen bonds with the Lycium $\beta$-sitosterol component, generating hydrophobic forces (Fig. 6A). ESR1 was able to form hydrogen bonds with the Lycium glycitein component, generating hydrophobic forces (Fig. 6B). EGFR was capable of forming hydrogen bonds with the Lycium quercetin component, generating hydrophobic forces (Fig. 6C). TNF could form hydrogen bonds with the Lycium quercetin component, generating hydrophobic forces (Fig. 6D). IL-6 was able to form hydrogen bonds with the Lycium quercetin component, generating hydrophobic forces (Fig. 6E). TP53 was capable of forming hydrogen bonds with the Lycium quercetin component, generating hydrophobic forces (Fig. 6F). IL-1$\beta$ was able to form hydrogen bonds with Lycium quercetin component, generating hydrophobic forces (Fig. 6G). Finally, TGFB1 could form hydrogen bonds with the Lycium quercetin component, generating hydrophobic forces (Fig. 6H).

Fig. 6.

Docking analysis of hub genes with core active components. (A) Mode of binding of Transforming Growth Factor Beta 1 (TGFB1) to $\beta$-sitosterol components in Lycium barbarum. (B) The binding mode of Estrogen Receptor 1 (ESR1) to the glycitein component in Lycium barbarum. (C) The binding mode of Epidermal Growth Factor Receptor (EGFR) to the quercetin component in Lycium barbarum. (D) Pattern of TNF binding to quercetin components in Lycium barbarum. (E) Binding mode of IL-6 to quercetin components in Lycium barbarum. (F) The binding mode of TP53 to quercetin components in Lycium barbarum. (G) Mode of binding of IL-1$\beta$ to quercetin components in Lycium barbarum. (H) Mode of binding of TGFB1 to quercetin components in Lycium barbarum. Interaction force line: blue dashed line indicates hydrogen bonding; black dashed line indicates hydrophobic interaction. Bright blue indicates the protein structure of the hub genes, yellow and orange highlight the active components, and dark blue shows the amino acid residues that interact with the active components.

We next investigated the therapeutic effects of LbGP on ORNJ. The RT-qPCR data showed that the mRNA levels of IL1B and IL-6 were significantly reduced in rat JBMMSCs treated with LbGP compared to the ORNJ group; however, the levels of these mRNAs remained higher than in the NC group. The mRNA levels of TP53 and TNF were also lower in the LbGP group, whereas those of EGFR and TGFB1 increased. The mRNA level of ESR1 showed no significant change (Fig. 7A). The ELISA results exhibited significantly decreased protein levels of TNF-$\alpha$, IL-1$\beta$, and IL-6 in the cell culture supernatant of the LbGP treatment group compared to the ORNJ group (Fig. 7B). Immunohistochemical analysis also showed significantly reduced positivity for P53, TGF$\beta$1, IL-6, and TNF in the Lycium barbarum peptide extract group compared to the ORNJ group. Positive staining for EGF was increased, while there were no obvious changes in positive staining for IL-1$\beta$ and ESR1 (Fig. 7C). Therefore, LbGP may alleviate ORNJ-induced tissue injury in rats by regulating the mRNA expression levels of TP53, EGFR, IL-6, and TNF and of the corresponding translated proteins.

Fig. 7.

RT-qPCR validation, ELISA validation, and immunohistochemical validation of paraffin-embedded sections of rat mandible tissue. (A) Relative mRNA expressions of hub genes in the normal, ORNJ, and treatment groups. (B) Protein expressions of tumor necrosis factor-alpha (TNF-$\alpha$), IL-1$\beta$, and IL-6. (C) Positive area of hub gene-related proteins. Bar = 30 µm. Data are presented as the mean $\pm$ standard deviation of three independent experiments (n = 3). Statistical significance was evaluated using ANOVA. *p 0.05; **p 0.01; ***p 0.001; ****p 0.0001.