, , and Are Related to Glucocorticoid-Induced Osteoporosis Occurrence According to Integrated Bioinformatics Analysis

Background

Glucocorticoids are potent anti-inflammatory agents frequently used to treat syndromes related to inflammation. Nonetheless, these widely used agents also lead to serious adverse side effects, including glucocorticoid-induced osteoporosis (GIOP), which shows an increasing incidence over the past few decades, along with a younger trend [1,2]. Compared with treatment for degenerative osteoporosis, GIOP treatment is mainly centered on the molecular biological mechanism and drug efficacy, which are closely associated with the promotion of bone remodeling and skeletal metabolism [3,4]. However, bisphosphonates, calcium supplementation, and additional typical anti-osteoporosis treatments do not yield satisfactory results among patients dependent on steroid agents. Further, there have been no large-scale clinical studies of sufficient duration to prove which drug treatment plan is most effective for GIOP. It is difficult to develop uniform guidelines to manage GIOP because diagnosis and treatment methods differ across diverse countries. At present, bisphosphonate treatment has been utilized as the standard for GIOP care; however, these agents have an undetermined therapeutic effect among patients receiving glucocorticoid treatment for longer than 2 years [5]. Consequently, understanding the precise molecular mechanism of GIOP pathogenesis is necessary for identifying more potent treatments to control osteoporosis occurrence and development. Although many factors, such as osteoclastogenesis [6], apoptosis [7], osteoblast autophagy [8], Wnt/β-catenin signaling pathway [9], and altered intestinal microbiota composition [10], have been found to be associated with GIOP, specific diagnostic and therapeutic biomarkers for GIOP have not yet been identified. As a result, further investigations are warranted to develop more diagnostic and prognostic biomarkers for GIOP.

Bioinformatics analysis is extensively used in screening and analyzing genes linked with the progression of diverse disorders, helping to overcome the limits of experimentation. Numerous gene profiles have been acquired based on public databases, such as Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA), expanding our knowledge of diseases. Due to different microarray platforms of datasets, sample sizes, and heterogeneities of species and tissues, limitations and inconsistent results may still exist, but integrated bioinformatics approaches may overcome these limitations. Weighted gene coexpression network analysis (WGCNA) has been developed as a novel approach for analyzing gene expression profiles among a variety of samples [11]. WGCNA can analyze gene sets with highly synergistic alterations as well as potential biomarker genes, and thus can be useful in identifying treatment targets according to the associations between gene sets as well as between phenotypes and gene sets. In recent years, WGCNA has been extensively utilized in genomic research, such as Parkinson disease [12], pancreatic cancer [13], glioblastoma [14,15], and non-small cell lung cancer [16].

In recent years, some collagen gene mutations are detected from probands with genetic disorders. Some of them show phenotypes that are difficult to distinguish from common diseases. For example, collagen type III alpha 1 chain (COL3A1) mutations are related to numerous human disorders, including osteoporosis [17], osteoarthritis [18], and aortic aneurysms [19]. The collagen alpha-3 (VI) chain, encoded by gene COL6A3, is a microfibrillar part of the extracellular matrix (ECM) [20]. Serpin family H member 1 (SERPINH1) has been reported to facilitate type I procollagen molecular stability, which exerts a vital role in maintaining normal bone mass level and quality and avoiding contracture and fracture [21]. To our knowledge, a functional link between the above 3 genes and GIOP has not been reported so far. Therefore, in this study, we used an integrated analysis of bioinformatics, which included the protein-protein interaction (PPI) network, WGCNA, and gene set enrichment analysis (GSEA), to analyze core genes and pathways involved in the pathogenesis of GIOP in endogenous Cushing syndrome (CS). The present results revealed that COL3A1, COL6A3, and SERPINH1 could serve as potential critical biomarkers associated with GIOP. As far as we know, this work is the first to reveal that COL1A3, COL6A3, and SERPINH1 are upregulated in endogenous CS patients with GIOP. Our current findings provide insights into the pathogenesis of GIOP and identify novel targets for the treatment of GIOP. The study design is shown in the flow chart in Figure 1.

Material and Methods

SEARCH STRATEGY AND ELIGIBILITY CRITERIA:

The mRNA expression data of patients with GIOP were download from the GEO database (http://www.ncbi.nlm.nih.gov/geo/) on June 30, 2020. The terms “corticosteroids” and “osteoporosis” were searched alone or in combination. The type of study was “expression profiling by array,” and the entry type was “datasets.” We included studies conforming to the following criteria in our analysis: (a) studies in which the datasets for gene expression covered the microarray chip technique, and (b) studies that compared gene expression profiles between non-GIOP and GIOP status. Articles not conforming to these criteria were excluded from our analysis. The database search was independently carried out by 2 researchers.

ESTABLISHMENT OF COEXPRESSION NETWORK AND ANALYSIS OF MODULE FUNCTIONS:

First, expression profiles of differentially expressed genes (DEGs) were examined for screening the appropriate genes and samples. Second, a coexpression network of DEGs was established using ‘WGCNA’ package in R language [23,24]. Afterwards, each pair-wise gene was functioned using Pearson’s correlation matrix. Third, the power function amn=|cmn|β(cmn=Pearson’s correlation between gene m and gene n; amn=adjacency between gene m and gene n) was applied to create a weighted adjacency matrix, while β was the soft threshold factor adopted to stress the strong associations across genes and to penalize the weak relationships. Fourth, topological overlap matrix (TOM) adjacency was converted to measure the gene network connectivity, which was deemed to be the total value of the adjacency to the remaining genes in generating the network. Mean linkage hierarchical clustering was created by the dissimilarity measure based on TOM, and the minimal size (gene group) was set at 50 for the gene dendrogram. Thereafter, genes that had similar expression patterns were clustered within the same gene module. Lastly, the module eigengene dissimilarity was determined. Then, the gene modules were used to perform functional enrichment analyses to identify the related modules affecting GIOP in endogenous CS patients.

IDENTIFICATION OF GIOP STATUS HUB MODULE:

For identifying modules showing significant associations with illness state traits (GIOP vs. non-GIOP), module eigengenes (which represent the first principal component in a module) [22] were associated with the external traits to identify correlations with the highest significance. Meanwhile, module membership (MM) indicated the relationship between gene expression patterns and module eigengenes. At the same time, gene significance (GS) measure was indicative of absolute value relationships of genes with the external traits. In this study, genes that showed the greatest GS and MM values in the modules of interest were identified as the natural candidates in later analysis [23–26].

HUB GENES VALIDATION:

A hub gene is substantially related to other genes within the module, which is suggested in previous studies to display functional significance. This study screened hub genes within the coexpression network from the GIOP phenotype-related module. Subsequently, each hub gene within the module was imported to the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) (https://string-db.org/) [27], and confidence >0.4 was chosen for creating the PPI network. Later, Cytoscape software (www.cytoscape.org/) [28] was used to visualize the PPI network. Nodes that had a great degree of connectivity were more important to maintaining the network stability. The Cytoscape plug-in CytoHubba was used for calculating every protein node degree. Any gene within the PPI network with a connectivity of ≥5 (node/edge) was screened to be a hub gene. Later, the common hub genes within the coexpression and PPI networks were deemed as “real” hub genes, and they were screened in later analysis. The Venn diagram was constructed using Venny 2.1.0 (https://bioinfogp.cnb.csic.es/tools/venny/index.html) for visualizing the common hub genes in the PPI and coexpression networks. Then, both the training and test sets were used for validation. For training set GSE30159 dataset and test set GSE129228 dataset, the real hub genes were compared between GIOP and normal controls. The paired-sample t-test was used for statistical analysis, and a difference of P<0.05 indicated statistical significance. Figures were plotted using GraphPad Prism (version, 8.0; GraphPad Software, Inc., La Jolla, Ca, USA).

GENE ONTOLOGY ANNOTATION AND KYOTO ENCYCLOPEDIA OF GENES AND GENOMES PATHWAY ENRICHMENT ANALYSES:

Gene Ontology (GO) functional annotation has been developed as an efficient way to carry out functional enrichment in a large scale. In addition, Kyoto Encyclopedia of Genes and Genomes (KEGG) is an extensively applied database that preserves extensive data on drugs, chemical substances, diseases, biological pathways, and genomes. The current work employed the Metascape software (http://metascape.org) [29] in GO as well as KEGG analysis on the DEGs. P<0.05 indicated statistical significance.

GENE SET ENRICHMENT ANALYSIS:

To further determine functions of the candidate hub genes, we conducted GSEA (https://software.broadinstitute.org/gsea/index.jsp) [30] to investigate the enrichment of previously determined biological processes within the DEGs-derived gene rank. Terms enriched in each gene were recognized with the thresholds for false discovery rate (FDR) q-value <0.25 along with nominal P-value <0.05.

Results

INCLUDED STUDY CHARACTERISTICS:

After a careful review, 2 microarray datasets (GSE30159 and GSE129228) were selected from the GEO database [31–33]. GSE30159 was based on platform GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array, and the GSE129228 dataset was based on platform GPL21103 Illumina HiSeq 4000 (Mus musculus). The present study utilized the GSE30159 dataset as the training set for constructing the PPI and coexpression networks for identifying “real” hub genes as well as related pathways. This dataset was obtained from bone biopsies from 18 patients with endogenous CS prior to surgery and at 3 months after surgery. In the study related to the GSE129228 dataset, high-dose dexamethasone was used to establish a mouse steroid-induced osteoporosis model. This dataset contained samples from 6 groups, including 2 normal samples, 2 experimental groups receiving 10 μM dexamethasone alone samples, and 2 GIOP model samples. The GSE129228 dataset was used as the test set for result validation in the current study.

DEGS IDENTIFICATION:

We utilized R language ‘limma’ package to identify DEGs within the GSE30159 microarrays. Using |logFC| ≥1 and P<0.05 as thresholds, 628 DEGs were detected, which included 303 upregulated genes and 325 downregulated genes. Benjamini-Hochberg correction was utilized for adjusted P-values. A volcano plot for DEGs from this microarray is shown in Figure 2.

CONSTRUCTION OF THE WEIGHTED COEXPRESSION NETWORK AND IDENTIFICATION OF HUB MODULES:

Gene coexpression networks were constructed using R language ‘WGCNA’ package. Then, altogether 23 321 genes were obtained. At first, the leading 5000 genes in terms of standard deviation (SD) values were chosen to perform the sample clustering based on phenotype by the use of the average linkage (Figure 3A). It was observed that 18 samples could be basically classified as 2 clusters. In addition, the Pearson’s correlation was also carried out. Power β=9 was chosen for guaranteeing a scale-free network (Figure 3B, 3C). Eight modules in total were mined, and the red, brown, and turquoise modules were the most tightly related to GIOP in endogenous CS patients (Figure 4A, 4B). Thereafter, the interactions among these 8 modules were also examined, followed by the plotting of a network heatmap (Figure 5A). According to these findings, every module served as an independent validation for one another, demonstrating the high level of independence across various modules, as well as the relative gene expression independence for every module. For purposes of exploring coexpression similarity among these 8 modules, eigengene connectivity was assessed, and then consensus correlation was subjected to clustering analysis (Figure 5B). A heatmap drawn on the basis of adjacencies also displayed similar findings (Figure 5C). In addition, intramodular analysis including MM (module significance) and GS (gene significance) was performed in those 8 modules. The turquoise module was then excavated to further explore the highly related genes. Figure 5D displays the scatter plots for GS regarding GIOP traits, together with illness state compared with MM within the turquoise module. In the case of GIOP, GS and MM (Figure 5D) were significantly positively correlated, which suggested that the turquoise module elements with the greatest importance (central) might show high correlation with such external traits.

IDENTIFICATION OF GIOP STATUS HUB GENES WITHIN HUB MODULES:

Subsequently, the DEGs-associated PPI network was used to identify 31 hub genes at the thresholds of connectivity ≥2 and confidence >0.4. The stricter factors were used in additional analyses, including module connectivity determined through absolute Pearson’s correlation coefficient (cor.geneModuleMembership >0.8), together with relationships of clinical characteristics determined based on absolute Pearson’s correlation coefficient (cor.geneTraitSignifcance >0.2). There were 262 highly connected genes identified in the turquoise module, of which, COL3A1, POSTN, COL6A3, COL14A1, SERPINH1, ASPN, OGN, THY1, NID2, and TNMD were detected in both the coexpression and PPI networks (Figure 6A). According to our results, each hub gene within endogenous CS patients with GIOP was upregulated. Therefore, the above 10 genes were identified to be real hub genes indicating GIOP status and were screened in later analyses (Figure 6B).

VALIDATION OF HUB GENES:

To investigate hub genes in endogenous CS with GIOP, the expression levels of COL3A1, POSTN, COL6A3, COL14A1, SERPINH1, ASPN, OGN, THY1, NID2, and TNMD were detected using the training set GSE30159 dataset and the test set GSE129228 dataset, respectively. In the training set, we found all hub genes except POSTN had statistically significant differences in endogenous CS (Figure 7A–7J). In the test set, SERPINH1 and POSTN were significantly upregulated in the GIOP model groups in comparison with the non-GIOP groups, which included CG (normal) and DEX (dexamethasone alone) as control groups (Figure 8A–8H). However, the expression levels of COL14A1 and TNMD were not found in the test set. After overlapping the results from the training set and test set, we found SERPINH1 was altered in the comparison between the GIOP and normal control samples.

GO ANNOTATION AND KEGG PATHWAY ENRICHMENT ANALYSES:

For better understanding of the gene functions within the turquoise module, Metascape software was used to perform GO enrichment analyses. Based on our results, “extracellular matrix” was the gene set with the highest significance (Figure 9A). The analysis also showed that GIOP was associated with blood vessel development, metallopeptidase activity, integrin binding, connective tissue development, regulation of neuron differentiation, response to growth factor, tissue regeneration, regulation of cell morphogenesis, negative regulation of cell migration, regulation of animal organ morphogenesis, syncytium formation by plasma membrane fusion, positive regulation of Ras protein signal transduction, intrinsic component of external side of plasma membrane, ossification, metallocarboxypeptidase activity, negative regulation of cell differentiation, cellular response to vitamin, cell-matrix adhesion, basement membrane, and regulation of the Wnt signaling pathway, planar cell polarity (Figure 9B). Meanwhile, based on KEGG analysis, DEGs were mostly enriched in the pathways in apoptosis, protein digestion and absorption, Rap1 signaling pathway, protein processing in endoplasmic reticulum, and adrenergic signaling in cardiomyocytes (Figure 9C).

GENE SET ENRICHMENT ANALYSIS:

For better elucidating the functions of possible hub genes, we carried out GSEA. The training set was divided into 2 groups according to the label GIOP vs. NON-GIOP. In the GSEA software, normalized enrichment score (NES) is used as a measure of the degree of enrichment, and P-value and FDR are used as measures of statistical significance. The degree of enrichment is scored. If it is positive, the pathway tends to be enriched in genes that are upregulated, and if it is negative, the pathway tends to be enriched in genes that are downregulated. Based on our observations, both gene sets showed correlations with GIOP, and “ECM receptor interaction” (Figure 10A) and “focal adhesion” were enriched (Figure 10B). Moreover, up-regulation of both COL1A3 and COL6A3 was significantly associated with ECM receptor interaction (Table 1) and focal adhesion (Table 2).

Discussion

In the present study in which an integrated bioinformatical study on GIOP was performed, an overlap method was employed to combine WGCNA, PPI network, and GSEA to identify pathway-related genes. As suggested by our results, the turquoise module was recognized to be of clinical significance by WGCNA. In later analyses, 10 genes (COL3A1, POSTN, COL6A3, COL14A1, SERPINH1, ASPN, OGN, THY1, NID2, and TNMD) between coexpression and PPI networks were identified to be the real hub genes. Subsequently, to investigate hub genes in GIOP, the expression levels of the 10 genes were detected using a training set and a test set, respectively. In the training set, we found that all hub genes except POSTN had statistically significant differences in GIOP. However, in the test set, there was no significant difference in the expression of COL3A1 and COL6A3 among the validation groups. We believe the reason for the inconsistent results was due to different microarray platforms of the 2 datasets, different sample sizes, and heterogeneities between the species and tissues. Consequently, to further identify the potential function of the hub genes in GIOP, GSEA was conducted to search KEGG pathways enriched in the highly expressed samples, as previously reported [34,35]. In our results, the focal adhesion and the ECM receptor interaction pathways notably showed significant enrichment within samples showing high COL1A3 and COL6A3 expression levels. Consistent with our findings, previous studies identified COL3A1 and COL6A3 as being significantly associated with osteoporosis [17,18,36,37]. Overall, we determined that COL3A1, COL6A3, and SERPINH1 exhibited a high correlation via the integrated analysis among the 10 hub genes, indicating the potentially vital roles of such genes during GIOP occurrence or development.

Over the last several years, several GIOP-related genes or gene pathways alterations have been discovered, and these have shed more light on endogenous CS with regard to its pathogenesis. Nonetheless, a large number of articles are single cohort studies, which report diverse findings. Lekva et al. [31] used global gene expression profiles based on bone biopsies of CS cases and identified a gene that coded the GC-induced leucine zipper (GILZ) as an extensively modulated gene within CS. In addition, the mRNA level of COL1A2 and circulating osteocalcin were also found to be related to GIOP. In subsequent experiments, Lekva et al. [32] used the same gene expression profiling and found another gene that encoded TXNIP to be an extensively modulated gene associated with GIOP treatment. In one recent study, Lu et al. [33] applied an RNA sequencing (RNA-seq) technique in combination with bioinformatic analysis and discovered that endothelial progenitor cell extracellular vesicles avoided GIOP in mice through inhibition of the ferroptotic pathway within osteoblasts. To better understand the pathogenesis of GIOP, potential biomarkers were identified using an integrated bioinformatics analysis combined with gene expression profile analysis in the present study.

Currently, bioinformatic analysis allows for identifying vital molecular networks based on gene expression data. In addition, integration of multiple-gene microarrays can help to identify the gene biomarkers with higher accuracy. Moreover, those integrated bioinformatic approaches contribute to overcoming hurdles in identifying GIOP-related hub genes. Consequently, the present work used the gene expression profile datasets of 2 cohorts in diverse groups in combination with bioinformatic approaches to analyze raw data and identify new disease pathogenesis as well as novel diagnostic and prognostic biomarkers. In accordance with Lekva et al. [31], GILZ was significantly downregulated in our study. Interestingly, COL1A2 was the same as the DEGs in the PPI network and was related to focal adhesion along with ECM receptor interaction based on our GSEA analysis. Consistent with our study, another GIOP-related study by Lekva et al. [32] showed that TXNIP in bone was significantly downregulated in CS patients, while osteocalcin was significantly upregulated. The current work aimed to detect genes or pathways involved in GIOP pathogenesis, thus narrowing the scope of targets and offering novel targets for later analysis. According to our results, COL3A1, COL6A3, and SERPINH1 are possible candidates and potential new therapeutic targets. In our opinion, the differences observed might be associated with the following reasons. First, sample heterogeneity in independent studies or research utilizing only a single cohort may lead to inconsistent results. Second, comparison between GIOP samples and non-GIOP samples may result in a potential bias due to disease heterogeneity. Third, different thresholds for DEGs may screen out a corresponding number of genes from microarray. Fourth, if the samples are relatively small while the number of dependent variables substantially increases, the analyses of gene microarrays may lead to “curse of dimensionality,” which results in an increase in statistical errors [38]. Fifth, different microarray platforms often lead to inconsistent results. It should be noted that a functional link between the above 3 genes and bone metabolism has not been previously shown [31–33,39].

To further investigate the effects of hub genes on modulating GIOP, we conducted a gene functional enrichment analysis. First of all, GO functional annotation and KEGG analyses were carried out. According to GO analysis, the hub genes that we identified showed significant enrichment in the ECM, which includes a complicated mixture of functional and structural macromolecules that play vital roles during organ and tissue morphogenesis, as well as in tissue and cell structural and functional maintenance. Such associations allow direct or indirect control of cellular events, such as proliferation, differentiation, migration, adhesion, and apoptosis. Additionally, integrins play roles as mechanoreceptors, which provide a physical link for transmitting force between the cytoskeleton and the ECM. According to the results of the KEGG pathway enrichment analysis, the apoptosis-related pathway was the most significantly enriched. Consequently, the effects of COL1A3 and COL6A3 might be associated with apoptosis. GSEA was conducted to test this hypothesis. It was noticeable that the ECM receptor interaction together with the focal adhesion pathway was significantly enriched within samples showing high COL1A3 and COL6A3 levels. Consequently, COL1A3 and COL6A3 could have potential roles within the apoptosis pathway, which is functionally correlated with apoptosis. As suggested in numerous works, apoptosis plays an important role during bone metabolism through actions such as bone absorption and formation [7,40–43].

Downregulation of COL3A1 has been reported to be associated with WNT signaling regulation, ECM component shifts, and the differentiation and proliferation of cells within the subchondral bone-articular cartilage unit in Frzb−/− mice [44]. Furthermore, Yuan et al. [45] confirmed that COL3A1 was related to the healing of fractures based on RT-PCR assay. In addition, Stéger et al. [46] identified that COL3A1 expression within ossified velvet antler was 10 to 30 times higher than that in the skeleton. These findings implicate COL3A1 as having a role human osteoporosis. However, studies on COL3A1 in the context of GIOP are rare, and more research is needed. In a previous study on bone and mineral research, COL6A3 was found to have a significant correlation with osteoclast-like cell expression [37]. Furthermore, Barik et al. [47] demonstrated that COL6A3 expression on superimposed nano-textured surfaces was related to track osseointegration on implant surfaces. A previous study reported that COL6A3 participated in chondrocyte hypertrophic differentiation as well as ECM integrity in cartilage and basement membrane [48]. Additionally, Chou et al. [49] indicated that COL6A3 expression showed a statistically significant correlation with articular cartilage degeneration severity and subchondral bone structural changes. However, the influence of COL6A3 expression on GIOP remains undefined; therefore, more data are required to verify the suggested effect. These findings indicate that the effects of collagen cross-linking on osteogenesis can serve as the candidate pathogenesis within skeleton disorders. These results are consistent with our findings above.

Similar to results obtained in our analysis, a previous study reported SERPINH1 leads to telomerase deficiency by engaging in ECM, which blocked the differentiation of bone marrow stem cells to osteoblasts, finally affecting their proliferation, commitment, maturation, and matrix mineralization [50]. Furthermore, a previous study reported that the specificity of SERPINH1 expression in skeletal tissues was related to the in vitro phenotype of human mesenchymal stem cells collected from osteoporosis cases [51]. It is notable that SERPINH1 is suggested to exert a vital part in osteogenesis imperfecta, which results in low bone mineral density and brittle bones [52–55]. According to a previous study on skeletal micromorphology and osteogenesis, BAPN remarkably hindered osteogenic gene expression and altered the bone microstructure by influencing the expression of SERPINH1 [56]. As far as we know, this work first revealed that SERPINH1 was upregulated in endogenous CS patients with GIOP. However, further research is required to confirm this finding.

In total, the integrated bioinformatic research was conducted on GIOP in the current work, and an overlap method was employed to combine WGCNA, PPI network, and GSEA pathway-related genes. Subsequently, the datasets were restricted to the matched GIOP with non-GIOP samples, and then the hub genes in every dataset were analyzed using the paired-sample t-test. Since rigorous screening approaches were utilized, our findings may have a high specificity for the detection of key molecules related to GIOP. Nonetheless, certain limitations should still be noted. First, this work was conducted based on microarray data collected via the GSE30159 dataset, which had a small sample size. Therefore, if the database has samples updates, more studies are warranted in the future. Second, different hub genes and enriched functions were found in our work, but the relationships need to be examined further. Because most of our recognized genes had not been previously associated with GIOP, further studies are needed to validate such gene expression within GIOP as well as healthy control tissue samples. Also, cells isolated from GIOP tissue samples should be cultured in vitro to determine the molecular mechanisms associated with the expression of these genes. Preclinical animal models using gene knockout may also identify the functions of the genes identified and assess their role in the progression of the GIOP.

Conclusions

In conclusion, this work applied an integrated method, including WGCNA, PPI network, and GSEA, in identifying and validating hub genes related to GIOP. Our findings explained, from a bioinformatics perspective, the potential key genes (COL3A1, COL6A3, and SERPINH1) and molecular mechanisms that play an important role in the GIOP process. Similar reports are still rare. These potential biomarkers could enhance the early diagnosis and treatment of GIOP. Unfortunately, independent validation experiments were not carried out in this study. Therefore, in our future studies, we will design and conduct more rigorous experiments to verify the above findings.