Identification of the occurrence and potential mechanisms of heterotopic ossification associated with 17-beta-estradiol targeting MKX by bioinformatics analysis and cellular experiments

Bioinformatics analysis of MKX and Structure prediction of MKX

We used data from the STITCH (http://stitch.embl.de) database to identify MKX-related compounds and proteins (Szklarczyk et al., 2016). We used data from the SymMap (http://www.symmap.org/) database to identify MKX-related traditional Chinese medicine (TCM) drugs and symptoms (Wu et al., 2019). We also applied a genomic and proteomic discovery model named OncoBinder to identify potential regulators of MKX (Van Coillie et al., 2016). The tool was used to predict the structure of MKX (Gasteiger, 2003).

Enrichment analysis

The GO and KEGG pathway enrichment analysis was performed as described in the citation and visualized by R software (Yu et al., 2012). The clusterProfiler package is used for enrichment analysis, while the org.Hs.eg.db package is used for ID conversion. Gene Set Enrichment Analysis (GSEA) analysis was performed using the clusterProfiler package of R software [version 3.14.3] (Subramanian et al., 2005). The reference gene collection used by GSEA is c2.cp.v7.2.symbols.gmt, and the gene set database used is MSigDB collections. Conditions were considered significantly enriched when the false discovery rate was < 0.25 and p.adjust was < 0.05.

Molecular docking analysis

Molecular docking analysis is an important method for predicting compound and protein binding (Friesner et al., 2004; He et al., 2019; Mehta, Miszta & Filipek, 2021). The structures of Pro-Tyr and 17-beta-estradiol were searched in the PubChem database. After the three-dimensional structure of MKX was predicted by the ExPASy tool, the AutoDockTools software was used to determine the active pockets of potential sites for molecular docking. The structure has been verified and structure minimization has been performed before the docking study. Finally, the AutoDock Vina software was used to perform molecular docking to verify the targeting of MKX with Pro-Tyr and 17-beta-estradiol. The lowest binding energy values were considered to visualize the results using the PyMOL software.

Molecule dynamics

The complex system of each compound with the target egg obtained by docking was used as the initial structure to perform all-atom molecular dynamics simulations (MD) separately. The charges of the small molecules were calculated using the antechamber module and the Hartree-Fock (HF) SCF/6-31G* of the gaussian 09 software. Afterwards, small molecules and proteins were described using force fields (Wang et al., 2004; Maier et al., 2015). The LEaP module was used to add hydrogen atoms to the system, a truncated tetrahedral TIP3P solvent cartridge was added at a distance of 10 Å from the system, and Na+/Cl- was added to the system for balancing the system charge.

Molecular dynamics simulations were performed using AMBER 18 software (Lee et al., 2020). After the energy optimization of the system was completed, a 200 ps warming of the system at a fixed volume and constant warming rate was used to slowly increase the system temperature from 0 K to 298.15 K. Finally, these composite systems are simulated with 30 ns NPT (isothermal isobaric) tethering under periodic boundary conditions, respectively.For the simulations, the truncation distance of the non-bond was set to 8 Å. The Particle mesh Ewald (PME) method was used to calculate the long-range electrostatic interaction, the SHAKE method was used to limit the bond length of hydrogen atoms, and the Langevin algorithm was used for temperature control (Larini, Mannella & Leporini, 2007; Mehta, Miszta & Filipek, 2021). Traces were saved at 10 ps intervals for subsequent analysis.

MMGBSA free energy calculation

The binding free energy between protein and ligand for all systems was calculated by the MM/GBSA method (Genheden & Ryde, 2015; Chen et al., 2020). The MD trajectory of 25–30 ns was used as the calculation in this study. And the specific equation was as follows:

$$\begin{array}{rl}\mathrm{\Delta }{G}_{bind}& =\mathrm{\Delta }{\mathrm{G}}_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{e}\mathrm{x}}-(\mathrm{\Delta }{\mathrm{G}}_{\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{e}\mathrm{p}\mathrm{t}\mathrm{o}\mathrm{r}}+\mathrm{\Delta }{\mathrm{G}}_{\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{a}\mathrm{n}\mathrm{d}})\\ & =\mathrm{\Delta }{\mathrm{E}}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}}+\mathrm{\Delta }{\mathrm{E}}_{\mathrm{V}\mathrm{D}\mathrm{W}}+\mathrm{\Delta }{\mathrm{E}}_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}}+\mathrm{\Delta }{\mathrm{G}}_{\mathrm{G}\mathrm{B}}+\mathrm{\Delta }{\mathrm{G}}_{\mathrm{S}\mathrm{A}}\end{array}$$

$\mathrm{\Delta }{\mathrm{E}}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}}$ represents internal energy, $\mathrm{\Delta }{\mathrm{E}}_{\mathrm{V}\mathrm{D}\mathrm{W}}$ represents van der Waals interaction, and $\mathrm{\Delta }{\mathrm{E}}_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}}$ represents electrostatic interaction. $\mathrm{\Delta }{\mathrm{G}}_{\mathrm{G}\mathrm{B}}$ and $\mathrm{\Delta }{\mathrm{G}}_{\mathrm{G}\mathrm{A}}$ were collectively referred to as the solvation free energy. Where GGB was the polar solvation free energy and GSA was the nonpolar solvation free energy. The GB model developed by Nguyen was used to calculate $\mathrm{\Delta }{\mathrm{G}}_{\mathrm{G}\mathrm{B}}$ (Nguyen, Roe & Simmerling, 2013).

Experimental cells and culture

The cells used in this study were human, primary bone marrow-derived mesenchymal stem cells (BM-MSCs, Normal, Human, ATCC PCS-500-012), purchased from the American Type Culture Collection (ATCC). BM-MSCs were used, which was supplemented with 7% Fetal Bovine Serum (FBS), 15 ng/mL insulin-like growth factor 1 (rhIGF-1), 125 pg/mL fibroblast growth factor-basic (rhFGF), 2.4 mM L-alanyl-L-glutamine, and 1% penicillin-streptomycin solution (P/S). All cells were cultured at 37 °C in a cell culture incubator with 5% CO₂, and the experiment was completed with five passages of cells. Passaged BM-MSCs (n = 1,000) were placed in each well of a 96-well plate, divided into three groups, and were exposed to different doses of 17-beta-estradiol (0, 0.001, 0.01 nmol/L). Accordingly, the groups were categorized as control, 0.001 nmol/L 17-beta-estradiol-treated, and 0.01 nmol/L 17-beta-estradiol-treated.

Real-time quantitative reverse transcription PCR

Real-time quantitative reverse tanscription polymerase chain reaction (RT-qPCR) was performed to determine the mRNA expression of runt-related transcription factor 2 (RUNX2). RT-PCR was performed to quantify the expression of genes and cytokines associated with RUNX2 at 1, 3, and 5 d after the induction. Total RNA was extracted from the BM-MSCs using the total RNA extraction reagent (TRIzol) and purified using the QIAquick PCR purification kit. The first-strand cDNA synthesis kit and random hexamer primers were used for the reverse transcription reaction to obtain the cDNA using DNAse I. The quantitative RT-PCR was performed using a fluorescent quantitative PCR kit (SYBR Green PCR Master Mix). β-actin was selected as the internal reference (Table S1). PCR cycling condition was adjusted as follows: 95 °C for 2 min, 94 °C for 30 s, 55 °C for 30 s, and then at 72 °C for 1 min with 40 cycles performed. 2-DDCt method was used to obtain fold expression results.

Cell viability assay

Cell viability and proliferation were measured using the Cell Proliferation Kit I (MTT; Roche, Basel, Switzerland) according to the manufacturer’s instructions and the method described in previous studies (Jiang et al., 2020; Zhao et al., 2021, p. 5). First, the transfected cells were stabilized by trypsin digestion. Following this, the cells were centrifuged, collected, resuspended in a single-cell suspension, and seeded in 96-well plates at the density of 1 × 10⁴ cells/well. The cell cultures were maintained. A multifunctional microplate reader was used to measure optical density (OD) values. Cell growth curves were plotted using OD as the vertical axis and time as the horizontal axis.

Statistics and analysis

Statistical and visual analyses were performed using R software [version 3.6.0]. The measurement data were expressed as x ± s and compared using the t-test; the count data were expressed as percentages and compared using the χ2 test. p < 0.05 was considered statistically significant.

Construction of a regulatory network for MKX, related genes and compounds based on the STITCH database

We used data from the STITCH database (http://stitch.embl.de) to identify MKX-related compounds and proteins and constructed MKX-related gene and compound regulatory networks (Fig. 1A). Compounds related to MKX are Pro-Try and arsenite. Proteins related to MKX include SCXA, PROX2, suppressor interacting 3 A/histone deacetylase 1 (SIN3A/HDAC1), PBX1, MYOD1, TCF15, and PBX4. The regulatory network has nine nodes and 13 edges. Based on these proteins, we completed the GO functional enrichment analysis. The signaling pathways in which the MKX-related proteins were found to be enriched are as Fig. 1B. These results suggest that MKX-related protein functions are enriched in biological processes of biological development and transcriptional regulation.

MKX-related symptom network constructed using SymMap

The SymMap database integrates TCM and MM data on drug treatments and their molecular mechanisms, as shown in Fig. 2. Herbs, ingredients of active compounds, and TCM and MM symptoms are included in this network. The compound 17-beta-estradiol, which can be derived from three herbs (Moschus, Cervi cornu Pantotrichum, and Armeniacae Semen Amarum), was found to be capable of targeting and binding to MKX. MKX was found to be potentially associated with nine diseases, which in turn were associated with 42 symptoms (Table 1). In summary, the network demonstrates both TCM and MM symptoms associated with MKX.

Molecular docking validation of the compounds Pro-Tyr and 17-beta-estradiol for target binding of MKX

The predicted three-dimensional structure of MKX based on the ExPASy tool is shown in Fig. 3A. The structures of Pro-Tyr and 17-beta-estradiol that were obtained from the PubChem database are shown in Fig. 3B. Molecular docking using the AutoDock Vina software revealed that Pro-Tyr and 17-beta-estradiol target binding of MKX were at different loci (Figs. 3C and 3D).

Stability of complexes and molecular dynamics analysis

Pro-Tyr as well as Estradiol fluctuated less during the kinetic simulations, indicating that they were both stably bound within the active pocket (Figs. 4A, 4B). RMSF showed that Pro-Tyr could better reduce the flexibility of the protein, which implied that Pro-Tyr was more stable with the protein system (Fig. 4C). As shown in Figs. 5A and 5B, both molecules were maintained at the same binding site after MD simulation analysis. The small molecule seems to be closer to the protein after MD simulation analysis, which indicates that both molecules bind stably to the target protein. The results of MM/GBSA calculations showed that Estradiol/protein binding free energy was −17.57 ± 2.97 kcal/mol and Pro-Tyr/Protein binding free energy was −13.79 ± 2.68 kcal/mol (Table 2). The binding energy of Estradiol/protein is mainly contributed by van der Waals interaction, and that of Pro-Tyr/Protein is mainly contributed by electrostatic interaction.

Identification and functional analysis of proteins that may regulate MKX

Analysis of the OncoBinder database revealed a significant regulatory correlation between MKX and a series of proteins (TBP, HSP90A, SIN3A, TBL1X, TBL1XR1, HSP90AA1, HDAC1, ELAVL1, GTF2B, and GTF2A1) (Fig. 6A). We performed the GO and KEGG pathway enrichment analysis to identify the functions of these proteins (Fig. 6B). The GO enrichment analysis revealed that the functions of these proteins were significantly enriched in protein deacylation, protein deacetylation, histone deacetylation, transcriptional repressor complex, histone deacetylase complex, transcription factor complex, repressing transcription factor binding, RNA polymerase II general transcription initiation factor activity, and RNA polymerase II repressing transcription factor binding. The KEGG enrichment analysis revealed that these proteins were significantly enriched in viral carcinogenesis, basal transcription factors, and the IL-17 signaling pathway. The results suggest that MKX and the proteins significantly associated with it play important roles in protein acetylation and transcriptional regulation.

By combining the results presented in Figs. 1–3, it may be stated that the compounds Pro-Tyr and 17-beta-estradiol may bind MKX and affect its interaction with SIN3A and HDAC1 (Fig. 7A). The integrated GSEA based on the MSigDB collections gene set database revealed that the low expression of MKX was significantly associated with osteoblast differentiation, signaling by NOTCH4 and Hedgehog, and the Hedgehog signaling pathway enrichment demonstrated significant correlation (Fig. 7B). We evaluated the proliferation of BM-MSCs after treatment with different doses of 17-beta-estradiol and observed an increase in the proliferation rate (Fig. 7C). This could be attributed to the upregulation of the HDAC1-SIN3A/bone morphogenetic protein (BMP) pathway-related RUNX2 mRNA expression by 17-beta-estradiol (Fig. 7D). The results of our study suggest that 17-beta-estradiol may upregulate the HDAC1-SIN3A/BMP pathway.