Prediction and Verification of Epimedium Flavonoids With Different Glycosylation Numbers in Reversing Glucocorticoid-Induced Bone Formation Inhibition by Molecular Docking and Zebrafish

Glucocorticoids have been detected in environmental waters, and their biological potency has raised concerns on their impact on aquatic vertebrates especially fish. Numerous researches showed that the continuous and direct contact of aquatic vertebrates with glucocorticoid contaminants in environmental water will cause bone formation inhibition. The aim of this study is to predict and verify the effect of icaritin (IT), icariin (ICA), and baohuside-I (BHG-I) in reversing glucocorticoid-induced bone formation inhibition (GIBFI) by molecular docking and zebrafish model. We contrasted their activity in reversing GIBFI from their affinity to bone metabolism proteins (OPG, RANKL, BMP-2, BMP-4, Runx-2) by molecular docking. Subsequently, zebrafish model was adopted to evaluate their reverse effects on GIBFI. Alizarin red staining coupled with image quantification were performed to evaluate the effects of ICA, IT, and BHG-I on skeleton stained area (SSA) and cumulative optical density (COD). Inductively coupled plasma-mass spectrometry was applied to determine the contents of bone mineral elements (CBME, Mg, K, Ca, Fe, Zn) in zebrafish bones. Docking results showed the receptors (BMP-2, BMP-4, and Runx2) all combined well to ICA, while BHG-I bound well to OPG, the affinity between IT and the above targets were the weakest. Fortunately, IT, ICA, and BHG-I significantly increased the SSA, COD, and the contents of Ca compared with the model group (p < 0.05) in the order of IT>ICA>BHG-I. In conclusion, the glycosyl groups increased the H-bond affinity between flavonoids and target sites, which weakened bone formation. IT, BHG-I, and ICA all alleviated GIBFI, but their intensity order was IT>ICA>BHG-I.


INTRODUCTION
Glucocorticoids are garnering research interests in recent years Hidasi et al., 2017;Willi et al., 2018) as their potential environmental impacts are increasingly being recognized. They are frequently prescribed in medicine due to their anti-inflammatory, anti-allergic, and immunosuppressive properties (Willi et al., 2019). These compounds that originate from a wide range of medical applications in humans and domestic animals may pose an ecotoxicological risk for aquatic organisms in contaminated waters (McNeil et al., 2016;Neale et al., 2020).
The main source of these compounds entering the aquatic environment is the excretion of pharmaceutical residues used by either humans or animals in their free or conjugate forms (Goh et al., 2016). Existing treatment processes in wastewater treatment plants (WWTPs) are not entirely effective in removing these contaminants (Pedrouzo et al., 2009;Zwart et al., 2020). This can lead to pollutants transferring to the rivers and oceans via effluent of WWTPs (Cavallin et al., 2021). A previous study reported that total glucocorticoid levels had been detected at >0.5-52 ng/L in receiving river waters and up to 390 ng/L at discharged sites in Beijing region, China. Glucocorticoids and their metabolites totaling about 47-96 ng/L have been reported in Swiss river waters. In France, pharmaceutical industrial waste effluent detected with dexamethasone and prednisone up to 23,000 and 300 ng/L, respectively, have been discharged into a river where its downstream water was sampled over multiple timepoints and estimated to contain 1-2,900 ng/L of dexamethasone and 50-1,260 ng/L of 6 α-methylprednisolone .
The continuous and direct contact of aquatic vertebrates with glucocorticoid contaminants in environmental water is a valid concern especially for fish (Leusch et al., 2017;Allijn et al., 2018). It has been shown that fishes that have long term exposure to hydrocortisone experience adverse effects in their locomotion, impeding aggressive behavior, a change in their immune response and possible modifications to their sexual behavior (Øverli et al., 2002;DiBattista et al., 2005). Also, it is demonstrated that in embryos, the glucocorticoid system is active after around 48 h post fertilization (Wilson et al., 2013). As it plays an important role in regulating hatching (Nesan et al., 2012), swimming activity (Castillo-Ramírez et al., 2019), circadian rhythm (Zhao et al., 2018), musculoskeletal development, growth and stress response (Apaydin et al., 2020), and exposure to exogenous glucocorticoids can be especially disruptive during early development. Relevant research reported prednisone induced decreased bone mineral density and bone loss (Schmid et al., 2020a;Schmid et al., 2020b). Furthermore, the long-term use of glucocorticoids will induce osteoporosis mediated by glucocorticoid receptor-dependent and -independent pathways (Zhao et al., 2016;Zhao et al., 2018). In our present study, we focused on the effects of glucocorticoids in aquatic environment on zebrafish bone formation and tried to find therapeutic drugs from traditional Chinese medicine.
Epimedii Folium (EF) has been used to repair fractures, bone and joint damage, and gonad dysfunctions for thousands of years . The major active flavonoids in EF were icariin (ICA), baohuside-I (BHG-I), BHG-II, and epimedin (A, B, C) (Li et al., 2015). There is little Icaritin (IT) in Epimedium herbs and extractions, while EF-flavonoids can be biologically converted into IT in vivo after hydrolysis of glycosidic bonds. IT was the structural parent nucleus of ICA and BHG-I, and the glycosylation number of IT, BHG-I, and ICA were 0, 1, 2, respectively ( Figure 1). Previous studies revealed that higher concentrations of IT were discovered in rat plasma, urine, and feces after oral administration of ICA or BHG-I (Sun et al., 2014;Sun et al., 2016). Pharmacological studies confirmed that ICA increased bone mass and prevent glucocorticoid-induced apoptosis in osteocytes (Feng et al., 2013;Wang et al., 2018). In addition, BHG-I activated EGFR-Akt-Nrf2 signaling and protected osteoblasts from Dexamethasone . Also, relevant studies have reported that IT exerted positive effects on human osteoblast proliferation and osteogenic function (Lim et al., 2017). Up to now, there are barely documents reported the activity difference of IT, BHG-I, and ICA in anti-bone formation inhibition.
The skeletons were in constant renewal, and this steady state of bone was mainly maintained by osteoclasts and osteoblasts. Considerable evidence indicated that the osteoprotegerin (OPG)/receptor activator of nuclear factor κB ligand (RANKL)/RANK system (OPG/RANKL/RANK) was relevant to osteoclast proliferation and differentiation (Liu and Zhang, 2015;Kushlinskii et al., 2017). Bone morphogenic protein (BMP), a member of the transforming growth factor-beta (TGF-b) super family (Aslani et al., 2019), was related to the differentiation of osteoblasts from mesenchymal stem cells or osteoprogenitor cells in bone marrow (Garg et al., 2017;Liu L. et al., 2018;Yang et al., 2018). It was demonstrated that pharmacologic blockade of RANKL in the OPG/RANKL/ RANK system was an effective treatment for osteoporosis. Similarly, antagonists, such as noggin (Hashimi, 2019), chordin (Tekari et al., 2017), and gremlin-1, also effectively restricted the differentiation of osteoblasts through binding to BMP-2 (Kišonaitė et al., 2016). All these results demonstrated that binding to bone metabolism-related protein targets was one important way for drugs to promote or inhibit bone turnover. The structural difference in the glycosyl groups of IT, BHG-I, and IT caused their discrepancy in affinity with related targets, which was one significant factor for their bioactivity in antibone formation inhibition.
In this paper, we primarily contrasted the activity difference of IT, BHG-I, and ICA on anti-bone formation inhibition and clarified the underlying mechanism. Molecular docking was carried out to elaborate the affinity between three flavonoid homologs (IT, BHG-I, and ICA) and protein targets. Immediately, the zebrafish models were used to verify the molecular docking results and evaluates their effects in antibone formation inhibition according to skeleton staining, quantification of the area of mineralized bones and cumulative optical density, and determination of bone mineral contents. Finally, the order of action intensity of these three flavonoids in anti-bone formation inhibition was obtained, and the mechanism of this difference was explained from the perspective of the interaction between chemical structure and bone metabolic protein targets.

Molecular Docking Simulation
The amino acid sequences of OPG (PDB ID: 3URF, CHAIN: Z, Length: 171 AA, Residues: 22-186), RANKL (PDB ID: 3URF CHAIN: A Length: 162 AA Residues: 162-317) and BMP-2 (PDB ID: 2QJB CHAIN: AB Length: 116 AA Residues 283-396) have structural data of Experimental Analytical proteins. BMP-4 and Runx2 amino acid sequences have no experimental analytical protein structure data, and 3D protein structure data was obtained by homology modeling. Data were calculated using Molecular Operating Environment 2015.10 software (Chemical Computing Group Inc., Quebec, Canada). Maestro was carried out to build the structure of ligands, and then they were transformed into 3D format. The docking module was used to dock the ligand, and the protein and docking scores were evaluated by Glide's scoring function.

Animals
Zebrafish farming was conducted according to the standard procedure (Varga et al., 2018). Adult zebrafish mate and spawn naturally under the condition of 10/14 h dark/light cycle. The zebrafish embryos and larvae were cultured in sixhole plates with blank E 3 medium solution (5.0 mM NaCl, 0.17 mM KCl, 0.33 mM CaCl 2 , and 0.33 mM MgSO 4 ) in Light incubator (PGX-150B) at 28.5°C and natural air circulation. Animal experiments were carried out in accordance with the Guidelines for Animal Experimentation of Jiangsu University (Zhenjiang, China), and the protocol was approved by the Animal Ethics Committee of this institution.

Treatment Regime
The zebrafish model by prednisolone (PNSL, 25 μM) was established according to our preliminary study (Jiang et al., 2018). Zebrafish larvae were randomly divided into seven groups (n = 15 larvae in each well with three replications) as follows: Control (CK, blank E 3 medium), DMSO (DMSO, 0. . From 6 to 10 DAF, larvae were fed with paramecia for 1 h every day. Following each feeding, the remaining paramecia were washed out, and the medium was replaced with blank E 3 medium-containing drugs or not (Jiang et al., 2018).

Behavior Analysis
After 4 days of drug deliveries, zebrafish locomotor activity was monitored by EthoVision behavior system. Prior to the start of tracking, the software needed to be calibrated. The video sampling rate was 25 frames per second (fps), based on the design recommendations (Wolter and Svoboda, 2020). First, under the "Trial List," one trial was selected. For the "Arena Settings," each well/arena was calibrated based on the diameter of the well. The diameter of the wells within the well plates used in this manuscript are 96 well plate (6.54 mm). For the "Detection Settings," dynamic subtraction was selected and the dark contrast and subject contour were adjusted to optimize tracking efficiency (Knafo and Wyart, 2018). Within the "Analysis Profiles," the selected dependent variables were distance moved, velocity, and time spent moving. These endpoints were based on the centerpoint activity of the larvae (Martella et al., 2016). The results were then exported to Excel and statistical analysis software suites. According to the methods described above, moving distance (MD), average speed (AS), travel frequency (TF), and hotspot are selected as the anti-GIBFI drug efficiency evaluation index in this model.

Whole-Mount Skeletal Staining
Alizarin red staining on zebrafish larvae was performed as previously described (Bruneel et al., 2015). Briefly, the zebrafish larvae were collected at 10 DAF, all the zebrafish were killed by 3-aminobenzoic acid ethyl ester methane sulfonate (MS-222, 100 mg/L) (Martini et al., 2021). After removal of MS-222 solution, zebrafish larvae were fixed in paraformaldehyde solution 4% in PBS (pH 7.4) and stained with Alizarin red staining (Vimalraj et al., 2020). Fresh prepared bleach with 1.5% H 2 O 2 and 1% KOH was added 12 h later after the staining solution was removed and washed., All samples were decolored with glycerol 1 h later, and the stained zebrafish was placed under a microscope to observe their stained bones (Lim et al., 2021). Images were acquired by microscope (Olympus IX71/IX81, Olympus Corporation, Japan). Digital images were analyzed using

Sample and Standard Solution Preparation for Inductively Coupled Plasma-Mass Spectrometry Detection
The mineral contents in zebrafish larvae were measured by inductively coupled plasma-mass spectrometry (ICP-MS) as reported in our previous study. The larvae were collected at 10 DAF (n = 15) and washed five times in double distilled water and then transferred to centrifuge tubes. Immediately, samples were digested with 70% HNO 3 in a microwave oven for 4 h. Determination of Ca, K, Mg, Zn, and Fe was carried out by the 7500cx ICP/MS system (Agilent Technologies, Santa Clara, CA, USA) equipped with a G3160B I-AS integrated autosampler. The standard solutions of Mg, K, Ca, Fe, and Zn with mass concentration at 10 mg/L were diluted into tuning solution (1 μg/ L) with 2% nitric acid. The series concentration of Ca standard solution was 0, 40, 80, 120, 160, 200 mg/L. The other series concentration of Mg, K, Zn, and Fe standard solution were 0, 10, 20, 40, 80, 160 μg/L.

Inductively Coupled Plasma-Mass Spectrometry Conditions
Agilent 7500 ICP-MS system was used for simultaneous determination of Ca, K, Mg, Zn, and Fe. The pressure of Ar and He were set to 700 and 40 KPa, respectively. Circulating water temperature was 20°C. The pressure was 230-400 KPa. The exhaust air volume was set to 5,000-7,000 L/min. Plasma power was 1.5 KW. Carrier gas flow rate was 0.9 L/min, and compensation gas flow rate was 0.25 L/min. Injection depth was 8 mm. Peristaltic pump speed was set to 0.1 r/s. Premix chamber temperature was 2°C.

Method Validation of Inductively Coupled Plasma-Mass Spectrometry
After the instrument was tuned, the series mixed standard working solution was injected into ICP-MS, and the signal response value of each element was measured. The mass concentration was taken as the abscissa, and the response signal value was taken as the ordinate. The correlation coefficient (r) of the standard curve were all higher than 0.9990. The detailed calibration curves, linear range, and the lower limit of quantitation (LLOQ) for mineralized elements are available in Table 1.

Statistical Analysis
All data were presented as the mean ± SD. Statistical analysis was performed using GraphPad Prism 5 (GraphPad software, USA). Differences were analyzed by one-way analysis of variance (Tukey, compare all pairs of columns). Differences with p < 0.05 were considered significant.

Molecular Docking Results
The important parameters obtained from molecular docking included docking scores, interacting residues, bonding energy, binding position and H-bonds interaction between receptors and ligands. The ICA-RANKL interaction had 4 H-bond residues as GLU292 (A), ASP174 (A), SER252 (A), and LYS205 (A) with a total binding energy of −12.8 kcal/mol. The interaction between BMP-2 and ICA also had 4 H-bond residues as GLU46 (A), GLY45 (A), CYS47 (A), and CYS79 (A) with a total binding energy of −11.5 kcal/mol. ICA bound to BMP-4 via H-bond interaction at LEU 386(B), pi-H interactions with TRP 325(B), and TRP 322(B) with a total binding energy of −3.6 kcal/mol. The ICA-RANK interacted at GLN122, LEU143, and GLN144 with a total binding energy of −3.7 kcal/mol. Therefore, the docking results indicated that ICA could bind to multiple protein targets including RANKL, BMP-2, BMP-4, RANK, and Runx2 with docking scores at −6.67224169, −6.69499111, −6.68158722, −6.81634808, and −6.87304401, respectively (Tables 2 and 3; Figure 2, Supplementary Figure S1).
Although BHG-I could interact with the above six protein targets, the combination between BHG-I and OPG via pi-H interactions with interacting residues TYR 71(Z) and VAL 72(Z), H-bond interaction with LYS 87(Z) had the maximum absolute docking score (−9.64769459, Table 2) and the lowest binding free energy (−3.5 kcal/mol, Table 3). Similarly, ligand IT was able to communicate with OPG, BMP-2/-4, and RANK, but their absolute docking score were relatively smaller than ICA and BHG-I, which showed that the combination of IT and the above targets were weak. The IT-OPG interaction had two residues as TYR71, VAL65 with a total binding energy of −1.5 kcal/mol. The IT connected with BMP-2 through GLU109 residues with a binding energy of −1.1 kcal/mol. The IT-RANK interaction had two residues as VAL91 and CYS92 with a total binding energy of −1.8 kcal/mol. The IT-BMP-4 interaction had four residues as GLU388, TRP325, TRP322, and TYR385 with a total binding energy of −7.3 kcal/mol. The IT connected with Runx2 through LEU168(A) residues with a binding energy of −5.7 kcal/mol. These results exhibited that the affinity between BHG-I and OPG was the strongest, and the combination of IT and above targets were all instability.

Effect of Icariin, Icaritin, and Baohuside-I on Bone Formation and Bone Mineralization
The skeleton stained area (SSA) and cumulative optical density (COD) were important indicators of bone mineralization. Alizarin Red staining was widely used to detect and quantify mineralized bones because it was capable of binding to calcium salts and showed red under the microscope (Figure 3). Compared with the DMSO group, the area of SSA and COD in Zebrafish larvae were decreased significantly under the treatment of 25 μM PNSL (MX, p < 0.01, Figure 4). These results indicated that PNSL reduced bone mineralization and inhibited osteogenic differentiation in zebrafish larvae. Compared with the MX group, the COD was significantly increased (p < 0.05, Figure 4A Figure 4B). These flavonoids showed a prospective concentration-dependent in reversing PNSL induced inhibition of bone formation.

Effect of Icariin, Icaritin and Baohuside on Zebrafish Behavior
The behavior analyzer-Ethovision-was used to track the movement of zebrafish. As shown in Figures 5 and 6, the moving speed (MS), moving distance (MD), and travel frequency (TF) of PNSL group was significantly lower than those of the DMSO, while for YTLSN and all three compounds, these parameters were close to the control. The results showed that different concentrations of IT, ICA, and BHG-I increased MD by 131.50%, 124.66%, and 107.38%, respectively, with partial significant (p-value < 0.001-0.01). TF for IT, ICA, and BHG-I were 169.39%, 137.66%, and 128.57%, respectively. Also, the hot plot revealed zebrafish activity degree (Supplementary Figure S2). These results with significant differences demonstrated that IT, ICA, and BHG-I could reverse GIBFI.

Icariin, Icaritin, and Baohuside Promote the Enrichment of the Elements Required for Bone Mineralization
Bone was composed of 69-80 wt% calcium phosphate and other trace elements. To confirm the therapeutic effect of ICA, IT, and BHG-I on osteogenesis, whole-body Ca, K, Mg, Zn, and Fe contents of larvae were measured by ICP-MS. Compared with the DMSO, the treatment of PNSL significantly decreased wholebody Ca, K, Mg, Zn, and Fe levels by 6.24, 1.67, 1.06, 0.11, and 0.75, respectively (p < 0.05). After the treatment of IT, ICA, and BHG-I, whole-body Ca, K, Mg, Zn, and Fe levels were significantly higher than PNSL group (p < 0.05, Supplementary Figure S3, Figures 5).

DISCUSSION
The OPG/RANKL/RANK signaling pathway plays a crucial role in regulating the bone remodeling process (Khosla, 2001). RANKL binds to RANK on osteoclasts to stimulate differentiation. OPG can also bind to RANKL to block this process and to control the remodeling process (Tyrovola,  2015; Kovács et al., 2019;Li et al., 2019). If ligands bound to OPG, blocking the binding of OPG to RANKL, which weakened the inhibitory effect of ligand to osteoclasts. On the other hand, if the ligands combined with RANKL which also inhibited the interaction of RANKL-RANK and attenuated the differentiation of osteoclasts, and finally showed bone resorption inhibition (Jiang et al., 2018). Antagonists, such as noggin (Chien et al., 2020), chordin (Huang et al., 2019) and gremlin-1 (Silvério de Barros et al., 2019), effectively restrict osteoblast differentiation through binding to BMP-2. Similarly, the combination between ligand and Runx2 blocked the function of Runx2 in promoting osteogenic differentiation. In summary, the ligands bind to BMP-2, BMP-4, or Runx2 inhibited the differentiation of osteoblasts, the stronger the binding, the more inhibition on the bone formation. In molecular docking, the lower bonding energy, the more stable binding between the ligands and the receptors. In addition, the larger the absolute value of docking scores, the more stable the combination. From the docking scores, the binding between BHG-I and OPG (−9.64769459) was the most stable compared with that of IT (−8.96370506) and ICA (Table 2). Meanwhile, the interaction between BHG-I and OPG had the lowest binding energy (−3.5 kcal/mol, Table 3, Figure 2). According to the docking results, the position of H-bonds was mostly located in glycosylation, which indicated that glycosylation affected the formation of H-bonds. For example, BHG-I bound to OPG via pi-H interactions with interacting residues TYR 71(Z), and VAL 72(Z), H-bond interaction with LYS 87(Z), while IT bound to OPG was pi-H interaction at VAL 65(Z) and H-bond interaction at TYR 71(Z). The glycosyl groups in BHG-I promoted the formation of H-bond and the interaction of ligand-receptor, which ultimately manifested as the weakening of anti-GIBFI. The formation of more H-bonds between ICA and multiple bone turnover targets was mainly attributed to the two glycosyl groups in its structure. Therefore, the binding between ICA and above protein targets was easier and more stable than the other two flavonoids because of the less glycosyl groups in IT and BHG-I (0 and 1, respectively). Based on molecular docking, the receptors (RANKL, BMP-2, BMP-4, and Runx2) all combined well to ICA, while BHG-I bound well to OPG, the combination of IT and the above targets was the weakest. It was predicted that the anti-GIBFI effect of IT was stronger than that of ICA and BHG-I. Zebrafish was an ideal animal model in vivo for studying bone deformations for its high skeletal and genetics similarity to human skeleton (Luo et al., 2016;Zhao et al., 2020). The zebrafish larvae contained the sufficient and necessary cells for both bone formation and resorption activity (Luo et al., 2016).
Alizarin red staining is a special stain used for bone staining. This stain can be used to determine the bone mineralization level based on the color, and the area of alizarin red staining can be used to determine efficacy of drugs affecting bone mineralization. The indicators observed in this study were the COD and SSA value, which directly reflects the differentiation and number of osteoblasts. Osteoblast differentiation and bone formation can be visualized in the zebrafish larvae by monitoring the changes of dyeing area. Our results indicated that PNSL exposure significantly inhibited osteogenic differentiation and bone mineralization in zebrafish larvae (p < 0.05, Figure 4). After the drug intervention, IT, ICA, and BHG-I all exerted positive effects on reversing PNSL-induced osteopenia in zebrafish (p < 0.05). According to the COD and SSA results, the reversal effect of IT and ICA on bone formation inhibition reached or even exceeded that of positive drugs, especially at 10 μM (Figure 4). This phenomenon was reproduced in the results of bone mineral element contents. As shown in Table 4 and Supplementary Figure S3, the levels of Mg, K, Ca, Fe, and Zn were achieved or even exceeded that of YTLSN group. BHG-I was weaker than IT and ICA in increasing SSA, COD, and bone mineralized elements. Behavioral changes of zebrafish have been linked to chemical exposure (Pitt et al., 2018;Nunes et al., 2020). Behavior analyzer-Ethovision XT-made it possible to examine numerous motor events and facilitates quantitative analysis of behavior (Zheng et al., 2021). The behavioral change in zebrafish is an important indicator to evaluate the anti-GIBFI effect of ICA, IT, and BHG-I. The MD and TF in PNSL group further supported that the construction of osteoporosis model was successful. MDs for IT, ICA, and BHG-I were 60.32%-131.50%, 71.59%-112.10%, 47.63%-97.38%, respectively. These results with significant differences demonstrated that ICA, IT, and BHG-I could improve dyskinesia of zebrafish to some extent. According to the above discussion, the docking results indicated the affinity difference was ICA>BHG-I>IT. The more stable the binding, the stronger the effect of inhibiting the formation of osteoblasts and the less the effect of anti-GIBFI. In consequence, the order of their intensity in reversing GIBFI should be IT>BHG-I>ICA. Theoretically, anti-GIBFI effect of BHG-I, a single-glucose-containing flavanone glycoside, should  have been stronger than ICA (the order should be IT>BHG-I>ICA), but the zebrafish experimental results were contrary to this (IT>ICA>BHG-I). ICA exhibited a more powerful potency against GIBFI than BHG-I with AMB values of 513,950.2, COD values of 171,177.8, and the significant increase in mineral element content. This was mainly because the combination between BHG-I and OPG was more stable, the inhibitory effect was stronger than ICA. Meanwhile, we noticed that the binding between ICA and RANKL (−6.67224169), BMP-2 (−6.69499111), BMP-4 (−6.81634808), Runx2 (−6.87304401) FIGURE 6 | The mechanism of affinity/bioactivity of EF-flavonoids and bone metabolism targets. ICA, icariin; BHG-I, baohuside I; IT, icaritin. In the figure, the orange square frames are in the symbol of the small molecule ligands ICA, BHG-I, and IT. The blue boxes represent the protein targets in OPG/RANKL/RANK system or BMP signaling pathway. The numbers in the orange square frames represent the affinity of the compound with bone metabolic protein targets. The number 1 represents the combination between compound and target is more stable. These three compounds are equivalent to antagonists of bone metabolic targets. The more stable the binding, the less the effect of anti-GIBFI. This graph clearly shows the combination between receptors and protein targets. EF, Epimedii Folium; GIBFI, glucocorticoidinduced bone formation inhibition. 0.85 ± 0.028*** 1.33 ± 0.10*** 36.79 ± 2.60*** 0.71 ± 0.086*** 0.12 ± 0.0073*** BHG-0.1 0.28 ± 0.033 0.20 ± 0.033 15.62 ± 1.08*** 0.058 ± 0.037 0.0094 ± 0.0062 BHG-1.0 0.39 ± 0.046 0.29 ± 0.076 18.88 ± 1.50*** 0.18 ± 0.037 0.041 ± 0.0032 BHG-10.0 0.47 ± 0.047* 0.53 ± 0.044*** 22.71 ± 1.64*** 0.30 ± 0.088 0.053 ± 0.0062* ICA-0.1 0.48 ± 0.075* 0.50 ± 0.074** 20.11 ± 1.00*** 0.14 ± 0.053 0.018 ± 0.0035 ICA-1.0 0.80 ± 0.18*** 1.12 ± 0.14*** 25.99 ± 0.94*** 0.28 ± 0.044 0.038 ± 0.0043 ICA-10.0 1.35 ± 0.049*** 1.48 ± 0.12*** 29.27 ± 1.67*** 0.58 ± 0.12*** 0.079 ± 0.0066*** IT-0.1 0.62 ± 0.077*** 0.66 ± 0.036*** 23.53 ± 0.80*** 0.52 ± 0.076*** 0.077 ± 0.011*** IT-1.0 1.21 ± 0.16*** 1.38 ± 0.11*** 30.69 ± 2.99*** 0.76 ± 0.044*** 0.11 ± 0.017*** IT-10.0 1.79 ± 0.067*** 1.87 ± 0.16*** 37.61 ± 2.33*** 1.01 ± 0.16*** 0.18 ± 0.025*** Frontiers in Environmental Science | www.frontiersin.org February 2022 | Volume 9 | Article 793527 was the most stable compared with that of BHG-I and IT (Table 2). However, the docking scores between BHG-I and the above protein targets were second only to icariin, with the value of −6.09861946, −5.87657118, −6.19196224, and −5.74639988, respectively, which meant BHG-I had good affinity to these bone metabolism proteins (RANKL, BMP-2, BMP-4, Runx-2). The more stable the binding, the less the effect of anti-GIBFI. This accounted for the reason why the anti-GIBFI effect of BHG-I was lower than ICA.
In conclusion, the ligands that bind to OPG, BMP-2, BMP-4, or Runx2 inhibited the differentiation of osteoblasts; the stronger the binding, the more inhibition on the bone formation. The receptors (BMP-2, BMP-4, and Runx2) all combined well to ICA, while BHG-I bound well to OPG; the combination of IT and the above targets was the weakest. IT, BHG-I, and ICA all alleviated bone formation inhibition induced by PNSL, but the order of their intensity in reversing GIBFI was IT>ICA>BHG-I; the most potential compound was found to be IT.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

ETHICS STATEMENT
The animal study was reviewed and approved by the Animal Ethics Committee of Jiangsu University.