Systematic analysis of long non-coding RNA and mRNA profiling using RNA sequencing in the femur and muscle of ovariectomized rats.

OBJECTIVE: To investigate the expression profile of lncRNAs in bone and skeletal muscle of ovariectomized (OVX) rats.
METHODS: Six-month-old female Sprague-Dawley (SD) rats were divided into OVX group (ovariectomized, n=12) and sham group (sham-operated, n=12). After 12 weeks, RNA-seq was used to analyze the differential expression of lncRNAs and mRNAs in femur and quadriceps between two groups. Dys-regulated expression of lncRNAs was confirmed by qRT-PCR. The cis and trans-regulatory functions were analyzed to determine their function and biological processes. Lastly, GO and KEGG analyses were performed to assess the biological relevance of genes in each profile.
RESULTS: A total of 17 lncRNAs and 440 mRNAs were differentially expressed in the femur. Thirteen lncRNAs and 292 mRNAs were differentially expressed in the quadriceps. qRT-PCR results were in consistent with the RNA-seq data. Among them, ENSRNOT00000090777 was found in both femur and quadriceps samples. Bioinformatics analysis found that LNC_004549 participated in the differentiation of skeletal and skeletal muscle.
CONCLUSIONS: The expression profile of lncRNAs was significantly altered in femur and quadriceps of OVX rat models, which may offer new insights into pathogenesis of osteoporosis and sarcopenia and potentially provide novel therapeutic targets.

Introduction

Osteoporosis (OP) is a major public health problem due to the association with fragility fractures[1]. Sarcopenia, a form of muscle loss associated with osteoporosis, which may elevate the risk of fracture by increasing falls risk[2,3]. It is also becoming an increasingly important health concern worldwide[4]. Unlike previous prevention of osteoporotic fractures merely focusing on bone, the effect of sarcopenia on the risk of osteoporotic fractures has received more attention[5]. Recently, experts have suggested integrating the pathogenesis of sarcopenia and OP and uniting them as a single therapeutic target[6].

Long noncoding RNA (lncRNA) is a class of transcripts that consist of more than 200 nucleotides in length and have limited protein-coding potential[7,8]. Many lncRNAs are differentially expressed in different tissues under various developmental and pathological conditions[9,10], which act as key regulators of biological and pathological processes[11]. Emerging evidence suggests that lncRNAs are involved in osteogenic[12]-[14] and osteoclast differentiation[15,16], skeletal muscle differentiation and regeneration[17]-[20]. However, these studies have focused on the physiological functions of lncRNAs at the cellular level, the pathological role of lncRNAs in bone and skeletal muscle, especially caused by postmenopausal OP, remains elusive. Therefore, the current research aims to investigate the expression profile and potential mechanism of lncRNA in bone and skeletal muscle in rat models with postmenopausal OP. In addition, OVX rat is the most commonly used model demonstrating bone loss, in which both sarcopenia and increased visceral fat occur simultaneously. Hence, OVX rat models were established and employed in this investigation.

In the present study, we analyzed the expression profile of lncRNAs and mRNAs in the femur and quadriceps of OVX rats. Gene Ontology (GO) annotation and KEGG pathway analysis were performed to identify potentially novel roles of lncRNAs in the pathophysiological process of OP and sarcopenia.

Materials and methods

Animals

Female Sprague-Dawley (SD) rats, aged 6 months, were obtained from the Laboratory Animal Center of Guangzhou University of Chinese Medicine (Guangzhou, China), and maintained in a 12h light-dark cycle at controlled temperature (22-24°C) and humidity (50-60%). Standard laboratory rodent diet and water were provided ad libitum. All rats were acclimatized for 7 days prior to formal experiments. Animals underwent dorsal ovariectomized (OVX, n=12) or sham operation (sham, n=12) according to FDA guidelines[21] and were left for 12 weeks to develop osteopenia in OVX. After 12 weeks, all 24 rat models were successfully established by testing bone mineral density (BMD). Then, the left femurs and quadriceps were removed and fixed with 10% neutral buffered formalin for 24 hours. Then, micro-computed tomography (micro-CT) and pathological staining of the left femurs were performed. The right femurs and quadriceps were sectioned and frozen in liquid nitrogen for RNA sequencing and qRT-PCR validation, etc. Euthanasia and disposal of carcasses were performed in accordance with the guidelines by the IAEC for animal experimentation. All experimental protocols were approved by the Ethics Committee of Guangzhou University of Chinese Medicine (No.2016046).

BMD

All rats were placed in a prone position after anesthesia, and BMD was measured by dual-energy X-ray absorptiometry with a small-animal high-resolution collimator (Discovery A/SL/W/C; Hologic, Bedford, MA, 01730 USA). BMD analysis was performed using small animal mode of the software (Hologic, Bedford, USA, v.13.2:3) and was calibrated at the beginning of each experiment. The volume of interest (VOI) was marked with respect to all lumbar vertebrae[22].

Micro-CT scan

The left femurs were measured by a cone-beam-type desktop micro-CT (µCT80, Scanco Medical, Brüttisellen, Zurich, Switzerland) and was evaluated by associated analysis software (µCT80 Evaluation Program v6.5-1, Scanco Medical, Switzerland)[22]. After micro-CT scan, cancellous bone of the distal femur at 1 mm above the growth plate was chosen as the VOI[23], which was restricted to an internal region of femur where trabecular and cortical bones were extracted by drawing freeform contours with the CT analyzer software.

The microstructure of cancellous bone was characterized using standardized techniques to determine the relative bone volume (BV/TV, %), the relative bone surface (BS/TV, 1/mm), trabecular thickness (Tb.Th, mm), trabecular number (Tb.N, 1/mm), trabecular separation (Tb.Sp, mm). The three-dimensional (3D) images were obtained through multi-planar reformation.

Pathological staining

The fixed left femurs were decalcified with ethylenediamine tetraacetic acid solution (EDTA) (Servicebio, G1105) for 4-6 weeks. The samples were dehydrated by standard graded alcohol solutions and embedded in paraffin. Bone tissues in paraffin were sectioned longitudinally to 4 µm sections. Slides were placed in Xylene to eliminate the paraffin at room temperature, followed by disposal of graded ethanol and distilled water. HE staining solution (Servicebio, G1005) was used for the HE staining. The stained femurs were observed and imaged under a light microscope (Nikon Eclipse E100, Tokyo, Japan).

RNA sequencing

Total RNA was isolated using Trizol reagent (GIBCO, Life Technologies, USA). The concentration of RNA was measured using Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA) and the integrity of RNA was assessed using the RNA Nano 6000 assay kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

A total amount of 3 µg RNA per sample was used as input material for the RNA sample preparations. Firstly, ribosomal RNA was removed by Epicentre Ribo-zero™ rRNA Removal Kit (Epicentre, USA), and rRNA free residue was cleaned up by ethanol precipitation. Subsequently, sequencing libraries were generated using the rRNAdepleted RNA by NEBNext® Ultra™ Directional RNA Library Prep kit for Illumina® (NEB, USA) according to the manufacturer’s recommendations. To select cDNA fragments of preferentially 150-200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). Then, 3 µl USER enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95°C before PCR. PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. At last, the obtained products were purified and library quality was assessed on the Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster kit v3-cBot-HS (Illumina) according to the manufacturer’s instructions. After cluster generation, the libraries were sequenced on an Illumina Hiseq 4000 platform and 150 bp paired-end reads were generated. The preparation of whole transcriptome libraries and deep sequencing were performed by Novogene Bioinformatics Technology Cooperation (Beijing, China).

RNA-Seq reads mapping and transcriptome reconstruction

The clean reads were mapped to the rat transcriptome sequence by bowtie2 v2.2.8, and the related algorithm in FPKM was used to estimate transcript abundance. Index of the reference genome was built using bowtie2 v2.2.8 and paired-end clean reads were aligned to the reference genome using HISAT2 v2.0.4. HISAT2 was run with ‘--rna-strandness RF’, other parameters were set as default.

The mapped reads of each sample were assembled by StringTie (v1.3.1)[24] in a reference-based approach. StringTie uses a novel network flow algorithm as well as an optional de novo assembly step to assemble and quantitate fulllength transcripts representing multiple splice variants for each gene locus.

Novel lncRNAs screening and differential expression identification

The following five steps were performed to acquire novel lncRNA: (1) The number of exons of the transcript ≥2; (2) The length of the transcript >200 bp; (3) Aligning the spliced transcript with the database using the Cuffcompare software, lncRNA overlapping the exon region of the spliced transcript was used as a database annotation lncRNA; (4) Calculating the level of FPKM (expected number of fragments per kilobase of transcript sequence per millions of base pairs sequenced) of each transcript by Cuffquant, and selecting a transcript of FPKM ≥0.5; (5) Using CNCI (Coding-Non-Coding-Index) (v2), CPC (Coding Potential Calculator) (cpc-0.9-r2) and PFAM (Pfam Scan (v1.3)) for coding potential prediction, transcripts with no coding potential in these software analysis results were used as lncRNA predicted by this analysis. Then, the expression of lncRNA and mRNA was calculated using the FPKM method, and the difference in the expression between samples was compared using Ballgown software. A P value of less than 0.05 was considered as statistical significance.

LncRNA expression verification with qRT-PCR

Differentially expressed lncRNAs obtained from RNA-seq were verified by qRT-PCR. All the primer sequences were shown in Table 1. Total RNA from the 20 rats (OVX 10, sham 10) was reverse transcribed to cDNA according to the manufacturer’s instructions (SuperScriptTM III Reverse Transcriptase, Invitrogen). qRT-RCR was performed on a ViiA 7 real-time PCR system (Applied Biosystems) using a 2X PCR master mix (Arraystar). The reaction conditions included the following three steps: pre-denaturation at 95°C, 10 minutes; 40 cycles of denaturation (95°C, 10 seconds) and annealing/extension (60°C, 60 seconds); melt curve established by 95°C, 10 seconds; 60°C, 60 seconds; 95°C, 15 seconds. All reactions were performed in triplicate. The relative fold change was calculated using the standard curve method normalized to β-actin.

Table 1

Sequences of primers used in the qRT-PCR.

Name	Primer	Product length (bp)
GAPDH(RAT)	F:5’ GCTCTCTGCTCCTCCCTGTTCTA3’	124
	R:5’ TGGTAACCAGGCGTCCGATA3’
ENSRNOT00000076147(F)	F:5’ CGGGCGAGGAAGAGCAAG 3’	119
	R:5’ TGGCACACCTGGGATAAAAGC 3’
LNC_005704(F)	F:5’ ATAGATGGTGTTGTCGGGGTC 3’	175
	R:5’ GGCAGTTAGAGGGAAGTGGAA 3’
LNC_004234(F)	F:5’ TTCCCATTCTGTTCCTTTCAC 3’	293
	R:5’ AAGCAAGCGACTTCAGTTCAT 3’
LNC_004120(F)	F:5’ TGGAGATGGAAGACAGAGACG 3’	159
	R:5’ TAAATGTGAGGACAAAGGGCA 3’
LNC_004549(F)	F:5’ AAATACAGAAACCAGCCTAAA 3’	119
	R:5’ GTCACACAGCAAGAGAATACG 3’
LNC_003963(F)	F:5’ ATTTCCCTGCTGGCTGAGGC 3’	172
	R:5’ AACTTCCCCCGCATTCTTTAA 3’
LNC_003965(F)	F:5’ CCTGCTCGCTGAGCATCCT 3’	243
	R:5’ TTCAGAGTTTGCGGGTTTGG 3’
ENSRNOT00000078855(M)	F:5’ ACCAGGTTGCCTAAGGAGGA 3’	91
	R:5’ ACAGGCGCGATCCCACTAT 3’
ENSRNOT00000088613(M)	F:5’ CAAAAGTTTATAGTGTGACTGGGAG 3’	88
	R:5’ CCTTGCCGTGATAGGAATGAG 3’
ENSRNOT00000089687(M)	F:5’ CCACTTCACATGAGCACCAT 3’	139
	R:5’ GTTCCTGAGAACTAAGCGACA 3’
ENSRNOT00000080982(M)	F:5’ ACCTCCAGTCAAATTCTCCTG 3’	86
	R:5’ AGCAACAGGGTCCTCTGAAAT 3’
LNC_002572(M)	F:5’ TGTTTCTGTATGTTGGTCCTG 3’	216
	R:5’ TGTATAAGCCACCTCATCTTT 3’
LNC_004788(M)	F:5’ TAAGCAGATACACTCAGTGGTCT 3’	138
	R:5’ AAGTTGGAAGGGTCAAGTAAG 3’
LNC_005429(M)	F:5’ CAGGTAGAGTCTGGTCGTCAC 3’	186
	R:5’ CTTCCTCCTACAAGACATCCAT 3’
LNC_004076(M)	F:5’ CATCTCTTGTTGATTTCACGTCT 3’	159
	R:5’ ATTCGGTGTCAGGAAGCATAG 3’

Abbreviations: qRT-PCR, quantitative real-time polymerase chain reaction; LNC, long noncoding; F, femur; M, muscle.

Gene Ontology (GO) and pathway analysis

First, the lncRNA target gene prediction was performed. Previous studies have shown that lncRNAs can regulate the expression of nearby protein-coding genes by cis-acting[10,25] and regulate distant genes by trans-acting[26]. Furthermore, the function of lncRNA can be predicted by its closely related cis or trans-regulated target protein-encoding genes[27]. Therefore, we selected protein-coding genes within 100 k upstream and downstream from lncRNA, and those with a Pearson’s correlation coefficient (PCC) greater than 0.95 as target genes to predict the biological function of lncRNA. Then, GO enrichment analysis was implemented by the GOseq R package, in which gene length bias was corrected. Pathway analysis was performed by mapping genes to the Kyoto Encyclopedia of Genes Genomes (KEGG) pathways. GO terms and pathway correlation with corrected P-value less than 0.05 were considered significantly enriched by differentially expressed genes.

Results

Establishment of rat models of OP induced by ovariectomy

The successfully established rat model of OP was verified by BMD, micro-CT scan and pathological staining. The values of BMD were significantly decreased in the OVX group as compared with those in the sham group (Figure 1A). As for bone microarchitecture parameters, the levels of BV/TV, BS/TV, Tb.Th and Tb.N were significantly decreased, while the Tb.Sp level was increased in the OVX group compared to that in the sham group (Figure 1A), suggesting the reduction of bone mass in OVX rats in contrast to the sham group. Moreover, HE staining of the distal femur in the OVX group revealed significantly reduced bone trabecula and a wider marrow space compared with those in the sham operation group, which was consistent with the micro-CT scan data (Figures 1B, C).

Figure 1

Establishment of the rat model of OP induced by ovariectomy. (A) The BMD level of lumbar vertebrae and BV/TV, BS/TV, Tb.Th, Tb.N and Tb.Sp were showed. (B) The 3-dimensional micro-CT images from different groups of distal femurs. (C) Representative photomicrographs of H&E-stained sections. The images are magnified 200×. The data are presented as the mean ± SD. *P<0.05, **P<0.01 vs the sham group.

Identification and features of lncRNAs and mRNAs

By RNA-seq, a total of 95 to 112 million and 81 to 116 million of clean reads were obtained from femurs and quadriceps femoris samples, respectively. A comparative genomic analysis of clean reads using HISAT2 showed that the percentage of reads aligned to the genome was higher than 92%. These aligned reads were spliced using StringTie to obtain the smallest possible transcript set. After that, the spliced transcripts were combined by cuffmerge software and those transcripts whose chain orientation was uncertain were removed. Finally, 301,462 and 261,138 transcripts were obtained from the femurs and quadriceps femoris samples for further analysis.

When a four-step analysis of the number of exons, the length of the transcript, the comparison with the known lncRNA database, and the expression of the transcript was performed on the transcripts obtained above, a total of 14130 and 11905 lncRNAs were screened in the femur and muscle, respectively. After that, in the fifth step, three tools including CNCI, CPC and PFAM were used to remove potential coding transcripts, and finally 7281 and 6267 putative non-coding transcripts were predicted as novel lncRNAs (Figures 2C, D). The screening results of each step are shown in Figures 2A, B.

Figure 2

Screening of lncRNA. Notes: charts showing the number of transcripts screened per step in femur (A) and quadriceps femoris (B). Venn diagrams depicting the common and unique numbers of noncoding transcripts predicted by four coding potential prediction software in femur (C) and quadriceps femoris (D). Abbreviations: CPC, coding potential calculator; PFAM, Pfam Scan (v1.3); CNCI, coding-non-coding-index;

In the present study, the structure and sequence conservation of lncRNAs and mRNAs were compared. Compared with mRNAs, the predicted lncRNAs (annotated lncRNAs and novel lncRNAs) from both femurs and quadriceps, were shorter in length, fewer in the number of exons, shorter in open reading frame (ORF) and lower in sequence conservation. All these show both the difference between lncRNA and mRNA, and also show that the predicted lncRNAs conform to the general characteristics of lncRNA. For novel lncRNAs from femurs, 70.5% were intergenic lncRNAs (lincRNA), 14.9% were antisense lncRNAs, and 14.7% were intronic lncRNAs. For novel lncRNAs from quadriceps femoris, 65.7% were intergenic lncRNAs (lincRNA), 14.1% were antisense lncRNAs and 20.1% were intronic lncRNAs, respectively.

Differential expression of lncRNAs and mRNAs

For femur samples, in the OVX group, 17 lncRNAs (9 up-regulated and 8 down-regulated) and 440 mRNAs (160 up-regulated and 280 down-regulated) were differentially expressed compared to the sham group. As shown in [Figures 3A and B], the differentially expressed mRNAs were widely distributed throughout almost each chromosome, and differentially expressed lncRNAs were mainly distributed on chromosome 1 and X.

Figure 3

Transcriptome profile of RNA-Seq data distinguishing two groups. Notes: (A) Volcano plot and chromosomal distribution map of differentially expressed lncRNAs in femur; (B) Volcano plot and chromosomal distribution map of differentially expressed mRNAs in femur; (C) Volcano plot and chromosomal distribution map of differentially expressed lncRNAs in quadriceps femoris; (D) Volcano plot and chromosomal distribution map of differentially expressed mRNAs in quadriceps femoris.

For quadriceps femoris samples, 13 lncRNAs (8 up-regulated and 5 down-regulated) and 292 mRNAs (101 up-regulated and 191 down-regulated) were differentially expressed in the OVX group compared to the sham group. As shown in [Figures 3C and D], the differentially expressed mRNAs were widely distributed throughout almost each chromosome, and differentially expressed lncRNAs were mainly distributed on chromosome 4 and 5.

Validation of RNA-Seq results by qRT-PCR

Since the number of differences in lncRNAs was small, qRT-PCR was performed to verify RNA-seq data for all differentially expressed lncRNAs. Among them, when designing the primers, ENSRNOT00000092808 and LNC_003121, LNC_000787, and LNC_000261, as well as ENSRNOT00000090777 and LNC_003964, LNC_003970, and LNC_003965 could not be distinguished, so these lncRNAs were removed. Finally, 8 lncRNAs from the femurs and 7 lncRNAs from the quadriceps were verified by qRT-PCR. qRT-PCR samples were obtained from 8 OVX rats and 8 sham-operated rats. As shown in Figure 4, the results of qRT-PCR were consistent with the RNA-seq data, and each lncRNA shared the same trend (up-regulated or down-regulated), indicating that the RNA-Seq analysis is reliable.

Figure 4

Validation of RNA-seq data by qRT-PCR. Notes: RNA-seq data of the femur and quadriceps muscles were verified by qRT-PCR. X-axis, the relative expression level of lncRNA; Y-axis, differentially expressed lncRNAs. Abbreviations: qRT-PCR, quantitative real-time polymerase chain reaction; RNA-seq, RNA sequencing.

GO and KEGG analyses

To predict the cis-regulatory function of lncRNA, protein-coding genes within 100 kb of the upstream and downstream of lncRNA were investigated. To predict the trans-regulatory function of lncRNA, Pearson correlation coefficient was calculated between the expression levels of lncRNAs and protein-coding genes. Table 2 shows the functional prediction results for the major dysregulated lncRNAs. To elucidate the functional categories of protein-coding genes associated with differentially expressed lncRNAs, GO enrichment analysis was performed on cis and trans-regulated target mRNAs. However, all the adjusted P-values were >0.05 and no significant enrichment was found. In addition, GO enrichment analysis was also performed on differentially expressed mRNAs. For the femur samples, as shown in Figure 5A, the significantly enriched terms mainly included biological process terms: organic substance metabolic process, primary metabolic process and cellular metabolic process; cellular component terms: cell part, intracellular and organelle; molecular function terms: protein binding, heterocyclic compound binding and organic cyclic compound binding. For quadriceps femoris samples, as shown in Figure 5B, the significantly enriched terms mainly included biological process terms: cellular metabolic process, carbohydrate metabolic process, single-organism carbohydrate metabolic process and generation of precursor metabolites and energy; cellular component terms: intracellular part, organelle and membrane-bounded organelle; molecular function terms: glyceraldehyde-3-phosphate dehydrogenase (NAD+) (phosphorylating) activity.

Table 2

Functional prediction of major dysregulated lncRNA.

Name	Tissue	log2 (fold change)	p-value	Target_prediction_method	Target_official_symbol
ENSRNOT00000078855	femur	-0.13	0.02	co_location	Neto2
ENSRNOT00000078855				co_location	Itfg1
ENSRNOT00000078855				co_location	Dnaja2
ENSRNOT00000078855				co_expression	Tax1bp1
LNC_005429		0.95	0.04	co_expression	Rgs5
LNC_005429				co_expression	Pcdhgb5
LNC_005429				co_expression	Atg16l2
LNC_005429				co_expression	Trmt1
LNC_005429				co_expression	Snrnp48
LNC_005429				co_expression	Celsr1
LNC_005429				co_expression	Prkcg
LNC_005429				co_expression	Rere
LNC_005429				co_expression	Safb2
LNC_005429				co_expression	Agbl5
LNC_005429				co_expression	Inppl1
LNC_005429				co_expression	Lrp6
LNC_005429				co_expression	Atp8b2
LNC_005429				co_expression	Fmn2
LNC_005429				co_expression	Rcbtb2
LNC_005429				co_expression	Klc2
LNC_005429				co_expression	Zfp335
LNC_005429				co_expression	Skil
LNC_005429				co_expression	Dgkd
LNC_005429				co_expression	Pbxip1
LNC_005429				co_expression	AABR07062799.2
LNC_005429				co_expression	Gprasp1
LNC_005429				co_expression	Car1
LNC_005429				co_expression	Krba1
LNC_005429				co_expression	Hivep3
LNC_005429				co_expression	Lgi4
LNC_005429				co_expression	Kmt2c
LNC_005429				co_expression	Dgka
LNC_005429				co_expression	Myb
LNC_005429				co_expression	Manea
LNC_005429				co_expression	Bcam
LNC_005429				co_expression	Map4
LNC_005429				co_expression	Rilpl1
LNC_005429				co_expression	Tp53bp1
LNC_005429				co_expression	Cacng7
LNC_005429				co_expression	Tub
LNC_005429				co_expression	Aida
LNC_005429				co_expression	Phf8
LNC_005429				co_expression	Scly
LNC_004234	quadriceps femoris	-0.43	0.04	co_location	AABR07062266.1
LNC_004234				co_location	AABR07062268.1
LNC_004234				co_expression	Xk
LNC_004549		4.05	0.02	co_location	Ndc1
LNC_004549				co_location	Hspb11
LNC_004549				co_location	Ldlrad1
LNC_004549				co_location	Lrrc42
LNC_004549				co_location	Dio1
LNC_004549				co_location	Yipf1
LNC_004549				co_expression	LOC100911516
LNC_004549				co_expression	Wdr36
LNC_004549				co_expression	Oxsm
LNC_004549				co_expression	Klf5
LNC_004549

Figure 5

The GO analysis of differentially expressed mRNAs. Notes: The histograms show the number of differentially expressed mRNAs enriched in each of the GO terms in the femur (A) and quadriceps (B), respectively.

To determine the most important biochemical metabolic pathways and signal transduction pathways involved in differentially expressed lncRNA-related protein-coding genes, KEGG pathway enrichment analysis was performed on cis and trans-regulated target mRNAs. In terms of femur samples, as shown in [Figures 6A, B and C], the major enrichment pathways of target genes that were cis-regulated by differentially express lncRNAs consisted of ribosome, endocytosis and mucin type O-glycan biosynthesis, the major enrichment pathways of target genes that were trans-regulated by differentially expressed lncRNAs included phosphatidylinositol signaling system, Fc gamma R-mediated phagocytosis and glycerolipid metabolism, the major enrichment pathways of differentially expressed mRNAs included protein processing in endoplasmic reticulum, spliceosome and Prion diseases. In terms of quadriceps femoris samples, as shown in Figure 6D, the significantly enriched pathways of target genes that were cis-regulated by differentially expressed lncRNAs were RNA transport, thyroid hormone synthesis and p53 signaling pathway, the major enrichment pathways of target genes that were trans-regulated by differentially expressed lncRNAs included biotin metabolism, fatty acid biosynthesis and primary immunodeficiency, the major enrichment pathways of differentially expressed mRNAs were HIF-1 signaling pathway, biosynthesis of amino acids and glycolysis/gluconeogenesis.

Figure 6

KEGG pathway analysis of dysregulated transcripts. Note: The top 20 significant pathways for co-localization and co-expression enrichment of lncRNAs in the femur are shown in A and B respectively and the top 20 significant pathways for differentially expressed mRNA in the femur are shown in C. The top 20 significant pathways for co-localization and co-expression enrichment of lncRNAs in the quadriceps femoris are shown in D and E respectively and the top 20 significant pathways for differentially expressed mRNA in quadriceps femoris are shown in F. The size of each circle stands for the number of significantly differentially expressed genes enriched in corresponding pathway. The rich factor was calculated using the number of enriched genes divided by the number of all background genes in corresponding pathway. Q value was calculated using the Benjamini–Hochberg correction. The pathway showing q value <0.05 are to be considered as statistically significantly over-represented.

Discussion

Sarcopenia is often accompanied by the development of osteoporosis[28]-[30]. Bone and muscle tissues are interconnected[31] and that both muscle mass and strength are positively associated with bone health[32]. The mechanical interactions between bone and muscle have been well established, and a growing body of research suggests a potential cross-talk between bone and muscle[33], including a sophisticated network of paracrine and endocrine signals[34]. OVX rat is the most commonly used model demonstrating bone loss[35,36], in which sarcopenia and increased visceral fat occur simultaneously[37]. Therefore, 6-month-old rats were employed in this study to establish an OVX model to simulate the physical condition of postmenopausal women. After the model was successfully established, the femurs and quadriceps were taken for RNA-seq. A total of 17 lncRNAs (9 up-regulated and 8 down-regulated) and 440 mRNAs (160 up-regulated and 280 down-regulated) were differentially expressed in femur samples and 13 lncRNAs (8 up-regulated and 5 down-regulated) and 292 mRNAs (101 up-regulated and 191 down-regulated) were differentially expressed in quadriceps femoris samples of OVX rats compared with sham rats. qRT-PCR results were consistent with the RNA-Seq results and showed the same trends of changes in expression levels, confirming the reliability of the RNA-Seq data.

Functional predictive analysis of differentially expressed lncRNAs revealed some interesting lncRNAs. In the femur sample, LNC_005429, up-regulated by 1.92-fold, could regulate the autophagy-related protein 16 and LDLR class B repeat by trans, which are mainly involved in autophagy and Wnt signaling pathways. In the quadriceps sample, LNC_004549, up-regulated by 16.56-fold, could cis-regulate nucleoporin protein Ndc1-Nup, which is mainly involved in RNA transcription, and trans-regulate Zinc finger, C2H2, which is mainly involved in skeletal and skeletal muscle differentiation. Since LNC_004549 can participate in the differentiation of skeletal and skeletal muscle at the same time, it will be the focus of our next research. The lncRNA, ENSRNOT00000090777, which was found in both femur and quadriceps samples will be another research focus to investigate the relationship between bones and muscles.

GO and KEGG pathway enrichment analyses were performed to identify the biological functions of differentially expressed lncRNAs and mRNAs. No significant enrichment was found for the target genes of the differentially expressed lncRNA, whether in the femoral sample or the quadriceps sample. GO enrichment analysis of these differentially expressed mRNAs revealed that these genes are primarily involved in basal metabolism and signal transduction. The KEGG pathway enrichment analysis revealed that the major enriched pathways in the femoral sample were the ribosome, phosphatidylinositol signaling system and protein processing in endoplasmic reticulum. The major enriched pathways in the quadriceps sample were RNA transport, biotin metabolism, AMPK signaling pathway and glycolysis/gluconeogenesis. Previous studies have found that AMPK may participate in the stabilization of β-catenin by phosphorylation of Ser552 and interact with the Wnt pathway[38]. Activation of AMPK could induce the expression of BMP-2 to gain bone mass and improve bone quality[39]. Therefore, AMPK is considered to be a potential new target for the treatment of OP[40]. This study found that the differentially expressed lncRNA target gene in the quadriceps of OVX rats is enriched in the AMPK pathway, and it can be assumed that lncRNA from muscle may affect bone through the AMPK pathway, which requires further research.

In this study, differentially expressed lncRNAs and mRNAs in femur and quadriceps femoris of OVX rats were identified. The cis- and trans-regulatory functions were analyzed to determine their function and biological processes. GO and KEGG pathway analyses were performed to annotate their potential functions. Moreover, novel lncRNAs coded as LNC_004549 were predicted to participate in the differentiation of skeletal and skeletal muscle simultaneously, and the annotated lncRNAs coded as ENSRNOT00000090777 was identified in both femur and quadriceps samples, both of which may play a potential role in the interaction between bones and muscles. Our findings can serve as a foundation for the identification of functional lncRNAs in bone and muscle of OP.

Funding

This study was jointly funded by the Natural Science Foundation of China (nos. 81674004, 81673786), Science and Technology Projects of Guangdong Province (no. 2016A020216024).