Thirty-one 8-week-old C57BL/6Jnifdc (C57BL/6) mice and thirty 8-week-old BALB/cCrSlcNifdc (BALB/c) mice (Experimental Animal Center of Guangzhou University of Chinese Medicine, BALB/c: SCXK 2008-0002, C57BL/6: SCXK 2008-0002) were given 14 days of adaptation period with similar feeding and housing conditions prior to treatment to prevent impact from environmental variance. Mice were randomly divided into 4 groups, including the control group (Ctrl, non-treatment), Low-dose (Dex, 0.2 mg/kg four times per week by i. p.), Mid-dose (Dex, 0.5 mg/kg four times per week by i. p.) and High-dose (Dex, 1.0 mg/kg four times per week by i. p.). The administration period was 60 days in total. Each group of mice was housed individually in a room with a 12 h:12 h light/dark cycle with free access to complete pellet diets and water. All experimental procedures and animal care were approved by Animal Use and Care Committee of Guangdong Pharmaceutical University.

Mesenchymal Stem Cells (MSCs) were cultured in Dulbecco’s Modified Eagle’s medium (DMEM) containing 10% fetal bovine serum (FBS), 100 IU/ml penicillin, and 100 μg/ml streptomycin at 37°C in 5% CO2. Cells were passed twice a week. Differentiation of MSCs towards osteoblasts were performed according to previous protocols (Hanna et al., 2018).

Human bone marrow MSCs of passage 8-10 were pre-conditioned with or without 1 μM μM concentration of Dex for 72 h followed by a subsequent 21-days osteogenic differentiation protocol by StemPro Osteogenesis Differentiation Kit according to manufacturer’s manual. For examining the enhancer marker changes instigated along osteogenesis, Dex- or vehicle-treated MSCs were subjected to the same 72 h pre-conditioning followed by a brief 24 h treatment in the osteogenic medium. These pro-osteogenic MSCs were then collected and subjected to Cleavage Under Targets and Tagmentation (CUT and Tag) chromatin immunoprecipitation profiling with mono-methylated H3K4 (H3K4me1) enhancer mark-specific antibody (Abcam-ab176877) along the standard protocol in accordance to the manufacturer’s manual enlisted in the hyperactive In-Situ ChIP library prep kit (Vazyme-TD902). Other antibodies used in the current research work are anti-GAPDH (#5174) and anti-NR3C1 (#3660) (Cell signaling technologies), anti-TAGLN3 (abcam-ab204902), anti-CD36 (abcam-ab133625) and anti-stabilin-1 (Santa Cruz sc-293254).

All experimental subjects were labelled with calcein under double subcutaneous injection regime of 10 mg/kg per day for two consecutive days, 13, 14, 2 and 3 days before animal sacrifice. The proximal tibia metaphysis (PTM) and middle part of tibia shaft (TX) were taken for bone histomorphometry, anatomized to remove soft tissue, and then embedded in methyl methacrylate thereafter sectioned by hard tissue Microtome (Leica 2155, Germany). The slices were stained with Masson-Goldner Trichrome bone strain method described as previously reported (Foot, 1933).

Measurements were performed with a digital system consisting of fluorescent microscope with bright-field channels. The system was coupled to a computer with a morphometry program “Bioquant OSTEO 2009” (Bioquant Corporation, United States). Static parameters of trabecular bone, including the trabecular bone area percentage (%Tb.Ar) and trabecular thickness (Tb.Th), number (Tb.N), separation (Tb.Sp). For the cortical bone, parameters including total tissue area, percentage cortical area (%Ct.Ar), percentage marrow area (%Ma.Ar), cortical bone width and cortical bone endosteal formation rate (E-BFR/BS, a dynamic parameters). Label escape correction was used for the calculation of BFR (Jia et al., 2020).

Osteoblasts isolated from calvarias of 1-day-old C57BL/6 and BALB/c neonatal mice (Experimental Animal Center, Sun Yat-Sen University: 44008500012901) were cultured in MEM with 10% FBS and 100 U/mL penicillin, respectively. Cells were randomly divided into 2 groups, including Control group (MEM with 10% FBS and 100 U/mL penicillin) and Dex-treated group (Dex, 5*10⁻⁷M, MEM with 10% FBS and 100 U/mL penicillin). Passage 3 of osteoblasts were treated with Dex for 72 h.

Total cellular RNA was extracted and performed RNA sequencing (RNAseq). RNA-seq was performed using the BGISEQ-500 platform (BGI). The quality of sequencing reads was checked using the FastQC program (version 0.11.5). After adaptor trimming and removal of poor-quality reads (Q < 15) using SOAPnuke (SOAPnuke v1.5.2, BGI, China), Bowtie (bowtie2-v2.2.5) were used to map RNA-seq reads to mouse reference genome version mm10. After assembling reads to the genome using RSEM (RSEM, v1.2.12), gene expression was normalized against total read count (Fragment Per Kilobase of transcript per Million mapped reads, FPKM).

Differentially expressed genes (DEGs) in Dex-treated cells were determined based on Pearson’s correlation calculated using cor function of RSEM. Selection criteria for DEGs in this study were adjusted p-value (reported as false discovery rate in RSEM) ≤ 0.001 and fold change ≥2.0. For biological interpretation of DEGs, KEGG (the Kyoto Encyclopedia of Genes and Genomes) pathway analysis and Gene Ontology (GO) enrichment analysis of DEGs were performed using DAVID (the Database for Annotation, Visualization and Integrated Discovery, https://david.ncifcrf.gov/). For the meta-analysis of RNA-seq data for PCA analysis adapted for CC murine calvaria osteoblast culture from Sabik et al. (GSE134081) and osteogenic stromal progenitor culture from Khayal et al. (GSE54461), they were publicly-available from Gene Expression Omnibus database at the corresponding accession numbers.

Statistical analyses were performed using the Statistical Package for Social Sciences (SPSS 13.0; SPSS, Inc., Chicago, IL). Differences among groups were detected using one-way ANOVA. Values were expressed as mean ± SD, and the difference between groups was defined as significant at p < 0.05.

We firstly demonstrated that various strains of inbred mice had differential sensitivity towards bone loss recapitulating features associated with GIO. Here, adult mouse from either C57BL/6 or BALB/c mice were administered with various dosages of Dex for 60 days, followed by subjecting their tibia to bone histomorphometric analysis.

Indeed, BALB/c mice were more sensitive to the debilitating effects of Dex on trabeculae bone loss, despite a low-dose of Dex could only exert insignificant reduction in cancellous bone mass (Figure 1A). In contrast, C57BL/6 are largely tolerable to various dosages of Dex with minimal changes in bone mass except for a slight trend of dose-dependent loss (Figure 1B). Furthermore, we noted a reduction of cancellous bone in BALB/c mice PTM largely due to increased trabecular separation with decreased total trabeculae number (Figures 1C,E).

In addition, there is almost no differential impact in Dex treatment towards the cortical bone components (Figure 2A). The total tissue area remained the same as revealed from TX (Figure 2B).% Ct. Ar and cortical bone width were decreased (p < 0.05) while % Ma. Ar was increased (p < 0.05) in C57BL/6 mice compared to the control group (Figures 2C,E). Dynamic parameter of cortical bone also demonstrated a decrease in bone formation in endosteum in both strains (p < 0.01) (Figure 2F). Finally, we also confirmed that Dex exerted no effect in chondrocyte differentiation in the growth plates, bearing similar histomorphological features between these two strains. (Supplemental Figure S1). Therefore, our in vivo data confirmed BALB/c and C57BL/6 mice had recapitulated the iatrogenic effect of GC despite with a differential sensitivity towards the prolonged drug treatment likely due to their diverse genetic background.

To explore the transcriptomic differences contributed to various response to Dex caused by genetic background, we compared the expression profile of calvaria cultures among the previously reported CC mice, also importantly COL3.1-expression gated osteogenic progenitors from C57BL/6 across time; together with our two common laboratory strains of inbred mice (BALB/c and C57BL/6) with or without Dex administration after day 10 of mineralizing treatment ex vivo. Of note, CC mice would serve as ideal comparator to our experimental groups since one of their initial founder parent strains was the C57BL/6 line.

With the advantage of sharing the tissue of origin, we performed principal component analysis (PCA) in the expression profiles via a dimension reduction approach to explain underlying variation amongst data-points. We observed the spread of expression profile of all CC mice with various genetic markup across the PC2 panel, despite only accounting intrinsically for 10% of total variance (Figure 3A). When locating all data-points as of day 10 for mineralization treatment (only data points indicated in circles), we noticed that all RNA-seq profile for mice model not being descendant of those CC breeds fall outside the spread but succumb to scatter across the spread in the PC1 panel, which accounted for 64% of the variation (Figure 3A).

From that PCA analysis, the spread of blue data-points from Khayal et al. were also included to account for variations along stromal cell differentiation from C57BL6 mice culture enriched exclusively with COL3.1 + osteogenic stromal progenitors. Intriguingly, a coherent top-to-bottom spread in expression profiles along PC2 was found like that reported in CC mice (Figure 3A). Notably, the spread in CC mice expression profile along PC2 is correspondingly reflective of inherent efficacy of stromal cell differentiation potential as it was previously reported from Sabik et al. that there is a tendency of increasing mineralization in these CC mice from top to down along PC2 (Figure 3B).

It will therefore be interesting to firstly look for gene ontologies of candidate genes accounting for the differential expression along the variance of PC2. It was found that a list of 50 genes were either co-expressed or silenced together in individual datasets (Figure 3B). Surprisingly, an overwhelmingly 35 of them were previously defined to fall within the module with eigengene defining immune response. There is a clear trend of increase in these immune-related genes in the stromal progenitors from day 2 towards day 18, where the protein-protein interactions among these genes marked up a co-expressed, interactive network with predictive value of the tendency of bone mineralization among CC mice calvaria osteoblasts (Figure 3C). Such inverse expression-mineral content relationship invites our further study in how such co-expression network drove the corresponding stromal cell differentiation.

In contrast to that of PC2, the spread along PC1 explained more variation among data-points (Supplementary Figure S2). Noticeably from the PCA, our C57BL/6 calvaria culture expression profile faithfully resembles that profile of Day 10 mineralizing stromal progenitors (Figure 3A), indicative of the predictive nature of PC1 in aligning expression profile from the same genetic background. Likewise, it can also be confirmed that the expression profile of BALB/c mice lie approximately midway in the genetic variation between the two datasets of mouse calvaria expression profiles. In parallel, there is an agreeably smaller spread in intrinsic background differences among those CC mice along PC1 where their data-points clustered together as expected. It should be reminded that despite C57/BL6 served as one of the 8 founder breeds for descendant CC pure breeds, none of these mice descendants carry an expression profile that blend in with data-points obtained from that C57/BL6 background.

We were therefore eager to study the underlying co-expressed gene network governing this PC1 spread (Figure 3D). One must also observe that the top-enriched 20 genes unmistakably fall within the module with eigengene defining osteogenesis. To sum up, PC1 explains a large baseline difference in intrinsic osteogenic potential of calvaria cultures while PC2 accounts for the subtle but significant co-expression changes along a panel of genes implicated in affecting efficacy in stromal cell differentiation, which invite further studies in the action of its underlying gene network.

Furthermore, we were also concerned with whether such variant baseline osteogenic potential shall give rise to differential sensitivities for GC towards osteogenesis. Since our drug treatment versus control remained at the closest proximity to each other, it was clearly indicated that intrinsic genetic differences defining the PC1 spread already lead a larger effect size that masks over the incremental horizontal vector gap under prolonged Dex administration. Hence, most of the discrepancies in expression profiles of isolated calvaria could again only be explained on the grounds of genetic background. Yet, it could nonetheless be confirmed that there was a significant larger vertical vector gap between data points along the PC2 in BALB/c compared to that in C57/BL6.

As mentioned above, it was confirmed that most of these PC2 genes belong to the same module defined by an immune-related eigengene based on weighted gene co-expression network analysis (WGCNA). These gene had also been subjected to gene-gene interaction network clustering, from which we had successfully constructed the stromal cell differentiation network (Figure 3E). Revisiting the gene expression of the entire network collectively, we can extrapolate a more significant expression changes in the BALB/c mice between the calvaria culture with or without treatment with Dex compared to that between the cultures isolated from C57BL/6 (Figure 3E). Of note, the top DEGs within this network are the murine model specific Fcrls and Cd36, both of which belongs to the scavenger receptor B family.

Our pair-wise comparison with or without drug treatment in both inbred mice strains pinpointed a collection of DEGs as a result in Dex pre-conditioning (Figure 4A). Compared to untreated BALB/c osteoblast, 949 DEGs (415 up-regulated and 534 down-regulated) were detected after Dex treatment; while for the C57BL/6 osteoblasts, there were 1119 DEGs (578 up-regulated and 541 down-regulated) were detected (Figure 4B). GO analysis demonstrated cell cycle, angiogenesis, cell adhesion, cell division, mitotic nuclear division, positive regulation of angiogenesis, positive regulation of cell migration and response to lipopolysaccharide commonly observed to be associated with Dex preconditioning in both BALB/c and C57BL/6 primary osteoblasts (Figures 4C,D). Most of these genes fall under PI3K-Akt signaling pathway in both BALB/c and C57BL/6 osteoblasts.

Nonetheless, the biological processes most commonly associated with the drug treatment differs between the two mice strain, with cell proliferation and extracellular matrix organization specifically regulated in C57BL/6 osteoblasts, while chemotaxis, immune system process and inflammatory response were top-hits in BALB/c osteoblasts instead (Figures 4C,D). In line with these processes, DEGs were clustered within focal adhesion and ECM-receptor interaction pathways specifically regulated by GC in C57BL/6 osteoblasts, while the TNF signaling pathway and other chemokine/cytokine signaling pathway, notably those belonging to the bone morphogenetic protein (BMP) superfamily, were specifically regulated by GC in BALB/c osteoblasts.

When taken together, we believe Dex regulates osteoblasts cell cycle, cell division and PI3K-Akt signaling pathway in both BALB/c and C57BL/6, while BMP signaling most enriched under the action of Dex in BALB/c osteoblasts which might partially explaining GC-induced bone loss in BALB/c mice (Figure 4E, Supplementary Table S1, S2). As such, we cross-validated known DEGs of BALB/c and C57BL/6 belonging to BMP signaling pathway implicated in osteogenesis for further investigation (Figure 4E and Supplementary Table S3). Several genes promote bone formation such as Smpd3 and Pdcd4 were increased in a higher fold-change by GC in C57BL/6 osteoblasts than in BALB/c osteoblasts, while bone formation inhibitors Scx and Sfrp4 were decreased in a higher fold-change by GC in C57BL/6 osteoblasts than in BALB/c osteoblasts. Compared with C57BL/6 primary osteoblasts, GC significantly upregulated the expressions of Nog, Fst, Gdf6, Gdf10 and Gpc3, which are well-known as osteogenesis inhibitors in BALB/c primary osteoblasts.

Our prolonged Dex treatment in animal model effectively recapitulated symptoms and expression profile changes associated with GIO. Nonetheless, natural variance across various mouse strains investigated during mineralizing treatment demonstrated substantial natural variation that masks over the subtle pathway changes along Dex treatment. Since collective increment of individual effect sizes underlying genetic predisposition between these mice were not sufficient to mark up a unified locus associated with GC susceptibility, we turned to an alternative approach in gathering epigenetic and pharmacological data to provide biological inference for GWAS.

To begin with, MSCs with osteogenic potential were allowed to pre-condition in either vehicle control, 10 nM or 1 μM Dex for 72 h before subjected to osteogenic differentiation (Figure 5A). There is a slight enhancement in osteogenic differentiation for MSC-derived osteoblast culture under a 20-days protocol when subject to Alizarin Red staining. However, micromolar levels of Dex did not promote strong osteogenic differentiation as expected. We had repeated the same set of experiment probing into the alkaline phosphatase (ALP) activity and had come up with the same conclusion (Figure 5B).

Micromolar levels of Dex is way beyond the fluctuating ranges of the normal hormonal level as in GC. Therefore, we speculated an overwhelmed GC signaling when MSCs were pre-conditioned with such high levels of GCs. By cross-referencing both MTT assay and CCK-8 assay, we observed significant reduction in proliferating MSCs during the 72 h-interval of Dex preconditioning (Figures 5C,D). Taken altogether, GC seemed to have imposed an overkill effect in osteogenic activity after pre-conditioning.

We were satisfied with our GC-disrupted osteogenic human MSCs model at such desensitizing Dex level, i.e., planned as the experimental set-up intended to simulate the iatrogenic effect of GC-induced osteoporosis. To confirm that GC prevented osteogenesis through the actions of its nuclear receptor NR3C1, therefore we performed lentiviral induction of pSico-shNR3C1 in MSCs followed by preconditioning in the same concentration (Figure 5E). The plasmid was adopted for the specific reversion of the hairpin-mediated knockdown effect of the nuclear receptor with transient Cre expression. Results clearly demonstrated that Cre-mediated silencing of NR3C1 rescued the osteogenic effect of MSCs even when subjected to Dex pre-conditioning (Figure 5E). Since we had demonstrated the involvement of NR3C1 in previous experiment, we had ruled out non-specific pharmacological inhibition of receptor proteins instigated by Dex.

The human stromal cell osteogenic model had been adopted as a reliable model in studying osteogenesis. In line with the findings summarized from our mouse model, we hypothesize that GC pre-conditioned MSCs shall result in modifications within the epigenetic landscape which eventually give rise to variable extent of bone mineralization. In fact, there is a global induction of enhancer accessibility by performing chromatin-immunoprecipitation of the H3K4me1 mark, a unified epigenetic mark reflecting the enhancer regions (Figure 6A).

Furthermore, we therefore summarize the differential peak calling in the epigenetic mark H3K4me1 which may result in programmable tissue-specific chromatin accessibility which indirectly affect transcription factor binding in enhancer sites (Figure 6A). Enlightened by our previous mice experiments, we understood that cell proliferation through PI3K-Akt pathway and the BMP superfamily signaling pathways are involved in GC-induced osteogenic changes. We therefore tested whether these called peaks were functionally relevant to our GC-induced changes.

By focusing on those belonging to the transgelin (TAGLN) family, most notably member 3 (TAGLN3) (BMP superfamily signaling pathway) and neuregulin 1 (NRG1) (PI3K-Akt pathway) as they were also previously reported to play a role in MSC-derived osteogenesis (Figure 6B). We confirmed significant relaxation in selected regions of chromatin spanning their corresponding promoter with enhanced RNA polymerase and NR3C1 occupancy normalized to housekeeping gene GAPDH as confirmed by ChIP-PCR (Figure 6C). We monitored the gene expression of TAGLN members and NRG1 along osteogenic differentiation. Interestingly, TAGLN members were induced along the time-course of osteogenesis while there was a spike in NRG1 expression after Dex pre-condition followed by a gradual reduction along osteogenesis. In contrast, mild expression changes in these genes were observed in vehicle-treated only MSCs without Dex pre-conditioning (Figure 6D). Besides, we had also carefully validated the changes in protein expression of NR3C1 and TAGLN3 under the influence of DEX pre-condition (Figure 6E).

Osteogenic transcription factors, including RUNX2, OTX and osteoblast marker COL1, ALP, OPN and OCN were examined in pre-conditioned MSCs with or without co-treatment with SB431542 (TGF-beta inhibitor). Surprisingly, co-treatment with TGF-beta inhibitor (TAGLN pathway antagonist) in pre-conditioned MSCs enhanced osteogenesis marker expression (Figure 6F). Co-culturing of pre-conditioned MSCs with either recombinant NRG-1/EGF (NRG-1 agonists) also effectively promoted proliferation of the pre-conditioned MSCs resulting in partial rescue of osteogenesis (Figures 6G,H).

Besides, by leveraging inputs from epigenetic information, potential disease-relevant genes could be inferred from suggestive significant loci to provide the interpretation of pathological mechanisms. These differential epigenetic information before and after drug treatment shall render corresponding changes in cis-regulatory DNA histone marker changes as an alternative cost-effective approach to infer and predict trait-associated loci, in particular expression quantitative traits loci (eQTLs) directly or indirectly from GWAS results where sampling size remains inadequate.

In line with the findings in the mouse model, several of these genes were interconnected in the co-expressed network to the scavenger receptor Cd36, which itself belongs to the reported osteogenic core module. We investigated the differential H3K4me1 peaks along the Cd36 locus (marked in green). There is generally a more active open chromatin at the site in untreated osteoblasts compared to those pre-treated with Dex. Moreover, we confirmed a differential expression in Cd36 expression inferred from both mouse calvaria RNA-seq and qPCR analysis in the differentiated human MSCs (Figures 6I,J). To further reinstate the differential expression of scavenger receptors under the influence of Dex, we probed the changes in protein expression of CD36 and STAB1 in MSCs differentiated in osteogenic medium (Figure 6K).

Summarizing from above, we believe analyzing the global changes in H3K4me1 mark along osteogenic differentiation with or without Dex conditioning may suggest previously unknown targets such as TAGLN family in the TGF-beta signaling pathway and NRG1 in the EGF/mitogenic pathway potentially implicated in excessive GC-repressed osteogenesis and proliferation. Furthermore, we also identified a core-like scavenger receptor Cd36 interaction network responsible for stromal cell osteogenic differentiation which would effectively be modified under the actions of GC, potentially serving as a proxy biomarker for the idiosyncrasy in GIO.