The Whole-transcriptome Landscape of Diabetes-related Sarcopenia Reveals the Specific Function of Novel lncRNA Gm20743

Model-identification of diabetes-related sarcopenia in db/db mouse

Fifteen-week-old db/db mice showed a phenotype of severe obesity (Fig. 1a, b) with obvious hyperglycaemia (>25 mmol/L) (Fig. 1c) compared to normal control db/m. Morphological changes in GAS were observed by haematoxylin and eosin (H&E), Oil-Red O, and Masson’s trichrome staining compared to db/m mice (Fig. 1d–f). The average cross-sectional area (CSA) of muscle fibres was markedly reduced (Fig. 1g), and lipid droplet infiltration and increased fibrosis were observed in GAS of db/db mice, indicating evident pathological changes of diabetes-induced sarcopenia. The forelimb grip strength of db/db mice was significantly decreased (Fig. 1h), and the muscle mass and size of lower limbs were decreased (Fig. 1i) compared with db/m mice. These results indicate that muscle strength and mass were significantly decreased in db/db mice. Decreased bone mineral density (BMD), lower lean mass, and increased fat mass were observed by dual‐energy X‐ray absorptiometry (DEXA, Hologic Discovery A, Hologic Inc) analysis in the db/db group (Fig. 1j–l). Serum biochemical parameters were indicative of the systematic metabolic disorder. Triglycerides, glucose, total cholesterol, low-density lipoprotein cholesterol, and high-density lipoprotein cholesterol were increased significantly in db/db mice compared to db/m control (Supplementary Table 1). Morphological abnormalities of mitochondria, swelling of the endoplasmic reticulum, reduced muscular glycogen, and sarcomere damage were observed by transmission electron microscopy (TEM) in GAS of db/db compared with the db/m group (Fig. 1m, n). Our data showed that muscular weakness and loss with subcellular structural dysfunction are evident in db/db mice.

Compared to db/m group, db/db mice extremely developed obesity (a, b), hyperglycaemia (c) (n = 10). H&E (d), Oil-red O staining (e), and Masson staining (f) in GAS of db/db and db/m group. Scale bar, 20 μM. Red arrow indicates adipocytes; orange arrow indicates atrophic myofibres; green arrows indicate lipid droplets; black arrows indicate muscle fibrosis. g Quantitative analysis of the average cross-sectional area (CSA) of myofibres. h The forelimb grasping strength and i muscle mass and size of lower limbs were decreased in db/db mice compared with db/m mice. The whole-body composition analysis of dual‐energy X‐ray absorptiometry (DEXA) (j) displayed more fat mass, less lean mass (k) and decreased bone mineral density (BMD) (l) in db/db vs. db/m (n = 5 replicates). Transmission electron microscopy (TEM) of GAS from db/db and db/m mice (m, n). The long red arrow indicates mitochondrial abnormalities and the short one indicates normal mitochondria; the blue arrow indicates myelinbody; the yellow arrows indicate muscular glycogen; the long green arrow indicates swelling endoplasmic reticulum and the short one indicates normal endoplasmic reticulum. Data were expressed as mean ± S. D. *P < 0.05 **P < 0.01. Note: H&E: haematoxylin and eosin; GAS: gastrocnemius; CSA: average cross-sectional area; DEXA: dual‐energy X‐ray absorptiometry; BMD: decreased bone mineral density; TEM: Transmission electron microscopy.

Identification of dif-mRNAs in GAS of db/db mice

To identify new molecular mechanisms during diabetes-accelerated sarcopenia, total RNA from GAS of db/db and db/m mice was used for whole transcriptome sequencing. Raw sequencing data were filtered, and clean reads were mapped to the reference genome (reference species: Mus_musculus. Reference Genome Version: GCF_000001635.26_GRCm38.p6). Filtered clean reads were spliced into putative transcripts (Supplementary Table 2). An absolute log₂ FC value of 0.585 was used as the standard to confirm the dif-mRNAs (false discovery rate, FDR < 0.05). Based on that standard, 607 upregulated mRNAs and 1332 downregulated mRNAs were identified in db/db vs. db/m (Supplementary Data 1–2). The heatmap and a volcano plot (Fig. 2a, b) showed differential mRNA expression between db/db vs. db/m (FDR < 0.05). The top 20 dysregulated mRNAs are shown in Supplementary Table 3. The top5 up-regulated mRNAs are Ttll7, Sorbs2, Sh3rf2, Prelp, and Ppp1r3c, while the top5 down-regulated genes include Xirp1, Vldlr, Vgll2, Ugp2, and Ucp3 in Supplementary Table 3.

a Hierarchical clustering of mRNA sequencing showed mRNA differential expression profiles (|log₂ FC | < 0.585, FDR < 0.05) of GAS between db/db and db/m (NC) group. Each row indicates a single gene, and each column indicates the GAS tissue of a mouse. b The scatter plot was performed for the alteration of dif-mRNAs in GAS of db/db vs. db/m (n = 3). Red plots represent up-regulated mRNAs (log₂ FC > 0.585, FDR < 0.05) and green plots represent down-regulated genes (log₂ FC < −0.585, FDR < 0.05). The significant GO terms (q-value < 0.05) of mRNA expression variation between two compared groups were listed (c, d). The top significant KEGG pathway terms (q-value < 0.05) of dif-mRNAs were listed (e, f). Note: NC: normal control group. GAS: gastrocnemius; dif-mRNAs: differentially expressed mRNAs; DEGs: differential expression of genes; FC: fold change; FDR: false discovery rate; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes.

Functional prediction of differentially expressed mRNAs (dif-mRNAs)

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed to identify potential biological functions and implications. Dif-mRNAs were exhibited in three modules of the GO pathway enrichment analysis in GAS of the db/db group. Upregulated mRNAs were mainly enriched in the extracellular matrix (ECM) organization (GO:0030198), negative regulation of the Wnt pathway (GO:0030178), retinoic acid catabolic process (GO:0034653), inflammatory response (GO:0006954), and collagen fibril organization (GO:0030199) in the biological process category of GO terms (Fig. 2c). Downregulated mRNAs were primarily enriched in terms such as immune system process (GO:0002376), and skeletal muscle contraction (GO:0006936) (Fig. 2d). The top 10 GO terms associated with dif-mRNAs of db/db mice are shown in Supplementary Table 4–5. For dif-mRNAs, enriched transcription factors (TFs), including Nr4a1, Mafa, Jun (also known as AP-1), ATF5, Trp63, and Tbx1, were differentially expressed in the db/db group (Supplementary Table 6). These data showed that differential TFs might be essential for the progression of diabetes-induced sarcopenia.

Then, 607 upregulated mRNAs and 1332 downregulated mRNAs were enriched in dysregulated pathways by KEGG analysis. Our data showed that 607 upregulated genes were enriched in 204 signalling pathways, which were mainly enriched in the pathway of neuroactive ligand-receptor interaction (PATH: map04080), cAMP signaling (PATH: map04024), ECM-receptor interaction (PATH: map04512), regulation lipolysis in adipocytes (PATH: map04923) and protein digestion and absorption (PATH: map04974) (Fig. 2e). KEGG enrichment showed that 1332 downregulated genes were enriched in 249 significant pathways, including the pathway of muscle contraction (PATH: map04260), cell adhesion molecules (PATH: map04514), TNF signaling (PATH: map04668), arginine and proline metabolism (PATH: map00330), and glutathione metabolism (PATH: map00480) (Fig. 2f). The top 10 KEGG pathways associated with dysregulated mRNAs of diabetes-related sarcopenia are shown in Supplementary Table 7–8, indicating that these significant changes in pathways may contribute to the development of diabetes-induced sarcopenia.

Based on shared biological function, location, or expression of genes, Gene set enrichment analysis (GSEA) has a broad scope of enrichment analysis and discovers potential mechanisms in disease. The dif-mRNAs were associated with significant pathways in GSEA, as shown in Supplementary Table 9 (p-value < 0.05). The GSEA results showed that the top 2 upregulated mRNAs were closely related to neuroactive ligand-receptor interactions (PATH: map04080) and the renin secretion pathway (PATH: map04924) (Supplementary Fig. 1a, b). In contrast, the top 2 downregulated mRNAs were highly correlated with the pathway of oxidative phosphorylation (PATH: map00190), arginine and proline metabolism (PATH: map00330) (Supplementary Fig. 1c, d).

Pathway-act-network

To investigate deeper interactions in significant pathways, 370 dif-mRNAs enriched into 10 mainly significant pathways of KEGG were subjected to a visual pathway-act-network analysis (Supplementary Fig. 2a). A pathway-act network was constructed with 233 dif-mRNAs enriched with 10 significant pathways of GSEA (Supplementary Fig. 2b). These data showed that significant metabolic pathway regulation might be involved in the development of diabetes-induced sarcopenia.

Identification of differentially expressed lncRNAs (dif-lncRNAs) in GAS of db/db vs. db/m

According to the screening criteria, 306 upregulated lncRNAs and 554 downregulated lncRNAs were identified in db/db vs. db/m (Supplementary Data 3–4). All dif-lncRNAs were widely found on all chromosomes, including the X and Y sex chromosomes (Fig. 3a). These dif-lncRNAs were divided into four categories based on different genic loci: exonic sense, intronic sense, antisense, and intergenic (Fig. 3b). Heatmap and volcano plot were constructed (Fig. 3c, d). We also identified the top 20 dysregulated lncRNAs in db/db mice (Supplementary Table 10).

a Circus plots showed dif-lncRNAs on mouse chromosomes for db/db vs. db/m. The outer circle is the chromosome, the second circle is all the differential dysregulated lncRNAs, and the inner circle is the significantly dif-lncRNAs (|log₂ FC | > 1 q-value < 0.05). b Pie charts showed the percentage of four types (exonic sense, intronic sense, antisense, intergenic sense) of dysregulated lncRNAs in db/db mice and the controls. c The cluster heatmaps and d scatter plots showing the dif-lncRNAs in db/db and db/m mice, respectively (n = 3). Note: dif-lncRNAs: differentially expressed lncRNAs.

Co-expression network construction and Real-Time Quantitative PCR (RT-qPCR) validation

Based on the co-expression relationship of lncRNA-mRNA, we generated diabetes-induced sarcopenia‐specific lncRNA‐mRNA co-expression networks. Through using correlation analysis of the expression, 24 lncRNA (log2 FC > 1, q-value < 0.05 and expression >2) and 860 mRNAs (Pearson correlation coefficient, PCC > 0.95, p-value < 0.05) were used to construct a highly correlated co-expression network. Several core dif-lncRNAs in the co-expression network were validated by RT-qPCR in db/db and db/m mice. The results showed that 1700047G03Rik and Gm31814 were upregulated, while Gm20743, Gm35438, Gm36131, and A330074k22rik were downregulated in GAS of the db/db group compared with db/m mice (Fig. 4a–f). Accumulated evidence showed that the elevated serum PA level in T2DM is involved in the development of insulin resistance and implicated in skeletal muscle inflammation, oxidative stress, and mitochondrial dysfunction by causing lipotoxicity, which leads to skeletal muscle loss^21,22. In our results, PA-treated C2C12 myotubes decreased the diameter of and inhibited the differentiation level of myotubes using fusion index analysis (Fig. 5a, b). Atrogin-1, MuRF-1, and Mstn, which serve as muscle atrophy markers, were higher relative to NC (Fig. 5c–e). These results suggested that PA efficiently induced muscle atrophy in C2C12 myotubes. The expression of these six dif-lncRNAs was further confirmed in a PA-induced muscle atrophy C2C12 cell model by RT-qPCR (Fig. 6a–f). The results found that the expression of most lncRNAs was consistent with animal models and RNA-seq data except for Gm38141 (Fig. 6e), which may be due to the differences in tissue and cellular components. Overall, the RT-qPCR validation of dysregulated lncRNAs of co-expression network in diabetes-associated sarcopenia mouse models and PA-treated C2C12 cells provide us with further indications to characterize their functions in the development of diabetes-associated sarcopenia. In addition, we found that 3 core dif-lncRNAs named Gm20743 (log2FC = −1.57, q-value < 0.001), Gm35438 (log2FC = −1.90, q-value < 0.001) and 1700047G03Rik (log2FC = 6.34, q-value < 0.001) had the highest number of interactions in the co-expression network and may contribute to the development of diabetes-related sarcopenia (Figs. 7a and  8a).

The relative expression of dif-lncRNAs (|log2 FC | > 1, expression >2 and q-value <0.01) was detected by RT-qPCR in GAS of db/db and db/m mice (a–f). The relative expression of closely related dif-mRNAs involved in significant pathways was detected in GAS of db/db and db/m mice (g–n) (n = 5). Data were expressed as mean ± S.D. *P < 0.05 **P < 0.01. Note: RT-qPCR: Real-Time Quantitative PCR.

The decreased CSA and numbers of myotubes were observed in PA-treated C2C12 cells through fluorescently labeled MyHC by confocal microscopy (a) and light microscopy (b). The higher expression of Fbox32 (c), Trim63 (d), Mstn (e) suggested that PA-induced muscle atrophy in C2C12 cells. Data were expressed as mean ± S.D. *P < 0.05 **P < 0.01. Note: PA: palmitic acid; MyHC: myosin heavy chain.

The relative expression of dif-lncRNAs (|log2 FC | > 1, expression >2 and q-value < 0.01) was detected in PA-induced C2C12 cells compared with NC group (a–f). The relative expression of highly correlated dif-mRNAs involved in the co-expression network was detected in PA-induced C2C12 cells and the control group (g–n). Data were expressed as mean ± S.D. *P < 0.05 **P < 0.01.

a The lncRNA Gm35438 and Gm20743 and 137 mRNAs were constructed in a co‐expression network with a PCC > 0.95. The red color represents upregulated genes, and the green color indicates down-regulation expression. Moreover, the top significant pathway terms of dysregulated genes associated with Gm35438 and Gm20743were displayed according to KEGG analysis (b, c), respectively. Note: PCC: Pearson correlation coefficient.

a The lncRNA 1700047G03Rik and 45 mRNAs were constructed in a co‐expression network with a PCC > 0.95. b The top significant pathway terms of dysregulated genes associated with 1700047G03Rik were displayed according to KEGG analysis. Note: PCC: Pearson correlation coefficient.

The potential function of candidate lncRNAs

We performed functional prediction of the co-expressed dif-mRNAs highly correlated with the three candidate lncRNAs by KEGG analysis, respectively. Functional analysis revealed that Gm20743 may be involved in mitochondrial functions like oxidative phosphorylation and glutathione pathway (Fig. 7b). Through RT-qPCR verification, the key genes (like ATP1B1, Gclc, Gclm, and GPX4) in the co-expression network which contributed to oxidative phosphorylation and glutathione pathway were downregulated in vivo (Fig. 4g–j) and in vitro (Fig. 6g–j) consistent with the RNA-seq results. These results suggested that Gm20743 may be involved in the regulation of mitochondrial function and redox homeostasis in skeletal muscle. KEGG analysis found that the dif-mRNAs co-expressed with Gm35438 are mainly enriched in muscle contraction and adrenergic signaling pathways (Fig. 7c). The co-expressed contractile proteins Tnnc1 and TPM3 in these differentially enriched pathways were downregulated in db/db mice (Fig. 4k, l) and PA-treated C2C12 cells (Fig. 6k, l) compared to their control groups, and these results suggested that Gm35438 may be related to skeletal muscle contraction and strength. Through co-expression network and bioinformatic analysis, we found the dif-mRNAs co-expressed with 1700047G03Rik (the top upregulated lncRNA in Supplementary Table 10) were mainly enriched in the neuroactive ligand-receptor interaction and glycolipid metabolism pathway (Fig. 8b). The RT-qPCR validation found that co-expressed key catalytic enzymes ACSL3 and DGAT2 in the glycolipid metabolism pathway were upregulated in both animal (Fig. 4m, n) and cell models of diabetic sarcopenia (Fig. 6m, n). It is suggested that 1700047G03Rik may be involved in the intramuscular lipid deposition of diabetes-related sarcopenia.

The role of Gm20743 in the development of diabetic sarcopenia

In this study, we specifically focused on the role of lncRNA Gm20743 on mitochondrial function and redox homeostasis in skeletal muscle cells. The key reasons are (1) First, based on the criteria |log2 FC | > 0.585, FDR < 0.05, and expression >2, we ranked lncRNAs in terms of FDR and found that lncRNA Gm20743 ranked in the top 20 differentially expressed lncRNAs. (2) Secondly, through a systematic literature review, many studies found that mitochondrial dysfunction and oxidative stress play an important role in sarcopenia; (3) Moreover, our previous study found that insulin resistance, mitochondrial dysfunction, and oxidative stress were involved in the development of muscle atrophy in PA-treated C2C12 cells^23,24; (4) In addition, through functional prediction of dif-mRNAs and lncRNA-mRNA co-expression network. These results showed that the down-regulated mitochondrial oxidative phosphorylation, citrate cycle, and glutathione metabolism were the most significantly enriched pathways based on GSEA (Supplementary Figure 1c, e, f). The KEGG analysis also revealed that the genes highly correlated with Gm20743 were mainly enriched in down-regulated oxidative phosphorylation and glutathione metabolism pathways (Fig. 7b). The above findings suggested that Gm20743 may play an important role in regulating mitochondrial dysfunction and oxidative stress in diabetic sarcopenia.

To verify the role of Gm20743 in diabetes-related sarcopenia, through fluorescence in suit hybridization (FISH), our results found that the distribution of Gm20743 was found in both the cytoplasm and nucleus of skeletal muscle C2C12 cells (Fig. 9a). Gm20743 was knocked down in C2C12 cells by siRNA transfection. The expression of Gm20743 was significantly downregulated by Gm20743si-1, Gm20743si-2, and Gm20743si-3 in C2C12 cells (Fig. 9b), and Gm20743 was overexpressed in C2C12 cells by lentiviral transfection and the expression of Lv-Gm20743 was verified by RT-qPCR (Fig. 9c). After knocking down Gm20743, we found that the extent of knockdown of Gm20743si-1, and Gm20743si-2 was lower in C2C12 cells. To avoid off-target situations, we used Gm20743si-1 and Gm20743si-2 to knock down Gm20743 in C2C12 cells, respectively. Through the detection of Mito-Sox and 2’, 7’-dichlorofluorescein-diacetate (DCFH-DA), our results found that knockdown Gm20743 increased mitochondrial reactive oxygen species (ROS) (Fig. 9d, e) and intracellular ROS (Fig. 9f, g) in C2C12 cells with Gm20743si-1 and Gm20743si-2, while overexpression of Gm20743 significantly reduced mitochondrial ROS (Fig. 9h, i) and intracellular ROS (Fig. 9j, k) in PA-induced C2C12 cells. These results suggest that Gm20743 can regulate mitochondrial function and redox homeostasis in C2C12 cells. Previous studies found that mitochondrial dysfunction and oxidative stress prevent the proliferation and differentiation in skeletal muscle myoblasts^25,26,27. To further assess whether Gm20743 is involved in the cell proliferation and myotube differentiation in skeletal muscle myoblasts, the detection of ethynyl-2′-deoxyuridine (EdU) and myosin heavy chain (MyHC) staining was determined. In our results, we found that knockdown Gm20743 decreased the number of EdU-positive cells in C2C12 myoblasts compared to the siNC group (Fig. 10a, b), and overexpressed Gm20743 increased the number of EdU-positive cells in PA-induced C2C12 myoblasts (Fig. 10c, d), suggesting Gm20743 improves cell proliferation in C2C12 cells. In addition, we observed that knockdown Gm20743 significantly decreased the number and diameter of myotube in C2C12 cells compared to the siNC group (Fig. 10e, f), while overexpressed Gm20743 significantly increased the number and diameter of myotube in PA-induced C2C12 myoblasts (Fig. 10g, h), compared to Lv-vector control. Together, these results indicate that Gm20743 may promote cell proliferation and myotube differentiation by preventing mitochondrial dysfunction and oxidative stress in skeletal muscle cells.

a Fluorescence in suit hybridization was used to observe the distribution of Gm20743 in C2C12 cells. The expression of Gm20743 in C2C12 cells transfected with Gm20743-siRNA (b) and Lv-Gm20743 (c) was determined by RT-qPCR. Mito-Sox staining was applied to detect mitochondrial ROS levels in C2C12 cells with Gm20743-siRNA (d, e), scale bar = 10 μm. DCFH-DA staining was performed to evaluate intracellular ROS levels in C2C12 cells with Gm20743-siRNA (f, g), scale bar = 100 μm. Mito-Sox staining was applied to detect mitochondrial ROS levels in C2C12 cells with Lv-Gm20743 (h, i), scale bar = 10 μm. DCFH-DA staining was performed to evaluate intracellular ROS levels in C2C12 cells with Lv-Gm20743 (j, k), scale bar = 100 μm. Data were expressed as mean ± S.D. Compared with the siNC or Lv-vector (NC) group, *P < 0.05 **P < 0.01. Note: si-1: Gm20743siRNA-1; si-2: Gm20743siRNA-2; si-3: Gm20743siRNA-3; DCFH-DA: 2’, 7’-dichlorofluorescein-diacetate; ROS: reactive oxygen species.

EdU staining was measured to evaluate cell proliferation in C2C12 cells with Gm20743-siRNA (a, b) and Lv-Gm20743 (c, d), scale bar = 100 μm. Immunofluorescence staining of MyHC was determined to test myotube differentiation with Gm20743-siRNA (e, f) and Lv-Gm20743 (g, h), scale bar = 100 μm. Data were expressed as mean ± S.D. Compared with the siNC or Lv-vector (NC) group, *P < 0.05 **P < 0.01. Note: si-1: Gm20743siRNA-1; si-2: Gm20743siRNA-2; si-3: Gm20743siRNA-3; EdU: 5-ethynyl-2’-deoxyuridine; MyHC: myosin heavy chain.