Identification and validation of ferroptosis-Related genes in sarcopenia

doi:10.21203/rs.3.rs-3122342/v1

Download PDF

Article

Identification and validation of ferroptosis-Related genes in sarcopenia

https://doi.org/10.21203/rs.3.rs-3122342/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: Ferroptosis, characterized by iron accumulation and lipid peroxidation, leads to cell death. Growing evidence suggests the involvement of ferroptosis in sarcopenia. However, the fundamental ferroptosis-related genes (FRGs) for sarcopenia diagnosis, prognosis, and therapy remain elusive. This study aimed to identify molecular biomarkers of ferroptosis in sarcopenia patients.

Methods: Gene expression profiles were obtained from the Gene Expression Omnibus (GEO) database. Differentially expressed genes (DEGs) between normal and sarcopenia samples were identified using the "limma" package in R software. FRGs were extracted from GeneCards and FerrDB databases. Functional enrichment analysis determined the roles of DEGs using the "clusterProfiler" package. A protein-protein network was constructed using Cytoscape software. Immune infiltration analysis and receiver operating characteristic (ROC) analysis were performed. mRNA-miRNA, mRNA-TF, and mRNA-drug interactions were predicted using ENCORI, hTFtarget, and CHIPBase databases. The network was visualized using Cytoscape.

Results: We identified 46 FRGs in sarcopenia. Functional enrichment analysis revealed their involvement in critical biological processes, including responses to steroid hormones and glucocorticoids. KEGG enrichment analysis implicated pathways such as carbon metabolism, ferroptosis, and glyoxylate in sarcopenia. Totally, 11 hub genes were identified, and ROC analysis demonstrated their potential as sensitive and specific markers for sarcopenia in both datasets. Additionally, differences in immune cell infiltration were observed between normal and sarcopenia samples.

Conclusion: The hub genes identified in this study are closely associated with ferroptosis in sarcopenia and can effectively differentiate sarcopenia from controls. CDKN1A, CS, DLD, FOXO1, HSPB1, LDHA, MDH2, and YWHAZ show high sensitivity and specificity for sarcopenia diagnosis.

Biological sciences/Immunology/Cell death and immune response

Biological sciences/Computational biology and bioinformatics/Genome informatics

biomarkers

ferroptosis

genes

sarcopenia

Sarcopenia is an ascensive and generalized disease, playing a certain role in the process of muscle loss and organs dysfunction. Sarcopenia is correlated with increased adverse outcomes containing falls, functional decline, frailty, as well as mortality¹. Nowadays, sarcopenia has become a primary consideration of medical research, and it causes an increasing of health and public expenditure for its high morbidity and mortality ². At present, some methods of evaluating sarcopenia include walking speed, calf circumference, bio-impedance analysis (BIA), handgrip strength, dual-energy X-ray absorptiometry and imaging methods, whereas none of them are highly sensitive or specific for the evaluation of sarcopenia ^3,4. Thus far, universally diagnostic methods for sarcopenia have been lacking. Besides, existing research has suggested that IFG-1 has low efficacy, testosterone or other anabolic steroids have adverse effects, others like myostatin, angiotensin converting enzyme inhibitors, MT-102-are being investigated². Thus, it might be possible to identify new drug targets for sarcopenia by screening gene network.

Ferroptosis has been explored recently and been mostly combined with a mass of iron accumulation and lipid peroxidation amid the procedure of cell death⁵. Ferroptosis involves a wide variety of biology processes (e.g., iron metabolism, lipid metabolism, oxidative stress, and biosynthesis of nicotinamide adenine dinucleotide phosphate (NADPH), glutathione (GSH) and coenzyme Q10 (CoQ10)⁶). Latest evidence has revealed that excessive iron accumulation in muscle cells could lead to loss of muscle mass. In other words, ferroptosis plays a certain role in the occurrence and progression of sarcopenia ⁷. Nevertheless, the latent role of ferroptosis-related genes (FRGs) with the applications in diagnosis, prognosis and therapy for sarcopenia has not been illustrated in detail.

Microarray and bioinformatics analyses have been conducted at the genomic level for screening genetic alterations. Notably, gene expression signatures, pathways, biological process and underlying targets for sarcopenia have been identified through microarray analysis ^8,9, and the above studies screened complex molecular significantly correlated with sarcopenia and metabolism, (e.g., tumor growth factor-β (TGF-β), interleukin-12 (IL-12), milk fat globule epidermal growth factor 8 (MFG-E8)). In addition, the shared genes (e.g., UPF3A, CSTB and PEA15) and pathways (e.g., AGE-RAGE signaling pathway in diabetic complications, Amoebiasis, Platelet activation) between sarcopenia and type 2 diabetes were identified¹⁰. Moreover, sarcopenia and 5-mRNA risk modules may take on a critical significance in clinical diagnosis as combined prognostic factors ¹¹. However, most studies are integrative or single gene analyses that have not focused on specific phenotypes and immunity. Moreover, factors probably producing bias (e.g., different microarray platforms, different statistical methods, and insufficient sample sizes) can all contribute to inconsistent consequences across the above studies. Thus, further integrated analyses are required.

In this study, the Sarcopenia was downloaded from Gene Expression Omnibus (GEO) ¹², and significant DEGs were determined. The FRGs obtained from GeneCards database and FerrDb database. The intersection genes between Sarcopenia and ferroptosis were selected using the Venn diagram. Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, and Gene set enrichment analysis (GSEA) were performed to provide more information regarding the molecular mechanism. Subsequently, the protein-protein interaction (PPI) network was contrasted, and the hub genes were screened. Next, we predicted transcription factor (TF), microRNA, and drug targets of hub genes. Lastly, the abundance of immune infiltrating cells was calculated using single-sample gene set Enrichment analysis (ssGSEA) algorithm, and the correlation between the infiltration degree of various immune cells and the expression level of FRGs in the samples was analyzed. The above hub genes screened in this study can play a vital role in clinical diagnosis and therapeutic monitoring of sarcopenia. Moreover, the PPI network, the microRNA network, the TF-network, and the drugs-network built in this study can provide an efficacious framework for illustrating the molecular mechanisms of sarcopenia at the systems biological level. In other word, this study provides novel insights for future research in sarcopenia.

Data collection and preprocessing

Two microarray datasets of sarcopenia GSE1428¹³ and GSE8479¹⁴ were downloaded from the GEO database using “GEOquery” R package¹⁵. The species selected was Homo sapiens. GSE1428 was based on the platforms of the GPL96[HG-U133A] Affymetrix Human Genome U133A Array, containing the overall gene expression profile of vastus lateralis muscle of 22 samples, including 10 young subjects (19–25 years old) and 12 older subjects (70–80 years old). All samples were contained in this study. GSE8479 was based on the platforms of the GPL2700Sentrix HumanRef-8 Expression BeadChip. A total of 65 samples were included, including 26 muscle biopsies from the lateralis femoris muscle of young subjects before exercise or immobilization, lateralis femoris muscle biopsy samples from 25 aged subjects before or after exercise as well as 14 after exercise or fixation. A total of 26 samples of muscle biopsy of lateral femoral muscle of young subjects before exercise or fixation and 25 samples of muscle biopsy of lateral femoral muscle of the old subjects before exercise or fixation were selected for inclusion. The “limma” R package was selected to merge and normalize the raw data¹⁶.

To evaluate the role played by ferroptosis in sarcopenia, we obtained 619 FRGs from GeneCard database¹⁷ using the keyword “ferroptosis”. Furthermore, 474 FRGs were downloaded from FerrDb database¹⁸ (http://www.zhounan.org/ferrdb/current/). Lastly, there were 883 FRGs obtained after merging and dereplication in all.

Identification of the differentially expressed genes

The “limma” package was applied to normalize the dataset GSE1428 and dataset GSE8479 at first to identify the potential mechanism and relevant biological characteristics with pathways of DEGs in sarcopenia. Subsequently, the “limma” R package was adopted to the gene expression profile to filter DEGs of sarcopenia. P-value < 0.05 and | logFC | > 0 were set as the cut-off criteria to determine the significant DEGs. The genes with logFC > 0 and P-value < 0.05 were up-regulated genes, while those with logFC < 0 and P-value < 0.05 were down-regulated genes.

To obtain the FRGs associated with sarcopenia, the intersection of DEGs was first drawn from GSE1428 dataset and GSE8479 dataset. Ferroptosis-related differentially expressed genes (FRDEGs) shown in the Venn diagram were selected. The results are expressed as volcano maps by using “ggplot2” R package and heatmap using “pheatmap” R package.

Function enrichment analysis

Gene Ontology (GO) analysis has been used as a general method for mass functional enrichment, including biological processes (BP), molecular functions (MF) and Cellular Component (CC)¹⁹. Kyoto Encyclopedia of Genes and Genomes (KEGG)²⁰ is a extensively used database storing information of Genomes, biological pathways, diseases and drugs. The “clusterProfiler” R package²¹ was adopted to perform GO and KEGG enrichment analysis on FRDEGs. P-value < 0.05 and FDR value (Q.value) < 0.25 were set as the screened criteria and achieved statistical significance. The benjamini-Hochberg method (BH) was employed to correct the P value.

GSEA

GSEA (Gene Set Enrichment Analysis)²² refers to an algorithmic method to analyze whether a specific gene set is statistically different between two biological states. It is usually emlployed to estimate changes in pathway and biological process in expression profile samples. In this study, “clusterProfiler” R package was employed for the enrichment analysis of all genes in GSE1428 and GSE8479 datasets. The parameters used in this GSEA enrichment analysis are presented as follows. The number of seed was 2020, and the calculation times were 1000. The respective gene set contained at least 10 genes and at most 500 genes. The Benjamini-Hochberg method (BH) was employed for P-value correction. The C2.cp.v7.2. symbols gene set was obtained from Molecular Signatures Database (MSigDB)²³. The screening criteria for significant enrichment included P-value < 0.05 and FDR value (Q.value) < 0.25.

Protein-protein interaction (PPI) network

PPI network is a STRING database including individual proteins interacting with each other²⁴. This network is beneficial to search known proteins and predict interactions between proteins. In this study, the database named “STRING” was adopted to build the PPI network associated with the FRDEGs. Subsequently, Cytoscape was adopted to visualize the PPI network model. The tightly connected local regions of the PPI network may represent molecular complexes with specific biological functions. Using DEGREE (DEGREE Correlation), MNC (Maximum Neighborhood Component), MCC (Maximal Clique Centrality), EPC (Edge Percolated Component), and CLONESS (Total Centrality) in cytoHubba²⁵ plugin, five algorithms selected the top 15 total FRDEGs as the key genes (hub genes, mRNA).

mRNA-miRNA, mRNA-TF, mRNA-drugs interaction network

ENCORI²⁶ database (https://Starbase.sysu.edu.cn/) is the 3.0 version of starBase database. In ENCORI database, miRNA-ncRNA, miRNA-mRNA, ncRNA-RNA, RNA-RNA, RBP-ncRNA and RBP-mRNA interactions are based on CLIP-seq and degradation sequencing data mining, providing a wide variety of visual interfaces to discuss microRNA targets. miRDB database²⁷ was employed for miRNA target gene prediction and functional annotation. In this study, miRNAs interacting with hub genes (mRNA) were forecasted in accordance with ENCORI and miRDB databases, and the mRNA-miRNA interaction network was plotted using the intersection of the two databases.

CHIPBase database²⁸ (Version 2.0) (https://rna.sysu.edu.cn/chipbase/) identified thousands of DNA binding motif matrix and its binding sites, and predicted millions of Transcription factor (Transcription factors, TF) and transcriptional regulatory relationships between genes from the DNA binding protein ChIP-seq data. hTFtarget database ²⁹(http://bioinfo.life.hust.edu.cn/hTFtarget) refers to a comprehensive database of human TF and target regulation. TF binding to hub genes (mRNA) was searched through CHIPBase database (Version 2.0) and hTFtarget database. Moreover, DGIdb (The Drug Gene Interaction Database)³⁰ (HTTPS://www.dgidb.org) was employed to predict the potential drugs or small molecule compounds interacting with hub genes (mRNA). Moreover, Cytoscape software was adopted to visualize the mRNA-miRNA, mRNA-TF, and mRNA-drugs interaction networks.

The ROC Curve

Receiver Operating characteristic Curve (ROC) is a coordinate schema analysis tool employed to select the best model, discard the second-best model, or set the optimal threshold within the same model. ROC curve is a comprehensive index indicating the sensitivity and specificity of continuous variables, as well as the correlation between sensitivity and specificity through the composition method. In general, the area value under the ROC curve ranges from 0.5 to 1. The closer AUC to 1, the better the diagnostic effect will be. The accuracy is relatively low as AUC ranged from 0.5 to 0.7, while that is medium as AUC ranged from 0.7 to 0.9. If the AUC exceeds 0.9, there is a high accuracy. The “pROC” R package was adopted to plot the ROC curves of hub genes against different subgroups (normal and sarcopenia), and the area under the curve (AUC) was obtained to evaluate the diagnostic effect of hub genes expression on disease.

Immune infiltration analysis (ssGSEA)

Single sample gene set Enrichment analysis (ssGSEA) is an extension of the GSEA method, allowing the definition of an enrichment score that represents the absolute degree of enrichment of the gene set in the respective sample in a given dataset. The single-sample gene set enrichment analysis (ssGSEA) algorithm³¹ was adopted to quantify the relative abundance of each immune cell infiltration. Moreover, the correlation between the infiltration degree of various immune cells and the expression level of FRDEGs in the samples were analyzed.

Statistical analysis

All analyses were based on R software (Version 4.1.0), and the continuous variables are expressed as mean ± standard error of mean (SEM). Wilcoxon rank sum test was performed for comparison between the two groups. Kruskal-wallis test was performed for comparison of three groups or more. The comparison of progression-free survival between the two groups was drawn through Kaplan-Meier analysis combined with the log-rank test. If not specifically specified, the results are all P -values less than 0.05, which is the standard of significant difference results.

Data normalization and FRDEGs analysis

In this study, bioinformatics methods were primarily employed to explore the biological characteristics of sarcopenia, and the flow chart showing the overall analysis process is illustrated in Fig. 1. In the GSE1428 and GSE8479 datasets, the subjects were divided into a sarcopenia group and a normal control group. Two groups of datasets were normalized, annotation probe and other data cleaning operations were performed, and a boxplot was drawn for data distribution before and after standardization (Fig. 2A-B).

Figure 1. Flow chart of overall analysis of sarcopenia biological characteristics explored by bioinformatics methods

Figure 2. Boxplot of GEO dataset before and after correction.

(A). Box diagram of gene expression distribution between samples of GSE1428 dataset before correction. (B). Box diagram of gene expression distribution among samples of GSE1428 dataset after correction. (C). Box diagram of gene expression distribution among samples of GSE8479 dataset before correction. (D). Box diagram of gene expression distribution among samples of GSE8479 dataset after correction.

A total of 619 FRGs were downloaded by searching the keyword "Ferroptosis" in GeneCards database. Moreover, 474 FRGs were downloaded from FerrDb database (http://www.zhounan.org/ferrdb/current/). A total of 883 FRGs were obtained after merging and dereplication.

To analyze the differential gene expression values in the Sarcopenia sample group (Sarcopenia group) in comparison with the Normal sample group (Normal group) in the dataset GSE1428 and GSE8479 respectively, “limma” R package was used to acquire genes that were differentially expressed in the two datasets. The results are expressed as follows: 12548 DEGs were obtained in Dataset GSE1428, 1149 of them meet threshold | logFC | > 0, P-value < 0.05. Under the threshold, there are 603 high expression genes in the disease group (low expression in the normal group, logFC is positive, up-regulated genes), 546 low expression genes in the disease group (high expression and logFC is negative in the normal group). Volcano maps were drawn from the differential analysis results of this dataset (Fig. 3A). Totally, 16742 DEGs were selected in the dataset GSE8479, and 5039 of them meet the threshold | logFC | > 0, P-value < 0.05. Under the threshold, there are 2485 high expression genes in the disease group ((low expression in the normal group, logFC is positive, up-regulated genes), 2554 low expression genes in the disease group (high expression in the normal group, logFC is negative). Volcano map was drawn based on the results of differential analysis of this dataset (Fig. 3B).

As to obtain FRDEGs, all the DEGs selected from the datasets GSE1428 and GSE8479 were intersected (| logFC | > 0, P-value < 0.05). Next, 591 common differentially expressed genes (Co-DEGs) were acquired, and a Venn diagram was generated (Fig. 3C). We also intersected the Co-DEGs and FRGs in the datasets, and obtained a total of 46 FRDEGs, and a Venn diagram (Fig. 3D) was generated. The 46 FRDEGs included DDIT4, FTL, FTH1, HMGB1, HSPB1, CISD1, TP63, CDKN1A, PLIN2, ZFP36, SREBF2, FZD7, BRD3, DECR1, CP, GOT1, EX1, MPC1, ABCC5, USP11, CS, ATG5, MTCH1, SLC38A1, DLD, CIRBP, STOML2, LDHA, MDH2, ANP32B, FKBP3, LPIN1, PGK1, ANXA1, YWHAZ, ACTC1, LRPPRC, MAP4, HNRNPA3, MYL6B, RUFY1, WBP11, FOXO1, PRPF8, RPS27A, MAGED2. In accordance with the results of Venn diagram, the differential expression of FRDEGs among different sample groups in GSE1428 dataset (Fig. 3E) and GSE8479 dataset (Fig. 3F) was analyzed. The heatmaps were generated using “pheatmap” R package to display the analysis results.

Figure 3. Differential expressed gene analysis of GEO datasets.

(A). Volcano map of DEGs for GSE1428 dataset. (B). Volcano map of DEGs for GSE8479 dataset. (C). Venn diagrams of DEGs in GSE1428 and GSE8479 datasets. (D). Venn diagram of co-DEGs and ferroptosis-related genes in the dataset. (E). Complex numerical heat map of differentially expressed genes associated with iron death for the GSE1428 dataset. (F). Complex numerical heat map of differentially expressed genes associated with iron death for the GSE8479 dataset.

Co-DEGs: Common differentially expressed genes.

Functional enrichment analysis underlying FRDEGs

The GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analyses of FRDEGs were conducted first to investigate the biological processes, molecular functions, cellular components, biological pathways, as well as their relationship with sarcopenia in 46 FRDEGs. The results of GO functional enrichment analysis and KEGG enrichment analysis were presented in bubble chart (Fig. 4A-B) and bar chart (Fig. 4C-D). In addition, the correlation between 46 FRDEGs and GO/KEGG enrichment analysis results was presented as a ring network diagram (Fig. 4E-F).

The results suggest that 46 FRDEGs are significantly enriched in response to steroid hormone, response to glucocorticoid, response to corticosteroid, metabolic process in the BP category, and organelle outer membrane, ubiquitin protein ligase binding in cellular component autophagosome, secondary lysosome in the CC category, and ubiquitin-like protein ligase binding, single-stranded DNA binding, coenzyme binding in the MF category (Table 1). A total of 46 FRDEGs were remarkably enriched in Carbon metabolism, Ferroptosis, Citrate cycle (TCA cycle), Glyoxylate and dicarboxylate metabolism, Pyruvate metabolism, Cysteine and methionine metabolism, Amino Acids, HIF-1 Signaling Pathway, 2-oxocarboxylic Acid metabolism, Glycolysis/Gluconeogenesis, Propanoate metabolism, as well as other pathways (Table 2).

Table 1

GO enrichment analysis results of Ferroptosis-related differentially expressed genes
ONTOLOGY	ID	Description	GeneRatio	BgRatio	pvalue	p.adjust	qvalue	Count
BP	GO:0048545	response to steroid hormone	7/46	385/18670	4.01818E-05	0.006702331	0.004754146	7
BP	GO:0051384	response to glucocorticoid	6/46	146/18670	1.4933E-06	0.002280405	0.001617553	6
BP	GO:0031960	response to corticosteroid	6/46	162/18670	2.7343E-06	0.002280405	0.001617553	6
BP	GO:0019362	pyridine nucleotide metabolic process	6/46	189/18670	6.64981E-06	0.002653115	0.001881927	6
CC	GO:0031968	organelle outer membrane	4/46	201/19717	0.001223747	0.061913794	0.050806452	4
CC	GO:0019867	outer membrane	4/46	203/19717	0.001269255	0.061913794	0.050806452	4
CC	GO:0005776	autophagosome	3/46	93/19717	0.001331479	0.061913794	0.050806452	3
CC	GO:0005767	secondary lysosome	2/46	14/19717	0.000475991	0.061913794	0.050806452	2
MF	GO:0031625	ubiquitin protein ligase binding	5/45	290/17697	0.000815692	0.022256728	0.017295115	5
MF	GO:0044389	ubiquitin-like protein ligase binding	5/45	308/17697	0.001067763	0.025492839	0.019809812	5
MF	GO:0003697	single-stranded DNA binding	4/45	113/17697	0.000191887	0.0091626	0.007120014	4
MF	GO:0050662	coenzyme binding	4/45	291/17697	0.006283042	0.085718645	0.066609694	4
GO: Gene Ontology;BP༚ biological process, BP༛ CC༚cellular component༛MF༚molecular function.

Table 2

KEGG enrichment analysis results of Ferroptosis-related differentially expressed genes
ONTOLOGY	ID	Description	GeneRatio	BgRatio	pvalue	p.adjust	qvalue	Count
KEGG	hsa01200	Carbon metabolism	5/30	118/8076	6.50825E-05	0.003514458	0.002945842	5
KEGG	hsa04216	Ferroptosis	4/30	41/8076	1.42431E-05	0.001538251	0.001289373	4
KEGG	hsa00020	Citrate cycle (TCA cycle)	3/30	30/8076	0.000175531	0.004739324	0.003972533	3
KEGG	hsa00630	Glyoxylate and dicarboxylate metabolism	3/30	30/8076	0.000175531	0.004739324	0.003972533	3
KEGG	hsa00620	Pyruvate metabolism	3/30	39/8076	0.000386305	0.008344188	0.006994153	3
KEGG	hsa00270	Cysteine and methionine metabolism	3/30	50/8076	0.000805981	0.014507664	0.01216042	3
KEGG	hsa00010	Glycolysis / Gluconeogenesis	3/30	67/8076	0.001887912	0.029127783	0.0244151	3
KEGG	hsa01230	Biosynthesis of amino acids	3/30	75/8076	0.002608454	0.031301446	0.026237079	3
KEGG	hsa04066	HIF-1 signaling pathway	3/30	109/8076	0.007450151	0.073146932	0.061312244	3
KEGG	hsa01210	2-Oxocarboxylic acid metabolism	2/30	19/8076	0.002193384	0.029610689	0.024819876	2
KEGG	hsa00640	Propanoate metabolism	2/30	34/8076	0.006951451	0.073146932	0.061312244	2
KEGG: Kyoto Encyclopedia of Genes and Genomes.

Figure 4. Ferroptosis-related DEGs functional enrichment analysis (GO) and pathway enrichment (KEGG) analysis.

(A-B) GO functional enrichment analysis of Ferroptosis-related DEGs (A) and KEGG pathway enrichment analysis (B) results Bubble diagram presentation. (C-D). Results of GO functional enrichment analysis (C) and KEGG pathway enrichment analysis (D) of Ferroptosis-related DEGs are shown in bar charts. (E-F). Annular network diagram of GO functional enrichment analysis (E) and KEGG pathway enrichment analysis (F) results of Ferroptosis-related DEGs.

Ferroptosis-related DEGs: Ferroptosis-related genes differentially expressed genes; GO: Gene Ontology; BP: Biological process; CC: Cellular component; MF: Molecular function; KEGG: Kyoto Encyclopedia of Genes and Genomes

GSEA enrichment analysis of the sarcopenia datasets

To determine differentially enriched pathways between sarcopenia and non-sarcopenia patients, the GSEA (Gene Set Enrichment Analysis) of gene expression profiles was performed. The results showed that in dataset GSE1428 (Table 3), the molecular pathways were significantly associated with oxidative phosphorylation (Fig. 5E) endoplasmic reticulum stress response in coronavirus infection (Fig. 5C), oxidative phosphorylation (Fig. 5B), glycolysis and gluconeogenesis (Fig. 5D) and other pathways (Fig. 5A-E). In dataset GSE8479 (Table 4), they were significantly enriched in reactome chaperone mediated autophagy (Fig. 5H), apoptosis-related network due to altered notch3 in ovarian cancer (Fig. 5I), mechano-regulation and pathology of YAPTAZ on the basis of hippo and non-hippo mechanisms (Fig. 5J), glycolysis gluconeogenesis (Fig. 5G) and other pathways (Fig. 5F-J).

Table 3

GSEA of dataset GSE1428
ID	setSize	enrichmentScore	NES	pvalue	p.adjust	qvalues
NABA_CORE_MATRISOME	192	0.359224489	1.591219408	0.001858736	0.168392046	0.16549529
PID_TAP63_PATHWAY	50	0.559261037	2.013902319	0.001876173	0.168392046	0.16549529
WP_MRNA_PROCESSING	112	0.409850575	1.696861139	0.001876173	0.168392046	0.16549529
PID_CMYB_PATHWAY	80	0.425781838	1.659216402	0.001886792	0.168392046	0.16549529
WP_ADIPOGENESIS	125	0.423679599	1.770682554	0.001908397	0.168392046	0.16549529
PID_P73PATHWAY	74	0.492236343	1.886005776	0.001912046	0.168392046	0.16549529
PID_P53_DOWNSTREAM_PATHWAY	120	0.382632438	1.593149301	0.001919386	0.168392046	0.16549529
PID_MYC_REPRESS_PATHWAY	58	0.497125351	1.833777637	0.001949318	0.168392046	0.16549529
REACTOME_DISEASES_ASSOCIATED_ WITH_GLYCOSAMINOGLYCAN_METABOLISM	35	0.581885739	1.932458778	0.001949318	0.168392046	0.16549529
REACTOME_PYRUVATE_METABOLISM	28	-0.668264848	-2.112287759	0.001956947	0.168392046	0.16549529
KEGG_CITRATE_CYCLE_TCA_CYCLE	30	-0.623445311	-2.003716563	0.001996008	0.168392046	0.16549529
WP_OXIDATIVE_PHOSPHORYLATION	30	-0.535266811	-1.720316051	0.001996008	0.168392046	0.16549529
WP_ENDOPLASMIC_RETICULUM_STRESS_ RESPONSE_IN_CORONAVIRUS_INFECTION	34	-0.501903728	-1.664692773	0.00204918	0.168392046	0.16549529
KEGG_OXIDATIVE_PHOSPHORYLATION	82	-0.548571102	-2.184397656	0.002145923	0.168392046	0.16549529
WP_GLYCOLYSIS_AND_GLUCONEOGENESIS	45	-0.534061743	-1.902402409	0.002150538	0.168392046	0.16549529
GSEA: Gene Set Enrichment Analysis.

Table 4

GSEA of dataset GSE8479
ID	setSize	enrichmentScore	NES	pvalue	p.adjust	qvalues
KEGG_SYSTEMIC_LUPUS_ERYTHEMATOSUS	57	0.717866394	2.230238825	0.001865672	0.04807219	0.040628643
REACTOME_COMPLEMENT_CASCADE	56	0.656205379	2.032523115	0.001879699	0.04807219	0.040628643
WP_PRIMARY_FOCAL_SEGMENTAL_ GLOMERULOSCLEROSIS_FSGS	70	0.635980629	2.024884952	0.001879699	0.04807219	0.040628643
REACTOME_INTERLEUKIN_4_AND_ INTERLEUKIN_13_SIGNALING	108	0.594235253	2.016724784	0.001808318	0.04807219	0.040628643
REACTOME_SMOOTH_MUSCLE_CONTRACTION	34	0.72233364	2.006980347	0.001949318	0.04807219	0.040628643
WP_COMPLEMENT_ACTIVATION	22	0.772084612	1.975263157	0.001941748	0.04807219	0.040628643
NABA_ECM_GLYCOPROTEINS	175	0.546853662	1.962865961	0.00173913	0.04807219	0.040628643
REACTOME_CHAPERONE_MEDIATED_AUTOPHAGY	22	0.764534312	1.955946841	0.001941748	0.04807219	0.040628643
PID_NOTCH_PATHWAY	59	0.626456076	1.954493131	0.001872659	0.04807219	0.040628643
REACTOME_RESPONSE_TO_METAL_IONS	14	0.854295521	1.943158541	0.001992032	0.04807219	0.040628643
BIOCARTA_CARDIACEGF_PATHWAY	18	0.79285177	1.941356478	0.001934236	0.04807219	0.040628643
BIOCARTA_COMP_PATHWAY	19	0.782625403	1.937007068	0.001934236	0.04807219	0.040628643
WP_APOPTOSISRELATED_NETWORK_ DUE_TO_ALTERED_NOTCH3_IN_OVARIAN_CANCER	53	0.573101748	1.75137624	0.001890359	0.04807219	0.040628643
WP_MECHANOREGULATION_AND_PATHOLOGY_OF_ YAPTAZ_VIA_HIPPO_AND_NONHIPPO_MECHANISMS	46	0.539252571	1.600002309	0.005671078	0.078188934	0.066082079
KEGG_GLYCOLYSIS_GLUCONEOGENESIS	61	-0.470240257	-1.490154504	0.027542373	0.184983731	0.15634066
GSEA: Gene Set Enrichment Analysis.

Figure 5. GSEA enrichment analysis of sarcopenia datasets.

(A). Four main biological characteristics of GSEA enrichment analysis for GSE1428 dataset. (B-E). Genes significantly enriched in WP_OXIDATIVE_PHOSPHORYLATION in GSE1428 dataset (Fig. 5E). WP_ENDOPLASMIC_RETICULUM_STRESS_RESPONSE_IN_CORONAVIRUS_INFECTION (Fig. 5C), KEGG_OXIDATIVE_PHOSPHORYLATION (Fig. 5B), WP_GLYCOLYSIS_AND_GLUCONEOGENESIS (Fig. 5D) pathway. (F). GSEA analysis of GSE8479 dataset showed four main biological characteristics. (G-J). Genes in GSE8479 dataset were significantly enriched in REACTOME_CHAPERONE_MEDIATED_AUTOPHAGY (Fig. 5H). "Wp_sisrelated_network_due_to_altered_notch3_in_ovarian_cancer" (Fig. 5I). WP_MECHANOREGULATION_AND_PATHOLOGY_OF_YAPTAZ_VIA_HIPPO_AND_NONHIPPO_MECHANISMS (Fig. 5J), KEGG_GLYCOLYSIS_GLUCONEOGENESIS (Fig. 5G) and other pathways (Fig. 5F-J).

GSEA: Gene Set Enrichment Analysis.

Protein-protein interaction network

A total of 46 FRDEGs (DDIT4, FTL, FTH1, HMGB1, HSPB1, CISD1, TP63, CDKN1A, PLIN2, ZFP36, SREBF2, FZD7, BRD3, DECR1, CP, GOT1, BEX1, MPC1, ABCC5, USP11, CS, ATG5, MTCH1, SLC38A1, DLD, CIRBP, STOML2, LDHA, MDH2, ANP32B, FKBP3, LPIN1, PGK1, ANXA1, YWHAZ, ACTC1, LRPPRC, MAP4, HNRNPA3, MYL6B, RUFY1, WBP11, FOXO1, PRPF8, RPS27A, MAGED2) were analyzed for protein-protein interaction using STRING database. The minimum required interaction score > 0.400 (medium confidence 0.400) was selected as the standard, constructing 46 FRDEGs of PPI network, using Cytoscape software to visualize the protein interaction relationship. The result indicates that only 31 FRDEGs in the PPI network are associated with other genes. Subsequently, the PPI network diagram (Fig. 6B) was plotted according to the degree scores of the above 31 FRDEGs in the PPI network. Next, cytoHubba plugin was adopted to analyze five algorithms (including MCC, MNC, EPC, Degree and Closeness). As a result, 11 common genes of the top 15 genes were selected as the hub genes (mRNA) (Fig. 6B-C).

Figure 6. Protein-protein Interaction Network (PPI Network).

(A). PPI network of Ferroptosis-related DEGs under Degree algorithm. (B). Common gene of the first 15 Ferroptosis-related DEGs under five algorithms including MCC, MNC, EPC, Degree, and Closeness selected by Wayne diagram. (C). PPI Network of 11 Hub genes in Degree algorithm.

Ferroptosis-related DEGs: Ferroptosis-related genes differentially expressed genes; PPI network: protein-protein interaction; Degree: Degree Correlation; MNC: Maximum Neighborhood Component; MCC: Maximal Clique Centrality; EPC: Edge Percolated Component; And you may say you may not so-called Centrality. The higher the Degree score, the larger the node area, the darker the color, and the more interactions there are.

Hub genes-miRNA interaction network; hub genes-transcription factor interaction network; hub genes-drug interaction network

The mRNA-miRNA data in ENCORI database and miRDB database were used to predict the correlation between 11 hub genes (mRNA) (LDHA, CS, MDH2, PGK1, DLD, YWHAZ, RPS27A, FOXO1, PLIN2, CDKN1A, HSPB1), and then the intersection part of the results from the two databases was selected to draw the mRNA-miRNA interaction network by Cytoscape software for visualization (FIG. 7A). According to the mRNA-miRNA interaction network, our mRNA-miRNA interaction network consists of 11 hub genes (mRNA) (LDHA, CS, MDH2, PGK1, DLD, YWHAZ, RPS27A, FOXO1, PLIN2, CDKN1A, HSPB1), including 189 miRNA molecules, composing 428 mRNA-miRNA interaction pairs.

We searched the CHIPBase database (Version 2.0) and hTFtarget database for transcription factors (TFS) that is binding to the hub genes. After downloading the interaction relationship found in the two databases, we took the intersection with 11 hub genes (mRNA), and finally obtained 11 hub genes (CCNB1, CEP55, KIF20A, MELK, MKI67, NCAPG, NUF2, RRM2, TTK, UBE2C, UHRF1) and 48 TFS were visualized by Cytoscape software (FIG.7B). In the mRNA-TF interaction network, hub gene YWHAZ has the most interaction relationships with TF, and YWHAZ has 40 pairs of mRNA-TF interaction relationships.

The DGIdb database was used to identify potential drugs or molecular compounds for 11 hub genes. As shown in the mRNA-Drugs interaction network (FIG.7C), we identified 56 potential drugs or molecular compounds corresponding to six hub genes (LDHA, MDH2, PGK1, FOXO1, CDKN1A, HSPB1). Among them, the drugs or molecular compounds that interact with hub gene HSPB1 are the most.

Figure 7. Hub Genes-miRNA interaction network; Hub genes-transcription factor interaction network; Hub Genes-Drug interaction network.

(A). Hub Genes-miRNA interaction network. (B). Hub genes-transcription factor interaction network. (C). Hub genes-drug interaction network. The green oval is mRNA; Yellow hexagon is miRNA; The pink quadrilateral represents transcription factors; The blue triangle is Drug.

PPI network: protein-protein interaction network; TF: Transcription factors.

Expression analysis and ROC verification of hub genes in two datasets

To further explore the differential expression of FRGs in the sarcopenia datasets, the correlation between expression levels and different groups of 11 hub genes (LDHA, CS, MDH2, PGK1, DLD, YWHAZ, RPS27A, FOXO1, PLIN2, CDKN1A, HSPB1) screened by PPI network in GSE1428 dataset (Fig. 8A) and GSE8479 dataset (Fig. 9A) was analyzed. Thus, ROC verification was performed. As depicted in Fig. 8, 11 hub genes were statistically significant among different groups in the GSE1428 dataset (Fig. 8A). Thereinto, the expression levels of hub genes FOXO1, PLIN2, RPS27A and YWHAZ among different groups were statistically significant (P < 0.05). The expression levels of hub genes CDKN1A, CS, DLD, HSPB1, LDHA, MDH2 and PGK1 among different groups were highly statistically significant (P < 0.01). Figure 8B-K demonstrated the ROC verification of hub genes in GSE1428 dataset. Hub genes CDKN1A (AUC = 0.892, Fig. 8B), CS (AUC = 0.883, Fig. 8C), DLD (AUC = 0.846, Fig. 8D), FOXO1 (AUC = 0.783, Fig. 8E), HSPB1 (AUC = 0.867, Fig. 8F), LDHA (AUC = 0.842, Fig. 8G), MDH2 (AUC = 0.854, Fig. 8H), PGK1 (AUC = 0.867, Fig. 8I), RPS27A (AUC = 0.812, Fig. 8J), YWHAZ (AUC = 0.783, Fig. 8K). The expression level in the GSE1428 dataset showed a certain correlation with the Normal sample group (Normal) and Sarcopenia sample group (Sarcopenia).

In the GSE8479 dataset (Fig. 9A), the expression levels of 10 hub genes (LDHA, CS, MDH2, PGK1, DLD, YWHAZ, RPS27A, FOXO1, CDKN1A, HSPB1) were statistically significant among different groups. Thereinto, the expression levels of hub genes LDHA, CS, MDH2, PGK1, DLD, YWHAZ, FOXO1 and CDKN1A in different groups were highly statistically significant (P < 0.001). The expression levels of HSPB1 and RPS27A in different groups were statistically significant (P < 0.05). Figure 9B-K shows the ROC verification of hub genes in the GSE8479 dataset. As depicted in the figure, the expression level of hub genes CDKN1A (AUC = 0.954, Fig. 9B), DLD (AUC = 0.915, Fig. 9D) and FOXO1 (AUC = 0.928, Fig. 9E) in the GSE8479 dataset are significantly correlated with the Normal sample group (Normal) and Sarcopenia sample group (Sarcopenia). The expression level of hub genes CS (AUC = 0.875, Fig. 9C), HSPB1 (AUC = 0.704, Fig. 9F), LDHA (AUC = 0.858, Fig. 9G), MDH2 (AUC = 0.859, Fig. 9H), PGK1 (AUC = 0.856, Fig. 9I) and YWHAZ (AUC = 0.852, Fig. 9K) in the GSE8479 dataset showed a certain correlation with Normal sample group (Normal) and Sarcopenia sample group (Sarcopenia). The expression level of hub gene RPS27A (AUC = 0.690, Fig. 9J) in the GSE8479 dataset showed low correlation with both Normal sample group (Normal) and Sarcopenia sample group (Sarcopenia).

Figure 8. Differential analysis of Hub genes expression in dataset GSE1428.

(A). Group comparison of the results of differential analysis of hub genes expression in GSE1428 dataset. (B-K). ROC curves of hub genes CDKN1A (B), CS (C), DLD (D), FOXO1 (E), HSPB1 (F), LDHA (G), MDH2 (H), PGK1 (I), RPS27A (J) and YWHAZ (K) in GSE1428 dataset.

AUC: Area Under Curve. The symbol NS is equal to P ≥ 0.05, which is not statistically significant. The symbol * is equivalent to P < 0.05, which is statistically significant; The symbol ** is equivalent to P < 0.01, which is highly statistically significant; The symbol *** is equivalent to P < 0.001, which is highly statistically significant.

Figure 9. Differential analysis of Hub genes expression in dataset GSE8479.

(A). Group comparison diagram of the results of differential analysis of hub genes expression in GSE8479 data set. (B-K). ROC curves of hub genes CDKN1A (B), CS (C), DLD (D), FOXO1 (E), HSPB1 (F), LDHA (G), MDH2 (H), PGK1 (I), RPS27A (J) and YWHAZ (K) in GSE8479 dataset.

Immune infiltration of hub genes

After disposing the expression profile data of sarcopenia dataset GSE1428 and GSE8479, ssGSEA algorithm was applied to count the infiltration of 29 kinds of immune cells. We plotted the infiltration of immune cells in the two datasets in the form of bar chart (Fig. 10A-B).

Afterwards, we calculated the correlation between the expression level of 11 hub genes (LDHA, CS, MDH2, PGK1, DLD, YWHAZ, RPS27A, FOXO1, PLIN2, CDKN1A) and the abundance of immune cell infiltration in dataset GSE1428 (Fig. 11A) and GSE8479 (Fig. 11I). For the combination of immune cells and hub genes with statistically significant positive and negative correlations, scatter plots of correlation between the expression level of hub genes and the abundance of immune cell infiltration in datasets GSE1482 and GSE8479 were drawn. As depicted in the figure, a moderate negative correlation was identified between the infiltration degree of DCs of immune cells and the expression level of hub gene DLD (Fig. 11B, r= -0.668, P < 0.001) in dataset GSE1482. A weak negative correlation was identified between the degree of DCs infiltration and the expression level of hub gene LDHA (Fig. 11C, r= -0.443, P-value = 0.039). A weak negative correlation was identified between the invasion degree of Para inflammation and the expression level of hub gene LDHA (Fig. 11D, r= -0.485, P-value = 0.022). A weak positive correlation was identified between the degree of Treg infiltration of immune cells and the expression level of hub gene LDHA (Fig. 11E, r = 0.451, P-value = 0.035), and a moderate negative correlation between the degree of DCs infiltration of immune cells and the expression level of hub gene PGK1 (Fig. 11F, r = -0.514, P-value = 0.015), a moderate negative correlation was identified between the infiltration degree of immune cells pDCs and the expression level of hub gene YWHAZ (Fig. 11G, r= -0.660, P-value = 0.001). A weak negative correlation was identified between the infiltration degree of immune cells Th1 cells and the expression level of hub gene YWHAZ (Fig. 11H, r= -0.459, P-value = 0.033).

In GSE8479 dataset, a moderate negative correlation was identified between the degree of DCs infiltration and the expression level of hub gene DLD (Fig. 11J, r= -0.622, P-value < 0.001). A weak negative correlation was identified between the degree of DCs infiltration and the expression level of hub gene LDHA (Fig. 11K, r= -0.469, P-value < 0.001). A weak negative correlation was identified between the invasion degree of Para inflammation and the expression level of hub gene LDHA (Fig. 11L, r= -0.318, P-value = 0.023). A weak positive correlation was identified between the degree of Treg infiltration of immune cells and the expression level of hub gene LDHA (Fig. 11M, r = 0.307, P-value = 0.028), and a weak negative correlation was identified between the degree of DCs infiltration of immune cells and the expression level of hub gene PGK1 (Fig. 11N, r = -0.405, P-value = 0.003). a weak negative correlation was identified between the infiltration degree of immune cells pDCs and the expression level of hub gene YWHAZ (Fig. 11O, r= -0.390, P-value = 0.005). A weak negative correlation was identified between the infiltration degree of immune cells Th1 cells and the expression level of hub gene YWHAZ (Fig. 11P, r= -0.379, P-value = 0.006).

Figure 10. Immune infiltration analysis in sarcopenia dataset GSE1428 and GSE8479.

(A). Histogram of immune infiltration results of immune cells in dataset GSE1428. (B). Histogram of immune infiltration results of immune cells in dataset GSE8479.

Figure 11. Correlation analysis between the expression of hub genes and the degree of immune cell infiltration.

(A). Results of correlation analysis between hub genes and the abundance of 29 types of immune cell infiltration in dataset GSE1428. (B-H). Results of correlation analysis between the infiltration degree of immune cell DCs and the expression level of hub gene DLD in dataset GSE1428 (B), results of correlation analysis between the infiltration degree of immune cell DCs and the expression level of Hub gene LDHA (C), results of correlation analysis on the degree of infiltration of immune cells Para inflammation and the expression level of hub gene LDHA (D), and results of correlation analysis on the degree of infiltration of immune cells Treg and the expression level of hub gene LDHA (E). Results of correlation analysis between the infiltration degree of immune cells DCs and the expression level of hub gene PGK1 (F). Results of correlation analysis between the infiltration degree of immune cells pDCs and the expression water of hub gene YWHAZ (G). Correlation analysis between the infiltration degree of immune cells Th1 cells and the expression level of hub gene YWHAZ (H). I. Results of correlation analysis between hub genes and the abundance of 29 types of immune cell infiltration in dataset GSE8479. (J-P). Results of correlation analysis between the infiltration degree of immune cell DCs and the expression level of hub gene DLD in dataset GSE8479 (J), and results of correlation analysis between the infiltration degree of immune cell DCs and the expression level of Hub gene LDHA (K). Results of correlation analysis on the degree of infiltration of immune cells Para inflammation and the expression level of hub gene LDHA (L), and results of correlation analysis on the degree of infiltration of immune cells Treg and the expression level of hub gene LDHA (M). Results of correlation analysis between the infiltration degree of DCs and the expression level of hub gene PGK1 (N). Results of correlation analysis between the infiltration degree of immune cells pDCs and the expression water of hub gene YWHAZ (O). Correlation analysis between the infiltration degree of immune cells Th1 and the expression level of hub gene YWHAZ (P).

The symbol NS is equal to P ≥ 0.05, which is not statistically significant. The symbol * is equivalent to P < 0.05, which is statistically significant; The symbol ** is equivalent to P < 0.01, which is highly statistically significant; The symbol *** is equivalent to P < 0.001, which is highly statistically significant. The correlation coefficient R is positive, indicating that the two variables may be positively correlated. The correlation coefficient is negative, indicating that there may be a negative correlation between the two variables. The absolute value of correlation coefficient above 0.8 is a strong correlation, and the absolute value between 0.5 and 0.8 is a moderate correlation. The absolute value between 0.3 and 0.5 indicates a general degree of correlation; An absolute value below 0.3 is weak or uncorrelated.

Sarcopenia refers to a generalized and systematic disease characterized by age-related degenerative amyotrophy accompanied by increasing risk of adverse outcomes. The development of sarcopenia is correlated with several biological processes (e.g., mitochondrial dysfunction, oxidative stress, and inflammation ³²). Genetic factors serve as a key point in the progression of sarcopenia ³³. One of the distinctive features of ferroptosis is that it’s an iron-dependent pattern of non-apoptotic cell death ³⁴. Recent research has suggested that overload of iron in muscle cells, especially as age increasing, can cause amyotrophy or dysfunction, providing evidences that sarcopenia has connection with iron accumulation^35–37. For instance, it’s been found iron overload can up-regulate the expression of P53, which whereafter repress the protein level of Slc7a11 (solute carrier family 7, member 11), a renowned FRG. The down-regulation of Slc7a11 then lead to the ferroptosis of muscle cells by accumulating lipid peroxidation products, which may be one of the pathogenesis of sarcopenia⁷. However, there is a lack of sensitive and specific biomarkers to identify the sarcopenia, and there are no potential pharmacological methods of sarcopenia in clinical application. Accordingly, the pathological and molecular mechanisms of sarcopenia for diagnosis and therapy in clinical should be comprehended. Considerable genes concerned with the development of sarcopenia have been screened (e.g., Estrogen-related receptor alpha (ERRα), Mitofusin 2 (MFN2), as well as Optic Atrophy 1 (Opa1)³⁸. ) However, brand-new classes biomarkers with great sensitivity, specificity, and efficiency are demanded for the diagnosis and prognosis of sarcopenia. In this study, an integrated analysis was conducted on two mRNA microarray datasets (GSE1428 and GSE8479). In general, 46 DEGs, playing probable roles in the pathologic processes related to ferroptosis in sarcopenia (FRDEGs), were identified.

In this study, 11 hub genes (including LDHA, CS, MDH2, PGK1, DLD, YWHAZ, RPS27A, FOXO1, PLIN2, CDKN1A, and HSPB1) were screened by the PPI network with the highest scores, which can serve as specific biomarkers of sarcopenia diagnosis. The consequences of functional enrichment suggested that the alterations of the above genes were primarily enriched in the response to steroid hormone, response to glucocorticoid, response to corticosteroid and pyridine nucleotide metabolic process. The pathways related to sarcopenia involve the carbon metabolism, ferroptosis, citrate cycle, glyoxylate and dicarboxylate metabolism, pyruvate metabolism, cysteine and methionine metabolism, glycolysis / gluconeogenesis, biosynthesis of amino acids, hif-1 signaling pathway, 2-oxocarboxylic acid metabolism and propanoate metabolism. Accordingly, the outcome posed the assumption that the modules and pathways mostly enriched have genetic effect on sarcopenia. Furthermore, the above hub genes predicted perhaps interact in a network and coregulate sarcopenia.

Among the above biological processes, steroid hormone is essential in skeletal muscle function which can limit inflammatory stress damage on skeletal muscle ³⁹. The imbalance of corticosteroid including glucocorticoid level can cause frailty, whose key convergent pathological effects include loss of muscle mass and strength ⁴⁰. Niacin is capable of curing systemic NAD deficiency and contributing to a better physical performance in mitochondrial myopathy ⁴¹. For pathways, Carbon metabolism includes biosynthesis (purines and thymidine), amino acid homeostasis (including glycine, serine, and methionine), epigenetic maintenance, and redox defense⁴². Amino acid hemostasis takes on a critical significance to protein synthesis in the body. In such a condition, an active mTORC1 complex plays a key role to sustain the balance. As amino acid close to exhuasted, mTORC1 becomes inactivated, which triggers autophagy, leading to degradation of cellular proteins as well as protein biosynthesis reduction⁴³. Oxidative stress is a common and pivotal process that involves in the evolution of age-related diseases by inducing structural and functional changes of various cellular components. A decrease of the redox defense is an underlying cause of sarcopenia⁴⁴.

The ROC results suggest that CDKN1A, CS, DLD, FOXO1, HSPB1, LDHA, MDH2, and YWHAZ have relatively high sensitivity and specificity in both datasets of sarcopenia. To be specific, CDKN1A is a usually used senescence and denervation markers (SDMs) related to myofiber atrophy⁴⁵. FOXO1 regulates amyotrophy-related procedures, containing metabolic changes and E3-ubiquitin ligase-mediated proteolysis ⁴⁶. HSPB1 serves as a main indicator of oxidative stress in muscle physiology⁴⁷, and its phosphorylated form increases with the development of age-related amyotrophy⁴⁸. Moreover, HSPB1 can down-regulate NF-κB in skeletal muscle, thus attenuating the amyotrophy afterwards⁴⁹. MDH2 is a vital prognostic biomarker in Duchenne muscular dystrophy⁵⁰. Based on the above hub genes, the mRNA-microRNA, mRNA-TF and mRNA-drugs network were also constructed, providing novel targets for diagnosis and treatment of sarcopenia.

Besides, an immune infiltration analysis was conducted on the FRDEGs of sarcopenia. The results suggest that the correlation between the extent of DCs infiltration of immune cells and the expression level of hub gene DLD and PGK1 is medium and negative. In early myositis, DCs may induce muscle injury through intercellular intercourse and inflammatory cytokine secretion, thus resulting in impaired muscle regeneration⁵¹. YWHAZ, a central hub protein for numerous signal transduction pathways and a significant role in tumor progression⁵², shows a medium negative correlation with pDCs. Together the above findings indicated that the FRDEGs DLD, PGK1 and YWHAZ contribute to recruiting and regulating immune infiltrating cells Dcs and pDCs in sarcopenia.

The correlation between sarcopenia and ferroptosis were discovered using several methods, which provides more insights. However, there were still certain limitations. First, there are no sufficient sample size relatively, the statistical power may be low. Microarray samples including patients with diverse grades of sarcopenia are required to comprehensively determine the latent molecular mechanisms of sarcopenia during occurrence and development. One major limitation was that there was no basic experiment in vitro or in vivo to verify and approve our findings. Another one was that there was a lack of sarcopenia-related clinical data (e.g., stage and conditions) to perform further exploration of the correlation between ferroptosis and sarcopenia. Furthermore, the miRNA, TF, and drug target we predicted should be further verified in practice.

In brief, ferroptosis, a phenotype of sarcopenia during its occurrence and progression, does play a significant role. The hub genes (including LDHA, CS, MDH2, PGK1, DLD, YWHAZ, RPS27A, FOXO1, PLIN2, CDKN1A, and HSPB1) exhibit relatively high sensitivity and specificity in the diagnosis of sarcopenia. The results also suggest that DLD, PGK1 and YWHAZ inactivate dendritic cells (DCs and pDCs). Moreover, the TF, miRNA, and drugs predicted in this study may provide evidence for pathological mechanisms and therapeutic targets. The above results somehow provide novel insights into the diagnostic and therapeutic tactics of sarcopenia.

Statements & Declarations

Funding

This study was supported by National Key R&D Program of China (2018YFC2000604).

Competing Interests

The authors have no relevant financial or non-financial interests to disclose.

Author Contributions

Yanzhong Chen conceptualized and designed the study. Yanzhong Chen and Yaonan Zhang collected and validated the data. Yanzhong Chen analyzed the data. Yanzhong Chen, Sihan Zhang and Hong Ren prepared the manuscript. Hong Ren and Yaonan Zhang revised the manuscript critically. Sihan Zhang searched the relevant references. All authors contributed to the article and approved the submitted version.

Ethics approval

There are no ethical issues involved in this study.

Data Availability Statement

The datasets generated during the current study are available in the NCBI (GEO) repository.

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE8479.

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE1428.

Cruz-Jentoft, A. J. & Sayer, A. A. Sarcopenia. Lancet 393, 2636–2646, doi:10.1016/s0140-6736(19)31138-9 (2019).
Dhillon, R. J. & Hasni, S. Pathogenesis and Management of Sarcopenia. Clin Geriatr Med 33, 17–26, doi:10.1016/j.cger.2016.08.002 (2017).
Cesari, M. et al. Biomarkers of sarcopenia in clinical trials-recommendations from the International Working Group on Sarcopenia. J Cachexia Sarcopenia Muscle 3, 181–190, doi:10.1007/s13539-012-0078-2 (2012).
Abellan van Kan, G. et al. Sarcopenia: biomarkers and imaging (International Conference on Sarcopenia research). J Nutr Health Aging 15, 834–846, doi:10.1007/s12603-011-0365-1 (2011).
Li, J. et al. Ferroptosis: past, present and future. Cell Death Dis 11, 88, doi:10.1038/s41419-020-2298-2 (2020).
Qiu, Y., Cao, Y., Cao, W., Jia, Y. & Lu, N. The Application of Ferroptosis in Diseases. Pharmacol Res 159, 104919, doi:10.1016/j.phrs.2020.104919 (2020).
Huang, Y. et al. Ferroptosis in a sarcopenia model of senescence accelerated mouse prone 8 (SAMP8). Int J Biol Sci 17, 151–162, doi:10.7150/ijbs.53126 (2021).
Chen, Y. Y. et al. Association Between Interleukin-12 and Sarcopenia. J Inflamm Res 14, 2019–2029, doi:10.2147/jir.S313085 (2021).
Li, H., Guan, K., Liu, D. & Liu, M. Identification of mitochondria-related hub genes in sarcopenia and functional regulation of MFG-E8 on ROS-mediated mitochondrial dysfunction and cell cycle arrest. Food Funct 13, 624–638, doi:10.1039/d1fo02610k (2022).
Huang, S., Xiang, C. & Song, Y. Identification of the shared gene signatures and pathways between sarcopenia and type 2 diabetes mellitus. PLoS One 17, e0265221, doi:10.1371/journal.pone.0265221 (2022).
Yang, H., Tian, W. & Zhou, B. Sarcopenia and a 5-mRNA risk module as a combined factor to predict prognosis for patients with stomach adenocarcinoma. Genomics 114, 361–377, doi:10.1016/j.ygeno.2021.12.011 (2022).
Pazit, L. et al. Safety and feasibility of high speed resistance training with and without balance exercises for knee osteoarthritis: A pilot randomised controlled trial. Phys Ther Sport 34, 154–163, doi:10.1016/j.ptsp.2018.10.001 (2018).
Giresi, P. G. et al. Identification of a molecular signature of sarcopenia. Physiol Genomics 21, 253–263, doi:10.1152/physiolgenomics.00249.2004 (2005).
Melov, S., Tarnopolsky, M. A., Beckman, K., Felkey, K. & Hubbard, A. Resistance exercise reverses aging in human skeletal muscle. PLoS One 2, e465, doi:10.1371/journal.pone.0000465 (2007).
Davis, S. & Meltzer, P. S. GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics 23, 1846–1847, doi:10.1093/bioinformatics/btm254 (2007).
Ritchie, M. E. et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 43, e47, doi:10.1093/nar/gkv007 (2015).
Stelzer, G. et al. The GeneCards Suite: From Gene Data Mining to Disease Genome Sequence Analyses. Curr Protoc Bioinformatics 54, 1.30.31–31.30.33, doi:10.1002/cpbi.5 (2016).
Zhou, N. & Bao, J. FerrDb: a manually curated resource for regulators and markers of ferroptosis and ferroptosis-disease associations. Database (Oxford) 2020, doi:10.1093/database/baaa021 (2020).
Yu, G. Gene Ontology Semantic Similarity Analysis Using GOSemSim. Methods Mol Biol 2117, 207–215, doi:10.1007/978-1-0716-0301-7_11 (2020).
Ogata, H. et al. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res 27, 29–34, doi:10.1093/nar/27.1.29 (1999).
Yu, G., Wang, L. G., Han, Y. & He, Q. Y. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics 16, 284–287, doi:10.1089/omi.2011.0118 (2012).
Subramanian, A. et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 102, 15545–15550, doi:10.1073/pnas.0506580102 (2005).
Liberzon, A. et al. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst 1, 417–425, doi:10.1016/j.cels.2015.12.004 (2015).
Szklarczyk, D. et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res 47, D607-d613, doi:10.1093/nar/gky1131 (2019).
Chin, C. H. et al. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol 8 Suppl 4, S11, doi:10.1186/1752-0509-8-s4-s11 (2014).
Li, J. H., Liu, S., Zhou, H., Qu, L. H. & Yang, J. H. starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res 42, D92-97, doi:10.1093/nar/gkt1248 (2014).
Chen, Y. & Wang, X. miRDB: an online database for prediction of functional microRNA targets. Nucleic Acids Res 48, D127-d131, doi:10.1093/nar/gkz757 (2020).
Zhou, K. R. et al. ChIPBase v2.0: decoding transcriptional regulatory networks of non-coding RNAs and protein-coding genes from ChIP-seq data. Nucleic Acids Res 45, D43-d50, doi:10.1093/nar/gkw965 (2017).
Zhang, Q. et al. hTFtarget: A Comprehensive Database for Regulations of Human Transcription Factors and Their Targets. Genomics Proteomics Bioinformatics 18, 120–128, doi:10.1016/j.gpb.2019.09.006 (2020).
Freshour, S. L. et al. Integration of the Drug-Gene Interaction Database (DGIdb 4.0) with open crowdsource efforts. Nucleic Acids Res 49, D1144-d1151, doi:10.1093/nar/gkaa1084 (2021).
Newman, A. M. et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 12, 453–457, doi:10.1038/nmeth.3337 (2015).
Marzetti, E. et al. Sarcopenia: an overview. Aging Clin Exp Res 29, 11–17, doi:10.1007/s40520-016-0704-5 (2017).
Pascual-Fernández, J. et al. Sarcopenia: Molecular Pathways and Potential Targets for Intervention. Int J Mol Sci 21, doi:10.3390/ijms21228844 (2020).
Dixon, S. J. et al. Ferroptosis: an iron-dependent form of nonapoptotic cell death. Cell 149, 1060–1072, doi:10.1016/j.cell.2012.03.042 (2012).
Altun, M. et al. Iron load and redox stress in skeletal muscle of aged rats. Muscle Nerve 36, 223–233, doi:10.1002/mus.20808 (2007).
DeRuisseau, K. C. et al. Aging-related changes in the iron status of skeletal muscle. Exp Gerontol 48, 1294–1302, doi:10.1016/j.exger.2013.08.011 (2013).
Aydemir, T. B. et al. Aging amplifies multiple phenotypic defects in mice with zinc transporter Zip14 (Slc39a14) deletion. Exp Gerontol 85, 88–94, doi:10.1016/j.exger.2016.09.013 (2016).
Migliavacca, E. et al. Mitochondrial oxidative capacity and NAD(+) biosynthesis are reduced in human sarcopenia across ethnicities. Nat Commun 10, 5808, doi:10.1038/s41467-019-13694-1 (2019).
Laurent, M. R. et al. Age-related bone loss and sarcopenia in men. Maturitas 122, 51–56, doi:10.1016/j.maturitas.2019.01.006 (2019).
Clegg, A. & Hassan-Smith, Z. Frailty and the endocrine system. Lancet Diabetes Endocrinol 6, 743–752, doi:10.1016/s2213-8587(18)30110-4 (2018).
Pirinen, E. et al. Niacin Cures Systemic NAD(+) Deficiency and Improves Muscle Performance in Adult-Onset Mitochondrial Myopathy. Cell Metab 31, 1078–1090.e1075, doi:10.1016/j.cmet.2020.04.008 (2020).
Ducker, G. S. & Rabinowitz, J. D. One-Carbon Metabolism in Health and Disease. Cell Metab 25, 27–42, doi:10.1016/j.cmet.2016.08.009 (2017).
Bröer, S. & Bröer, A. Amino acid homeostasis and signalling in mammalian cells and organisms. Biochem J 474, 1935–1963, doi:10.1042/bcj20160822 (2017).
Fougere, B., van Kan, G. A., Vellas, B. & Cesari, M. Redox Systems, Antioxidants and Sarcopenia. Curr Protein Pept Sci 19, 643–648, doi:10.2174/1389203718666170317120040 (2018).
Solovyeva, E. M. et al. New insights into molecular changes in skeletal muscle aging and disease: Differential alternative splicing and senescence. Mech Ageing Dev 197, 111510, doi:10.1016/j.mad.2021.111510 (2021).
Liu, L. et al. Identification of a KLF5-dependent program and drug development for skeletal muscle atrophy. Proc Natl Acad Sci U S A 118, doi:10.1073/pnas.2102895118 (2021).
Poussard, S., Pires-Alves, A., Diallo, R., Dupuy, J. W. & Dargelos, E. A natural antioxidant pine bark extract, Oligopin®, regulates the stress chaperone HSPB1 in human skeletal muscle cells: a proteomics approach. Phytother res 27, 1529–1535 (2013).
Yamauchi, J. et al. (-)-epigallocatechin gallate inhibits prostaglandin D2-stimulated HSP27 induction via suppression of the p44/p42 MAP kinase pathway in osteoblasts. Prostaglandins Leukot Essent Fatty Acids 77, 173–179, doi:10.1016/j.plefa.2007.09.001 (2007).
Dodd, S. L., Hain, B., Senf, S. M. & Judge, A. R. Hsp27 inhibits IKKbeta-induced NF-kappaB activity and skeletal muscle atrophy. Faseb j 23, 3415–3423, doi:10.1096/fj.08-124602 (2009).
Signorelli, M. et al. Longitudinal serum biomarker screening identifies malate dehydrogenase 2 as candidate prognostic biomarker for Duchenne muscular dystrophy. J Cachexia Sarcopenia Muscle 11, 505–517, doi:10.1002/jcsm.12517 (2020).
Ladislau, L. et al. Activated dendritic cells modulate proliferation and differentiation of human myoblasts. Cell Death Dis 9, 551, doi:10.1038/s41419-018-0426-z (2018).
Gan, Y., Ye, F. & He, X. X. The role of YWHAZ in cancer: A maze of opportunities and challenges. J Cancer 11, 2252–2264, doi:10.7150/jca.41316 (2020).

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

Identification and validation of ferroptosis-Related genes in sarcopenia

Status:

Version 1

Abstract

Figures

Introduction

Materials and methods

Data collection and preprocessing

Identification of the differentially expressed genes

Function enrichment analysis

GSEA

Protein-protein interaction (PPI) network

mRNA-miRNA, mRNA-TF, mRNA-drugs interaction network

The ROC Curve

Immune infiltration analysis (ssGSEA)

Statistical analysis

Results

Data normalization and FRDEGs analysis

Functional enrichment analysis underlying FRDEGs

GSEA enrichment analysis of the sarcopenia datasets

Protein-protein interaction network

Hub genes-miRNA interaction network; hub genes-transcription factor interaction network; hub genes-drug interaction network

Expression analysis and ROC verification of hub genes in two datasets

Immune infiltration of hub genes

Discussion

Conclusions

Declarations

References

Additional Declarations

Status:

Version 1