LncRNA‐mRNA expression profiles and functional networks in osteoclast differentiation

Abstract Human osteoclasts are differentiated from CD14+ monocytes and are responsible for bone resorption. Long non‐coding RNAs (lncRNAs) have been proved to be significantly involved in multiple biologic processes, especially in cell differentiation. However, the effect of lncRNAs in osteoclast differentiation is less appreciated. In our study, RNA sequencing (RNA‐seq) was used to identify the expression profiles of lncRNAs and mRNAs in osteoclast differentiation. The results demonstrated that expressions of 1117 lncRNAs and 296 mRNAs were significantly altered after osteoclast differentiation. qRT‐PCR assays were performed to confirm the expression profiles, and the results were almost consistent with the RNA‐seq data. GO and KEGG analyses were used to predict the functions of these differentially expressed mRNA and lncRNAs. The Path‐net analysis demonstrated that MAPK pathway, PI3K‐AKT pathway and NF‐kappa B pathway played important roles in osteoclast differentiation. Co‐expression networks and competing endogenous RNA networks indicated that ENSG00000257764.2‐miR‐106a‐5p‐TIMP2 may play a central role in osteoclast differentiation. Our study provides a foundation to further understand the role and underlying mechanism of lncRNAs in osteoclast differentiation, in which many of them could be potential targets for bone metabolic disease.


| INTRODUC TI ON
Human osteoclasts are differentiated from CD14 + monocytes under the induction of macrophage colony stimulating factor (M-CSF) and receptor activator of NF-κB ligand (RANKL). 1,2 Tartrateresistant acid phosphatase (TRAP) and cathepsin K (CTSK) are the osteoclast differentiation markers and are in charge of bone resorption. 3 During osteoclast differentiation, mononuclear cells fuse together and transform into multinucleated osteoclasts. 4,5 The actin cytoskeleton of osteoclasts is rearranged and a tight junction between the bone surface and basal membrane arises to form a sealing zone. 3,6 This sealing zone is acidified by hydrogen ions, and then, TRAP and CTSK are secreted to resorb bone and an absorption pit is formed. 3,7 In normal conditions, osteoclasts and osteoblasts maintain the balance of bone metabolism. 8 Unfortunately, abnormal differentiation of osteoclasts occurs in many skeletal diseases such as giant cell tumours of the bone, multiple myeloma, osteoporosis and inflammatory rheumatic diseases. [9][10][11][12] Over the past decades, increasing numbers of regulating factors involved in the differentiation of osteoclasts have been discovered, including cytokines, molecules of signalling pathways, natural compounds and others. 13,14 However, the epigenetic regulation of osteoclast differentiation that leads to the dysregulation of osteoclasts is less well known.
Epigenetics analyse the heritable changes in gene expression when the nucleotide sequence of a gene does not change. 15,16 Long non-coding RNAs (lncRNAs) belong to a kind of epigenetic regulation. Transcripts of lncRNAs are usually >200 nt and they do not encode proteins. [17][18][19] Much of the non-coding genome has historically been regarded as junk DNA. 20 However, with advances in sequencing technology, an in-depth examination of the non-coding genome has been achieved; this revealed that the rate of non-coding genome is more than 98%. 20 LncRNAs are now known to participate in many biological processes, such as cell proliferation, differentiation and apoptosis. [21][22][23][24] However, little is known about the roles of lncRNAs playing in osteoclast differentiation.
LncRNAs carry out their functions through diverse mechanisms involving miRNAs, mRNAs and proteins. 25,26 miRNAs have emerged as important regulators of osteoclast differentiation. 27,28 Recent studies demonstrate that lncRNAs serve as competing endogenous RNAs (ceRNAs) by competing for binding to miRNAs. [29][30][31] Studying this novel RNA crosstalk will provide significant insight into gene regulatory networks of osteoclast differentiation.
In this study, we explored the lncRNA-mRNA expression profiles in human osteoclast differentiation by RNA-seq and confirmed the expression profiles by qRT-PCR assays. Moreover, GO and KEGG analyses were used to predict potential cellular functions of these differentially expressed mRNA and lncRNAs. The Path-net analysis, co-expression networks and ceRNA networks were applied to predict gene regulatory networks of osteoclast differentiation. Our study provides a new perspective on the regulation of osteoclast differentiation by lncRNAs.

| Cell isolation
Peripheral blood was acquired from healthy donors with the approval of the Ethics Committee of Sun Yat-sen Memorial Hospital.

| Tartrate-resistant acid phosphatase staining
CD14 + monocytes (2.5 × 10 5 /cm 2 ) were plated onto 24-plate wells and cultured with 25 ng/mL M-CSF and 50 ng/mL RANKL to induce osteoclasts. On day 9, cells were stained for tartrate-resistant acid phosphatase (TRAP) activity with the leucocyte acid phosphatase kit (Sigma, St. Louis, MO, USA) according to the manufacturer's protocol. TRAP-positive cells which had at least three nuclei were considered osteoclasts.

| F-actin assay
To detect podosomes in osteoclasts, we carry out F-actin assays.
After fixing with 4% paraformaldehyde for 5 minutes, cells were washed three times with PBS. Then, cells were stained with FITCconjugated phalloidin (Sigma) and 4′,6-diamidino-2-phenylindole (DAPI) for 40 minutes. In the end, cells were washed with PBS thoroughly and observed with an Axio Observer fluorescence microscope (Carl Zeiss, Jena, Germany).

| Real-time quantitative reverse transcriptionpolymerase chain reaction (qRT-PCR)
Cells were treated with TRIzol (Invitrogen, Carlsbad, CA, USA) and total RNA was extracted. Then, 1 µg of RNA was transcribed to cDNA using a PrimeScript RT reagent kit (TaKaRa, Dalian, China) according to the manufacturer's protocol. The qRT-PCR assay was carried out on a LightCycler ® 480 PCR system (Roche, Basel, Switzerland) using SYBR Premix Ex Taq (TaKaRa). The procedure for the quantitative-PCR assay was 95°C for 30 seconds, 40 cycles at 95°C for 5 seconds and 60°C for 20 seconds. Each sample was divided into triplicate and the mean expression of each gene was calculated using the 2 −ΔΔCt method. Primers for all the genes are available in Table S1.

| Western Blot analysis
After lysing in RIPA buffer (Sigma), cells were centrifuged at 10 000 g for 30 minutes. Then, proteins in the supernatant were quantified with a BCA assay kit (Sigma), separated by 10% sodium dodecyl sulphate-polyacrylamide gel electrophoresis and elec-

| RNA-seq analysis
Total RNA was extracted from cells after adding TRIzol (Invitrogen), and ribosomal RNA was removed using the Ribo-Zero™ kit

| Co-expression network construction
The co-expression network was established by calculating the Pearson correlation coefficient and P-value between multiple genes.
In this study, the transcripts were filtered using a COR of >0.85 and a P-value of <0.05. The co-expression network was illustrated using Cytoscape software (available online: https://cytos cape.org).

| Competing endogenous RNA network construction
The lncRNAs and mRNAs were selected to predict miRNA targets using the miRbase. Then the miRNAs obtained from the predictions were screened with the miRanda and TargetScan programs.
Afterwards, lncRNAs and mRNAs possessing miRNA recognition elements (MREs) for the targeted miRNAs were predicted using RNA22. The competitive endogenous RNA (ceRNA) network was established and illustrated using Cytoscape software.

| Statistical analysis
Statistical analysis was conducted with SPSS 22.0 software (Chicago, IL, USA). Student's t test was used for the comparisons between two groups. A P-value of <0.05 was considered significant.

| Osteoclast differentiation and identification
To induce osteoclast differentiation, we cultured human CD14 + monocytes with 50 ng/mL RANKL and 25 ng/mL M-CSF. TRAP staining was performed on day 9 and TRAP + osteoclast-like cells were found after culture with M-CSF and RANKL ( Figure 1A). F-actin assays indicated that these osteoclast-like cells can cause bone resorption, and this conclusion is further confirmed by pit formation assays ( Figure 1A). As shown in Figure 1B,C, osteoclasts expressed higher levels of TRAP and CTSK compared to monocytes, as determined by qRT-PCR and Western Blot analyses.

| Expression profiles of lncRNAs and mRNAs in osteoclast differentiation
To identify possible lncRNAs and mRNAs participated in the dif-  Table 1 and Table 2, respectively.

| Validation of the expression profiles by qRT-PCR
To confirm the data of RNA-seq, 10 lncRNAs and mRNAs from above tables were selected for the qRT-PCR analysis. Their expressions in undifferentiated or differentiated osteoclasts are shown in Figure 3A

| GO analysis
GO analysis showed that a total of 2442 GO terms in Biological Process were differentially expressed in osteoclast differentiation, and the three most prevalent GO terms were cellular process, single-organism process and single-organism cellular process. A total of 280 GO terms in Cellular Component were differentially expressed in osteoclast differentiation, and the three most prevalent GO terms were cell, cell part and organelle. A total of 374 GO terms in Molecular Function were differentially expressed in osteoclast differentiation, and the three most prevalent GO terms were binding, protein binding and heterocyclic compound binding. The top 10 significantly enriched GO terms are presented in Figure 4A.

| Pathway analysis
The KEGG pathway analysis showed that 147 pathways achieved statistically significant enrichment. The most enriched pathways included the MAPK signalling pathway and metabolic pathways ( Figure 4B), suggesting that these signalling pathways may play important roles in osteoclast differentiation.
The interaction of all significantly enriched pathways was done as a Path-net ( Figure 5). We found that the MAPK signalling pathway obtained the highest number of interactions with other pathways, highlighting its vital role in osteoclast differentiation ( Figure 5). Moreover, the PI3K-AKT signalling pathway and the NF-κB signalling pathway also played a significant role in osteoclast differentiation ( Figure 5). To illustrate the correlation of monocytes and mature osteoclasts, gene set enrichment analysis (GSEA) was performed.
GSEA plots showed that genes differentially expressed between monocytes and osteoclasts were involved in osteoclast differentiation ( Figure S1A). GSEA also revealed that these differentially expressed genes were associated with the MAPK pathway ( Figure S1B), supporting the previous results of GO and KEGG analyses.

| Construction of the lncRNA-mRNA coexpression network
To explore the potential interactions between mRNAs and lncRNAs, a lncRNA-mRNA co-expression network was established. 79 lncRNAs were selected from the differentially expressed lncRNAs and further built as 140 pairs of co-expression lncRNA-mRNA ( Figure 6). In the co-expression network, the one that obtained the highest number of interactions was lnRNA ENSG00000257764.2 ( Figure 6). Also, we found that ENSG00000257764.2 interacted with TIMP2 ( Figure 6).
Top 5 in the lncRNA-mRNA interaction networks were showed in Table S2. The expressions of lncRNAs and mRNAs were confirmed by qRT-PCR, and their correlations were analysed by Pearson correlation coefficient. The results showed that lncRNA PVT1 had a strong negative correlation with FCN1, while the other four had a strong positive correlation ( Figure S2). These results were almost consistent with the data from the lncRNA-mRNA interaction networks.

| Construction of the competing endogenous RNA network
A ceRNA network was constructed to display the interactions between miRNAs, mRNAs and lncRNAs (Figure 7). We found that hsa-miR-106a-5p, hsa-miR-548c-3p and hsa-miR-6133 were en-

| D ISCUSS I ON
In this study, through RNA-seq analysis, a number of lncRNAs and mRNAs differentially expressed in osteoclast differentiation were identified. Then the cellular events and biological pathways in osteoclast differentiation were identified with GO, KEGG and Pathnet analysis. Furthermore, co-expression networks presented the interactions between lncRNAs and mRNAs, and the core regulatory factors in osteoclast differentiation. Finally, a ceRNA network was established to display the interactions between lncRNAs and miR-NAs. We found that ENSG00000257764.2-miR-106a-5p-TIMP2 may play a central role in osteoclast differentiation. Our study provides a basis to further understand the role and underlying mechanism of lncRNAs in osteoclast differentiation and many of them could be potential targets for bone metabolic disease.   to find the significant pathways involved in osteoclast differentiation. Recently, some researchers have found that non-canonical pathways also take part in osteoclast differentiation. 53 He et al showed that the p38 pathway in myeloma cells regulated osteoclast differentiation. 54 In fact, p38, together with ERK and JNK, belong to the MAPK signalling pathway. 55 In our study, as shown in Regarding the regulatory functions of lncRNAs, we constructed the lncRNA-mRNA co-expression network and found that lncRNA ENSG00000257764.2 obtained the highest number of interactions and interacted with TIMP2. TIMP2 are inhibitors of matrix metalloproteinases (MMPs), which are a group of endopeptidases that regulate the formation and function of osteoclast. 56 In the qRT-PCR results mentioned above, we found that ENSG00000257764.2 was significantly down-regulated in the differentiation of osteoclast.
These results indicated that ENSG00000257764.2 may play a central role in regulating osteoclast differentiation by regulating the expression of TIMP2. However, its significant roles have not been studied in osteoclast differentiation and further experiments are required to test this hypothesis.
LncRNAs carry out their functions through a variety of mechanisms including the ceRNA hypothesis. 29,57 The ceRNA hypothesis suggests that transcripts with common miRNA binding sites compete for the post-transcriptional regulation of genes. 29,57 This hypothesis has attracted extensive attention these years as an important functional mechanism of lncRNAs. However, there are no studies that link miRNA and lncRNA in osteoclast differentiation. In our study, a ceRNA network was established for the differentially expressed lncRNAs. We found that hsa-miR-106a-5p was enriched in the ceRNA network, and interacted not only with ENSG00000257764.2, but also TIMP2. Therefore, we hypothesized that ENSG00000257764.2-miR-106a-5p-TIMP2 may be a crucial axis in the regulation network of osteoclast differentiation since it interacts with the greatest number of differentially expressed mRNAs and lncRNAs. However, the role of ENSG00000257764.2-miR-106a-5p-TIMP2 in osteoclast differentiation remains unknown and further studies are necessary.
In conclusion, we investigated the lncRNA-mRNA expression profiles in osteoclast differentiation through RNA-seq and identified the potential regulatory mechanisms via bioinformatic analyses. We aim to reveal the role of lncRNAs in osteoclast differentiation and osteoclast-related bone disorders. Our results may help illuminate the mechanism of osteoclast differentiation and provide potential targets for osteoclast-related bone disorders. However, our study also has some limitations, such as a lack of functional research of these differentially expressed lncRNAs. More studies are needed to explore the roles of these differentially expressed lncRNAs in osteoclast differentiation.

F I G U R E 6
Construction of the lncRNA-mRNA co-expression network. A node with a yellow ring indicates lncRNA, and a node without a yellow ring indicates mRNA. Up-regulated lncRNAs and mRNAs are shown in blue and down-regulated mRNAs and lncRNAs are shown in red

ACK N OWLED G EM ENTS
Our study was financially supported by the National Natural Science

CO N FLI C T O F I NTE R E S T
The authors confirm that there are no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author on request.