A Comprehensive Bioinformatics Analysis of the Molecular Mechanisms Common to Mitophagy and Atherosclerosis

1. Introduction

Atherosclerosis (AS) is a chronic inflammatory disease of the blood vessels. Cell death, vascular endothelial damage, and vascular inflammatory infiltration are the key pathological mechanisms of AS [1]. Plaque rupture, thrombosis, and stenotic blood flow obstruction in vascular lumens may lead to coronary artery disease, ischemic cerebral infarction, and peripheral artery disease. Hence, AS is a major cause of high mortality and disability worldwide [2]. AS is most commonly treated with statins which reduce low-density lipoprotein (LDL) levels but not the incidence of major adverse cardiovascular events (MACEs) [3,4]. Thus, the complex mechanisms of AS onset and progression must be elucidated to develop efficacious prophylactic and therapeutic measures for it.

Mitophagy is selective mitochondrial degradation. It is induced by various factors that are deleterious to cells and clears damaged and dysfunctional mitochondria via a unique molecular mechanism [5]. Mitophagy regulates cardiovascular activity through different pathways and is implicated in the pathophysiology of various cardiovascular diseases (CVDs) such as hypertension, arteriosclerosis, myocardial ischemia-reperfusion injury, and heart failure (HF) [6,7,8]. Mitochondrial dysfunction and AS are associated mainly with oxidative stress, abnormal ATP biosynthesis, extensive cell injury, and apoptosis [9,10]. Unstable atherosclerotic plaques have a typical autophagic ultrastructure. Nevertheless, the mechanisms connecting mitophagy and AS are unknown. Therefore, it is necessary to assess the feasibility of mitophagy-related genes (MRGs) as diagnostic and therapeutic biomarkers of AS [11].

Single-cell RNA sequencing (scRNA-Seq) techniques reveal similarities and differences between and among individual cells [12]. ScRNA-Seq is now widely used in the field of CVD research as it generates detailed transcriptional profiles of the various cell types related to AS [13,14]. ScRNA-Seq detects and identifies rare and easily overlooked cell populations and helps clarify the molecular mechanisms of various diseases [13]. Therefore, the present study used single-cell datasets to investigate the cellular-level expression of AS biomarkers.

Mendelian randomization (MR) uses genetic variants as instrumental variables (IVs) to analyze the causality between exposure factors and disease outcomes [12]. As genotypes are determined at conception, they are at a low risk of confounding and reverse causality. Hence, causality between exposure factors and disease outcomes can be effectively analyzed [15]. Prior studies showed that the genetic predisposition to arteriosclerosis may be evaluated through causal reasoning [16,17,18]. Here, we applied MR analysis to (a) screen genes associated with AS causality and (b) facilitate precision research.

To the best of our knowledge, variations in the MRGs related to AS have seldom been reported. For this reason, the present work aimed to explore the cause/effect association between specific mitophagy-related genes and AS onset and progression. Bioinformatics and MR analyses were implemented to detect biomarkers associated with mitophagy in AS. Single-cell analysis identified the key cell types connected to AS and measured their AS biomarker expression levels. It is believed that the discoveries made herein will facilitate clinical AS prophylaxis and therapy.

2. Results

2.1. A total of 4,001 module genes and 102 DEGs1 were obtained in AS

The Wilcoxon signed-rank test revealed that the MRG ssGSEA score was significantly lower for the AS than the normal tissue group (Figure 1a). The clustering analysis indicated that there was no obvious outlier in the GSE43292 dataset (Figure S1). Seven was selected as the soft threshold. The gene interactions were consistent with the scale-free network distribution (Figure 1b). The black (correlation coefficient (COR) = 0.56; p < 0.05) and blue (COR = -0.51; p < 0.05) modules were the most strongly positively and negatively correlated, respectively, with the MRG ssGSEA score (Figure 1c,d) and accounted for 4,001 module genes. The differential expression analysis generated 102 DEGs1 of which 49 were upregulated and 53 were downregulated in AS (Figure 1e,f).

2.2. ScRNA-Seq analysis identified 818 DEGs2

Before QC, 1,000 genes were detected in the cells, the total expression level for all genes was < 25,000, the proportion of mitochondrial genes was < 20%, and the proportion of erythrocytic genes was ~0 (Figure S2). After QC, we obtained 51,721 cells and 22,615 genes (Figure S3). The top 2,000 highly variable genes were determined and the top 10 of these were displayed (Figure 2a). The cells were then grouped into 19 taxa based on an optimal latitude value of 20 (Figure 2b-c). The annotation revealed eight distinct cell types including B cells, endothelial cells, monocytes, neutrophils, NK cells, smooth muscle cells, T cells, and tissue stem cells. However, T cells were the majority in both the AS and normal tissue groups. Endothelial cells and monocytes accounted for the second largest proportions in the normal and AS tissue groups, respectively (Figure 2d-e). The expression levels of monocytes, NK cells, and T cells were higher in the AS tissue group than in the normal tissue group. The opposite was true for neutrophils (Figure 2f). We obtained 818 DEGs2 after the repeats were removed.

2.3. IL1RN, MME, ACP5 and CD36 detected and identified by MR analysis

Twenty-two genes had significant causal relationships with AS (Table 1). Acid phosphatase 5, tartrate resistant (ACP5), CD36, IL1RN, and MME were obtained by intersecting all 22 genes with the 818 DEGs2 and the 4,001 module genes (Figure 3a). The positive and negative slopes of the scatterplot show genes that are risk and protective factors for AS, respectively. Here, IL1RN and MME are protective factors while ACP5 and CD36 are risk factors for AS (Figure S4). The forest plot confirmed the results of the scatterplot (Figure S5). The funnel plot revealed that while ACP5 and CD36 conformed to Mendel’s second law, IL1RN and MME did not. A sensitivity analysis validated the reliability of the MR analysis (Figure S6). The Q p-values were all > 0.05 for ACP5, CD36, and IL1RN. Hence, there was no heterogeneity in the exposure factors or the outcome datasets. However, the Q p-value was < 0.05 for MME. Thus, the exposure factors and outcome datasets were heterogeneous (Table 1). As the IVW method was used, however, any potential heterogeneity would have only a minimal impact on the results. The horizontal pleiotropy test disclosed that all four candidate genes had p-values > 0.05. Thus, there were no confounders in the present study (Table 2). The LOO test indicated that there were no major bias points. Therefore, the results were reliable and stable (Figure 3b).

2.4. The ANN model had greater diagnostic power for AS than any other biomarker

The LASSO regression analysis identified CD36, IL1RN, and MME as candidate key genes (Figure 4a), and their expression was verified in both the GSE43292 and GSE100927 datasets. All three genes were significantly upregulated in the AS groups of both datasets (Figure 4b-c). For this reason, they were designated AS biomarkers. All three genes had AUC > 0.7 in both datasets. Hence, the biomarkers had adequate diagnostic efficacy (Figure 4d-e).

We then constructed an ANN network based on the three foregoing biomarkers. A hidden layer containing two neurons was detected (Figure 5a). The ROC curves indicated that the diagnostic power of the ANN model was greater than those of the biomarkers as the AUC values of the former were 0.842 and 1 in the GSE43292 and GSE100927 datasets, respectively (Figure 5b-c). The subcellular localization analysis revealed that all three genes were distributed primarily in the extracellular region and, to a lesser extent, in the cytoplasm, endoplasmic reticulum (ER), mitochondria, and nucleus (Figure 5d).

2.5. Significant strong positive correlations between the biomarkers and monocytes

The GSEA showed that CD36, IL1RN, and MME were significantly enriched in 63, 58, and 65 pathways, respectively. The top five pathways were displayed and they disclosed that CD36, IL1RN, and MME were enriched in lysosomes, hematopoietic cell lineage, the Toll-like receptor signaling pathway, and cytokine-cytokine receptor interaction (Figure 6a-c). The degrees of infiltration of the 22 immune cells are illustrated in Supplementary Figure S7. Of the ten immune cells obtained by filtration, five differed between the AS and normal tissue groups. Naive B cells, CD8 T cells, activated NK cells, and monocytes were significantly upregulated in the normal tissue group whereas M0 macrophages were significantly upregulated in the AS group (Figure 6d). The monocytes were obtained by intersecting five different immune cells while the nine immune cells were obtained by LASSO analysis (Figures S8 and S9). The ROC curve for the monocytes had AUC = 0.675. Hence, monocytes had relatively superior diagnostic efficacy for AS (Figure 6e). The correlation analysis showed significant correlations between the biomarkers and the monocytes (|R| > 0.4; p < 0.05) (Figure 6f).

2.6. SNHG16 regulated MME through hsa-miR-216a-3p, hsa-miR-3681-3p and hsa-miR-128-3p

We predicted three miRNAs, 11 lncRNAs, and 17 transcription factors (TFs) related to MME. We then established a ceRNA network comprising one mRNA, three miRNAs, and 11 lncRNAs as well as a TF-mRNA-miRNA network comprising one mRNA, three miRNAs, and 17 TFs (Figure 7a,b). We found that SNHG16 regulates MME via hsa-miR-216a-3p, hsa-miR-3681-3p, and hsa-miR-128-3p and predicted 12 drugs targeting CD36, IL1RN, and MME. ABT-510 targeted CD36; diacerein, haloperidol, and methotrexate targeted IL1RN; and candoxatril, LCZ696, pepinemab, sampatrilat, SLV-334, gallopamil, ilepatril, and sacubitril targeted MME (Figure 7c).

2.7. Cellular- and transcriptome-level verification of CD36, IL1RN, and MME expression

Biomarker expression was verified at the cellular level and a box plot showed that the expression score was higher for monocytes than for neutrophils, NK cells, or T cells (Figure 8a). The CD36, IL1RN, and MME expression levels were higher in the monocytes than in the other key cells (Figure 8b).

3. Discussion

Recent studies have shown that mitophagy may impede the progression of AS by scavenging damaged mitochondria, promoting mitochondrial regeneration, hindering ROS accumulation, mitigating vascular endothelial damage, inhibiting macrophage apoptosis, and suppressing the activation of thermoprotein domain-related protein 3 (NLRP3) inflammasome in the nucleotide oligomerization domain (NOD)-like receptor [9,10,19,20]. ScRNA-Seq generates detailed transcriptional profiles of the various cell types associated with AS. This technique has been widely used to elucidate the molecular mechanisms of various CVDs. MR analysis clarifies the causal relationships between the exposure factors and the outcomes of CVDs. Investigating the modes of action of MRGs in AS may facilitate the clinical diagnosis and treatment of this disorder. Here, scRNA-Seq, transcriptome analysis, and AS-related GEO datasets identified three mitophagy-related biomarkers in AS, namely, CD36, IL1RN, and MME.

The human CD36 gene is located on chromosome 7q11.2 and regulates long-chain fatty acid (LCFA) uptake [21]. We discovered that CD36 is a risk factor in AS, is upregulated in AS cells, is internalized by foam cells, and exacerbates AS progression [22]. CD36 inhibits macrophage migration in oxLDL tissue, activates the NLRP3 inflammasome, and induces an inflammatory response [23].

Most IL-1 family cytokines are secreted extracellularly. However, IL1RN encodes a signaling peptide that secretes cytokines via the ER and the Golgi apparatus [24]. Here, we observed that IL1RN exerted a protective effect against AS. IL1RN represses AS pathogenesis as it is upregulated when M2 macrophages are transformed into M1 macrophages [25].

The transmembrane glycoprotein membrane metalloendopeptidase (MME) degrades various substrates. It might be implicated in tumor onset and progression as its mRNA and protein levels are significantly lower in tumors than in adjacent healthy/normal tissues [26]. To the best of our knowledge, the present work is one of the first to find and report that MME is involved in the development of AS.

We observed that the lysosome, hematopoietic lineage, Toll-like receptor signaling, and cytokine-cytokine receptor interaction pathways were co-expressed by all three biomarkers. The aforementioned pathways are associated with AS [27,28,29,30].

Lysosomes participate in degradation, the innate and adaptive immune systems, and nutrient sensing [31]. They protect blood vessels against oxidative damage and improve vascular proliferation. Lysosomal dysfunction may contribute to AS pathology by modulating autophagy, inflammasomes, apoptosis, and lysosomal biosynthesis [30,32,33].

Reorganization of the morphology of the hematopoietic lineage reduces the numbers and alters the phenotypes of macrophages. The rapid and efficient clearance of apoptotic cells by macrophages is vital to the amelioration of AS. Nevertheless, macrophage dysfunction may accelerate the onset and progression of this condition [34,35].

Toll-like receptors (TLRs) are expressed on dendritic cells (DCs) and macrophages. They sense pathogenic factors, induce adaptive immunity, and regulate certain acute and chronic inflammatory diseases [36,37]. AS has been categorized as the latter. Endothelial dysfunction and certain immune cell interactions initiate TLRs. The activation of other inflammatory pathways also initiates and advances AS [38,39].

The biomarkers detected herein were closely associated with the AS immune microenvironment. Naive B cells, CD8 T cells, activated NK cells, and monocytes were upregulated in the normal tissue group whereas M0 macrophages were upregulated in the AS tissue group. Naive B cells are major components of adaptive immunity [40]. In response to IL-5 stimulation, they may enhance oxLDL uptake and the formation of foam cells secreting oxLDL-specific native immunoglobulin M (IgM) antibodies. Foam cells promote the clearance of apoptotic cells and, by extension, confer protection against AS [41]. T cells are key drivers and modifiers in AS pathogenesis [42]. CD8+ T cells regulate monocyte and macrophage accumulation early in AS progression [29,43]. CD8+ T cells regulate monocyte and macrophage accumulation early in AS progression [42]. In AS plaques, CD8+ T cells are cytotoxic, cause macrophage death, and form necrotic core [43]. NK cells are activated in high-fat environments, migrate toward AS lesions, and promote the death of diseased cells that abound in vulnerable unstable plaques. Hence, the killing function of NK cells is critical in plaque rupture and AS onset [44]. Monocytes are components of the innate immune system and are activated by low-density lipoproteins (LDL) and other stimulants. Activated monocytes release inflammatory factors and are considered proinflammatory cells [41]. During AS, monocytes differentiate into macrophages and the latter generate clusters of foam cells that become necrotic cores and contribute to the formation of AS plaques [45]. Within the latter, inflammatory cell infiltration results from the accumulation of substances secreted by macrophages. When M0 macrophages transform into M1 macrophages, proinflammatory cytokines are secreted and plaque inflammation is promoted. Inflammatory macrophages and foam cells also secrete proinflammatory cytokines and matrix metalloproteinases (MMPs) which destabilize and promote the rupture of plaque [46].

Here, we predicted that LCZ696 and sampatrilat target AS-related biomarkers and could ameliorate HF. Though MME biomarkers are potential therapeutic/drug targets, their specific targeting mechanisms remain to be elucidated [47,48].

4. Materials and Methods

4.1. Data extraction

The GSE43292, GSE100927, and GSE159677 datasets were acquired from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/) and their sequencing platforms were GPL6244, GPL17077, and GPL18573, respectively. GSE43292 was a training set and included 32 AS and 32 adjacent normal carotid artery tissue samples. GSE100927 was a validation set and included 29 AS and 12 adjacent normal carotid artery tissue samples. The single-cell dataset GSE159677 included three AS and three adjacent normal carotid artery tissue samples. Furthermore, 29 mitophagy-related genes (MRGs) were selected from the Reactome Pathway Database (https://reactome.org/download-data) by searching for keywords in the Molecular Signatures Database (MsigDB) (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) [49].

4.2. Weighted Gene Co-expression Network analysis (WGCNA)

The score for the single-sample gene set enrichment analysis (ssGSEA) of the MRGs calculated for all the samples in the GSE43292 dataset served as a trait. WGCNA v. 1.71 [50] was used to construct a co-expression network that would identify the module genes most strongly correlated with the ssGSEA score. The ssGSEA enrichment score was applied to all MRGs based on all the training set genes via GSVA v. 1.46.0 [51]. The Wilcoxon signed-rank test was used to compare the differences in the MRG scores between the AS and normal tissue samples. Clustering analysis was conducted on all the samples in the GSE43292 dataset to detect and eliminate any potential outliers. A soft threshold was screened to conform the gene interactions to the scale-free distribution. A dynamic tree cut (DTC) algorithm was used to identify various modules according to the predetermined soft threshold. The screening criteria were: (a) minimum module size = 100, (b) deepsplit value = 4, and (c) mergcutheight value = 0.25. The module most strongly correlated with the ssGSEA score was selected and the genes therein were designated module genes.

4.3. Single-cell RNA-seq (scRNA-Seq) analysis

Limma v. 3.50.1 [52]was used to screen differentially expressed genes (DEGs) between the AS and normal tissue groups in the GSE43292 dataset. The foregoing genes were designated DEGs1 and the screening parameters were (a) p.adjust < 0.05 and (b) |log2 fold change (FC)| > 1. A volcano plot was generated with ggplot2 v. 3.4.1 [53]to visualize the DEGs. The top ten upregulated and top ten downregulated genes were labeled according to their |log2FC| values. A DEG heatmap was plotted with ComplexHeatmap v. 2.14.0 [54].

4.4. Differential expression analysis

The GSE159677 dataset (single-cell transcriptome sequencing data) was subjected to quality control (QC) with Seurat v. 4.1.0 [55]. The screening criteria were: (a) 200–2,500 genes/cell, (b) < 10% of all genes are mitochondrial, and (c) < 3% of all genes are erythrocytic. Highly variable genes were identified using the FindVariableFeatures function in Seurat and were based on large coefficients of variation (CV) between cells. The top 2,000 highly variable genes were selected for further analysis. A principal components analysis (PCA) was performed on them to reduce their linear dimensionality and determine the optimal latitude value for cell clustering. Unsupervised clustering analysis was conducted on the filtered cells with the FindNeighbors (https://satijalab.org/seurat/reference/findneighbors) and FindClusters (https://satijalab.org/seurat/reference/findclusters) functions in Seurat v. 4.1.0 [55]. The FindAllMarkers function (https://satijalab.org/seurat/reference/findallmarkers) was employed to identify important marker genes associated with various clusters. The cell types were named based on the expression patterns of the marker genes in each cluster. SingleR v. 1.8.0 [56] was applied to notate the clustering and identify the different cell types. The cell clustering plot was visualized by uniform manifold approximation and projection (UMAP) (https://github.com/lmcinnes/umap). Scale plots were utilized to determine the abundance of each cell type in the AS and normal tissue samples. The Kruskal-Wallis test was used to screen significant differences in the cell types between the AS and normal tissues. A differential expression analysis was performed using the FindMarkers function (https://satijalab.org/seurat/reference/findmarkers) in Seurat v. 4.1.0 to identify the differentially expressed cells between the AS and normal tissue [55] as well as the differentially expressed genes (DEGs) designated DEGs2. The screening thresholds were (a) |log2FC| > 1 and (b) adjust.p < 0.05.

4.5. MR analysis between DEGs1 and AS

The eQTL (expression quantitative trait loci) dataset for DEGs1 and the ukb-d-I9_CORATHER (coronary atherosclerosis) dataset associated with AS were obtained from the Integrative Epidemiology Unit Genome-Wide Association Study (IEU OpenGWAS) Database (https://gwas.mrcieu.ac.uk/) and applied as exposure factors and outcomes, respectively. The ukb-d-I9_CORATHER dataset included 14,334 AS samples, 346,860 control samples, and 13,586,589 single-nucleotide polymorphisms (SNPs). Exposure factors and outcomes were determined and IVs were screened with TwoSampleMR(27). The SNPs that were significantly associated with the exposure factors (p < 5 × 10-6; clump = TRUE) and unrelated to the outcomes were selected as the IVs. The linkage disequilibria (LD) of the IVs were removed by setting R2 = 0.001 and kb = 10. The Harmonise_data function (https://mrcieu.github.io/TwoSampleMR/articles/harmonise.html) in TwoSampleMR v. 0.5.6 [57] unified the effect allele and effect size while the MR function was applied to combine MR Egger, weighted median, inverse variance weighted (IVW), simple mode, and weighted mode algorithms for the MR analysis [58]. IVW was the main method used (p < 0.05). The assumptions made about the IVs for the MR analysis were that they (a) are closely related to exposure, (b) are unaffected by confounding factors associated with exposure and outcome, and (c) affect the results only through exposure and not via any other biological pathways. VennDiagram v. 1.7.1 (https://cran.r-project.org/web/packages/VennDiagram/index.html) identified the candidate genes by intersecting the genes obtained by the MR analysis, the module genes, and the DEGs. The results of the MR analysis of the candidate genes were merged via scatter, forest, and funnel plots. The reliability of the MR analysis was assessed by sensitivity analysis comprising a heterogeneity test, a horizontal pleiotropy test, and a Leave-One-Out (LOO) analysis performed by the mr_heterogeneity, mr_pleiotropy_test, and mr_leaveoneout functions, respectively, in TwoSampleMR v. 0.5.6 [57].

4.6. Biomarker identification

Candidate genes in the GSE43292 dataset were identified via the least absolute shrinkage and selection operator (LASSO) regression algorithm in glmnet v. 4.1-4 [59]. Candidate key genes with consistent expression trends in both the GSE43292 and GSE100927 datasets were designated biomarkers. Receiver operating characteristic (ROC) curves were plotted to evaluate the diagnostic power of the biomarkers in the foregoing datasets.

4.7. Artificial neural network (ANN) and subcellular localization analysis

ANNs were constructed for the GSE43292 and GSE100927 datasets to assess the overall diagnostic performance of the biomarkers. The ROC curve was plotted to evaluate the diagnostic power of each ANN model. Area under the curve (AUC) values > 0.8 for the ANN model in the GSE43292 and GSE100927 datasets indicated strong diagnostic power. FASTA files were obtained from the GSE43292 dataset to clarify the mechanisms of the biomarkers in AS. The subcellular localizations of the biomarkers were predicted with mRNALocater (http://bio-bigdata.cn/mRNALocater/) [60], and their scores across different sites were calculated.

4.8. Biomarker enrichment and immune infiltration analyses

Based on the median biomarker expression values, the GSE43292 dataset was divided into high- and low-expression groups for the differential expression analysis. The |log2FC| values were calculated and arranged in descending order. A gene set enrichment analysis (GSEA) was performed using clusterProfiler v. 4.2.2 [61]. The Kyoto Encyclopedia of Genes and Genomes (KEGG) c2.cp.keg.v7.4.symbols.gmt gene set from the MSigDB Database (https://www.gsea-msigdb.org/gsea/msigdb) was the background. The CIBERSORT algorithm (https://cibersortx.stanford.edu/) was used to calculate the infiltration of 22 immune cells in all samples from the GSE43292 dataset and explore immune infiltration in AS. The immune cells were filtered based on zero expression of 22 immune cells in 75% of all samples. A box plot was then constructed to detect significant (p < 0.05) differences in the expression levels of the remaining immune cells between the AS and normal tissue samples. Then the 22 immune cells were screened by LASSO analysis. Candidate immune cells were identified from the intersection of both parts. ROC curves were then plotted to evaluate the diagnostic power of the candidate immune cells. Spearman’s rank correlation analysis was then conducted on the GSE43292 dataset to assess the relationship between the biomarkers and the candidate immune cells.

4.9. Regulatory network construction and drug prediction for biomarkers

Competing endogenous RNAs (ceRNA) and TF-mRNA-miRNA networks were constructed to elucidate the mechanism by which the biomarkers are regulated. The targetscan (http://www.targetscan.org/vert 80) and mirdb (http://www.mirdb.org) databases in multiMiR v. 1.20.0 [62] were used to predict the upstream biomarker miRNAs. The key miRNA was identified by intersecting the predicted miRNA results. The starbase database (http://starbase.sysu.edu.cn/index.php) was employed to predict the long non-coding RNAs (lncRNAs) targeting the key miRNA. The hTFtarget database (https://ngdc.ac.cn/databasecommons/database/id/6946) was applied to predict possible biomarker transcription factors (TFs). Cytoscape v. 3.9.1 [63] was utilized to visualize the ceRNA and TF-mRNA-miRNA networks. The biomarkers were imported into the Drug-Gene Interaction Database (DGIDB) (https://dgidb.org) to find small-molecule drugs that could potentially treat AS.

4.10. Validation of biomarker expression

Biomarker expression was verified in the GSE159677 dataset. AUCell v. 1.28.1 [64]analyzed the biomarker expression scores for the differentially expressed cells. Key cells were those with the highest biomarker expression scores. A box plot was constructed to represent the biomarker expression in the key cells.

4.11. Statistical analysis

R software (version 4.2.1) was used to process and analyze the data. Correlation analysis was performed using corrplot package (version 0.92). The p value less than 0.05 was considered statistically significantR v. 4.2.1 (https://cran.r-project.org/bin/windows/base/old/4.2.1/) was used to process and analyze all data. Corrplot v. 0.92 (https://github.com/taiyun/corrplot) was used to perform a correlation analysis. P < 0.05 indicated a statistically significant difference.

5. Conclusions

The present work identified three mitophagy-related AS biomarkers that could facilitate the diagnosis and treatment of AS and improve its prognosis. Nevertheless, this study had certain limitations, and the modes of action and the clinical applications of the aforementioned biomarkers must be clarified in future research.