The Skeletal Muscle Transcriptome Profile of Elderly Men with Metabolic Syndrome Based on Weighted Gene Co-Expression Network Analysis Open Access

Metabolic syndrome (MS) is a multicomponent metabolic disorder characterized by abdominal obesity, insulin resistance, hypertension, and hyperlipidemia [1]. MS increases the risk of type 2 diabetes [2] and cardiovascular diseases [3]. With the spread of the Western lifestyle across the globe, the prevalence of MS is rapidly increasing worldwide, and MS has become a genuinely global problem [1].

Skeletal muscle, the primary site of insulin-stimulated glucose metabolism and fatty acid metabolism, plays a vital role in developing insulin resistance and MS [4, 5]. Sarcopenia and deleterious metabolism are common in the elderly (EL) [6]. The skeletal muscle’s overall functional, structural, and biochemical alterations have been extensively studied for chronological aging, but the molecular mechanisms implicated remain specified [7]. In addition, skeletal muscle mass (SMM) is involved in the development of MS. One cohort study shows that an increase in relative SMM overtime prevented the development of MS [7]. The morphological and metabolic adaptations in the skeletal muscle were also present in high fat-feeding [8]. Therefore, it is necessary to identify the transcriptional characteristic underlying skeletal muscle of EL men with MS to explore the abnormal function and molecular mechanism of skeletal muscle in MS.

High-throughput sequencing and bioinformatics provide new ideas for studying molecular mechanisms and therapeutic targets of diseases [9]. In human and mouse studies, the transcriptome of skeletal muscle in MS was analyzed by bioinformatics methods [10, 12]. However, none of those studies had used weighted gene co-expression network analysis (WGCNA), which was used to describe the correlation patterns between expressed genes to find highly correlated gene clusters, known as modules, and identify the characteristics of gene modules and the hub genes in the modules [13]. WGCNA has been successfully applied in cancer [14], ischemic stroke [15], Alzheimer’s disease [16], Parkinson’s syndrome [17], and other diseases [18] to find the molecular mechanisms involved in the occurrence and development of diseases and to screen biomarkers or potential therapeutic targets.

Therefore, the present study analyzed the GSE136344 dataset to get more insight into the transcriptome characteristics of skeletal muscle of EL men with MS. GO and KEGG enrichment analyses, gene interaction network analyses, and WGCNA were performed on the differentially expressed genes to profiled functional and molecular mechanisms of skeletal muscle in the occurrence and development of MS.

The data discussed in this publication have been deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO Series accession number GSE136344 (https://www.ncbi.nlm.nih.gov/geo) [12]. The dataset included 12 healthy young (YO) adult men, 11 healthy EL men, and 7 EL men diagnosed with MS (SX) for at least 10 years. As evidenced by physical and clinical examinations, YO and EL subjects were all in good health. EL and SX men were selected from the PROgnostic indicator OF cardiovascular and cerebrovascular events cohort [19, 20]. PROgnostic indicator OF cardiovascular and cerebrovascular events is a unique cohort of EL subjects perfectly well characterized since 2001 for many clinical, behavioral, and biological criteria [19]. MS was diagnosed when 3 of 5 components occurred: waist circumference ≥102 cm, triglycerides ≥1.7 mmol/L, HDL cholesterol <1.04 mmol/L, blood pressure ≥130/85 mm Hg, and/or fasting glucose ≥5.6 mmol/L [12]. For each subject, a single-needle biopsy was taken under local anesthesia from the superficial portion of the left vastus lateralis using a percutaneous technique [21]. GSE146869 [22], GSE145412 [23], and GSE23561 [24] datasets were analyzed to discover the differentially expressed genes in blood samples of subjects with MS.

The limma package of R language was applied to identify the differential expression genes [25]. The mRNA meeting |FC|> 1.5 and p value < 0.05 were the required data. The heat map and volcano plot of the differential expression genes were plotted using pheatmap and ggplot2 packages.

GO and KEGG pathways were the biological sequence analysis methods that can effectively cluster functional genes into different biological processes (BPs), mainly used to study DNA and protein-related issues [26]. Next, the DAVID database (https://david.ncifcrf.gov/tools.jsp) was taken advantage of by performing GO and KEGG analysis on differential expression genes. These analyses were mapped with bioinformatics (http://www.bioinformatics.com.cn/). The p value and false discovery rate were controlled at the 0.05 threshold.

The gene interaction network analyses used the STRING database version 11.0 (http://string-db.org) [27]. The interaction network of differential expression genes was visualized by Cytoscape 3.6.0 software (http://www.cytoscape.org/).

The WGCNA package of R language was applied to perform WGCNA on the differential expression genes. This study included 65 differentially expressed genes in the weighted gene co-expression network model. The optimal soft threshold (power) was 7, and the min module size was 5. The interaction network of intramodular genes was visualized by Cytoscape 3.6.0 software.

R 3.5.1 was used for statistical analysis. Values were means ± SE. p values < 0.05 were considered statistically significant.

The GSE136344 dataset included three groups of men: healthy YO, healthy EL men, and EL men with MS (SX). Compared with the YO group, 2,170 differentially expressed genes were identified in the EL group, including 892 upregulated genes and 1,278 downregulated genes (Fig. 1a). In the SX group, 765 differentially expressed genes were identified compared with the EL group, including 447 upregulated genes and 318 downregulated genes (Fig. 1b). In total, 2,131 differentially expressed genes were identified in the SX group compared with the YO group, including 996 upregulated genes and 1,135 downregulated genes (Fig. 1c). There were 36 co-upregulated genes (Fig. 1d) and 29 co-downregulated genes (Fig. 1e) to represent the confect of aging factor and metabolic disorder factor on gene expression. The 65 co-differentially expressed genes were visualized by a heat map (Fig. 1f). With this 65 co-differentially expressed genes, functional and pathway enrichment analyses, gene interaction network analyses, and WGCNA were further performed. In addition, among the 65 co-differentially expressed genes, 9 differentially expressed genes were found in both skeletal muscle and blood samples of subjects with MS, including F13A1, GDA, INCA1, KBTBD2, LEPR, LGMN, SLPI, SMAP1, and TUBA8 (online suppl. Fig. S1; for all online suppl. material, see www.karger.com/doi/10.1159/000530216).

GO analyses and KEGG analyses were performed to summarize the functional and pathway enrichment of the co-differentially expressed genes. In GO analyses, upregulated genes were enriched in 23 BP terms (Fig. 2a). Among the 23 BP terms, metabolism-related BPs were enriched, such as response to reactive oxygen species, purine nucleotide catabolic process, superoxide anion generation, response to fatty acid, superoxide metabolic process, reactive oxygen species metabolic process, response to peptide hormone, and cellular response to glucose stimulus. Immune-related BPs, including positive regulation of viral entry into host cells and response to cytokine, were also enriched. The function of muscle was affected during MS by affecting BPs including positive regulation of stress fiber assembly, aging, positive regulation of ERK1 and ERK2 cascade, positive regulation of I-kappaB kinase/NF-kappaB signaling, cellular response to transforming growth factor beta stimulus, cellular response to cAMP, and apoptotic process. The co-downregulated genes were enriched in 2 BP terms, including positive regulation of cytokinesis and fatty acid metabolic process. Enrichment analysis of the KEGG pathway showed that co-differentially expressed genes were mainly related to 3 pathways (Fig. 2b), including AMPK signaling pathway, carbon metabolism, and cytokine-cytokine receptor interaction.

In the interaction network analyses, the 65 co-differentially expressed genes were filtered into the one interaction network complexes, containing 29 nodes and 29 edges (Fig. 2c). The CAT, CYBB, IFI30, TIMP1, and S100A10, interacted with at least four co-differentially expressed genes, were defined as the hub genes. Among the 65 co-differentially expressed genes, CAT, MDH2, and PFKP were enriched into carbon metabolism pathway. PFKP, LEPR, and IRS1 were enriched into cytokine-cytokine receptor interaction pathway.

Next, WGCNA was further performed to filter the 65 co-differentially expressed genes into different modules. First, there was no notable outlier according to the sample dendrogram and trait heat map (Fig. 3a). Therefore, all the patients were included in the weighted gene co-expression network model. In this study, the optimal soft threshold (power) was selected as 7 (online suppl. Fig. S2A, B) and the min module size was 5. By average linkage hierarchical clustering, five modules were identified and represented in different colors, including yellow module, blue module, green module, turquoise module, and brown module (Fig. 3b). The blue module and yellow module were positively correlated, while both two were negatively correlated with the other three modules (online suppl. Fig. S2C).

Based on the WGCNA results, the interaction network of intramodular genes is presented in Figure 3c–g. In the blue module, CYBB, RNASE1, and LGMN were the hub genes, by the most interaction (Fig. 3c). In turquoise, CTNNBIP1 and ST3GAL3 were the hub genes (Fig. 3d). NNMT and NPNT were the hub genes of the brown module (Fig. 3e). CXorf57 and SLPI were the hub genes of the yellow module (Fig. 3f). CTGF and MAP10 were the hub genes of the green module (Fig. 3f). Therefore, 15 hub genes were found in this study (Fig. 3h), including 11 hub genes found in WGCNA, CAT, CYBB, IFI30, TIMP1, and S100A10, the hub genes found in gene interaction network analyses. The functions of hub genes in skeletal muscle of the MS are summarized in online supplementary Table 1.

In addition, the differentially expressed genes of SX compared with EL, which respect the differential expression in MS independent from aging, were also performed by WGCNA (online suppl. Fig. S3A, B). Seven modules were identified, including turquoise module, yellow module, blue module, red module, green module, black module, and brown module (online suppl. Fig. S3B, C). Turquoise module, yellow module, red module, blue module, brown module, and green module were correlated with MS. Black module, yellow module, and brown module were correlated with aging (online suppl. Fig. S3D).

In the present study, the GSE136344 dataset was analyzed to get more insight into the transcriptome characteristics of skeletal muscle of EL men with MS. Among the healthy YO man group, healthy EL men group, and EL men with MS group, 65 co-differentially expressed genes, including 36 co-upregulated genes and 29 co-downregulated genes, were found maybe regulated by age factor and MS factor. Those co-differentially expressed genes were enriched into 25 BP terms and 3 KEGG pathways. WGCNA was performed to filter the 65 co-differentially expressed genes into 5 modules. Fifteen hub genes, found in the interaction network analyses and WGCNA, may play an essential role in regulating the function of skeletal muscle of EL men with MS.

Current therapeutic options for MS are limited to individual treatments for hypertension, hyperglycemia, and hypertriglyceridemia, as well as dietary control measures and regular exercise [28]. Regular exercise and an increase in relative SMM prevented the development of MS [29]. During the development of aging and MS, 65 genes, including 36 co-upregulated genes and 29 co-downregulated genes, were differentially expressed in skeletal muscle. Similar to other studies, immune and metabolic functions were involved in developing MS [28, 30, 31]. In the present study, two immune functions, for example, positive regulation of viral entry into the host cell and response to cytokine BP, may involve in the developing of MS. Six metabolic functions, including cellular response to glucose stimulus, fatty acid metabolic process, purine nucleotide catabolic process, reactive oxygen species metabolic process, response to fatty acid, and superoxide metabolic process, may change in skeletal muscle of EL men with MS. A healthy skeletal muscle is characterized by the ability to switch easily between glucose and fat oxidation and will promote the oxidation of fat acid as a source for energy [10]. The changes in metabolic functions in skeletal muscle probably were the critical feature of MS. In addition, upregulated genes, LGALS1, ADAMTSL4, TNFRSF10B, and ZFP36L1, regulated apoptotic process; upregulated genes, CAT, TIMP1, and CTGF, regulated aging; upregulated genes, NOX4, CTGF, and S100A10, regulated a positive regulation of stress fiber assembly. Skeletal muscle is the primary site of insulin-stimulated glucose metabolism and fatty acid metabolism [4, 5]. Changes in BPs may lead to changes in muscle function and induce the development of aging and MS.

Skeletal muscle activity requires energy through anabolism and catabolism of glycogen, carbohydrates, and fat, all of which are important for energy storage and supply [32, 33]. Recent research reported that specific mRNA isoforms that changed significantly with age in skeletal muscle were enriched for proteins involved in adipogenesis [34]. In the present study, we found regulation of fat cell differentiation as one of the BP terms enriched in GO analyses of differentially expressed genes in EL men with MS.

In KEGG enrichment analysis, carbon metabolism was enriched, indicating abnormal metabolic functions of skeletal muscle in MS. Furthermore, the cytokine-cytokine receptor interaction pathway was enriched in KEGG enrichment analysis, indicating that inflammatory genes played essential roles in developing MS. In addition, the AMPK signaling pathway, enriched in our study, was reported as a function to alleviate MS and aging effect [35, 36].

Bioinformatics analysis mainly focuses on strong-effect genes and genes with known functions, but it is challenging to find weak-effect genes. However, WGCNA is to observe the function of genes with similar expression trends, which is an excellent supplement to weak-effect gene analysis [37]. A total of five modules, including blue module, turquoise module, brown module, yellow module, and green module, were found in the present study. Therefore, we can take these five modules as a whole as the characteristic targets of MS and further explore the molecular mechanism of metabolic disorder induced by the modules in skeletal muscle.

Among the fifteen hub genes, CAT[38, 40], CYBB[41, 45], IFI30[46, 47], TIMP1[48], NNMT[49], and CTGF[50] were reported regulated functions of skeletal muscle in MS. However, S100A10, LGMN, RNASE1, CTNNBIP1, ST3GAL3, NPNT, CXorf57, SLPI, and MAP10 were novel discovered genes that may also play an important role in the occurrence and development of MS. Future research should focus on the potential mechanisms of those fifteen hub genes in C2C12 cells, C. elegans, or mice models, to identify therapeutic targets for reducing MS risk.

In conclusion, 65 differentially expressed genes and 5 modules may regulate the function of skeletal muscle of EL men with MS, among which fifteen hub genes may play an essential role in the occurrence and development of MS. Our results contribute to a better understanding of the molecular events of the skeletal muscle in MS. Future research should focus on the potential mechanisms of those fifteen hub genes to identify therapeutic targets for reducing MS risk.

We thank the GEO database for sharing the data.

No ethical approval nor informed consent was required in this study due to the public availability of data in the GEO database.

The authors declare no conflict of interest.

This study was supported by the National Natural Science Foundation of China (No. 81872647), Natural Science Research of Jiangsu Higher Education Institutions of China (19KJB330007), Scientific Research Projects for Outstanding Teachers of Xuzhou Medical University (D2019011), and the Post-graduate Research and Practice Innovation Program of Jiangsu Province (No. KYCX20_2445).

Xing Ge was the principal investigator responsible for the study design and data definitions. Wang Qingqing conducted the study and provided the other data. Li-Chun Xu provided theoretical guidance and reviewed the manuscript. Other authors (Chaoran Zhang, Jiafei Xie, and Tingting Yao) participated in the literature collection, analysis, and summary. All authors have read and agreed to the published version of the manuscript.

The data discussed in this publication have been deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO Series accession number GSE136344 (https://www.ncbi.nlm.nih.gov/geo). Further inquiries can be directed to the corresponding author.

Xing Ge and Qingqing Wang are co-first authors.