Identification of potential pathogenic genes associated with osteoporosis.

OBJECTIVES: Osteoporosis is a chronic disease. The aim of this study was to identify key genes in osteoporosis.
METHODS: Microarray data sets GSE56815 and GSE56814, comprising 67 osteoporosis blood samples and 62 control blood samples, were obtained from the Gene Expression Omnibus database. Differentially expressed genes (DEGs) were identified in osteoporosis using Limma package (3.2.1) and Meta-MA packages. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses were performed to identify biological functions. Furthermore, the transcriptional regulatory network was established between the top 20 DEGs and transcriptional factors using the UCSC ENCODE Genome Browser. Receiver operating characteristic (ROC) analysis was applied to investigate the diagnostic value of several DEGs.
RESULTS: A total of 1320 DEGs were obtained, of which 855 were up-regulated and 465 were down-regulated. These differentially expressed genes were enriched in Gene Ontology terms and Kyoto Encyclopedia of Genes and Genomes pathways, mainly associated with gene expression and osteoclast differentiation. In the transcriptional regulatory network, there were 6038 interactions pairs involving 88 transcriptional factors. In addition, the quantitative reverse transcriptase-polymerase chain reaction result validated the expression of several genes (VPS35, FCGR2A, TBCA, HIRA, TYROBP, and JUND). Finally, ROC analyses showed that VPS35, HIRA, PHF20 and NFKB2 had a significant diagnostic value for osteoporosis.
CONCLUSION: Genes such as VPS35, FCGR2A, TBCA, HIRA, TYROBP, JUND, PHF20, NFKB2, RPL35A and BICD2 may be considered to be potential pathogenic genes of osteoporosis and may be useful for further study of the mechanisms underlying osteoporosis.Cite this article: B. Xia, Y. Li, J. Zhou, B. Tian, L. Feng. Identification of potential pathogenic genes associated with osteoporosis. Bone Joint Res 2017;6:640-648. DOI: 10.1302/2046-3758.612.BJR-2017-0102.R1.

The objective of this study was to identify key genes in osteoporosis.

Genes including VPS35, FCGR2A, TBCA, HIRA, TYROBP, JUND, PHF20, NFKB2, RPL35A and BICD2 were considered to play potential roles in the pathology of osteoporosis.

Bioinformatics analysis was used to identify and study the biological function of differentially expressed genes.

The validation sample for quantitative reverse transcriptase-polymerase chain reaction was small and it is necessary to collect larger samples for further validation.

Introduction

Osteoporosis, characterised by the impairment of bone microarchitecture and the loss of bone mass and strength, has become an important clinical problem in ageing populations.[1,2] The spine is the most common site of osteoporotic fractures, followed by the hip, forearm and proximal humerus.[3] Osteoporosis is characterised by deterioration of the microstructure of bone, specifically at trabecular sites, which leads to pain, deformity, disability and possibly death.[4,5] A variety of factors contribute to the development of osteoporosis, such as genetic variants, gender, steroid production, age, lifestyle and environment.[6-8] Additionally, low calcium intake, cigarette-smoking and intake of excessive alcohol may be secondary causes.[4,9] Generally, osteoporosis is considered a silent disease because it is asymptomatic until a fracture occurs. Recently, the treatment of osteoporosis has been mainly pharmaceutical, but treatment may not be satisfactory due to its time-consuming nature and high cost, as well as the side effects of the drugs.

As Sims et al[10] have reported, a number of studies have explored the mechanism of osteoporosis. The members of the Wnt signalling pathway, such as Wnt3a, low-density lipoprotein receptor-related protein 5, secreted frizzled-related protein 1 and sclerostin, have been suggested to be related to variation in bone mineral density (BMD).[10] It is known that osteoporosis is defined clinically by measuring the BMD with heritability estimates of 0.5 to 0.9.[11] This further demonstrated that BMD is an important clinical marker in osteoporosis. The underlying aetiology of osteoporosis is still not fully understood and the identification of novel therapeutic target genes for osteoporosis is needed.

It is worth mentioning that microarray data analysis is available to identify vital genes and gene regulatory networks associated with disease.[12] In this study, we downloaded the data sets GSE56815 and GSE56814 from the Gene Expression Omnibus database, and identified the differentially expressed genes (DEGs) in blood samples obtained from osteoporosis patients, followed by Gene Ontology (GO)[13] and Kyoto Encyclopedia of Genes and Genomes (KEGG)[14] enrichment analyses and interaction network construction between transcription factors (TFs) and DEGs. We aimed to find involvement of key genes in osteoporosis which may be potential modulators in the pathology of osteoporosis.

Methods

Microarray data

The two microarray data sets (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE56815 and https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE56814) were downloaded using the GEOquery package with R version (The R Foundation for Statistical Computing, Vienna, Austria) from the Gene Expression Omnibus (GEO) database. These two microarray data sets were from the blood samples of 67 osteoporosis patients (based on World Health Organization (WHO) criteria) and 62 control groups. It is reported that the WHO recommends the use of BMD of the spine and proximal femur measured by double energy X-ray absorption as the benchmark to diagnose osteoporosis and its severity.[15] Therefore, diagnosis of osteoporosis depended on the BMD in these two datasets; these osteoporosis patients all had low BMD. In addition, the patients with osteoporosis did not receive any anti-osteoporosis medication or other medication that may affect bone metabolism. Detailed information of two data sets is shown in Table I.

Table I.

Two Gene Expression Omnibus (GEO) datasets of osteoporosis

Gene type	GEO ID	Platform	Sample count (control group: patients with osteoporosis)	Notes
mRNA	GSE56815	GPL96 [HG-U133A] Affymetrix Human Genome U133A Array	80 (40:40)	Liu YZ, USA, 2016[16]
mRNA	GSE56814	GPL5175 [HuEx-1_0-st] Affymetrix Human Exon 1.0 ST Array [transcript (gene) version]	49 (22:27)	Liu YZ, USA, 2016[16]

mRNA, messenger ribonucleic acid

Screening of DEGs

Normalised microarray data were downloaded and gene expression value was calculated as the mean value of its corresponding probe values. Limma package[17] and Meta-MA[18] package were used to identify the DEGs between osteoporosis and control group; p-values and false discovery rates (FDR) were further calculated. FDR < 0.01 was selected as the threshold for screening DEGs.

Functional analyses of DEGs

GO and KEGG enrichment analyses were carried out for the identified DEGs using GeneCodis3 (http://www.genecodis.cnb.csic.es/analysis).[19-21] Significant GO terms and KEGG pathways were identified according to the threshold of p < 0.05.

Network construction of DEGs and TFs

To gain deeper insight into the molecular functions of DEGs, the gene regulatory relationships between DEGs and TFs were selected based on human TF binding sites data and genetic coordinate position information, which were available at the UCSC ENCODE Genome Browser.[22] The identified TFs were considered to be associated with DEGs and the regulatory network between DEGs and TFs was visualised using Cytoscape software (Cytoscape Corporation, San Diego, California).

Quantitative reverse transcriptase-polymerase chain reaction (qRT-PCR) in vitro

Four women diagnosed with osteoporosis with low BMD were enrolled in this study, as were three women with high BMD but no diagnosis of osteoporosis. The clinical information of osteoporosis patients is shown in Table II. All blood samples were collected for further qRT-PCR experimentation. All participating individuals provided informed consent with the approval of the ethics committee.

Table II.

The clinical information of osteoporosis patients

Patient number	Gender	Age (yrs)	Weight (kg)	Height (cm)	Smoking history	Drinking history	Family history	Menopause status	BMD	Diagnostic method	Symptom
1	Female	82	50	160	No	No	No	Yes	-2.57	DXA	Elderly osteoporosis
2	Female	81	45	155	No	No	No	Yes	-2.66	DXA	Elderly osteoporosis
3	Female	72	55	163	No	No	No	Yes	-2.93	DXA	Elderly osteoporosis
4	Female	66	60	165	No	No	No	Yes	-2.82	DXA	Elderly osteoporosis
5	Female	69	53	158	No	No	No	Yes	-0.49	DXA	Normal
6	Female	79	46	162	No	No	No	Yes	-0.58	DXA	Normal
7	Female	73	57	157	No	No	No	Yes	-0.43	DXA	Normal

BMD, bone mineral density; DXA, double energy X-ray absorption

Total ribonucleic acid (RNA) of the blood samples was extracted using TRIzol Reagent (Invitrogen, Carlsbad, California), in accordance with the manufacturer’s protocols. Then 1 ug RNA was applied to synthesise DNA by SuperScript III Reverse Transcriptase (Invitrogen) and qRT-PCR was performed in an ABI 7500 Real-time PCR system (Invitrogen) with SYBR Green PCR Master Mix (Invitrogen). Glyceraldehyde 3-phosphate dehydrogenase was used as internal control and all reactions were performed in triplicate. Relative gene expressions were analysed by the 2-△△Ct method.

Receiver operating characteristic (ROC) analyses

By using the pROC package[23] in R language, we performed the ROC analyses to assess the diagnostic value of DEGs. The area under the curve (AUC) under the binomial exact confidence interval was calculated and the ROC curve was generated.

Results

DEGs identification

In this study, we identified 1320 DEGs with significantly altered expression in osteoporosis blood samples, of which 855 were up-regulated and 465 were down-regulated compared with the control groups. The top 20 (ten up-regulated and ten down-regulated) DEGs abbreviations are defined in Table III and presented in Table IV. The results indicated that the expression pattern of these DEGs could observably distinguish the osteoporosis samples from control groups.

Table III.

Differentially expressed gene abbreviations

Abbreviation	Full name
ALKBH1	alkB homolog 1, histone H2A dioxygenase
BICD2	BICD cargo adaptor 2
CDKN2D	cyclin dependent kinase inhibitor 2D
DPP8	dipeptidyl peptidase 8
FCGR2A	Fc fragment of IgG receptor IIa
HIRA	histone cell cycle regulation
IER2	immediate early response 2
JUND	JunD proto-oncogene, AP-1 transcription factor subunit
METTL4	methyltransferase like 4
NBEAL2	neurobeachin like 2
NFKB2	nuclear factor kappa B subunit 2
NIF3L1	NGG1 interacting factor 3 like 1
NIT2	nitrilase family member 2
PAF1	PAF1 homolog, Paf1/RNA polymerase II complex component
PHF20	PHD finger protein 20
POGLUT1	protein O-glucosyltransferase 1
RPL35A	ribosomal protein L35a
SAP130	Sin3A associated protein 130
SH3GLB2	SH3 domain containing GRB2 like, endophilin B2
TBCA	tubulin folding cofactor A
TYROBP	TYRO protein tyrosine kinase binding protein
VPS35	VPS35, retromer complex component

Table IV.

The top 20 differentially expressed genes (DEGs) in osteoporosis

ID	Symbol	Combined ES	p-value	FDR	Regulation
55737	VPS35	1.41E+00	5.58E-13	4.37E-09	Up
56983	POGLUT1	1.40E+00	7.55E-13	4.37E-09	Up
54878	DPP8	1.39E+00	1.32E-12	5.09E-09	Up
60491	NIF3L1	1.31E+00	1.32E-11	3.81E-08	Up
51230	PHF20	1.25E+00	9.79E-11	2.03E-07	Up
56954	NIT2	1.26E+00	1.05E-10	2.03E-07	Up
64863	METTL4	1.21E+00	1.77E-10	2.92E-07	Up
79595	SAP130	1.21E+00	3.95E-10	5.71E-07	Up
2212	FCGR2A	1.16E+00	1.32E-09	1.70E-06	Up
8846	ALKBH1	1.17E+00	2.41E-09	2.79E-06	Up
4791	NFKB2	-1.06E+00	1.97E-08	6.92E-06	Down
6165	RPL35A	-1.07E+00	2.51E-08	7.98E-06	Down
23299	BICD2	-1.05E+00	2.55E-08	7.98E-06	Down
6902	TBCA	-1.07E+00	2.88E-08	8.35E-06	Down
56904	SH3GLB2	-1.03E+00	3.04E-08	8.58E-06	Down
7290	HIRA	-1.02E+00	3.50E-08	9.22E-06	Down
54623	PAF1	-1.03E+00	3.76E-08	9.68E-06	Down
23218	NBEAL2	-1.05E+00	4.25E-08	1.07E-05	Down
9592	IER2	-1.02E+00	5.34E-08	1.19E-05	Down
1032	CDKN2D	-1.07E+00	6.00E-08	1.24E-05	Down

ES, effect size; FDR, false discovery rate

GO and KEGG analyses of DEGs

Among 1320 DEGs, 1251 genes were recognised and significantly involved in different GO terms and KEGG pathways. The top 15 enriched GO terms, such as biological process, molecular function and cellular component of the identified DEGs are shown in Table V. The results showed that these DEGs were mainly involved in the GO terms associated with apoptosis processes, gene expression and signal transduction. On the other hand, the results of KEGG analysis revealed that a total of 15 pathways were enriched; for example, Ubiquitin mediated proteolysis and osteoclast differentiation (Table VI).

Table V.

Top 15 enriched Gene Oncology (GO) terms in osteoporosis

GO ID	GO term	Genes (n)	FDR
	Biological process
GO:0006915	Apoptotic process	75	2.26E-17
GO:0010467	Gene expression	56	2.31E-14
GO:0007165	Signal transduction	105	2.99E-14
GO:0015031	Protein transport	50	9.06E-11
GO:0006355	Regulation of transcription, DNA-dependent	119	1.16E-10
GO:0006810	Transport	61	4.45E-10
GO:0048011	Nerve growth factor receptor signalling pathway	33	1.19E-09
GO:0016070	RNA metabolic process	36	1.73E-09
GO:0016032	Viral reproduction	41	2.23E-09
GO:0044419	Interspecies interaction between organisms	41	2.24E-09
GO:0008380	RNA splicing	36	3.86E-09
GO:0042981	Regulation of apoptotic process	31	4.62E-09
GO:0000278	Mitotic cell cycle	38	8.32E-09
GO:0016071	mRNA metabolic process	32	8.35E-09
GO:0051437	Positive regulation of ubiquitin-protein ligase activity involved in mitotic cell cycle	18	9.82E-09
	Molecular function
GO:0005515	Protein binding	446	2.10E-93
GO:0000166	Nucleotide binding	180	7.43E-24
GO:0003677	DNA binding	142	4.86E-16
GO:0005524	ATP binding	126	5.21E-16
GO:0046872	Metal ion binding	193	1.71E-14
GO:0003723	RNA binding	69	5.81E-14
GO:0016787	Hydrolase activity	90	1.13E-13
GO:0016740	Transferase activity	59	2.82E-09
GO:0008270	Zinc ion binding	129	3.94E-09
GO:0003676	Nucleic acid binding	68	2.21E-08
GO:0003700	Sequence-specific DNA binding transcription factor activity	71	1.04E-07
GO:0016874	Ligase activity	40	1.10E-07
GO:0003824	Catalytic activity	37	3.27E-06
GO:0008233	Peptidase activity	44	4.47E-06
GO:0042803	Protein homodimerisation activity	45	5.19E-06
	Cellular component
GO:0005737	Cytoplasm	464	3.22E-78
GO:0005634	Nucleus	457	4.04E-71
GO:0005829	Cytosol	246	6.44E-58
GO:0005739	Mitochondrion	155	1.62E-33
GO:0005654	Nucleoplasm	116	6.90E-31
GO:0016020	Membrane	295	3.54E-30
GO:0005730	Nucleolus	141	3.18E-24
GO:0005622	Intracellular	156	2.02E-18
GO:0005794	Golgi apparatus	97	7.96E-18
GO:0005783	Endoplasmic reticulum	99	1.43E-17
GO:0005789	Endoplasmic reticulum membrane	73	4.91E-16
GO:0016021	Integral to membrane	259	8.22E-14
GO:0005743	Mitochondrial inner membrane	43	2.07E-12
GO:0005625	Soluble fraction	44	1.45E-09
GO:0000139	Golgi membrane	46	1.55E-09

FDR, false discovery rate; RNA, ribonucleic acid; mRNA, messenger RNA; ATP, adenosine triphosphate

Table VI.

Top 15 enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) terms in osteoporosis

KEGG ID	KEGG term	Count	FDR	Genes
Hsa04120	Ubiquitin mediated proteolysis	24	2.83E-08	UBE2A,MDM2,HERC1,UBE2Z,UBE4B,FZR1,DET1,CDC23,FBXO4,ANAPC2,ANAPC13,CDC27,UBE2W,UBE2G1,NEDD4,SOCS3,UBE2D3,UBE2E3,UBE4A,VHL,CUL1,UBE2D4,UBE2N,HERC2
Hsa04380	Osteoclast differentiation	22	1.18E-07	TYROBP,NFKB2,JUND,IL1R1,TYK2,NFKBIA,CREB1,TNFRSF1A,PPP3CC,LILRB3,SOCS3,FOSB,CYBB,FCGR2A,NFKB1,LILRA2,JUN,MITF,MAPK14,JUNB,PIK3CG,IFNAR1
Hsa05016	Huntington's disease	26	1.32E-07	POLR2J,NDUFA8,ATP5G1,COX6A2,POLR2H,POLR2B,NDUFC1,CREB1,COX5B,NRF1,GRIN2B,PLCB3,GNAQ,TBP,UQCRC2,NDUFA7,UQCR10,NDUFA13,BAX,NDUFB1,SDHD,PPARGC1A,NDUFB5,NDUFS3,SDHA,ATP5B
Hsa00190	Oxidative phosphorylation	22	1.45E-07	NDUFA8,ATP5G1,COX6A2,NDUFC1,COX10,COX5B,ATP5L,ATP5I,ATP6V1E1,ATP6V0E2,UQCRC2,NDUFA7,UQCR10,NDUFA13,NDUFB1,SDHD,NDUFB5,ATP6V1B2,NDUFS3,SDHA,LHPP,ATP5B
Hsa04141	Protein processing in endoplasmic reticulum	24	1.76E-07	SEC63,SSR2,UBE4B,CANX,ATF6,RPN2,SSR4,PPP1R15A,STT3A,UBE2G1,UBE2D3,UBE2E3,RRBP1,PDIA6,CUL1,BAX,UBE2D4,HERPUD1,SEL1L,DNAJC10,EIF2S1,LMAN2,MBTPS2,SIL1
Hsa04142	Lysosome	20	5.87E-07	NPC2,SCARB2,CTSC,GLA,LAPTM4A,MANBA,CLN5,FUCA1,AGA,LIPA,ASAH1,NPC1,PPT1,ARSB,HEXB,DNASE2,CTSO,IDS,IGF2R,LAPTM5
Hsa05010	Alzheimer's disease	23	9.66E-07	NDUFA8,ATP5G1,COX6A2,NDUFC1,IDE,ATF6,TNFRSF1A,COX5B,PPP3CC,GRIN2B,ITPR2,PLCB3,GNAQ,UQCRC2,NDUFA7,UQCR10,NDUFA13,NDUFB1,SDHD,NDUFB5,NDUFS3,SDHA,ATP5B
Hsa00020	Citrate cycle (TCA cycle)	10	1.50E-06	CS,DLAT,IDH3A,SDHD,MDH1,DLST,ACO2,SDHA,PDHB,SUCLG2
Hsa03040	Spliceosome	19	1.54E-06	SF3A3,WBP11,PRPF6,RBM17,DHX38,PPIE,TXNL4A,PRPF38B,EIF4A3,PRPF3,RBM25,SART1,SRSF3,CDC5L,CRNKL1,ACIN1,LSM3,SNRNP40,SNRPB2
Hsa03050	Proteasome	11	2.77E-06	SHFM1,PSMC2,PSMD11,PSMC4,PSME3,PSMB3,PSMD1,PSMD2,PSME1,PSMD8,PSMA4
Hsa00970	Aminoacyl-tRNA biosynthesis	10	2.27E-05	IARS,DARS2,RARS,WARS,DARS,GARS,NARS2,NARS,WARS2,EPRS
Hsa05200	Pathways in cancer	31	2.43E-05	VEGFB,MDM2,SOS1,NFKB2,ARNT,HRAS,CCNE1,NFKBIA,E2F1,CCND1,PTEN,TPM3,RALB,BCR,RALBP1,CTBP1,CTBP2,FZD2,NFKB1,TFG,JUN,VHL,FLT3,MAP2K2,MITF,BAX,VEGFC,SOS2,PIK3CG,MLH1,DVL1
Hsa03013	RNA transport	19	2.53E-05	EIF2S3,RPP30,EIF3H,NUP98,PRMT5,POM121,EIF3D,NUP43,NUP155,NXT2,EIF3G,RPP38,EIF4E2,EIF4A3,SNUPN,POP4,EIF2S1,ACIN1,EIF3E
Hsa05220	Chronic myeloid leukemia	13	3.34E-05	MDM2,SOS1,HRAS,NFKBIA,E2F1,CCND1,BCR,CTBP1,CTBP2,NFKB1,MAP2K2,SOS2,PIK3CG
Hsa04144	Endocytosis	22	3.68E-05	MDM2,CXCR4,SH3GLB2,ADRB2,RAB11FIP5,HRAS,FOLR2,ADRBK1,RAB11FIP3,GRK6,TSG101,ARAP1,NEDD4,HGS,PARD6A,SRC,RAB5C,RAB5A,ARFGAP1,CHMP2A,PDCD6IP,RNF41

FDR, false discovery rate; TCA, tricarboxylic acid; tRNA, trnaser ribonucleic acid; RNA, ribonucleic acid

Transcriptional regulatory relationships between DEGs and TFs

Regulatory relationships were predicted between top 20 DEGs (ten up-regulated and ten down-regulated) and TFs, and the regulatory network was established and visualised by Cytoscape software (Cytoscape Corporation) (Fig. 1). In this network, there were 6038 interactions pairs involving 88 TFs. The top seven TFs covering the most downstream genes were FOXD3, Nkx2-5, Pax-4, Oct-1, HNF-4, Pax-6 and COMP1.

Fig. 1

Transcription factors (TFs)-differentially expressed genes interaction networks. Rectangles and ellipses represent the TFs and target genes, respectively. The red and green colors represent up-regulation and down-regulation, respectively.

qRT-PCR

Among the identified DEGs, VPS35, FCGR2A, TBCA, HIRA, TYROBP and JUND were selected to verify the integrated result. The qRT-PCR results showed that TYROBP and JUND were up-regulated, while VPS35, FCGR2A, TBCA and HIRA were down-regulated. The expression of TBCA, HIRA and TYROBP were consistent with integrated analyses except VPS35, FCGR2A and JUND. The qRT-PCT results are shown in Figure 2.

Fig. 2

Validation of differentially expressed genes (DEGs) in the osteoporosis blood by quantitative reverse transcriptase-polymerase chain reaction. *p < 0.05, †p < 0.01.

ROC curve analysis

We performed ROC curve analyses and calculated the AUC to assess the discriminatory ability of DEGs (VPS35, HIRA, PHF20 and NFKB2) in data set GSE56815. The AUC of four DEGs including VPS35 (0.789), HIRA (0.77), PHF20 (0.851) and NFKB2 (0.741) was > 0.7 (Fig. 3). PHF20 had the largest AUC among these four DEGs. For the diagnosis of osteoporosis, the sensitivity (proportion of true positives) and 1-specificity (proportion of false positives) of VPS35 was 52.5% and 92.5%, respectively; the sensitivity and 1-specificity of HIRA was 60% and 87.5%, respectively; the sensitivity and 1-specificity of PHF20 was 77.5% and 85%, respectively; and the sensitivity and 1-specificity of NFKB2 was 67.5% and 75%, respectively.

Fig. 3

The receiver operating characteristic (ROC) curves of a) VPS35, b) HIRA, c) PHF20 and d) NFKB2 between osteoporosis patients and healthy controls. The ROC curves were used to show the diagnostic ability of these selected differentially expressed genes (DEGs) with 1-specificity (x-axis; the proportion of false positives) and sensitivity (y-axis; the proportion of true positives).

Discussion

Osteoporosis is a complex disease that is characterised by reduced bone mass and the deterioration of bone-tissue microarchitecture, ultimately leading to increased risk of fractures.[24] In the present study, we identified 855 up-regulated and 465 down-regulated genes in blood samples of osteoporosis patients compared with control groups. It is well known that osteoclast formation and activation are critical events in maintaining the normal bone mass and structure. Herein, GO and KEGG analyses indicated that these DEGs were significantly involved in osteoclast differentiation, which demonstrated the important role in osteoclast development of osteoporosis. Additionally, the TF-DEGs regulatory network was constructed involving top 20 (ten up-regulated and ten down-regulated) genes and 88 TFs, which further illustrated the role of these DEGs under the regulation of TFs in osteoporosis. Finally, qRT-PCR in vitro validated the expression patterns of several genes including VPS35, FCGR2A, TBCA, HIRA, TYROBP and JUND. Some results have been inconsistent with the microarray analyses, probably because of the heterogeneity of the studies, including different inclusion criteria and the small number of patients in the validation set. In summary, ten genes including VPS35, FCGR2A, TBCA, HIRA, TYROBP, JUND, PHF20, NFKB2, RPL35A and BICD2 were considered to play a potential role in the pathology of osteoporosis.

It has been reported that VPS35 is highly expressed in osteoclasts as well as osteoblasts and loss of function will increase hyper-resorptive osteoclast formation.[25] Specific VPS35 knockdown in the osteoblast lineage resulted in mildly lowered bone mass in the primary spongiosa.[26] Our study found up-regulated expression of VPS35, which indicated that it may function in the regulation of osteoclast and osteoblast activity in osteoporosis. In addition, VPS35 had significant diagnostic value for osteoporosis, which may serve as a diagnostic biomarker of osteoporosis. FCGR2A has been reported to be related to rheumatoid arthritis and is a target gene of drugs in rheumatoid arthritis treatment.[27,28] It has been indicated that FCGR2A is also associated with ankylosing spondylitis and axial spondyloarthritis.[29] In the current study, we found an increased expression of FCGR2A in osteoporosis, suggesting that FCGR2A may be associated with the pathology of osteoporosis.

In human cell lines, TBCA knockdown will decrease a mass of a- and b-tubulin levels, subtle changes in the microtubule cytoskeleton and cell death.[30] Herein, it is down-regulated in osteoporosis blood samples compared with controls. Therefore, we suggest that TBCA may participate in the formation process of cytoskeleton in osteoporosis. HIRA has been recorded associated with H3.3-containing nucleosomes in transcriptional active genomic loci in bone tumours.[31] In this study, we found that HIRA was down-regulated in osteoporosis blood samples, which suggested that HIRA may play a significant role in histone modification in bone development in osteoporosis. In addition, HIRA has diagnostic value for osteoporosis, suggesting that HIRA may be associated with the pathology of osteoporosis and may serve as a diagnostic biomarker of osteoporosis.

TYROBP is a protein involved in osteoclast differentiation and function, such as the generation of the actin cytoskeleton, which is important for bone resorption.[32] In the current study, the expression of TYROBP was up-regulated, indicating that it may play a vital role in the regulation of osteoclast differentiation in osteoporosis. JUND has proven to control bone formation, osteoblast proliferation, and differentiation.[33-36] Furthermore, it is a DEG in mesenchymal stem cells in osteoporosis.[37] Our results showed that JUND was expressed differentially in osteoporosis, which demonstrated the essential role of JUND in bone formation and osteoporosis.

It has been reported that PHF20-null mice show delay in bone formation, defects in skeletal composition and haematopoiesis.[38] We found that PHF20 was up-regulated, suggesting that PHF20 may play a role in bone development through different mechanisms or signal pathways in human osteoporosis. It is noted that PHF20 had a significant diagnostic value for osteoporosis, which may be considered as a biomarker in osteoporosis diagnosis. NFKB2 is reported up-regulated under estrogen treatment in needle bone biopsies of osteoporosis.[39] Additionally, it is a key gene in the TRAIL pathway for osteoporosis fractures.[40] It is noted that bone resorption marker TRAP5b was undetectable in Nfkb2+/– and Nfkb2–/– mice, and was slightly but significantly increased by TNF injection in Nfkb2+/– mice,[41] which suggested the relationship between NFKB2 and bone resorption. Moreover, our study revealed that NFKB2 was significantly enriched in osteoblast differentiation, which further indicated that it was an essential molecule in bone metabolism. Additionally, we found that NFKB2 was remarkably associated with diagnosis and may play a valuable role in the clinical and laboratory diagnosis of osteoporosis.

RPL35A participates in cytoplasmic ribosomal protein pathways in osteoarthritis chondrocytes.[42] In addition, it is linked to inherited bone marrow failure syndromes.[43] Herein, we found that RPL35A was down-regulated in osteoporosis blood samples compared with normal controls, which provided another pathogenic role in bone disease. BICD2 mutation is involved in spinal muscular atrophy.[44-47] Moreover, it has been reported that mutations in BICD2 will cause early onset non-length dependent lower-limb predominant weakness and contractures.[48] In the present study, we found decreased expression of BICD2 in osteoporosis, suggesting the vital role in bone formation and metabolism.

There are limitations to our study. First, the sample size in the qRT-PCR data set was small and larger numbers of blood samples of osteoporosis patients are needed for further research. Secondly, the deregulated DEGs in osteoporosis were identified and the definite biological function was not investigated in our study. In vivo and in vitro experiments are essential for elucidation of the biological roles of DEGs in osteoporosis in future work.

In summary, we identified several key genes including VPS35, FCGR2A, TBCA, HIRA, TYROBP, JUND, PHF20, NFKB2, RPL35A and BICD2 involved in the regulation of bone formation and metabolism under the regulation of TFs (FOXD3, Nkx2-5, Pax-4, Oct-1, HNF-4, Pax-6 and COMP1) in osteoporosis. It is noted that VPS35, HIRA, PHF20 and NFKB2 had significant diagnostic value for osteoporosis and may serve as diagnostic biomarkers of osteoporosis. Our results may provide important information for studying the pathogenic mechanisms and consequences of osteoporosis.