Identification of Aberrantly Expressed Long Non-Coding RNAs and Nearby Targeted Genes in Male Osteoporosis

Introduction

Primary osteoporosis is characterized by decreased bone mass, low bone mineral density (BMD) and loss of microarchitectural integrity with consequent increases the risk of fracture.¹ Primary osteoporosis mainly occurs in postmenopausal women and senile men. Like postmenopausal osteoporosis, male osteoporosis is becoming a serious health problem worldwide with aging.² The loss of BMD is approximately 1% per year in senile men and 20% of senile man aged over 50 years suffer an osteoporosis-related fracture.³ One third of men with osteoporosis suffer a hip fracture, a major risk factor for osteoporotic fracture-related death.⁴

In the last decades, there has been considerable effort to uncover the mechanisms underlying osteoporosis, especially the postmenopausal osteoporosis. Bone loss of osteoporosis was induced by remodeling imbalance between bone resorption and bone formation.⁴ Calcium,⁵ parathyroid hormone,⁶ estrogen⁷ and vitamin D⁸ are well-recognized regulators in osteoporosis. Multiple cytokines such as interleukin-1 (IL-1), IL-6, macrophage colony-stimulating factor (M-CSF), tumor necrosis factor-α (TNF-α) and receptor activator of nuclear factor-kB ligand (RANKL) were found to be involved in osteoporosis.^9,10 Immune system plays a crucial role in osteoporosis as well.^4,8 However, the molecular mechanisms of osteoporosis in men are not fully known.

Recently, accumulated evidence has indicated that long non-coding RNAs (lncRNAs) are involved in various biological progress by regulating the gene expression at transcriptional, epigenetic, or post-transcriptional levels.^11–13 It is reported that lncRNAs can alter the modification and location of transcription factors, act as precursor for siRNAs or miRNAs, regulate the alternative splicing of pre-mRNAs through SR complex, regulate histone modification by interacting with modification factors and bind to chromatin modification complexes to regulate chromatin remodeling and structure.¹⁴ Previous studies have reported that several lncRNAs such as DNACR and MEG3 are associated with postmenopausal osteoporosis.^12,13,15

In this study, we investigated different expression profiles of lncRNAs and mRNAs between male osteoporosis and normal control by high throughput RNA sequencing. Differentially expressed lncRNAs (DElncRNAs) and differentially expressed genes (DEGs) in male osteoporosis were further identified. Since lncRNAs can regulate the expression of their nearby gene by cis-regulatory effects,¹¹ the nearby targeted DEGs of DElncRNAs in male osteoporosis were identified as well. To further study the functions of DEGs in male osteoporosis, functional annotation was carried out.

Materials and Methods

Patients and Samples

We collected full blood samples from three men with osteoporosis and two healthy men. Detailed full blood collection was as follows: (1) the full blood was collected with BD tube (2.5 mL per tube) at room temperature (18–25°C); (2) the full blood was slowly reversed and mixed 8 to 10 times to mix well with the anticoagulant in the BD tube; (3) the full blood was placed vertically at room temperature for 12–24 h to fully contacted with the RNA stabilizer in the BD tube; (4) the full blood was transferred to the refrigerator at −20°C for 24 h, then transferred to the refrigerator at −80°C for long-term preservation.

The inclusion criteria of male osteoporosis were as follows: (1) according to world health organization (WHO) standards, patients whose T-score of the femoral neck or lumbar vertebrae or total hip is ≤-2.5 were diagnosed as osteoporosis; (2) patients had a new vertebral compression fracture with no trauma history; (3) patients had no prior history of drug use that might affect bone metabolisms such as glucocorticoids, thyroid hormone, parathyroid hormone, fluoride, calcitonin, thiazides, barbiturates, antiepileptic drugs, vitamin D and calcium. The exclusion criteria of male osteoporosis were as follows: (1) patients had diseases that affect bone metabolic such as kidney disease, liver disease, thyroid disease, diabetes, hyperprolactinemia, rheumatoid arthritis, ankylosing spondylitis, chronic abdominal bay absorption adverse symptoms, malignant tumor, blood disease and abnormal reproductive system and related medical history; (2) patients had traumatic fractures. In addition, the inclusion criteria of normal controls (men with degenerative osteoarthritis of the lumbar spine) were as follows: (1) individuals were matched by gender, ageand frequency of the case group; (2) T-scores of the femoral neck and lumbar vertebrae (L_1-4) and total hip are all ≥-1. The exclusion criteria of normal controls were as follows: (1) People with a previous history of orthopedic system disease. The detailed information of participants is shown in Table 1. All the participants submitted the signed informed consent. The protocol for this study has been approved by the Ethics Committee of Beijing Friendship Hospital, Capital Medical University (2017-P2-084-01).

Table 1 Patients Characteristics in RNA Sequencing

RNA Sequencing

We isolated total RNA from these five blood samples through the TRIZOL reagent (Invitrogen, Carlsbad, CA, USA). Nanodrop ND-2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA) was used to test the concentration and purity of total RNA. 2% agarose gel was used to check the integrity of RNA. We also obtained the RNA integrity number (RIN) value by an Agilent 2100 Bioanalyzer. The threshold for cDNA library construction was as follows: (1) the total amount of RNA was more than 5μg; (2) the concentration of RNA was more than 200ng/mL; (3) the value of OD260/280 was 1.8–2.2.

We removed ribosomal RNA by Ribo-Zero Magnetic kit (EpiCentre, Madison, WI, USA). RNA was fragmented into about 200 base pairs randomly. By using Truseq RNA sample Prep Kit (Illumina, Inc., San Diego, CA, USA), cDNA library was constructed. Firstly, the first cDNA strand was synthesized with RNA fragments primed with random hexamer primers. Secondly, we synthesized the second cDNA strand with dUTP instead of dTTP. Thirdly, end repair, 3ʹend adenylation and adapter ligation of double-stranded cDNAs were performed with End Repair Mix. Lastly, the second cDNA was digested by UNG enzyme (Illumina, Inc., San Diego, CA, USA). After the construction of the cDNA library, we amplified the library by performing polymerase chain reaction (PCR) for 15 cycles. Then, the amplified cDNA library was purified and quantified by Certified Low Range Ultra Agarose (Bio-Rad) and Picogreen (Molecular probes) on TBS380 (Turner Biosystems), respectively. The bridge PCR on cBot was performed by Truseq PE Cluster Kit v3-cBot-HS. In the end, sequencing of each library was performed by Illumina HiSeq 2000.

Data Processing

The FASTQ sequence data derived from RNA-seq results were obtained through Base Calling. After removed the low-quality reads by cutadapt (http://cufflinks.cbcb.umd.edu/), we obtained the high-quality clean data. Fasta and gff files of human reference genome (Ensemble GRCh38.p7) were downloaded from NCBI. TopHat (http://tophat.cbcb.umd.edu/) and Ensemble gene annotation were used to align clean reads to Ensemble GRCh38.p7. Cuffquant (http://cufflinks.cbcb.umd.edu/) was used to determine the expression of both mRNAs and lncRNAs and outputted the normalized data. The DEGs and DElncRNAs in male osteoporosis compared to normal control were obtained by using an R package, DEGseq. P-value <0.05 was the threshold of DEGs and DElncRNAs in male osteoporosis. The heat-maps of top 50 up- and down-regulated DEGs and DElncRNAs expression profile were performed with cluster in R language.

Protein-Protein Interaction (PPI) Network of DEGs

Based on Biological General Repository for Interaction Datasets (BioGRID) (http://thebiogrid.org/), we explored the proteins interacted with proteins encoded by the top 10 up-regulated and down-regulated DEGs in male osteoporosis. After removing these proteins that were not belonging to DEGs, we constructed the PPI network of top 10 up- and down-regulated DEGs in male osteoporosis.

Nearby Targeted DEGs of DElncRNAs

In order to identify the target DEGs of DElncRNAs, the DEGs transcribed within a 100 kbp window upstream or down-stream of DElncRNAs were researched and considered as nearby cis targets of DElncRNAs. The network of DElncRNA-near targeted DEGs was constructed by using Cytoscape.

Functional Annotation of All DEGs and Nearby Targeted DEGs of DElncRNAs

To further study the biological functions of identified all of DEGs and nearby targeted DEGs of DElncRNAs, we performed the Gene Ontology (GO) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of DEGs by online GeneCoDis3 software. False discovery rate (FDR) <0.05 was regarded as the cut-off for significant GO terms and KEGG pathways.

DElncRNAs-Nearby Targeted DEGs Interaction Network

Firstly, we constructed the DElncRNAs-nearby targeted DEGs interaction network by using all of the DElncRNAs and DEGs in male osteoporosis compared to normal control. Then, this DElncRNAs-nearby targeted DEGs interaction network and the male osteoporosis-specific PPI network were merged with the top 10 up- and down-regulated DEGs.

In vitro Validation

One DElncRNA-nearby targeted DEGs interaction pair of LINC02009-CCR2 was selected randomly for validation. A total of 36 full blood samples (16 cases) from male osteoporosis patients and healthy individuals (20 normal controls) were used for validation of CCR2 and LINC02009. The clinical information of these individuals was shown in the supplementary Table 1. We performed QRT-PCR in an ABI 7500 real-time PCR system. Relative gene expression was analyzed using 2^−ΔΔCt method. The human GAPDH gene was used as an endogenous control for mRNA expression in the analysis. All the participants submitted the signed informed consent. The protocol for this study has been approved by the Ethics Committee of Beijing Friendship Hospital, Capital Medical University (2017-P2-084-01).

Results

LncRNAs and mRNAs Expression Profile

A total of 3296 DEGs (1428 up-regulated DEGs and 1868 down-regulated DEGs) were obtained. A total of 204 DElncRNAs including 57 up-regulated DElncRNAs and 147 down-regulated DElncRNAs were obtained. The top 10 up- and down-regulated DEGs and top 10 up- and down-regulated DElncRNAs were shown in Tables 2 and 3, respectively. The heat-maps of the top 50 up- and down-regulated DElncRNAs and DEGs expression profile are illustrated in Figure 1.

Table 2 Top 10 Up- and Down-Regulated DEGs in Male Osteoporosis

Table 3 Top 10 Up- and Down-Regulated DElncRNAs in Male Osteoporosis

Figure 1 The heat-maps of top 50 up- and down-regulated DEGs and DElncRNAs in male osteoporosis. (A) the heat-map of top 50 up- and down-regulated DEGs in male osteoporosis. Row and column were used to represent DEGs and samples, respectively. The color scale was used to represent the expression levels of DEGs. (B) the heat-maps of top 50 up- and down-regulated DElncRNAs in male osteoporosis compared to normal control. Row and column were used to represent DElncRNAs and samples, respectively. The color scale was used to represent the expression levels of DElncRNAs.

PPI Network of DEGs in Male Osteoporosis

The PPI network consisted of 101 nodes and 105 edges (Figure 2). Based on this network, HLA-C (degree=18), KRT1 (degree=17) and IFIT1 (degree=15) were three hub proteins in male osteoporosis.

Figure 2 The male osteoporosis-specific PPI network. The PPI network was consisted of 101 nodes and 105 edges. Ellipses (nodes) were the proteins encoded by DEGs and ellipses with black border represented the proteins encoded by DEGs which belonged to the top 10 up- and down-regulated DEGs in male osteoporosis. Red and blue colors were used to represent up- and down-regulation in male osteoporosis, respectively.

Nearby Targeted DEGs of DElncRNAs

After searching the DEGs transcribed within a 100 kbp window upstream or down-stream of DElncRNAs, we obtained 168 DElncRNAs-nearby targeted DEGs pairs which consisted of 116 DElncRNAs and 162 DEGs. The DElncRNAs-nearby targeted DEGs network is illustrated in Figure 3. LOC101926887 (degree=5), LOC100131626 (degree=5) and TMX2-CTNND1 (degree=5) were the top three DElncRNAs covered most nearby targeted DEGs.

Figure 3 The DElncRNAs-nearby targeted DEGs interaction network in male osteoporosis. The DElncRNAs-nearby targeted DEGs interaction network was consisted of 278 nodes and 168 edges. Ellipses and rhombuses were used to represent the DEGs and DElncRNAs in male osteoporosis, respectively. Red and blue colors were used to represent up- and down-regulation in male osteoporosis, respectively.

Functional Annotation of All DEGs and Nearby Targeted DEGs of DElncRNAs

Firstly, we performed the functional enrichment analysis of all DEGs. In the GO enrichment analysis (Figure 4A–C), apoptotic process (FDR=4.78E-14), regulation of apoptotic process (FDR=8.32E-10), protein binding (FDR=3.42E-86), nucleotide binding (FDR=1.69E-22), nucleus (FDR=2.15E-73), cytoplasm (FDR=3.48E-72) were significantly enriched GO terms in male osteoporosis. Oxidative phosphorylation (FDR=3.52E-11), pyrimidine metabolism (FDR=8.20E-09), osteoclast differentiation (FDR=8.54E-09) and chemokine signaling pathway (FDR=5.61E-15) were significantly enriched pathways in male osteoporosis (Figure 4D). The detailed signaling pathway of osteoclast differentiation is shown in Figure 5. Then, we performed the functional annotation of nearby targeted DEGs of DElncRNAs. Based on the GO enrichment analysis (Figure 6), type I interferon-mediated signaling pathway (FDR=2.66E-05), cytokine-mediated signaling pathway (FDR=3.78E-05), protein binding (FDR=4.34E-09), zinc ion binding (FDR=9.54E-06), cytoplasm (FDR=1.26E-09) and cytosol (1.30E-06) were significantly enriched GO terms of nearby targeted DEGs of DElncRNAs in male osteoporosis. However, we did not obtain an enriched pathway for nearby targeted DEGs of DElncRNAs in male osteoporosis based on the KEGG enrichment analysis.

Figure 4 Functional annotation of all DEGs in male osteoporosis. (A) top 15 most significantly enriched biological process; (B) top 15 most significantly enriched cellular components; (C) top 15 most significantly enriched molecular functions; (D) top 15 most significantly enriched KEGG pathways. The x-axis shows -log FDR and y-axis shows the terms of GO/KEGG terms.

Figure 5 Pathway of osteoclast differentiation that enriched in DEGs in male osteoporosis. The red and green rectangles were represented the particles which regulated by the up- and down-regulated DEGs in male osteoporosis, respectively.

Figure 6 Functional annotations of nearby targeted DEGs of DElncRNAs in male osteoporosis. (A) top 15 most significantly enriched biological process; (B) top 15 most significantly enriched cellular components; (C) top 15 most significantly enriched molecular functions. The x-axis shows -log FDR and y-axis shows the terms of GO/KEGG terms.

DElncRNAs-Nearby Targeted DEGs Interaction Network

After merged the DElncRNAs-nearby targeted DEGs interaction network and the PPI network (involved top 10 up- and down-regulated DEGs), we obtained key DElncRNA-nearby targeted DEGs (top 10 up- and down-regulated DEGs) interaction network (Figure 7).

Figure 7 The DElncRNAs-nearby targeted DEGs interaction network in male osteoporosis. The DElncRNAs-nearby targeted DEGs interaction network was consisted of 26 nodes and 22 edges. Ellipses and rhombuses were used to represent the DEGs and DElncRNAs in male osteoporosis, respectively. Ellipses with black or blue border represented the DEGs which belonged to the top 10 up- and down-regulated DEGs in male osteoporosis compared to normal control. Red and blue colors were used to represent up- and down-regulation in male osteoporosis, respectively.

QRT-PCR Confirmation of DEGs

To manifest the expression of high throughput RNA sequencing for male osteoporosis, we compared the results of qRT-PCR and RNA sequencing of DElncRNA-nearby targeted DEGs interaction pair LINC02009-CCR2. In Figure 8, expressions of LINC02009 and CCR2 were all down-regulated and consistent with our RNA sequencing results.

Figure 8 QRT-PCR verification results. The X-axis was indicated different samples and the Y-axis was indicated the relative level of CCR2 mRNA (A) and LINC02009 lncRNA (B) for qRT-PCR results.

Discussion

Accumulated evidence has emphasized the importance of lncRNAs in various cancers and other diseases, such as osteoporosis.^11–13 Herein, we obtained the lncRNA and mRNA expression profiles between male osteoporosis and normal control by high throughput RNA sequencing. Based on functional annotation, we identified several DElncRNA-DEGs interaction pairs, which may be associated with the process of male osteoporosis.

Compared to normal control, HBB (hemoglobin subunit beta) was the most significantly up-regulated gene, while HBA1 was a significantly down-regulated gene in male osteoporosis. HBB and HBA1 encode the β and α 1 subunit of hemoglobin, respectively. As a major component of hemoglobin, iron is an important factor of bone metabolism. Iron can disorder the process of bone formation, make toxic effects on osteoblasts and reduce the recruitment of osteoblastic-lineage cells.¹⁶ It is pointed out that hemoglobin level is related to BMD and bone mass negatively and independently.¹⁷ Moreover, higher bone loss was found in men with lower hemoglobin level than women with lower hemoglobin level.¹⁷ Moreover, we found an interaction of “SNHG5-SYNCRIP-HBA1-HBB” based on the DElncRNA-nearby targeted DEGs interaction network. SYNCRIP (synaptotagmin binding cytoplasmic RNA interacting protein) was a down-regulated gene in male osteoporosis compared to normal control. SYNCRIP has been reported to be involved in bone morphogenesis protein signaling pathway.¹⁸ These three DEGs (HBA1, HBB and SYNCRIP) may be involved in the pathogenesis of male osteoporosis by their interactions. Small nucleolar RNA host gene 5 (SNHG5) is a member of the non-protein-coding multiple snoRNA host gene family which was significantly down-regulated in male osteoporosis in our study. SNHG5 has been reported to be associated with various cancers.^19,20 SNHG5 is involved in gastric cancer cell proliferation and migration.¹⁹ SNHG5 can be served as a tumor marker of malignant melanoma²⁰ and promotes tumor cell survival in colorectal cancer.¹⁹ However, no association between SNHG5 and male osteoporosis has been reported. Our study suggested that down-regulated SNHG5 may play an important role in male osteoporosis through regulating the expression of SYNCRIP, HBB and HBA1.

Based on the male osteoporosis-specific PPI network, major histocompatibility complex, class I, C (HLA-C) was a hub gene, which was up-regulated in male osteoporosis compared to normal control. HLA-C is a member of the HLA class I heavy chain paralogues. HLA class I molecule consists of a heavy chain and a light chain which plays a key role in the immune system. Accumulated evidence has demonstrated the relationship between the immune system and bone.²¹ T cells play a crucial role in the formation and activity of osteoclast and osteoblast in multiple skeletal diseases, such as osteoporosis.²² The production of T-lymphocyte, T-helper and T-cytotoxic cells is decreased with aging.²³ Additionally, other osteoporosis-related factors including estrogen, calcium, parathyroid hormone, and vitamin D play a regulatory in immune system.²³ The role of immune system in the pathogenesis of osteoporosis is mainly induced by decreased estrogen and secondary hyperparathyroidism.²³ In this study, HLA-C was a targeted nearby cis-DEG of HCG27 (HLA complex group 27). In this study, HCG27 was a significantly up-regulated lncRNA in male osteoporosis. Our finding emphasized the importance of HCG27 in the process of male osteoporosis.

Interferon-induced protein with tetratricopeptide repeats 1 (IFIT1) was another hub gene based on PPI network which was up-regulated in male osteoporosis. Depending on the DElncRNAs-nearby targeted DEGs interaction network, IFIT1 could be interacted with three other up-regulated DEGs including IFIT2, IFIT3 and IFIT5. These four genes were interferon-induced proteins that contain tetratricopeptide repeats.²⁴ As an innate immune bottleneck, IFIT1 is expressed during the response of osteoblasts to immune cytokine interferon-β.²⁵ IFIT2 is a regulator of innate immunity as well. Both IFIT1 and IFIT2 play a role in the fracture healing process in osteoporosis,²⁶ which provides new clues for osteoporotic fracture repair therapies. LOC101926887 can interact with the IFIT1-IFIT2, IFIT1-IFIT3 and IFIT1-IFIT5 interaction pairs. In addition, LOC101926887 was one of the up-regulated lncRNA covered most nearby targeted DEGs. Through regulating the expression of targeted nearby DEGs (IFIT1, IFIT2, IFIT3, IFIT15 and LIPA), LOC101926887 may be an important regulator of male osteoporosis. LOC100131626 and TMX2-CTNND1 were the other two DElncRNAs that covered most nearby targeted DEGs in male osteoporosis. LOC100131626 was a down-regulated lncRNA in male osteoporosis containing five targeted nearby DEGs (KMT2A, ATP5L, UBE4A, CD3D and CD3E). TMX2-CTNND1 was an up-regulated lncRNA in male osteoporosis containing five targeted nearby DEGs (CTNND1, TMX2, ZDHHC5, CLP1 and YPEL4). These two DElncRNAs may play crucial roles in male osteoporosis through regulating their nearby targeted DEGs as well.

Compared to normal control, Loc105372801 and KCNQ1OT1 was the most significantly up-regulated and down-regulated DElncRNA in male osteoporosis, respectively. No previous studies have reported functions of Loc105372801 and KCNQ1OT1 in male osteoporosis. KCNQ1 opposite strand/antisense transcript 1 (KCNQ1OT1) is only expressed on the paternal allele.²⁷ KCNQ1OT1 is closely related to the pathogenesis of familial Beckwith-Wiedemann syndrome.²⁸ Loss of imprinting of KCNQ1OT1 locus by epigenetic disruption is related to colorectal carcinogenesis.^29–31 A novel tetranucleotide repeat polymorphism (rs35622507) within KCNQ1OT1 confers risk for hepatocellular carcinoma.³⁰ KCNQ1OT1 involves in breast cancer through the epigenetic suppression of the cyclin-dependent kinase inhibitor 1C.³¹ The precise function of KCNQ1OT1 in male osteoporosis needs further research.

In our study, chemokine signaling pathway was a significantly enriched pathway in male osteoporosis. As one of the DEGs enriched in chemokine signaling pathway, C-C motif chemokine receptor 2 (CCR2) is a receptor for ligands like monocyte chemoattractant protein-1 (MCP-1). CCR2 is served as a regulator of bone loss, which is expressed on preosteoclasts and osteoclasts.³² The MCP-1/CCR2 axis has been demonstrated to be involved in fusion and maturation of osteoclast and bone resorption.^33,34 In addition, gene variants of MCP-1 and CCR2 are reported to be risk factors for osteoporosis.³⁵ In young adult mice, knockout of the Ccr2 results in delay in fracture healing.^36,37 In the present study, CCR2 and LINC02009 were all significantly down-regulated in male osteoporosis, which was validated by qRT-PCR. CCR2 was a nearby targeted DEG of a down-regulated DElncRNA, LINC02009. Our result suggested that LINC02009 may be associated with the process of osteoporosis by regulating the expression of CCR2 and other nearby targeted DEGs including lactotransferrin (LTF) and leucine-rich repeat containing 2 (LRRC2).

Osteoclast differentiation was another significantly enriched pathway for DEGs in male osteoporosis which was demonstrated to be crucial pathway of osteoporosis. These DEGs enriched in osteoclast differentiation may be important regulators in male osteoporosis as well.

Taken together, identified several DElncRNA-nearby targeted DEGs interaction pairs such as SNHG5-SYNCRIP-HBA1-HBB, HCG27-HLA-C, LINC02009-CCR2, LOC101926887-IFIT1-IFIT2, LOC101926887-IFIT1-IFIT3 and LOC101926887-IFIT1-IFIT5 may play an important role in the pathogenesis of male osteoporosis. The validated DElncRNA-nearby targeted DEGs interaction pair of LINC02009-CCR2 may be important regulators in the process of male osteoporosis as well. Our findings may contribute to underlying the function of lncRNAs and mechanism of male osteoporosis and provide new clues for novel diagnostic and therapeutic strategies for male osteoporosis. However, there are limitations to our study. Firstly, the sample size in the RNA sequencing is small, larger numbers of blood samples or bone samples are further needed for further validation of identified lncRNAs and genes. Secondly, we did not validate the expression of identified other DElncRNAs and DEGs in male osteoporosis or any cell lines. Thirdly, we did not perform the deeper mechanism study of male osteoporosis, some in vivo or in vitro experiments are further needed.