Key circular RNAs identified in male osteoporosis patients by whole transcriptome sequencing

Background Osteoporosis (OP) is a systemic disease with bone loss and microstructural deterioration. Numerous noncoding RNAs (ncRNAs) have been proved to participate in various diseases, especially circular RNAs (circRNAs). However, the expression profile and mechanisms underlying circRNAs in male osteoporosis have not yet been explored. Methods The whole transcriptome expression profile and differences in mRNAs, circRNAs, and microRNAs (miRNAs) were investigated in peripheral blood samples of patients with osteoporosis and healthy controls consisting of males ≥ 60-years-old. Results A total of 398 circRNAs, 51 miRNAs, and 642 mRNAs were significantly and differentially expressed in osteoporosis compared to healthy controls. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis showed that the host genes of significantly differentially expressed circRNAs were mainly enriched in the regulation of cell cycle process: biological process (BP), organelle part cellular components (CC), protein binding molecular function (MF), Toll-like receptor signaling pathway, tumor necrosis factor (TNF) signaling pathway, and thyroid hormone signaling pathway. circRNA-miRNA-mRNA regulatory network was constructed using the differentially expressed RNAs. Moreover, key circRNAs (hsa_circ_0042409) in osteoporosis were discovered and validated by qPCR. Conclusions The key cicrRNAs plays a major role in the pathogenesis of osteoporosis and could be used as potential biomarkers or targets in the diagnosis and treatment of osteoporosis.


INTRODUCTION
Osteoporosis (OP) is a systemic disease with osteopenia and microstructural deterioration that increases the risk of fracture susceptibility, especially in the spine, buttocks, and wrists (Kuo & Chen, 2017). Male osteoporosis is a common age-related degenerative disease, characterized by impaired bone formation and low bone turnover. Postmenopausal osteoporosis is related to excessive bone resorption caused by estrogen deficiency. According to the report of the World Health Organization (WHO), the number of OP patients is increasing rapidly (Cabral et al., 2016), and accounts for about 6.6% of the total population in China (Zhang et al., 2018a). Men have a higher mortality rate after fracture than women (Bliuc & Center, 2016). Previous studies have identified some biochemical indexes of OP, which can be used for the diagnosis and monitoring of patients with OP (Wang et al., 2017), such as serum collagen type I N-terminal pre-peptide (PINP) and crosslinked C-terminal peptide (CTX) (Yang et al., 2019). However, conventional biochemical markers are not effective in determining the possible secondary causes of osteoporosis in men (Fink et al., 2016), while currently available biochemical markers cannot detect all risk factors for fractures (Cheng et al., 2017). Therefore, it is identifying new biomarkers to improve the diagnosis and treatment of osteoporosis in men is imperative.
Accumulating evidence shows that ncRNAs are associated with various diseases through indirect or direct regulation of the corresponding gene expression (Yang et al., 2018). ncRNAs, including miRNA and circRNA, play a crucial role in the occurrence, development, and progression of cancer (Dou et al., 2016). Importantly, previous studies have suggested that circulating miRNAs may be used as a critical biomarker for osteoporosis (Fu et al., 2018;Materozzi et al., 2018).
circRNAs are a class of endogenous, abundant, non-polyadenylated RNAs with a covalently closed, continuous loop structure (Zhang, Yang & Xiao, 2018b). These RNAs are associated with various biological processes, and their dysregulated expression are implicated in human diseases, including diabetes, Alzheimer's disease, tumors, and cardiovascular disease because of their high stability and prevention from RNA exonuclease degradation (Qiao et al., 2018;Zhao et al., 2017). Interestingly, the complex regulatory interactions between different types of ncRNAs have fundamental roles in the development of multiple diseases (Peng et al., 2018;Zhong et al., 2018). circRNAs are wellknown as miRNA sponge in inhibiting the function of miRNA via competing endogenous RNA (ceRNA) network (Zhong et al., 2018). For instance, circRNA-ZNF609, containing multiple binding sites for miR-150-5p, regulates ATK3 expression in Hirschsprung's disease through ceRNA network (Fu et al., 2018). It has also been proved that circRNA MYLK binds competitively to miRNA29a-3p, thereby increasing the expression of the target genes VEGFA, DNMT3B, and ITGB1, involved in the progression of bladder cancer (Fu et al., 2017). Together, these findings suggested that mRNAs, miRNAs, and circRNAs play a major role in various human diseases, such as osteoporosis (Mandourah et al., 2018).
In the present study, whole transcriptome sequencing was carried out on monocytes from male healthy controls and osteoporosis patients. Key circRNAs involved in the pathogenesis of osteoporosis were identified by bioinformatics analysis. Thus, our findings provide a basis for further in-depth study of pathogenic genes and the rapid, simple diagnosis, and treatment of osteoporosis in men.

Patients and samples
Study participants were collected from October 2016 to November 2017, and bone mineral density (BMD) was examined in the Second Affiliated Hospital of Harbin Medical University. All peripheral blood samples, including healthy controls and patients with osteoporosis, were collected. Healthy controls were defined by spine BMD T-score ≥−1.0 SD, while osteoporosis was defined by spine BMD T-score ≤−2.5 SD. All participants were males, aged ≥ 60-years-old. The detailed characteristics of the study samples are shown in Table 1. This study was approved by the Ethics Committee of The Second Affiliated Hospital of Harbin Medical University (#KY2016-198). All patient samples were obtained at the time of diagnosis, and informed consent was signed at the The Second Affiliated Hospital of Harbin Medical University.

RNA isolation and RNA sequencing
Total RNA was isolated from mononuclear cells of 6 peripheral blood samples (3 OP and 3 healthy controls) using TRIzol reagent (Sigma, St. Louis, USA), following manufacture's protocol. An equivalent of 5 µg RNA was utilized as input material for the RNA sample preparations. Libraries were constructed utilizing rRNA depleted and RNase R digested RNAs or NEBNext Multiplex Small RNA Library Prep Set for Illumina (NEB, USA), according to manufacturer's instructions. After cluster generation on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot or TruSeq SR Cluster Kit v3-cBot-HS (Illumina), the library preparations were sequenced on an Illumina HiSeq 2500/4000 platform. The flowchart was as follows (Fig. 1).

Bioinformatics analysis
circRNAs were predicted using find_circ (Memczak et al., 2013) and CIRI2 (Zeng et al., 2017) to reduce false positives. The predicted circRNA results of the two software were intersected based on the position of circRNAs on chromosome. Stringent filter criteria were applied to select candidate circRNAs as follows: at least junction reads ≥5 in one samples or junction reads ≥2 in all samples of one group. The gene expression level was quantified using TPM (readCount ×1,000,000)/libsize. Deseq2 (Love, Huber & Anders, 2014) was employed to perform differentially expressed gene analysis with the cutoff fold-change > 1 and adjust p-value < 0.05. GOseq and KOBAS (Xie et al., 2011) were used to carry out Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, respectively. miRTarBase was used to predict the miRNAs that target the differentially expressed mRNAs, while miRanda was utilized to predict the binding sites of miRNA and circRNA. Cytoscape was employed to construct the miRNA-circRNA-mRNA regulatory network. The flowchart was as follows (Fig. 2).

Statistical analysis
qRT-PCR experiment was repeated three times. SPSS was utilized for statistical analysis with independent t -test. P < 0.05 was considered as statistically significant.

Differentially expressed circRNAs, miRNAs and mRNAs
Raw data(raw reads) of fastq format were processed using in-house Perlscripts. Clean data (clean reads) were obtained by removing the reads containing adapter, reads containing ploy-N, and low quality reads from raw data. Differential expression analysis of the two groups was performed using the DESeq R package (1.10.1), which further determined the differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting P-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted P-value identified by DESeq were termed as differentially expressed. A total of 12,839 circRNAs, including 5,682 circRNAs, were novel. 398 circRNAs were differentially expressed between OP and healthy control. Of those, 195 circRNAs were upregulated while 203 circRNAs were downregulated ( Fig. 3 and Table 2) with the cutoff fold-change > 1 and adjusted P-value (padj) < 0.05. We also identified 642 differentially expressed mRNAs (305 upregulated and 337 downregulated) and 51 miRNAs (28 upregulated and 23 downregulated), respectively.

Functional enrichment analysis of differentially expressed circRNAs
GO and KEGG pathway enrichment analysis were carried out using the host genes of significantly differentially expressed circRNAs. These circRNAs are mainly enriched in biological processes (BP), including metabolic process (GO:0008152), cellular metabolic process (GO:0044237), biological regulation (GO:0065007) and regulation of cellular process (GO:0050794); cellular component (CC) like intracellular organelle (GO:0043229)  Table 3). The differentially expressed circRNAs were also enriched in viral carcinogenesis, Toll-like receptor signaling pathway, tumor necrosis factor (TNF) signaling pathway, and thyroid hormone signaling pathway ( Fig. 4B and Table 4).

qPCR experiment validation
The circrNA-mirNA-mrna regulatory network was constructed, and several key circRNAs related to miRNA and mRNA. Among these, hsa_circ_0042409 was linked to 7 miRNAs and 26 mRNAs, and KLC1 expression was regulated by inhibiting miRNA hsa-mir-195-5p. qPCR showed that the expression level of circRNA hsa_circ_0042409 and KLC1 mRNA was significantly increased in male osteoporosis patients, while that of hsa-miR-195-5p was significantly decreased with the cutoff of P-value < 0.05 (Fig. 7) compared to healthy controls.

Figure 3 (A) Volcano plot (A) of differentially expressed circRNAs between OP and healthy control.
Red indicates upregulated while green represents downregulated. (B) Heatmap of differentially expressed circRNAs between OP and healthy control with red denoting hig expression and blue signifying low expression.
Full-size DOI: 10.7717/peerj.11420/ fig-3 DISCUSSION Osteoporosis (OP) is a systemic disease with reduction in bone mass and deterioration of microstructure augmenting the risk of fragility and susceptibility to fracture, especially in the spine, hip, and wrists. Increasing evidence has revealed that ncRNAs participated in various diseases by directly or indirectly regulating the corresponding gene expression. Furthermore, ncRNAs are associated with various diseases through indirect or direct regulation of the corresponding gene expression (Yang et al., 2018). ncRNAs, including miRNA, long lncRNA, and circRNA, play a crucial role in the occurrence, development, and progression of cancer (Dou et al., 2016). Importantly, previous studies have suggested that circmiRNAs may be used as vital biomarkers for osteoporosis (Fu et al., 2018;Materozzi et al., 2018). Wei et al. (2012) found that miR-34s inhibit osteoblast proliferation and differentiation in mice by targeting SATB2. Xia et al. (2016) also discovered that miR-31-5p and miR-424-5p were downregulated in cartilage-derived mesenchymal stem cells (CMSCs) from the degraded cartilage. Moreover, Yin et al. (2018) reported previously that circRUNX2 regulated RUNX2 to prevent osteoporosis via hsa-miR-203. Herein, we performed whole transcriptome sequencing and investigated the regulatory mechanisms and functions of non-coding RNAs, especially circRNAs in the pathogenesis of osteoporosis. A total of 398 circRNAs, 51 miRNAs, and 642 mRNAs were identified to be significantly differentially expressed in osteoporosis compared to healthy controls. Moreover, GO and KEGG enrichment analysis illustrated that the host genes of significantly  differentially expressed circRNAs were mainly enriched in the regulation of cell cycle processes, such as BP, organelle part CC, protein binding MF, Toll-like receptor signaling pathway, TNF signaling pathway, and thyroid hormone signaling pathway. Based on the circRNA-miRNA-mRNA regulatory network, some key circRNAs in osteoporosis were  discovered further, such as hsa_circ_0042409, hsa_circ_0001924, hsa_circ_0003990, and hsa_circ_0000983. hsa_circ_0042409 was linked to 8 miRNAs and 26 mRNAs in the ceRNA regulatory network, and was upregulated in osteoporosis. In addition, it regulated the expression of KLC1, RNH1, CPEN1, and STXBP2 by inhibiting hsa-miR-195-5p, hsa-miR-30b-5p, hsa-miR-32b-5p, hsa-miR-378d, hsa-miR-424b-5p, and hsa-miR-6763-5p. Qu, Chu & Wang (2017) demonstrated that miR-195-5p suppresses osteosarcoma cell proliferation and invasion by suppressing naked cuticle homolog 1. Zhang et al. (2019) discovered that DOC2B promoted insulin sensitivity in mice via a novel KLC1-dependent mechanism in skeletal muscle. Importantly, we found that circRNA hsa_circ_0042409 and KLC1 mRNA were significantly increased while hsa-miR-195-5p was significantly decreased in male osteoporosis patients by whole transcriptome sequencing and further validated by qPCR. These findings coincidentally suggested that key circRNA, hsa_circ_0042409, is associated with the development of osteoporosis. The increased expression of hsa_circ_0042409 regulated the expression level of KLC1 by spongy hsa-mir-195-5P, thus promoting the pathogenesis of osteoporosis.

CONCLUSIONS
Although further studies might be needed to support these findings, key cicrRNAs (hsa_circ_0042409 et al.) play a major role in the pathogenesis of osteoporosis and could be used as potential biomarkers or targets in the diagnosis and treatment of osteoporosis.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
The authors received no funding for this work.