Identification and Characterization of MicroRNAs by High Through-Put Sequencing in Mesenchymal Stem Cells and Bone Tissue from Mice of Age-Related Osteoporosis

The functional deficiencies of bone marrow-derived mesenchymal stem cells (MSCs) may contribute to the aging process and age-related diseases, such as osteoporosis. Although it has been reported that microRNAs (miRNAs) played an important role in mechanisms of gene regulation of aging, and their expression profiles in MSCs osteogenic differentiation were established in recent years, but it is still elusive for the dynamic patterns of miRNAs in aging process. Importantly, the miRNAs in aged bone tissue had not been yet reported so far. Here, we combined high through-put sequencing with computational techniques to detect miRNAs dynamics in MSCs and bone tissue of age-related osteoporosis. Among the detected miRNAs, 59 identified miRNAs in MSCs and 159 in bone showed significantly differential expressions. And more importantly, there existed 8 up-regulated and 30 down-regulated miRNAs in both MSCs and bone during the aging process, with the majority having a trend of down-regulation. Furthermore, after target prediction and KEGG pathway analysis, we found that their targeted genes were significantly enriched in pathways in cancer, which are complex genetic networks, comprise of a number of age-related pathways. These results strongly suggest that these analyzed miRNAs may be negatively involved in age-related osteoporosis, given that most of them showed a decreased expression, which could lay a good foundation for further functional analysis of these miRNAs in age-related osteoporosis.


Introduction
Aging increases the risk to develop several diseases, such as osteoporosis, which is regarded as highly age-related [1], featuring the bone loss and susceptibility to fragility fractures [2,3]. Theoretically, as an attractive therapeutic candidates for several diseases and degenerative applications, mesenchymal stem cells (MSCs) have also been evaluated as a vital factor in the development of osteoporosis, since it keeps the ability of selfrenewal and differentiation into multiple cell types, including osteocytes, chondrocytes, and adipocytes [4,5,6,7,8]. Indeed, recently, it has been reported that the functional deficiencies of MSCs could lead to the declining of bone integrity and function in the elderly [9,10]. Moreover, other studies showed that MSCs were involved in continuous maintenance and repair of bone during aging [9,11,12]. Environmentally, the neighboring aged tissues should also be taken into account, for they could exert potential impact on MSCs, further causing the decrease in regenerative potential of bone [13,14,15]. Furthermore, epigenetically, microRNAs (miRNA) could influence many basic cellular and pathological processes by regulating gene expression post-transcriptionally through binding to complementary sequences in the 39 untranslated region (39 UTR) of target mRNAs [16]. And recent studies showed that miRNAs participated in regulation of aging and a variety of age-associated pathways [17]. For example, Inukai et al. [18] screened miRNA expression by using brains from mice of different ages and found that there was a global downward trend of miRNA expression during aging. Additionally, Mori et al. [19] observed an attenuated miRNA processing in adipose tissues during aging in mice, worm, and human, suggesting that it could be a conserved feature of aging. Hence, as for osteoporosis, being an age-related disease, it may be closely linked to miRNA expressions. However, to date, only limited types of tissues or organs have been analyzed in the mouse and human [15]. Yet little is known about the miRNA expression in the bone of mouse at different ages.
Here, in this study, to acquire a better understanding of potential contributions of miRNAs to age-related osteoporosis, we used high through-put deep sequencing to perform a comprehensive survey of miRNA expressions in MSCs and bone tissue based on a mouse model of osteoporosis. This study could provide reliable evidence for further uncovering the development mechanisms of osteoporosis in human during aging.

Construction of a Small RNA Library [20] by Deep Sequencing
MSCs. In order to detect the miRNAs dynamics in agerelated osteoporosis, a small RNA library of MSCs samples in three groups were obtained by deep sequencing. Overall, 14,565,370 clean reads were detected in young group (2 m), 13,834,317 clean reads were detected in adult group (8 m), and 13,001,553 clean reads were detected in old group (25 m). After alignment to the mouse genome (mm9), the results indicated that 11,843,150 (81.31%) reads in young MSCs, 11,194,206 (80.92%) reads in adult MSCs and 10,666,211(82.04%) reads in old MSCs were matched to mm9. The small RNAs were classified into different categories according to their biogenesis and annotation. Among the total reads, the ratio of noncoding RNA, including rRNA, tRNA, snRNA and snoRNA accounted for 1.99%, 1.3% and 1.03% in young, adult and old group, respectively (Table S1). The remaining small RNAs were retained for further analysis.
Bone tissue. Total RNA of the bone tissue samples of femurs and tibiae from three groups were extracted, subsequently, we used high through-put sequencing to identify miRNA profiles in bone. 14,979,697 clean reads, 14,597,780 and 10,899,301 clean reads were detected in young, adult and old group, respectively. But unlike the MSCs, there had a much lower genome-matched rate in three bone groups. The ratio was 54.75% (8,201,928), 68.63% (10,017,983) and 64.00% (6,975,893) in young, adult and old bone, respectively. At the same time, the rates of rRNA, tRNA, snRNA and snoRNA were increased in bone; they were 32.8%, 28.8% and 18.9% in young, adult and old group, respectively (Table S1).

MiRNAs Profiles in Age-related Osteoporosis
To investigate miRNA expressions in age-related osteoporosis, genome-aligned reads correspond to known miRNAs were compared with miRBase database (miRBase18) to obtain the miRNA count. The total 1157 miRNAs were detected; however, only 587, 578, and 574 miRNAs were expressed in 2 m-MSCs, 8 m-MSCs and 25 m-MSCs, respectively (Table S2). For bone tissue samples, there were 676, 657 and 603 miRNAs expressions in 2 m-bone, 8 m-bone and 25 m-bone, respectively (Table S3).

Differential Expressions of Known miRNA
To analyze the variance of miRNAs that is differentially regulated during bone aging process of mice, we compared the known miRNA expressions between two samples to find out the differentially expressed miRNA. The expression of miRNA in two samples was shown by plotting Log2-ratio figure and scatter plot. For MSCs, the expression levels of miRNAs were shown in Table  Table S5 and Figure 1D-F. To generate a convinced list of miRNAs variation in age-related osteoporosis, we filtered these data step by step. Firstly, selecting two-fold altered expression between the 2 m and 25 m samples, and ensuring Pvalue ,0.05. The expression data were in Table S6. In MSCs, there were 59 miRNAs differently expressed, the miRNA expression changes were shown in Figure 2A. Comparing the numbers between down-regulated and up-regulated miRNAs (30 versus 29), we found that it reached a new homeostasis for the whole miRNA pattern in aged MSCs. In other words, there existed the transition of steady state of miRNA dynamics in different aged MSCs. However, from Figure 2B, we could find the changes of miRNA expressions in 2 m and 25 m bone, in 159 changed miRNAs, majority of them were down-regulated in old bone (114 versus 45). Interestingly, some miRNAs changes trends were opposite in MSCs and bone (red cells in Table S6). So, secondly, considering the miRNA expressions showed a consistent changing trend either decreasingly or increasingly among the three samples of MSCs and bone. Finally, screening the common miRNAs existed in the MSCs and bone, and ensuring the P-value ,0.05 (Table S7). The results showed that there were 8 upregulated and 30 down-regulated known miRNAs during bone aging (Table 1). Biological functions of corresponding miRNAs are involved in aging and osteoblastic differentiation was shown in Table 2.
In order to validate the sequencing data, we selected several miRNAs from Table 2 for additional qRT-PCR validation, which the minimum normalized read count of miRNAs was 5 in young, adult and old groups, including miR-210 [29], miR-22 [31,32], miR-31 [36,37], and miR-10b [16] (Figure 3). Except for miR-10b, other miRNAs showed similar expression trends in qPCR and deep sequencing data. This result suggests that high through-put sequencing is an effective method for identifying mature miRNAs, and shows our prospect on the further research on the function of single miRNA.

Computational Target Prediction of Changed miRNAs and Pathway Analysis
To identify the potential targets of changed miRNAs, we computationally identified mRNA targets of miRNAs using the software Mireap (software developed by BGI). To annotate the target genes corresponding to certain biological function during bone aging, we used GO (Gene Ontology) enrichment analysis to reveal the functions significantly related with predicted target gene candidates of changed miRNAs. The significantly enriched GO terms of biological process in genes repressed by changed miRNAs were shown in Table 3 and Table 4. To further understand the biological functions of changed miRNAs, KEGG pathway analysis was carried out among the predicted target genes of miRNAs. The KEGG analysis could reveal the main pathways which the target gene candidates are involved in. For up-regulated miRNAs, we found that cell cycle pathway (ko04110, P-value = 0.01621647), pathways in cancer (ko05200, P-value = 0.01038872) and calcium signaling pathway (ko04020, P-value = 0.006787428) were significantly targeted KEGG pathways ( Figure S1A-C). These pathways were correlated with cell proliferation, osteogenic differentiation, and aging. And for down-regulated miRNAs, pathways in cancer (ko05200, Pvalue = 0.01555212) was significantly targeted KEGG pathway ( Figure S1D).

KEGG Pathways Analysis of Target Genes of miR-31
Due to the fact that each miRNA has hundreds of target genes, which function through multiple pathways, in either a positive or negative manner, it is rather difficult to perform a general evaluation of the effects of all the down-regulated miRNAs on aging and osteogenic differentiation in our study. Hence, here only miR-31 was chosen as representative, for it has been shown to be associated with aging and osteogenic differentiation in some other studies. After KEGG analysis, its target genes were indeed enriched in the p53 and Wnt signaling pathways, which are important parts of ko05200 ( Figure 4, Table S8). Notably, based on the literature, p53 signaling pathway is associated with aging and Wnt pathway is associated with cell differentiation.

Animals
Male C57BL/6J mice (Experimental Animal Center of Fourth Military Medical University, China), 2 months-old (young group), 8 months-old (adult group) and 25 months-old (old group, osteoporosis) were used, each group had three samples.

Ethics Statement
All animal experiments were performed under an animal study protocol approved by the ethics committee of Fourth Military Medical University.

Cell Culture and Bone Collect
Femurs and tibiae were isolated from nine mice. The bone marrow was flushed out with a-MEM (Gibco, New York, NY, USA) containing 20% FBS (GIBCO, New York, NY, USA) and seeded in 90-mm dishes, cultures were refeeded every 3,4 days, maintained about 15 days, then cells were passaged. The femurs and tibiae without bone marrow were stored in liquid nitrogen.

RNA Preparation
Total cellular RNA was extracted from the second passage MSCs using Trizol Reagent (Invitrogen, Carlsbad, CA, USA), according to the manufacturer's instructions. And the total RNA of femurs and tibiae was extracted as described by the office of shared research facilities in the University of Chicago (http://fgf. uchicago.edu/protocols_microarray.php), with slight modifications. In brief, guanidine buffer (containing 7.5 uLb-ME per mL) was replaced by Trizol reagent, and phenol/chloroform/ isoamyl alcohol was replaced by chloroform/isoamyl alcohol.

High Through-put Deep Sequencing
High through-put deep sequencing was carried out using Hiseq2000 Sequencer (Illumina); the sequencing procedure was described by Chen et al. [38,39]. Then the small RNA digitalization analysis based on the the image files generated by the sequencer were processed. The 50 nt sequence tags from HiSeq sequencing go through the data cleaning first, which included getting rid of the low quality tags and several kinds of contaminants from the 50 nt tags. Clean reads that correspond to known miRNAs were compared with a miRBase database (release 18.0) and the total copy number of each sample was normalized to 1,000,000. The data was processed using software developed by BGI. The high through-put sequencing data for this study have been submitted to the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/Traces/sra) under accession number SRA072825 and SRA074415.

Known MiRNA Expression Analysis
Considering the diversity of miRNA expression changes with age, we studied samples (MSCs and bone) from 2, 8 and 25 months-old mice to generate a convinced list of miRNAs variation in age-related osteoporosis. Primarily, we compared the known miRNA expressions between two samples, including 2 m versus 8 m, 8 m versus 25 m and 2 m versus 25 m. Then, the data were screened step by step. MiRNA expressions between the 2 m and 25 m samples were analyzed in detail, the miRNA with two-fold changed expressions (P-value,0.05) were selected. After that, we analyzed the data of 8 m samples, selected the miRNAs with consistent changing trend during aging process for further study. At last, we confirmed the convinced list of miRNAs variation by screening the common miRNAs existed in the MSCs and bone.

MiRNA Target Prediction and Functional Analysis
Targets for known miRNAs were predicted using software Mireap (software developed by BGI). GO enrichment analysis of target gene candidates was carried out according to the GO terms in the database (http://www.geneontology.org/). KEGG pathway analysis is also used for the target gene candidates. Genes with FDR,0.05 were considered as significantly enriched in target gene candidates.

Quantitative RT-PCR Assay
For real time PCR validation, miRNA cDNA was obtained according to specification of the One step PrimeScript miRNA cDNA synthesis kit (TaKaRa, Japan), Equal amounts of cDNA were used in duplicate and amplified with SYBR Premix Ex Taq

Statistics
Differences were statistically evaluated by Student's t-test. Unless otherwise stated, P-value,0.05 was considered to be statistically significant.

Discussion
Global expression profiling of miRNAs in various aged tissues or organs of different species have been performed. However, few data on miRNAs pattern are available so far in terms of agerelated osteoporosis. In present results, we provided miRNAs dynamics in MSCs and bone tissue of age-related osteoporosis and further performed target prediction and KEGG pathway analysis, which strongly suggested that multiple miRNAs may be involved in this complex process.
As multipotent adult stem cells, MSCs are capable of differentiating into osteoblast cells, however, many studies proved that this ability declined with age [10,12,40,41]. On the molecular level, an increasing number of miRNAs have been identified to regulate osteogenic differentiation of MSCs [8,42,43,44]. Therefore, miRNAs are crucial in the bone aging  Table 3. Functional annotation clusters of enriched GO biological processes predicted to be suppressed by upregulated miRNAs in age-related osteoporosis.  process. However, aside from Wagner et al. [5] revealed an upregulation of mir-371, mir-369-5p, mir-29c, mir-499 and let-7f upon in human MSCs replicative senescence, there was no reports for miRNA expression profiling of MSCs in aging process. In our results, compared with the young group, 30 down-regulated and 29 up-regulated miRNAs were identified in the old MSCs group. The overall level of changes was no significant difference. Due to the fact that one cell type failed to represent the overall miRNA alteration in one organism [15], so, next we compared the miRNAs expression changes between the young and old bone tissues. Among which in total 114 miRNAs and 45 miRNAs showed significant decreasing and increasing expressions, respectively. Notably, the majority of miRNAs were down-regulated (fold change $2) in bone, meaning that aged tissue environment might play important roles in bone aging process. After carefully screening step by step, 8 up-regulated and 30 down-regulated miRNAs in both MSCs and bone were confirmed. That is to say, during aging process, the bone exhibits predominant miRNA down-regulation, which was in agreement with other miRNA studies on aging tissues and organisms [16,45]. However, compared with other studies, some known SA-mi-RNAs [46], like miR-34a [47,48], miR-499 [5,49], miR-605 [50], miR-17 [16,51] and so on, there was no significant difference in our study. One possibility was the difference of tissue specific, or because of in our study only common miRNAs in MSCs and bone were considered. Possibly, the descend trend could be associated with Dicer, since recently study has provided evidence that the Dicer was declined with age [17], this might explain the decreases of multiple miRNAs in aging process. Additionally, other groups demonstrated that Dicer deficiency in osteoblast dominantly suppressed the bone formation [52,53], and further showed that Dicer generated miRNAs are essential for the bone formation. Therefore, it is speculated that the down-regulated miRNAs may play protective roles in the bone aging.
Large numbers of studies have reported that miRNAs regulate osteoblast differentiation and bone formation [42,54]. So, the next interesting point for us is how these changed miRNAs functioned in the bone aging process. When referring to literature, the mechanism of miRNAs worked by inhibiting gene expression, including positive and negative transcriptional factor and cytokine regulators, which are target genes of miRNAs [55], and regulate bone phenotype through certain biological processes [54]. We firstly used GO enrichment analysis to reveal the biological processes related with predicted target gene candidates of changed miRNAs. Results showed that the GO term transport (GO:0006810) significantly was enriched in the predicted target genes of both up-regulated and down-regulated miRNAs. Especially for the up-regulated miRNA related target genes, they enriched in GO terms regulation of response to stimulus (GO: 0048583) and cell differentiation (GO: 0030154). This result suggests us the up-regulated miRNAs might suppress some genes associated with the ability to respond on external stimulus. For example, the response of the cells to higher ROS (reactive oxygen species) levels, the MSCs from old mice might reduce ROS defense for the low expression level of related genes, further resulting in bone aging [13]. On the other hand, osteogenic differentiation of MSCs decline with age is clear, so, the differentiation-related gene was inhibited by the up-regulated miRNAs might be one of the reasons that caused age-related osteoporosis. For example, miR-214, which displayed up-regulation, could target ATF4, a positive transcription factor regulated osteoblast function, to inhibit bone formation [56].
Some studies reported that miRNA can affect pathways involved in aging, including p53 pathway, IGF pathway, mTOR pathway and sirtuins [15]. In our study, we analyzed the KEGG pathways of target genes that were related with significant alterations in miRNA expressions. Cell cycle pathway, calcium signaling pathway and pathways in cancer were significantly targeted KEGG pathways of up-regulated miRNAs, thereinto, pathways in cancer were affect by down-regulated miRNAs as well. These results suggest that the up-regulated miRNAs inhibited their target genes, which might mediate cell cycle control and regulate calcium ion homeostasis, which were important factors influencing the MSCs aging and mineralisation. Moreover, cancer is an age-related disease as well, so the miRNA regulates the pathways in cancer was reasonable in our study. However, according to the Figure S1, almost all of genes in cancer pathways (ko05200) were targeted by up-and down-regulated miRNAs, especially for most down-regulated miRNAs, it was hard for us to make appropriate judgments on which pathways play more important roles during bone aging process. So, in this study, we selected one of validated miRNAs to study its enriched pathways of likely targets, and analyzed the genes with known roles in aging or osteogenesis.
MiR-31, which was down-regulated in our study, and some published studies reported the function of miR-31 in aging, radioresistance, cancer and cell differentiation. For example, Lynam-Lennon et al [57] found that miR-31 directly targeted 13 genes involved in DNA repair to decrease cellular defense against DNA damage. Goff et al [58] reported that miR-31 might be an inhibitor of the osteocyte differentiation pathway. Another study also showed that miR-31 was under-expressed in osteo-differentiated MSCs [43]. Rokah et al [59] found that miR-31 regulated cell cycle related genes, like cyclinD1. In our result, miR-31 was down-regulated during bone aging process, after KEGG enrichment analysis, we found that several negative regulator of Wnt signal pathway, such as DKK1, GSK-3b and CCND1 were targeted by miR-31. It has been recognized that Wnt signaling plays an important role in the regulation of bone formation [60] and aging [61]. Moreover, we also found cyclinD1 was a target gene of miR-31, which was a crosstalk of Wnt pathway and p53 pathway. It is noteworthy that p53 is well known as an age-related protein [62], and has ability to decrease bone mass [63]. As stated above, miR-31 might play an important role in regulating the agerelated osteoporosis through its target genes, which were regulators in Wnt and/or p53 signal pathways.
In conclusion, age-related osteoporosis is a complex process, so, it was impossible to specify individual miRNA or signal pathway has decisive influence on bone aging. It was important to know how the networks of miRNAs and signal pathways operate temporally; the research on the networks and roles of miRNAs and/or their targeting genes in age-related osteoporosis was our further object.    Table S4 Comparison of known miRNAs expression levels in MSCs between two samples determined by deep sequencing. *: greater than 2.0-fold changes and Pvalue,0.05; **: greater than 2.0-fold changes and P-value,0.01. (XLSX) Table S5 Comparison of known miRNAs expression levels in bone between two samples determined by deep sequencing. *: greater than 2.0-fold changes and P-value,0.05; **: greater than 2.0-fold changes and P-value,0.01. (XLSX)