Identification of DEGs and Functional enrichment analysis
The bioinformatics analysis flow chart of this study is illustrated in Fig. 1. The patients baseline features of each muscle sample were also collected from the GEO database and described in Additional file 2. The genes with a |log2FC|>1 and p < 0.05 between SP and NSP were defined as DEGs, and we identified 808 DEGs (207 downregulated and 601 upregulated), which is illustrated by a volcano plot and heatmap (Fig. 2a, 2b). To examine the different biological phenotypes and the potential mechanisms of the DEGs, GO and KEGG enrichment analyses were performed. “Epidermis development” for BP (GeneRatio = 0.163, p < 0.0001, count = 71), “cornified envelope” for CC (GeneRatio = 0.0513, p < 0.0001, count = 23) and “structural constituent of skin epidermis” for MF (GeneRatio = 0.014, p < 0.0001, count = 6) were the most significant GO items (Fig. 2c). And KEGG enrichment analysis showed “Viral protein interaction with cytokine and cytokine” (GeneRatio = 0.0497, p = 0.0003, count = 9) was the most significant item (Fig. 2d). Besides, several immune response processes, such as “humoral immune response”, “cytokine activity”, “interleukin-1 receptor binding”, “IL-17 signaling pathway” were also significantly different in GO analysis and KEGG analysis (Additional file 3), indicating that immune-related mechanisms were involved in the sarcopenia.
Identification of potential signaling pathways in sarcopenia
To uncover the deeper mechanism underlying sarcopenia, ORA enrichment analysis was performed. We explored the alteration of sarcopenia in 8 major collections of MSigDB and analyzed the differentially expressed gene sets systematically between the NSP and SP. By using ORA, we identified 288 differentially expressed gene sets (Additional file 4) and the bubble plot represents the top 2 or 4 differentially expressed genes sets for each collection (Fig. 3a). Furthermore, many immune related gene sets, such as “ZHOU_INFLAMMATORY_RESPONSE_LPS_UP” (GeneRatio = 0.0303, p < 0.0001) in C2, and “GO_ANTIMICROBIAL_HUMORAL_RESPONSE” (GeneRatio = 0.0167, p < 0.0001) in C5, were also significantly up-regulated in sarcopenia (Additional file 4).
In order to avoid the disadvantages of artificially setting thresholds of DEGs in KEGG and ORA, GSVA was also performed to explore the significantly altered Hallmark pathways, and the expressional levels of gene sets for each sample was shown in the heatmap (Fig. 3b). Through GSVA, we identified 16 DEPs, including 14 upregulated and 2 downregulated pathways (Fig. 3c). “OXIDATIVE_PHOSPHORYLATION” (ES = 2.1680) and “G2M_CHECKPOINT” (ES=-1.8655) were mostly downregulated and upregulated in sarcopenia, respectively. Besides, levels of the “ESTROGEN_RESPONSE”, “APOPTOSIS” and “GLYCOLYSIS” were also significantly elevated in sarcopenia (Fig. 3c).
Identification of DETFs, keyDEGs and DEPs in sarcopenia
Sequence data for 318 TFs was retrieved from the Cistrome database (http://cistrome.org/). The DETFs were screened with p < 0.05 and |log2FC|>0.7 by using EdgeR (Fig. 3d). We identified 4 DETFs, including 3 upregulated and 1 downregulated and the expression of DETFs in two groups were displayed in the volcano plot in Fig. 3e. Among the four DETFs, only the PAX5 exhibited significantly decreased expression while TP73, BCL11A and TFAP2C were increased in sarcopenia.
Then, the Pearson correlation analysis was performed to screen the DETFs related key DEGs and final DEPs. The DEGs, DETFs and DEPs, whose correlation coefficient > 0.4 and p < 0.01 with each other, were selected for further analysis. Finally, 9 key DEGs, including 1 downregulated and 8 upregulated, and 11 final DEPs were included, and the results of DETFs and key DEGs were illustrated by the heatmap in Fig. 3f. The immune-related pathways “INFLAMMATORY_RESPONSE”, metabolic abnormalities-related pathways “PI3K_AKT_MTOR_SIGNALING”, “XENOBIOTIC_METABOLISM”, “HYPOXIA”, and apoptosis-related pathways “APOPTOSIS” were the most relevant among those 11 finals DEPs (Additional file 5).
Construction of novel transcriptional regulatory network with immune signature in sarcopenia
In order to explore the distinct involved immune response events or immune cells, the ssGSEA algorithm is applied to calculate expressional levels of immune gene sets of each skeletal muscle biopsy. In this study, we analyzed 29 immune-associated gene sets, and the expression for each sample was displayed by the heatmap in Fig. 4a. Finally, 1208 differentially expressed immune response related events or immune cells were identified with |correlation|>0.4 and p < 0.01 (Additional file 6).
To determine the immune-related transcriptional regulatory network in skeletal muscle with sarcopenia, Pearson correlation analysis was also performed to identify the relationship between immune gene sets, DETFs, key DEGs and DEPs. The interaction pairs of immune gene sets with correlation coefficients > 0.4 and p < 0.0001 were included in the subsequent analysis. The interaction pairs of immune gene sets with correlation coefficients > 0.4 and p < 0.0001 were included in the subsequent analysis. We identified 6 aberrantly expressed immune gene sets in sarcopenia from 29 immune gene sets, including 5 types immune cells and 1 immune response. The 6 immune gene sets correlated with sarcopenia were iDCs (correlation coefficient = 0.698, p < 0.001, positive), mast cells (correlation coefficient=-0.621, p < 0.001, negative), Th2 cells (correlation coefficient = 0.488, p = 0.001, positive), TIL (correlation coefficient = 0.460, p = 0.003, positive), B cells (correlation coefficient = 0.437, p = 0.005, positive) and type II IFN response (correlation coefficient = 0.511, p < 0.001, positive). Finally, the co-analysis result for DETFs, key DEGs, final DEPs and immune gene sets was shown by the co-expression heatmap in Fig. 4b. In order to show our results more clearly, a network was plotted by Cytoscape 3.7.1. Finally, sarcopenia-related hypothesis built on the bioinformatics was displayed by signaling diagram Fig. 4c.
Identification of specific inhibitors by C-Map
To provided potential treatment strategy for sarcopenia, the C-Map was utilized to identify potential inhibitors of network. We screened 1309 revealed bioactive compounds in C-Map database (Additional file 7) and 18 potential specific inhibitors targeting the network was identified with p < 0.01 (Fig. 5a). Among those compounds, “tanespimycin” (enrichment=-0.403, p < 0.0001), “trichostatin A” (enrichment=-0.365, p < 0.0001), “vorinostat” (enrichment=-0.568, p = 0.00038), “thioridazine” (enrichment = 0.396, p = 0.00242) were the most potential inhibitors in sarcopenia according to the results of C-Map. Based on the clue database, the detail information of Thioridazine (Fig. 5b), Trichostatin A (Fig. 5c), Vorinostat (Fig. 5d) and Tanespimycin (Fig. 5e) were further confirmed. Finally, the compound of Tanespimycin was mostly potential medicine for sarcopenia according to the results of literature review and C-Map.
Clinical characterization of patients with and without sarcopenia
In order to support our hypothesis, the histological and morphological analyses of skeletal muscle between SP and NSP were performed. The biopsy specimens were obtained from a prospective cohort study, and 12 participants were diagnosed with sarcopenia, and 12 were non-sarcopenia according to the 2019 AWGS (Additional file 8). The muscle mass and function were significantly difference as grip strength, lean body mass and ALMi for NSP and SP were 34.2 ± 3.6 kg vs 23.0 ± 2.8 kg, 43.5 ± 2.7 kg vs 38.8 ± 1.6 kg, and 7.6 ± 0.3 kg/m² vs 6.3 ± 0.3 kg/m², respectively. The laboratory test results revealed that the sarcopenia patients had higher levels of inflammation and worse function reserve as the C-reaction protein, leukocyte, hemoglobin, albumin for NSP and SP were 24.3 ± 7.0 g/dL vs 7.2 ± 10.7 g/dL, 6.8 ± 1.2* 109 /L vs 9.0 ± 1.9*109 /L, 137.0 ± 10.5 g/L vs 122.7 ± 9.9 g/L, 41.8 ± 2.3 vs 38.9 ± 1.8*109 /L, respectively. There was no difference in the age, BMI, NRS2002 score, platelets, and the content of fat mass between two groups, and all patients’ demographics information was summarized in Table 1.
Table 1
Clinical characteristics of patients with or without sarcopenia.
Variables | Non-sarcopenia (n = 12) | Sarcopenia (n = 12) | p-value (Control vs Sarcopenia) |
Age (years,Mean ± SD) | 74.8 ± 4.7 | 75. 4 ± 4.1 | N.S. |
BMI (kg/m²,Mean ± SD) | 21.9 ± 1.4 | 21.8 ± 1.9 | N.S. |
Grip strength (kg,Mean ± SD) | 34.2 ± 3.6 | 23.0 ± 2.8 | < 0.001 |
Lean body mass(kg,Mean ± SD) | 43.5 ± 2.7 | 38.8 ± 1.6 | < 0.001 |
Fat Mass(kg,Mean ± SD) | 15.5 ± 1.1 | 15.9 ± 1.1 | N.S. |
ALMi(kg/m²,Mean ± SD) | 7.6 ± 0.3 | 6.3 ± 0.3 | < 0.001 |
NRS2002 score(Media,IQR) | 2,1.0 | 2,1.8 | |
History of tobacco use(n) | 3 | 4 | |
History of alcohol use(n) | 2 | 3 | |
Drugs (n) | | | |
Antihypertensive drugs | 2 | 2 | |
Antihyperglycemic drugs | 1 | 2 | |
Antihyperlipidemic drugs | 1 | 1 | |
Laboratory data | | | |
C-reaction protein(g/dL) | 24.3 ± 7.0 | 37.2 ± 10.7 | 0.002 |
Leukocyte(10^9 /L) | 6.8 ± 1.2 | 9.0 ± 1.9 | 0.003 |
Hemoglobin(g/L) | 137.0 ± 10.5 | 122.7 ± 9.9 | 0.002 |
Platelets(10^9 /L) | 195.8 ± 68.5 | 222.2 ± 71.5 | N.S. |
Albumin(g/L) | 41.8 ± 2.3 | 38.9 ± 1.8 | 0.003 |
Identified of muscle fiber size and ultrastructural change in sarcopenia
The H&E analysis demonstrated significantly changes in the diameter of muscle fibers between two groups (Fig. 6a). In cross-sections, muscle fiber atrophy was observed, and the CSA was significantly decreased in SP (1026 ± 402µm2) when compared with the NSP (1687 ± 198µm2) (Fig. 6b).
TEM further revealed the aberrantly ultrastructure change of myofibrils (Fig. 6c). In NSP, the myofibrils were plump, well-organized, and attached closely with each other. However, in SP, myofibrils appeared atrophied mildly and the interspace between myofibrils (black asterisk) enlarged. Many myofilaments and sarcomeres became ruptured, fused, or disappeared (dark red arrows) in SP, and the number of lipid droplets (black arrows) between myofibrils also significantly increased. Besides, the size and shape of mitochondria varied greatly, and the percentage of abnormal mitochondria (white arrows) was also significantly increased in SP. All these anomalies in quantity, morphology, and distribution of myofibrils, mitochondria, and lipid droplets indicated the occurrence of muscle disfunction in SP.
Identification of DETFs and key DEGs in the study population by RT-qPCR
In order to re-predict the transcriptional regulatory pattern of DETFs and key DEGs, the RNA expressional levels of biopsies were detected by RT-qPCR. The RNA expression result of 4 DETFs showed that TP73 and TFAP2C were significantly upregulated, and the PAX5 was obviously downregulated in SP, while the BCL11A had no significant change between two groups (Fig. 6d). Then, we tested the RNA expressional levels of 9 key DEGs. As the results showed in Fig. 6e, the expression of KRT10, ACTG2, DSP, PERP, KRT2 and KRT1, which were associated with TP73, BCL11A and TFAP2C, increased significantly, while the expression level of SERPINA5, which might be positive regulated by PAX5, was significantly decreased. At the same time, the expression of PKP1, SFN, might be positive regulated by P73, BCL11A and TFAP2C, had no significant change between SP and NSP.
External validation of the regulatory mechanism with multiple online databases
For further validation of the regulatory mechanism between the DETFs and key DEGs, a comprehensive retrieval of public databases was performed to search for the presence of DETFs putative binding sites located within the transcriptional regulatory region of the key DEGs. Firstly, we found the promoters of key DEGs with NCBI, and made sure that the potential direct transcription factor binding motifs (TFBM) of PAX5, TFAP2C and TP73 were located in promoters of SERPINA5, ACTG2 and DSP respectively, and the results were illustrated in Fig. 6f. Next, JASPAR in-silico analysis was performed to figured out the binding site sequences of DETFs. The precise binding sites and the relative score were illustrated in Table 2 with a relative profile score threshold higher than 80%. In the promoter region of SERPINA5, ACTG2 and DSP, 4 putative binding sites for PAX5 and TP73, and 5 putative binding sites for TFAP2C were detected with the relative score from 0.80 to 0.85. All the key DEGs reported in the Table 2 show more than one putative binding site for DETFs, thus suggesting that all of them may be potentially bound by DETFs within their regulatory region.
Table 2
JASPAR analysis results for DETFs binding sites within the promoter of key DEGs.
DEG Name | TF Matrix ID | TF Name | Score | Relative score | Sequence ID | Start | End | Strand | Predicted sequence |
Serpina5 | MA0014.3 | MA0014.3.PAX5 | 8.354057 | 0.835880295 | NC_000014.9:94579426–94581426 | 1427 | 1438 | + | ATGCATGACTTT |
MA0014.3 | MA0014.3.PAX5 | 7.273882 | 0.813180180 | NC_000014.9:94579426–94581426 | 171 | 1182 | - | GTGCATGGCCGT |
MA0014.3 | MA0014.3.PAX5 | 6.896367 | 0.805246625 | NC_000014.9:94579426–94581426 | 228 | 239 | - | ATCCGTGTCCTA |
MA0014.4 | MA0014.3.PAX6 | 6.66627 | 0.800411082 | NC_000014.9:94579426–94581426 | 115 | 1176 | - | GGCCGTGATCAT |
ACTG2 | MA0815.1 | MA0815.1.TFAP2C | 6.9714537 | 0.833060697 | NC_000002.12:73891008–73893008 | 1861 | 1873 | + | CGCCTTGGGGGTA |
MA0815.1 | MA0815.1.TFAP2C | 6.731798 | 0.829739041 | NC_000002.12:73891008–73893008 | 1861 | 1873 | - | TACCCCCAAGGCG |
MA0815.1 | MA0815.1.TFAP2C | 5.5529466 | 0.813400009 | NC_000002.12:73891008–73893008 | 1601 | 1613 | + | AGACCCAGGGGCC |
MA0815.1 | MA0815.1.TFAP2C | 5.440579 | 0.811842579 | NC_000002.12:73891008–73893008 | 1601 | 1613 | - | GGCCCCTGGGTCT |
MA0815.1 | MA0815.1.TFAP2C | 4.90462 | 0.804414123 | NC_000002.12:73891008–73893008 | 1862 | 1874 | - | TCGCCTTGGGGGT |
DSP | MA0861.1 | MA0861.1.TP73 | 10.096636 | 0.850966062 | NC_000006.12:7539671–7541671 | 1861 | 1873 | + | CGCCTTGGGGGTA |
MA0861.1 | MA0861.1.TP73 | 8.66093 | 0.835814809 | NC_000006.12:7539671–7541671 | 1861 | 1873 | - | TACCCCCAAGGCG |
MA0861.1 | MA0861.1.TP73 | 8.087384 | 0.829762085 | NC_000006.12:7539671–7541671 | 1601 | 1613 | + | AGACCCAGGGGCC |
MA0861.1 | MA0861.1.TP73 | 6.944413 | 0.817700116 | NC_000006.12:7539671–7541671 | 1862 | 1874 | - | TCGCCTTGGGGGT |