The present study was approved by the Ethics Committee of the Second Affiliated Hospital of Xi’an Jiaotong University (Ethical Approval number No. 2019035). All tissue donors provided written informed consent for this study. Finally, subchondral bone samples from six patients who underwent primary total hip arthroplasty (three ONFH patients and three patients in the control group with femoral neck fracture) were obtained from the Second Affiliated Hospital of Xi’an Jiaotong University. All the samples were snap-frozen and stored in liquid nitrogen until further analysis.

Frozen subchondral bone samples were rapidly ground in liquid nitrogen. Total RNA was extracted by TRIzol reagent (Invitrogen) following the manufacturer’s instructions. The NEBNext Poly(A) mRNA Magnetic Isolation Module (New England Biolabs), RiboZero Magnetic Gold Kit (Illumina), and KAPA Stranded RNA-Seq Library Prep Kit (Illumina) were used to achieve RNA enrichment and sequencing library generation. An Agilent Bioanalyzer 2100 system was used to conduct the quality control analysis. Finally, high-throughput RNA sequencing was performed based on the Illumina HiSeq 6000 sequencing platform (Illumina) using the TruSeq SR Cluster Kit (Illumina).

After trimmed and reads filtering (Solexa pipeline program, Cutadapt software, and FastQC software), further lncRNA and mRNA alignment were both conducted by Hisat2 according to human reference genome indexing (hg38). Then we performed principal component analysis to calculate the heterogeneity between samples and remove the outliers. The differentially expressed lncRNAs (DELs) and differentially expressed mRNAs (DEGs) were defined as lncRNAs or mRNAs with |log2FoldChange| ≥ 1 and a p value < 0.05 [calculated by “edgeR” package (Robinson et al., 2009)].

Gene Ontology (GO) function enrichment analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathway enrichment analysis, and gene set enrichment analysis were conducted by the “clusterProfiler” package with the threshold set at p value < 0.05.

Human BMSCs were purchased from Procell (CP-H166, China) and cultured in minimum essential medium-α (α-MEM, Gibco, United States) supplemented with 10% fetal bovine serum (FBS, Gibco, United States) at 37°C in a humidified environment containing 5% CO₂. When the BMSCs reached 70% confluence in six-well plates, osteogenic induction was performed using osteogenic medium (100 nM dexamethasone, 50 μg/ml ascorbic acid, and 10 mM β-glycerophosphate, Sigma-Aldrich).

We performed alkaline phosphatase (ALP) staining (Yan et al., 2020) on days 3 and 7 after the induction of osteogenic differentiation to measure the osteogenic differentiation phenotype of BMSCs using a BCIP/NBT ALP staining kit (C3206, Beyotime, China). Briefly, the plates were washed with phosphate-buffered saline, fixed with 4% paraformaldehyde, washed with ddH₂O, ALP stained with BCIP/NBT working solution at room temperature, and washed again with ddH₂O. Finally, ALP staining image analysis was conducted using ImageJ software.

In this study, we performed reverse transcription using StarScript II First-strand cDNA Synthesis Mix (A223-02, Genestar, China). Subsequently, 2× Universal SYBR Green Fast qPCR mix (RK21203, Abclonal, China) was used to conduct the quantitative real-time PCR (RT-qPCR) according to the manufacturer’s instructions. The results were analyzed by using the 2−ΔΔCt method, and GAPDH was used as an endogenous reference gene to normalize the gene expression data. The sequences of primers are as follows:

Human GAPDH (forward primer GAC​AGT​CAG​CCG​CAT​CTT​CT, reverse primer GCG​CCC​AAT​ACG​ACC​AAA​TC), human GAS5 (forward primer GTT​GTG​TCC​CCA​AGG​AAG​GAT​GAG, reverse primer TGT​CTA​ATG​CCT​GTG​TGC​CAA​TGG), Human RUNX2 (forward primer AGC​AGC​ACT​CCA​TAT​CTC​TAC​TAT, reverse primer CAT​CAG​CGT​CAA​CAC​CAT​CAT), rat GAS5 (forward primer GCA​AGC​TCC​ACA​CAA​GGT​CCT​TC, reverse primer TGT​TCA​AGC​ATC​CAT​CCA​GTC​ACC), rat GAPDH primer (purchased from Sangon Biotech, China, B661204).

All animal experiments included in this study were approved by the Laboratory Animal Care Committee of Xi’an Jiaotong University. A total of 24 adult SPF Sprague–Dawley rats weighing 350–400 g were purchased from the Medical School of Xi’an Jiaotong University Animal Center. The rats were divided into two groups: the ONFH group and the control group. The rats were anesthetized with sodium pentobarbital (50 mg/kg, intraperitoneal injection), the fur was removed, and the skin was sterilized around the surgical area. An incision was made in the hip joint (right side) and then separated by layer to expose the femoral head. The hip joint was dislocated, the round ligament was removed, blood supply to the femoral head was stopped, the wound was disinfected, and then the incision layer was replaced. All the rats were killed after 9 weeks. Six femoral heads on the right side were collected for further micro-CT analysis, and 18 subchondral bone samples from the femoral heads were collected for RNA extraction and quantitative real-time reverse transcription polymerase chain reaction (RT-qPCR).

The main bone tissue parameters of the femoral heads were evaluated using high-resolution micro-CT scanning and reconstruction. Finally, the values of bone volume/tissue volume (BV/TV, %) were calculated to reflect the bone mass condition in the femoral heads, and the trabecular thickness (Tb.Th, mm) and trabecular space (Tb.Sp, mm) were calculated to reflect the trabecular structure.

Independent sample t-tests and Mann–Whitney U tests were conducted for the difference analysis according to the distribution of independent variables. IBM SPSS Statistics 22 software was used to perform all the statistical analyses, and p < 0.05 was considered statistically significant.

We performed lncRNA and mRNA annotation based on RNA-seq data, and the expression levels of 2966 lncRNAs and 13925 mRNAs were detected in clinical subchondral bone samples in this study. One sample was recognized as an outlier sample and removed through principal component analysis (PCA, shown in Supplementary Files S1, S2). Then, according to the thresholds as |log2FoldChange| ≥ 1 and p value < 0.05, a total of 126 differentially expressed lncRNAs (DELs) and 959 differentially expressed mRNA (DEGs) were identified (Figures 1A,B). Additionally, our results in Figure 1C showed that the lncRNA GAS5 expression level of ONFH patients was significantly different from that of the control group patients (log2FoldChange < −1 and p value = 0.003).

GO function enrichment analysis indicated that these DELs and DEGs were highly associated with biological processes like extracellular matrix organization, ossification, and bone mineralization; molecular functions like extracellular matrix structural constituent and signaling receptor activator activity; and cell components like the collagen-containing extracellular matrix and the external side of the plasma membrane. In addition, KEGG signaling pathway enrichment analysis suggested that these DELs and DEGs may play roles in signaling pathways like the PI3K-Akt signaling pathway, the HIF-1 signaling pathway, and cytokine–cytokine receptor interaction (Figures 2A,B). In order to improve the accuracy of enrichment analysis, we further performed GSEA analysis for all these DELs and DEGs. The results showed that they were involved in ossification processes, receptor regulator activity, and PPAR signaling pathways (Figures 2C–F). Therefore, it is necessary to investigate the association between lncRNA GAS5 expression and BMSC osteogenic differentiation.

lncRNA and mRNA expression profiles during BMSC osteogenic differentiation periods (0 h, 4 h, 1 day, 3 days, 7 days, 14 days) in GSE113253 datasets were downloaded from the GEO database (Rauch et al., 2019). We further calculated the co-expression correlation between the novel biomarker lncRNA GAS5 and mRNAs. In total, 683 co-expression genes were identified based on the thresholds of |Pearson correlation coefficient| > 0.9 and p < 0.05 (Supplementary File S3). Then, an osteogenesis-associated gene list consisting of 399 genes (including ossification, osteoblast differentiation, and osteoblast proliferation) was obtained from the MsigDB database (Supplementary File S4). The results of intersection analysis in Table 1 showed that 11 osteogenesis-associated genes were co-expressed with lncRNA GAS5 (PTPN11, HNRNPC, CCDC47, DHX9, HSPE1, FBL, RDH14, TWIST2, PHB, RSL1D1, DHX36). In addition, we performed lncRNA–mRNA interaction analysis and lncRNA–protein interaction analysis based on the RNAInter database (http://www.rnainter.org/) and the AnnoLnc2 database (http://annolnc.gao-lab.org/) to investigate the potential functional route of GAS5. In total, 1870 lncRNA GAS5–mRNA interaction pairs and 523 lncRNA GAS5–protein interaction pairs were obtained (Supplementary Files S5, S6).

The expression data of lncRNA GAS5 in GSE113253 revealed its expression trend during the osteogenic differentiation of BMSCs. However, analysis of the differences in lncRNA GAS5 expression levels between osteogenic differentiation induction groups and normal proliferation medium groups can make these results more convincing. Therefore, in this study, we conducted an ALP staining assay and RT-qPCR on days 3 and 7 after the induction of BMSC osteogenic differentiation to verify the association between lncRNA GAS5 and the development of the BMSC osteogenic phenotype. ALP activity increased gradually from day 3 to day 7, with the ALP activity in the osteogenic differentiation induction groups being significantly higher than that in the normal proliferation medium groups (Figures 3A,B). The RT-qPCR results in Figure 3C indicated that the expression of osteogenic marker RUNX2 is significantly upregulated in osteogenic differentiation group on day 7 (p = 0.002). Finally, lncRNA GAS5 expression levels in the osteogenic differentiation induction and normal control groups were tracked at these two separate osteogenic differentiation points. The results showed that the expression level of lncRNA GAS5 was significantly upregulated (p < 0.05) in the osteogenic differentiation induction groups compared to the normal proliferation medium groups (Figures 3D, E).

ONFH rat models were constructed, and the main bone tissue parameters of the femoral heads were measured using micro-CT 9 weeks after induction. The subchondral trabecular bone in the weight-bearing area of the femoral head was thinner and sparser according to the results of micro-CT and hematoxylin–eosin staining (Figures 4A–F). In addition, the BV/TV (%) and Tb.Th (mm) values were significantly lower in the ONFH group than in the control group. The value of Tb.Sp (mm) significantly increased in the ONFH group (Figures 4G–I). These parameters showed the changes in the bone mass level and trabecular structure, indicating the development of ONFH. RT-qPCR was performed to evaluate the expression of lncRNA GAS5 in the femoral head subchondral bone tissue. The results showed that lncRNA GAS5 expression was significantly decreased in the ONFH group compared to the control group (Figure 4J), consistent with our findings in clinical femoral head subchondral bone tissue. Therefore, lncRNA GAS5 can act as an osteogenic biomarker for ONFH development.