Integrated analysis of long non-coding RNAs (lncRNAs) and mRNA expression profiles identifies lncRNA PRKG1-AS1 playing important roles in skeletal muscle aging

Overview of RNA-sequencing and identification of mRNA and lncRNAs in old group

The RNA-Sequencing data from 6 subjects were analyzed and a total of 583,406,044 raw reads were obtained. After quality control, 582,678,196 clean reads were left. The base average error rate of clean reads was 0.024%, and the average Q20 and Q30 values were 98.48 and 95.24%, respectively. The average GC content was 47.9% (Table 1).

A total of 424 DEGs were identified, including 271 up-regulated and 153 down-regulated genes (Figure 1A, Supplementary Table 1). In addition, 152 DElncRNAs including 76 up-regulated and 76 down-regulated genes were obtained (Figure 1B and Supplementary Table 2).

Figure 1. Volcano plot of differentially expressed genes (DEGs) (A) and differentially expressed lncRNAs (DElncRNAs) (B). The red dots represent upregulated genes or lncRNAs and blue dots represent downregulated genes or lncRNAs.

Functional enrichment analyses of DEGs

Based on the threshold of FDR < 0.05, we obtained 772 GO terms for the DEGs. The top significant GO terms were related with immune response, such as “immune system process” (gene count: 133; FDR = 0) and “regulation of immune system process” (gene count: 84; FDR = 0) and “immune response” (gene count: 94; FDR = 0) (Figure 2, Supplementary Table 3). Simultaneously, we obtained 15 significant KEGG pathways for the DEGs, including “Complement and coagulation cascades” (gene count = 21, FDR = 0), “Phagosome” (gene count = 15; FDR =0.0016), and hematopoietic cell lineage (gene count = 11; FDR = 0.0037) (Figure 3 and Supplementary Table 4).

Figure 2. Gene ontology enrichment analyses of DEGs. ‘Round’ represents biological process term and ‘tri-angle’ represents cell component term.

Figure 3. KEGG enrichment analyses of DEGs. ‘Round’ represents cellular processes term, ‘tri-angle’ represents environmental information processing, ‘square’ represents human diseases, ‘plus sign’ represents metabolism term and ‘cross within square’ represents organismal systems term.

Co-expression network of DElncRNAs and DEGs

Based on the Pearson correlation coefficient > 0.9, a co-expression network containing 7 DElncRNAs, 33 DEGs and 51 edges was built (Figure 4). Five upregulated lncRNAs, including AC004797.1 (degree = 17), PRKG1-AS1 (Protein kinase CGMP-dependent 1-antisense 1; degree = 16), MAPT-AS1 (microtubule associated protein tau –antisense 1; degree = 7), AC012254.3 (degree = 4) and CASC19 (degree = 3) and two downregulated lncRNAs, including AC022148.1 (degree = 2) and AC103740.1 (degree = 2) were included in this network. The DEGs of ITK (IL2 inducible T cell kinase) and DSC2 (desmocollin 2) could be positively regulated by the lncRNAs AC004797.1 and PRKG1-AS1.

Figure 4. Co-expression network of DEGs and DElncRNAs. The red ellipse represent upregulated gene and the green ellipse represent downregulated gene. The yellow rhombus represents upregulated lncRNA and blue rhombus represents downregulated lncRNA. The lines between the genes and lncRNAs indicate that there is a co-expression relationship between the two.

Validation of key genes by qRT-PCR

To further validate the reliability of RNA sequencing, we performed qRT-PCR on five DEGs, including SERPINE1 (Serpin family E member 1), OPRD1 (Opioid receptor delta 1), ITK, TXNRD1 (Thioredoxin reductase 1) and TDGF1 (Teratocarcinoma-derived growth factor 1), as well as five DElncRNAs, including CASC19, AC103740.1, AC004797.1, PRKG1-AS1 and GPRC5D (G protein-coupled receptor class c group 5 member D)-AS1. The RNA sequencing data showed that ITK, TXNRD1, CASC19, AC004797.1, and PRKG1-AS1 were upregulated, while SERPINE1, OPRD1, TDGF1, AC103740.1 and GPRC5D-AS1 were downregulated in the muscle of old group. As shown in Figure 5, the qRT-PCR results of ITK, TXNRD1, AC004797.1, PRKG1-AS1, SERPINE1, OPRD1, TDGF1, and GPRC5D-AS1 were in line with the RNA sequencing data, while no significant difference on CASC19 and AC103740.1 were detected between muscle of young and old group by qRT-PCR.

Figure 5. Validation of key DEGs and DElncRNAs by quantitative real-time polymerase chain reaction (qRT-PCR) in skeletal muscle of young group and old group. Difference between young group and old group was analyzed by students’ t test. * P < 0.05, ** P < 0.01.

Establishment of a dexamethasone-induced muscle atrophy cell model

To further investigate the roles of critical DElncRNAs, we established a muscle atrophy cell model by dexamethasone. The expression of myoblast determination protein 1 (MyoD), a skeletal muscle-specific bHLH transcription factor, was gradually decreased in a dose-dependent manner. The expression of MyoD at 10 nM and 15 nM dexamethasone treatment groups was significantly different from that of control group (P < 0.001), indicating the muscle atrophy cell model was successfully established (Figure 6A). The expression levels of three DElncRNAs, including PRKG1-AS1, AC004797.1 and GPRC5D-AS1 were determined under dexamethasone treatment. As expect, the expression levels of PRKG1-AS1 and AC004797.1 was increased in a dose-dependent manner, while GPRC5D-AS1 was decreased in a dose-dependent manner (P < 0.05, Figure 6B).

Figure 6. Validation of MyoD (A) and key DElncRNAs (B) by qRT-PCR in a dexamethasone-induced muscle atrophy cell model. Different concentrations of dexamethasone (Dex, 5 mM, 10 mM and 15 mM) were added in human skeletal muscle myoblasts and incubated for 48 h. Difference among groups was analyzed by ANOVA with Dunnett’s multiple comparison test. * P < 0.05, ** P < 0.01, ***P < 0.001, compared with control group.

Knockdown of PRKG1-AS1 increased cell viability and decreased cell apoptosis

We selected PRKG1-AS1 for further research. The expression of PRKG1-AS1 was knocked down by siRNA. qRT-PCR suggested that the expression of PRKG1-AS1 in siRNA1, siRNA2 and siRNA3 was all decreased and the decrease of PRKG1-AS1 in siRNA3 group was the most significant (P < 0.01, Figure 7A). Therefore, we selected siRNA3 for further experiment. CCK-8 assay showed that knock-down of PRKG1-AS1 could significantly increase cell viability in a time-dependent manner (P < 0.05, Figure 7B). Besides, flow cytometry demonstrated that cell apoptosis was significantly decreased by knocking-down of PRKG1-AS1 compared with the model group (P < 0.01, Figure 7C, 7D).

Figure 7. The effect of knock-down of PRKG1-AS1 on cell viability and cell apoptosis. (A) PRKG1-AS1 was knocked down by small interference RNA (siRNA) and the efficiency was detected by qRT-PCR. Difference among groups was analyzed by ANOVA with Bonferroni’s multiple comparison test. * P < 0.05, ** P < 0.01 compared with blank group; # P < 0.05, ## P < 0.01 compared with siNC (negative control) group. (B) Cell viability was tested by CCK-8 assay. Dexamethasone (15 mM) was added in human skeletal muscle myoblasts to establish atrophy cell model. Si-PRKG1-AS1 or siNC was transfected into human skeletal muscle myoblasts and incubated for 24 h, 48 h and 72 h. Cell viability was tested by CCK-8 assay. ** P < 0.01 compared with control group; # P < 0.05, ## P < 0.01 compared with model group; & P < 0.05, && P < 0.01 compared with siNC group. (C) Cell apoptosis was tested by flow cytometry. Dexamethasone (15 mM) was added in human skeletal muscle myoblasts to establish atrophy cell model. Si-PRKG1-AS1 or siNC was transfected into human skeletal muscle myoblasts and incubated for 48 h. (D) Quantitative analysis of cell apoptosis. * P < 0.05, ** P < 0.01 compared with control group; # P < 0.05 compared with model group; && P < 0.01 compared with siNC group.

Knockdown of PRKG1-AS1 affected mRNA and protein expression of muscle regulatory factors

In order to validate the effect of PKG1-AS1 at molecular level, we detected the mRNA and protein expression of muscle regulatory factors, including MyoD, MyoG, Myf2c and Myf5. As shown in Figure 8, these four factors were significantly decreased in the dexamethasone-induced muscle atrophy cell model (P < 0.05). After transfection with si-PRKG1-AS1, their expression was significantly upregulated both at mRNA level (P < 0.05). Western blot analysis showed consistent results with qRT-PCR, except for myf5, which showed no significant difference among groups.

Figure 8. The effect of knock-down of PRKG1-AS1 on muscle regulatory factors. (A) qRT-PCR analyzes gene expression of MyoD, MyoG, Myf2c and Myf5. Dexamethasone (15 mM) was added in human skeletal muscle myoblasts to establish atrophy cell model. Si-PRKG1-AS1 or siNC was transfected into human skeletal muscle myoblasts and incubated for 48 h. * P < 0.05, ** P < 0.01 compared with control group; ## P < 0.01 compared with model group; & P < 0.05, && P < 0.01 compared with siNC group. (B) Protein expression of MyoD, MyoG, Myf2c and Myf5 detected by western blot. Dexamethasone (15 mM) was added in human skeletal muscle myoblasts to establish atrophy cell model. Si-PRKG1-AS1 or siNC was transfected into human skeletal muscle myoblasts and incubated for 48 h.

Ethics statement

The study was performed according to protocols approved by the Ethics Committee of First Hospital of Jilin University (Changchun, China). Written informed consent for participating in this study has been received from all subjects.

Sample collection

The old and young skeletal muscle from 6 subjects (age, 17–81 years) was collected during surgery (3 samples in each group), and stored at -80° C. The baseline characteristics of the subjects are shown in Table 2. The inclusion criteria for enrollment were as follows: 1) not participated in exercise training before 1 week of the surgery; 2) without any disease directly affecting skeletal muscle tissue morphology and/or function; 3) overall healthy. Total RNA was isolated from each individual sample with TRIzol reagent (Invitrogen, USA). The concentration and purity of total RNA were measured by Nanodrop2000 (Thermo Fisher, Waltham, MA, USA), and the integrity was detected using agarose gel electrophoresis.

Library preparation for sequencing

Five μg RNA per sample was cleared for rRNA by beads Ribo-Zero Magnetic Kit (EpiCentre, Madison, WI, USA). The mRNA was randomly fragmented and single-strand cDNA was synthesized using random hexamer primer, the second strand cDNA synthesis was subsequently performed, and dTTP were replaced by dUTP. Then Illumina adaptor sequences were ligated to the end-repaired DNA fragments. The libraries were sequenced using the Hiseq2000 Truseq SBS Kit v3-HS (200 cycles).

Quality control analysis of original sequencing data

Raw reads of fastq format were cleaned to remove empty reads, adapter sequences and fragments smaller than 25 bp, non-unique oligonucleotide (AGCT) reads at the 5' end, reads with over 10% N sequences, and low quality reads, in which the number of bases with a quality value Q ≤ 10 was > 50%. In addition, Q20%, Q30% and GC% of the clean data were calculated.

Sequence alignment to reference genome and library quality assessment

Reads were mapped with Tophat (v2.0.9) to the human genome sequence (Ensemble GRCh38), and the Mapped Reads were analyzed with hisat2 (v 2.1.0, https://ccb.jhu.edu/software/hisat2/index.shtml). In addition, the sequence duplication, the non-uniform read distribution, the saturation for gene expression, the discreteness of insert sequence was evaluated with RSeQC (v2.6.4, http://rseqc.sourceforge.net/).

Evaluation of mRNA expression

The Stringtie v1.3.3 package (http://ccb.jhu.edu/software/stringtie/) normally used to process read alignments and the reference annotation was applied to estimate gene expression level per million mapped reads (FPKM) score.

Pearson's Correlation Coefficient (r2) used as an indicator of correlations between two independent biological replicates was calculated by plot_cor_exp (v1.1.0). The closer r2 is to 1, the stronger the correlation between the two replicate samples.

Analysis of DEGs

In the process of DEGs detection by edger (v3.24, http://www.bioconductor.org/packages/release/bioc/html/edgeR.html), |log2FC (fold change)| >1 and p value < 0.05 were used as screening criteria, and multi-test adjustment method (False Discovery Rate, FDR) was used to correct the p values. The clustering analysis of expression patterns was performed on DEGs using the distance calculation algorithm by plot_cluster_exp v1.1.0.

Gene ontology (GO) analysis of DEGs

GO (Gene Ontology, http://www.geneontology.org/) database was used to classify genes according to the biological processes, cellular components and molecular functions. The p value was corrected by Bonferroni, Holm, Sidak and FDR approach. When the corrected p-value (FDR) < 0.05, it is considered that there is a significant enrichment of this GO function.

KEGG pathway annotation of DEGs

To determine the most important biochemical metabolic pathways and signal transduction pathways involved in DEGs, the kegg_enrichment v2.1.0 was used with the Fisher's exact test. In order to control the false positive rate, the Benjamini–Hochberg (BH) procedure was used for multiple tests. KEGG pathway with corrected p value < 0.05 was defined as significantly enriched.

Identification of DElncRNAs

The lncRNA expression level of FPKM score was calculated by Stringtie based on the gene annotation information of the lncRNA. In addition, the Pearson's r2 score between samples was calculated by plot_cor_exp. DElncRNAs with |log2FC (fold change)| > 1 and p value < 0.05 were screened, and the FDR controlled multi-test correction based on BH method was used to correct p values. The clustering analysis of expression patterns was performed on DElncRNAs using the distance calculation algorithm by plot_cluster_exp v1.1.0.

DELncRNA-DEGs co-expression association analysis

Gene co-expression analysis could reveal the mechanism of transcriptional regulation. The interaction relationship was clarified via analyzing the PCC between DEGs and DElncRNAs in different samples. The interactions with PCC > 0.9 were filtered to construct a co-expression network using cytoscape software version 3.7.1.

Cell culture and treatment

Human skeletal muscle myoblasts were purchased from LONZA Pharma and Biotech (Tokyo, Japan) and were cultured in DMEM containing 1% penicillin/streptomycin and 10% fetal bovine serum (Thermo Fisher Scientific, Waltham, MA, USA) at 37° C and 5% CO₂. The medium was changed to high-glucose DMEM supplemented with 2% horse serum, 2% glutamine and 1% penicillin/streptomycin to induce myotubes. The medium was replaced every 48 h. On the fourth day after incubation with high-glucose DMEM, different concentrations of dexamethasone (Dex, 5 mM, 10 mM and 15 mM) were added and incubated for another 48 h. The expression of MyoD was determined by qRT-PCR to confirm the success of model establishment.

Cell transfection

Skeletal muscle myoblasts were divided into four groups: control group (cells without treatment), model group (cells treated with Dex), si-NC (negative control) group (cells treated with Dex and transfected with negative control siRNA) and si-PRKG1-AS1 group (cells treated with Dex and transfected with siRNA-PRKG1-AS1). Si-NC and si-PRKG1-AS1 were synthesized by Guangzhou Ribobio biotechnology Co., Ltd (Guangdong, China). Skeletal muscle myoblasts were seeded into 24-well plate (1 × 10⁵/well). Si-NC (20 pM) or si-PRKG1-AS1 (20 pM) were transfected into cells at 70% confluence by 1 μL Lipofectamine 2000 (Thermo Fisher).

Cell counting kit- 8 (CCK-8) assay

Skeletal muscle myoblasts were seeded on 96-well plate. 10 μ1 10% CCK-8 reagent was added to each well and incubated in the dark for 2 h. The absorbance at 450 nm was determined by a microplate reader (MK3, Thermo fisher).

Flow cytometry analysis

Cell apoptosis was determined by flow cytometry (FACSCalibur, BD, San Jose, CA, USA) as described in previous study [50]. Cell apoptosis was determined by an Annexin V-FITC/Propidium Iodide (PI) apoptosis kit (BD, USA). Briefly, skeletal muscle myoblasts at 48 h post-transfection were centrifuged at 200 ×g for 5 min and re-suspended in 1× Binding buffer. Cell suspension (100 μL) was transferred into test tube, then FITC-Annexin V (5 μL) and PI (5 μL) were added and incubated at room temperature (25° C) for 15 min in dark. Cell apoptosis was tested on flow cytometry within 1h.

qRT-PCR analysis for DElncRNA and DEGs

Tissue samples (50 mg) or skeletal muscle myoblasts of each group were extracted for total RNA by TRIzol reagent according to the manufacturer's instructions (TaKaRa, Dalian, China, Product code: 9109). Then reverse transcription reaction was performed for cDNA synthesis with PrimeScript™RT Master Mix (Perfect Real Time) (TaKaRa, Product code: RR036A). qRT-PCR was conducted to amplify genes using the following conditions: 50.0° C for 3min, 95.0° C for 3min, followed by 40 cycles of denaturation at 95.0° C for 10s and annealing-extension at 60.0° C for 30s. After reaction, the melting curve was evaluated by heating from 60° C to 95° C with 0.5° C for 10s increments. The primer sequences were shown in Table 3.

Western blot

Skeletal muscle myoblasts at 48 h post-transfection were added with RIPA lysis (Beyotime, Shanghai, China) to isolate total protein. Protein concentration was determined by BCA method (Thermo Fisher, USA), followed by separated on SDS-PAGE. Then, protein was transferred on to PVDF membrane (Millipore, USA), blocked with 5% skim milk and incubated with primary antibodies of anti-MyoD (Cal. No. 18943-1- AP, Proteintech, Rosemont, IL, USA; 1:1000), anti-MyoG (Cal. No. ab77232, Abcam, Cambridge, MA, USA; 1:1000), anti-Mef2c (Cal. No. 10056-1-AP, Proteintech, Rosemont, IL, USA; 1:1000), anti-Myf5 (Cal. No. ab125078, Abcam, Cambridge, MA, USA; 1:1000) and anti-GAPDH (Cal. No. 10494-1-AP, Proteintech, Rosemont, IL, USA; 1:1000) overnight at 4° C. On the second day, horseradish Peroxidase conjugated goat anti-rabbit IgG (H+L) (Cal. No. 111-035-003, Jackson ImmunoResearch, West Grove, PA) was added and incubated at 37° C for 2h. Chemiluminescence was developed by ECL system (Millipore, USA).

Statistical analysis

All experiments were repeated for three times. Experimental data were expressed as mean ± standard deviation (SD) and were processed in Graphpad prism 5(Graphpad Software, San Diego, CA, USA). P < 0.05 represented statistical significance.