Abstract

   Introduction: Moyamoya disease (MMD) is a chronic cerebrovascular
   disease that can lead to ischemia and hemorrhagic stroke. The
   relationship between oxidative phosphorylation (OXPHOS) and MMD
   pathogenesis remains unknown.

   Methods: The gene expression data of 60 participants were acquired from
   three Gene Expression Omnibus (GEO) datasets, including 36 and 24 in
   the MMD and control groups. Differentially expressed genes (DEGs)
   between MMD patients MMD and control groups were identified. Machine
   learning was used to select the key OXPHOS-related genes associated
   with MMD from the intersection of DEGs and OXPHOS-related gene sets.
   Gene ontology (GO), Kyoto encyclopedia of genes and genomes (KEGG),
   gene set enrichment analysis (GSEA), Immune infiltration and
   microenvironments analysis were used to analyze the function of key
   genes. Machine learning selected four key OXPHOS-related genes
   associated with MMD: CSK, NARS2, PTPN6 and SMAD2 (PTPN6 was upregulated
   and the other three were downregulated).

   Results: Functional enrichment analysis showed that these genes were
   mainly enriched in the Notch signaling pathway, GAP junction, and RNA
   degradation, which are related to several biological processes,
   including angiogenesis, proliferation of vascular smooth muscle and
   endothelial cells, and cytoskeleton regulation. Immune analysis
   revealed immune infiltration and microenvironment in these MMD samples
   and their relationships with four key OXPHOS-related genes. APC
   co-inhibition (p = 0.032), HLA (p = 0.001), MHC I (p = 0.013), T
   cellco- inhibition (p = 0.032) and Type I IFN responses (p < 0.001)
   were significantly higher in the MMD groups than those in the control
   groups. The CSK positively correlated with APC co-inhibition and T
   cell-co-inhibition. The NARS2 negatively correlated with Type I IFN
   response. The SMAD2 negatively correlated with APC co-inhibition and
   Type I IFN response. The PTPN6 positively correlated with HLA, MHC I
   and Type I IFN responses.

   Discussion: This study provides a comprehensive understanding of the
   role of OXPHOS in MMD and will contribute to the development of new
   treatment methods and exploration of MMD pathogenesis.

   Keywords: moyamoya disease, oxidative phosphorylation, pathogenesis,
   RNA sequencing, immunity

1 Introduction

   Moyamoya Disease (MMD) is a cerebrovascular disease characterized by
   chronic progressive occlusion of the distal end of the internal carotid
   artery or the main branches of the Willis circle. This occlusion causes
   the compensatory proliferation of small blood vessels at the brain
   base, presented as a ‘smoke-like’ vascular network on digital
   subtraction angiography (DSA) ([42]Gonzalez et al., 2023). MMD may have
   catastrophic consequences such as hemorrhagic or ischemic stroke. In
   addition, it can cause serious symptoms such as motor dysfunction,
   sensory dysfunction, cognitive disorders, and epilepsy ([43]Liu et al.,
   2023; [44]Fang et al., 2024). However, its pathogenesis remains
   unclear. Currently, intracranial revascularization surgery is the
   primary treatment for MMD ([45]Ihara et al., 2022). Although several
   targets related to MMD pathogenesis, such as RNF213 and DIAPH1, have
   been discovered recently, corresponding drugs have not been developed
   based on these targets ([46]He et al., 2023a; [47]Fang et al., 2024).
   Few effective drugs for treating MMD are available because of the
   unknown molecular mechanisms underlying its occurrence ([48]Ihara et
   al., 2022). Therefore, the mechanism of MMD should be urgently
   clarified to improve treatments.

   Oxidative phosphorylation (OXPHOS), a crucial biochemical process in
   eukaryotic cells, is the final stage of aerobic respiration, followed
   by glycolysis and citric acid cycle ([49]Wilson, 2017). It reduces
   oxygen to generate high-energy phosphate bonds in the form of adenosine
   triphosphate (ATP) and is the most efficient stage in the electron
   transport chain to generate ATP ([50]Nolfi-Donegan et al., 2020). The
   ATP generated through OXPHOS is the primary energy source for cells and
   plays a crucial role in various physiological and pathological cellular
   processes ([51]Lei et al., 2017). Recently, numerous studies have
   investigated the role of OXPHOS in vascular diseases. A previous study
   has reported that CHCHD4 orchestrates mitochondrial OXPHOS and
   antagonizes the aberrant growth and migration of pulmonary artery
   smooth muscle cells in pulmonary arterial hypertension (PAH) ([52]Wang
   et al., 2023). The glutamine antagonist 6-diazo-5-oxo-L-norleucine
   (DON) inhibits vascular smooth muscle cell (VSMC) proliferation and
   migration by attenuating OXPHOS in atherosclerosis ([53]Park et al.,
   2021). OXPHOS is essential for purinergic receptor-mediated angiogenic
   responses in vascular endothelial cells (ECs) ([54]Lapel et al., 2017).
   Pathological mutations in the OXPHOS system subunits can increase
   vascular endothelial growth factor (VEGF) production and pathological
   angiogenesis ([55]Bayona-Bafaluy et al., 2019). Thus, OXPHOS plays an
   important role in vascular disease pathogenesis.

   Previous studies have demonstrated that OXPHOS is related to abnormal
   EC, VSMC, and angiogenesis, which might play important roles in
   vascular disease pathogenesis. As a cerebral vascular disease, the
   pathological characteristics of MMD also includes abnormal EC
   proliferation, abnormal VSMC migration, pathological angiogenesis, and
   vascular remodeling ([56]Takagi et al., 2016; [57]He et al., 2023b;
   [58]He et al., 2023c). This generates evidence that OXPHOS disorder or
   regulation is closely related to MMD mechanism. However, studies on the
   relationship between MMD and OXPHOS remain scarce. Although a previous
   study has mentioned their relationship ([59]Kanamori et al., 2021),
   further exploration is needed to study the role of OXPHOS in MMD, such
   as signaling pathways, related genes, and the immune response.
   Therefore, investigating the association between OXPHOS and MMD is
   urgently needed.

   In this study, we collected MMD-related gene expression datasets and
   OXPHOS-related gene sets from Gene Expression Omnibus (GEO) and
   Genecards database. Machine learning algorithms of Least absolute
   shrinkage and selection operator (LASSO) and support vector
   machine-recursive feature elimination (SVM-RFE) were utilized to
   selected the key OXPHOS-genes related to MMD. We finally identified 4
   key OXPHOS-related genes, CSK, NARS2, PTPN6, and SMAD2, and analyzed
   their function by KEGG and GO enrichment analysis, which showed that
   these genes were most relevant to MMD pathogenesis. GSEA analysis were
   also performed to provide the profile of immune infiltration and
   microenvironment about these genes. Ultimately, we identified the
   associations between these 4 key OXPHOS-related genes and the
   pathogenesis of MMD.

2 Materials and methods

2.1 Data acquisition

   Gene expression data were retrieved from RNA sequencing data from the
   GEO database ([60]https://www.ncbi.nlm.nih.gov/geo/info/datasets.html).
   The series matrix data files [61]GSE189993, [62]GSE157628, and
   [63]GSE141024 enrolled 32, 20, and 8 patients, respectively, including
   11, 9, and 4 in the control group, respectively, and 21, 9, and 4 in
   the MMD group, respectively. All the series matrix files were annotated
   using the [64]GPL16699 file. OXPHOS-related gene sets were downloaded
   in the GeneCards database ([65]https://www.genecards.org/). Those
   OXPHOS-related genes include protein-coding genes, functional element
   genes and RNA genes. The relevance score of GeneCards is a
   comprehensive score that measures the correlation between genes and
   specific keywords or topics based on multiple data sources and
   algorithms. Choosing a higher reliability score threshold can reduce
   the interference of noise and unrelated genes, enhance the reliability
   of analysis results, and ensure that the number of genes is not too
   large. Referring to several previous studies about screening genes with
   specific functions in GeneCards dataset ([66]Li et al., 2022; [67]Liu
   et al., 2022), we set the threshold of relevance score >15 to ensure
   that the selected genes have sufficient correlation and potential
   biological significance of oxidative phosphorylation.

2.2 Gene differential expression analysis

   After merging the three datasets, the correction of the microarray data
   was performed using a surrogate variable analysis (SVA) package to
   reduce batch effects among them. Subsequently, we used principal
   component analysis (PCA) to draw a sample distribution map before and
   after correction, visually demonstrating the changes in batch effects
   before and after correction. The R package ‘limma’ was utilized to
   identify the differentially expressed genes (DEGs) between MMD and
   control samples by analyzing the gene expression data from the GEO
   cohort. Due to the small sample size in this study, p-value <0.05 was
   set as the criteria for filtering the DEGs to reduce the omissions of
   some potential biologically significant genes, rather using padj or FDR
   value. Meanwhile, the |log2 Fold Change (FC)| >1 was also set as the
   criteria for DEGs identification combined with p < 0.05 to ensure that
   the screened genes have significant differences in expression levels.
   Therefore, the final criteria of filtering DEGs was set as p < 0.05 and
   |log2 FC| >1, which was also be utilized in previous studies ([68]Liu
   et al., 2020; [69]Wang et al., 2022).

2.3 GO and KEGG functional enrichment analysis

   GO and KEGG enrichment analyses were performed to evaluating the
   biological function, biological process, molecular function, and
   cellular component of the intersection genes between DEGs and
   OXPHOS-related gene set. When performing the enrichment analysis, we
   utilized the R package ‘clusterProfiler’ to perform the GO and KEGG
   enrichment analysis. The default threshold for ‘clusterProfiler’
   package to do enrichment analysis is a p-value <0.05 and a Q value
   <0.2. In this study, to make the enrichment results more accurate, we
   further limited the Q value as less than 0.05, which is a relative
   strict standard. Finally, we set the cutoffs of GO and KEGG enrichment
   analysis as a p-value <0.05 and a Q value <0.05, referring to several
   previous studies ([70]Liu et al., 2023; [71]Li et al., 2024).

2.4 Key genes and possible disease diagnostic indicator selection

   LASSO logistic regression and SVM-RFE algorithm were used to select the
   potential disease diagnostic markers. LASSO logistic regression was
   conducted using ‘glmnet’ package. When performing LASSO logistic
   regression, we used the ‘glmnet’ package to fit gene expression levels
   based on the glmnet function, with the family set as gaussian
   distribution, and automatically selects the optimal lambda value. By
   using the ‘cv.glmnet’ function for 10 fold cross validation, the
   optimal regularization parameter lambda was selected, and the
   performance curves of the model were plotted for each lambda value
   during the cross validation process. Finally, we selected the optimal
   lambda value by observing the lowest point of the curve. SVM-RFE
   machine learning model was constructed using ‘e1071’ package to further
   identify the significance of indicators for disease diagnosis. When
   performing SVM-RFE, we first used the SVM-RFE method to perform feature
   selection on DEGs to obtain the ranking of each gene. SVM model was
   established using the “e1071”R package, and the SVM model is
   constructed based on the top 10 genes to evaluate the lowest error rate
   of the model. The gene sets with the error rate reaching the lowest
   point and the genes obtained by the LASSO were intersected to obtain
   the key genes.

2.5 Gene set enrichment pathway analysis

   Gene set enrichment analysis (GSEA) was used to investigate the DEGs
   between the high- and low-expression groups based on the expression
   profiles of patients with MvMD. GSEA was performed using a predefined
   gene set ‘c2. cp. Kegg’ to rank them according to their differential
   expression levels in the two types of samples. The preset gene sets
   were examined to determine whether they were enriched at the top or
   bottom of the ranking table. This study used GSEA software
   ([72]http://www.broadinstitute.org/gsea) to compare the differences in
   KEGG signaling pathways between high expression and low expression
   groups, and explored the molecular mechanisms of key genes in the two
   groups of patients. The times of substitutions was set as 1,000 and the
   type of substitution was set as ‘phenotype.’ Pathways with a p-value
   <0.05 were considered significantly enriched.

2.6 Immune cell infiltration analysis

   Single-sample GSEA (ssGSEA) was performed to evaluate the immune cells
   types in the microenvironment using ‘GSVA’ R package. We performed
   ssGSEA to quantify the immune cells in the expression profile and
   calculated the relative proportions of 19 types of infiltrating immune
   cells, including activated Dendritic Cells (aDCs), Antigen-Presenting
   Cells co-inhibition (APC_co_inhibition), Antigen-Presenting Cells
   co-stimulation (APC_co_stimulation), B lymphocytes (B_cells), C-C
   Chemokine Receptor (CCR), Checkpoint, Human Leukocyte Antigen (HLA),
   Inflammation-promoting, Macrophages, Major Histocompatibility Complex
   class I (MHC_class_I), Neutrophils, Parainflammation, plasmacytoid
   Dendritic Cells (pDCs), T cell co-inhibition, T follicular helper cells
   (Tfh), Tumor-Infiltrating Lymphocytes (TIL), Regulatory T cells (Treg),
   Type I Interferon Response (Type_I_IFN_Reponse) and Type II Interferon
   Response (Type_II_IFN_Reponse). Spearman’s correlation analysis was
   used to analyze the gene expression and immune cell content. p < 0.05
   was considered statistically significant.

2.7 Statistical analysis

   All statistical analyses were performed using the R software (version
   4.2.2) and the figures were plotted using the ‘ggplot’ R package
   (version 3.3.6). Statistical significance was set at p < 0.05.

3 Results

3.1 DEG identification and functional enrichment analysis in the MMD cohort

   The MMD-related datasets [73]GSE189993, [74]GSE157628, and
   [75]GSE141024 were downloaded from the GEO database and used in this
   study, including 60 samples from patients with MMD (n = 36) and control
   participants (n = 24). The detailed clinical information for the
   participants was listed in [76]Supplementary Table S1 in
   [77]Supplementary Material. Correction was performed using the SVA
   algorithm to reduce the batch effect between microarrays. The PCA plot
   shows that this process reduced the batch effect ([78]Figure 1A).

FIGURE 1.

   [79]FIGURE 1
   [80]Open in a new tab

   Screening of the phosphorylation-associated genes for MMD and
   functional enrichment. (A) PCA plot for gene expression data of the
   participants. (B) Volcano plot of DEGs differentially expressed between
   MMD patients and control participants in the GEO dataset. (C) Venn plot
   of the overlapping DEGs and phosphorylation-associated genes. Finally,
   27 intersecting genes were selected. (D) GO enrichment analysis for the
   27 intersecting genes. E, KEGG enrichment analysis for the 27
   intersecting genes. MMD, Moyamoya disease; PCA, Principal component
   analysis; GEO, Gene Expression Omnibus; DEGs, Differentially expressed
   genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and
   Genomes.

   A total of 1,548 genes were differentially expressed between patients
   with MMD and control participants (850 were upregulated and 698 were
   downregulated in the MMD group), according to p < 0.05 and |logFC|>1.
   The volcano plot displays an overview of the DEGs between patients with
   MMD and control participants ([81]Figure 1B). The detailed information
   of these 1,548 DEGs were listed in [82]Supplementary Table S2 in
   [83]Supplementary Material. A total of 389 OXPHOS-related genes with
   >15 relevance score was screened from OXPHOS-associated gene sets in
   the GeneCards database ([84]https://www.genecards.org/). The detailed
   information of these 389 OXPHOS-related genes were listed in
   [85]Supplementary Table S3 in [86]Supplementary Material. Finally, 27
   intersecting genes were selected from the screened OXPHOS-associated
   genes and genes that were differentially expressed between the patients
   with MMD and control participants ([87]Figure 1C). The 27 intersecting
   genes were also listed in [88]Supplementary Table S3 in
   [89]Supplementary Material.

   GO biological process enrichment analysis revealed that these 27 genes
   were mainly enriched in peptidyl tyrosine phosphorylation, peptidyl
   tyrosine modification, and gland development. GO molecular function
   enrichment analysis revealed that these genes were mainly enriched in
   protein tyrosine kinase activity, phosphotyrosine residue binding and
   phosphatase binding ([90]Figure 1D). The GO cellular component
   enrichment was not significant. KEGG enrichment analysis showed that
   these genes were mainly enriched in the apelin signaling pathway,
   PI3K–Akt signaling pathway, and endocrine resistance ([91]Figure 1E).

3.2 Identification of key phosphorylation-associated DEGs by machine learning

   To further identify the key genes affecting MMD, LASSO logistic
   regression and SVM algorithms were performed based on the previously
   obtained intersecting genes. The results identified 13 genes as MMD
   characteristic genes using LASSO regression ([92]Figures 2A,B). We also
   used the SVM algorithm to identify characteristic genes for MMD and
   screened the top seven genes with the lowest error rates ([93]Figure
   2C). The genes screened by LASSO regression and the SVM algorithm were
   then intersected, and four genes were selected, including three
   downregulated genes, CSK, NARS2, and PTPN6, and one upregulated gene,
   SMAD2 ([94]Figure 2D). These four genes were considered key genes in
   MMD and further investigated.

FIGURE 2.

   [95]FIGURE 2
   [96]Open in a new tab

   Identification for the key phosphorylation-associated genes for MMD by
   machine learning algorithm and functional enrichment. (A, B)
   Identification for characteristic genes for MMD by LASSO regression
   algorithm. Finally, 13 intersecting genes were selected. (C)
   Identification for characteristic genes for MMD by SVM algorithm.
   Finally, 7 intersecting genes were selected. (D) Venn plot of the
   overlapping genes identified by LASSO regression and SVM algorithm.
   Finally, 4 intersecting genes were selected as key
   phosphorylation-associated genes for MMD. (E) GSEA analysis for gene
   CSK. (F) GSEA analysis for gene NARS2. (G) GSEA analysis for gene
   PTPN6. (H) GSEA analysis for gene SMAD2. MMD, Moyamoya disease; LASSO,
   Least Absolute Shrinkage and Selection Operator; SVM, Support Vector
   Machine; GSEA, Gene set enrichment analysis.

3.3 Functional enrichment analysis for the key genes in MMD

   Next, we investigated the signaling pathways enriched in these four
   genes and explored the potential molecular mechanisms affecting MMD
   progression. GSEA results showed that CSK was enriched in the
   adipocytokine signaling pathway, neurotrophin signaling pathway,
   proteasome, neuroactive ligand–receptor interaction, olfactory
   transduction, and RNA degradation ([97]Figure 2E). NARS2 was enriched
   in aminoacyl-tRNA biosynthesis, inositol phosphate metabolism,
   neuroactive ligand–receptor interaction, olfactory transduction, PPAR
   signaling pathway, and RNA degradation ([98]Figure 2F). PTPN6 was
   enriched in arachidonic acid metabolism, the B cell receptor signaling
   pathway, glycosylphosphatidylinositol (GPI) anchor biosynthesis,
   neuroactive ligand–receptor interaction, primary immunodeficiency, and
   proteasomes ([99]Figure 2G). SMAD2 was enriched in the cytosolic NDA
   sensing pathway, GAP junction, glycerolipid metabolism, Notch signaling
   pathway, pentose phosphate pathway, and toll-like receptor signaling
   pathway ([100]Figure 2H).

3.4 Analysis of immune infiltration and microenvironment

   Next, the correlation between the four key genes and immune cell
   infiltration was analyzed. The heatmap shows the proportion of 19
   immune cell types in each sample ([101]Figure 3A). A thermodynamic
   chart displays the relationship with each immune cell type ([102]Figure
   3B). The results showed that APC co-inhibition (p = 0.032), CCR (p =
   0.019), checkpoint (p = 0.001), HLA (p = 0.001), MHC I (p = 0.013), T
   cell-co-inhibition (p = 0.032), TFH cells (p < 0.001), TIL cells (p =
   0.001), Treg cells (p = 0.012), and Type I IFN responses (p < 0.001)
   were significantly higher in the MMD groups than those in the control
   groups ([103]Figure 3D).

FIGURE 3.

   [104]FIGURE 3
   [105]Open in a new tab

   Analysis of immune infiltration and microenvironment of the four key
   phosphorylation-associated genes for MMD. (A) Heatmap for the
   proportion of 19 kinds of immune cells in each sample. (B)
   Thermodynamic plot of the relationship between each kind of immune
   cells. (C) The relationship between the gene CSK (upper of the figure)
   and NARS2 (lower of the figure), and immune cells. (D) Violin plot for
   the difference of immune infiltration between MMD and control group.
   MMD = Moyamoya disease.

   Next, we analyzed the relationships between the four key genes and
   immune cells. The results showed that CSK significantly positively
   correlated with T cell co-inhibition and APC co-inhibition, and
   significantly negatively correlated with B cells and inflammation
   promotion ([106]Figure 3C). NARS2 significantly negatively correlated
   with pDCs and the Type I IFN response ([107]Figure 3C). SMAD2
   positively correlated with macrophages, but significantly negatively
   correlated with TFH cells, APC co-inhibition, pDCs, and the Type I IFN
   response ([108]Supplementary Figure S1F). PTPN6 positively correlated
   with TIL, CCR, macrophages, HLA, neutrophil MHC I, and inflammation
   ([109]Supplementary Figure S1F). Bubble plots showed the relationship
   between the four key genes and the immune microenvironment, including
   chemokine, immunoinhibitor, receptor, immunostimulatory, and MHC
   ([110]Supplementary Figures S1A–E).

3.5 In depth analysis of the correlation between the key genes and MMD

   To further explore the potential molecular mechanisms of the four key
   genes in MMD, we analyzed other related genes. We first analyzed
   differences in the expression of MMD-related genes obtained from the
   GeneCards database between the MMD and control groups. The results
   showed that APOE, DNM2, IL10, TGFβ1, and VWF expression levels in MMD
   group were significantly higher than that in control group ([111]Figure
   4A). Next, we analyzed the correlation between MMD-related genes and
   these four key genes. Significant correlations were observed among
   these variables. Among these, PTPN6 was significantly and positively
   correlated with DNM2, whereas CSK was significantly and negatively
   correlated with APOE ([112]Figure 4B).

FIGURE 4.

   [113]FIGURE 4
   [114]Open in a new tab

   Analysis of the correlation between the key genes and MMD. (A) Box plot
   for the expression levels of DEGs for MMD related to the key
   phosphorylation-associated genes. (B) Bubble plot for the relationship
   between DEGs for MMD and the key phosphorylation-associated genes. MMD
   = Moyamoya disease; DEGs = Differentially expressed genes.

4 Discussion

   A considerable proportion of patients with MMD may suffer from serious
   cerebral hemorrhage or ischemic stroke; thus, timely intervention is
   required. Currently, intracranial revascularization surgery is the main
   treatment for treating MMD ([115]Ihara et al., 2022). Due to possible
   postoperative complications, such as stroke or infection, surgery
   cannot achieve a satisfactory prognosis in some patients with MMD
   ([116]Kuroda and Houkin, 2008; [117]Liu et al., 2023). Therefore, novel
   treatment methods should be urgently developed. However, owing to the
   unclear mechanism of MMD, there is a lack of effective drugs to prevent
   its progression, as well as effective animal models ([118]Rallo et al.,
   2021; [119]Ihara et al., 2022). Although some biological targets have
   recently been identified by genomics, proteomics, or other methods,
   animal modeling and drug research based on these targets have achieved
   unsatisfactory results ([120]Rallo et al., 2021; [121]Ihara et al.,
   2022). This suggests that other mechanisms may be involved in MMD
   development. Several studies have reported the role of OXPHOS in
   vascular diseases ([122]Lapel et al., 2017; [123]Bayona-Bafaluy et al.,
   2019; [124]Park et al., 2021; [125]Wang et al., 2023), indicating that
   OXPHOS may also be involved in MMD pathogenesis. However, few studies
   have examined the relationship between OXPHOS and MMD. Therefore, in
   this study, we systematically investigated the role of OXPHOS-related
   genes, signaling pathways, and immune responses in the pathogenesis of
   MMD using bioinformatics. Four key genes were selected: CSK, NARS2,
   PTPN6, and SMAD2.

   SMAD2 is a signal transducer and transcriptional modulator that
   mediates multiple signaling pathways ([126]Nomura and Li, 1998). This
   protein mediates transforming growth factor beta (TGFβ) signal, and
   thus regulates multiple cellular processes, such as cell proliferation,
   apoptosis, and differentiation ([127]Kamato and Little, 2020). In our
   study, we found a significantly downregulation of SMAD2 expression. As
   for previous studies of SMAD2, Petersen et al. found that SMAD2
   knockdown developed an aggressive phenotype in breast cancer cells and
   enhanced angiogenesis ([128]Petersen et al., 2010). Cao et al. found
   that activin receptor-like kinase 7 promotes VSMC apoptosis by
   activating the Smad2/3 signaling in diabetic atherosclerosis ([129]Cao
   et al., 2022). These studies provided evidence that SMAD2
   downregulation in patients with MMD might promote angiogenesis,
   proliferation and phenotype switching of VSMCs, which are the important
   characteristics of MMD mentioned before ([130]Ma et al., 2022; [131]He
   et al., 2023c; [132]He et al., 2023d). As for the results in our study,
   pathway enrichment analysis also showed that SMAD2 was enriched in the
   Notch, GAP junction, and toll-like receptor signaling pathways, which
   are also associated with angiogenesis, vascular development,
   endothelial cell and VSMC growth in previous studies ([133]Beyer and
   Berthoud, 2018; [134]Shafeghat et al., 2022; [135]Cuervo et al., 2023).
   In summary, these results generated the evidence that SMAD2 may play an
   important role in MMD pathogenesis and vascular changes.

   CSK is a negative regulator of Src family kinases (SFKs) recruited to
   the plasma membrane by binding to transmembrane or adapter proteins to
   suppress signaling via various surface receptors by phosphorylating and
   maintaining several inactive positive effectors, such as FYN or LCK
   ([136]D and CJB, 1999; [137]Sun and Ayrapetov, 2023). It plays
   important roles in cell growth, migration, differentiation, and immune
   response ([138]Zhu et al., 2023). There were several studies about the
   CSK on endothelial cells. Alejandro et al. found that CSK activation
   renders the endothelial barrier more susceptible to TNF-α through p38
   activation, which leads to endothelial tight junction loss ([139]Adam
   et al., 2016). Ulf et al. reported that CSK is involved in endothelial
   cell proliferation inhibition ([140]Baumeister et al., 2005). Previous
   studies showed that the CSK has an inhibitory effect on the
   proliferation of endothelial cells. In this study, we found a
   downregulated of CSK expression in the MMD group. Accordingly, we
   inferred that may promote the angiogenesis through stimulating the
   endothelial cell proliferation, which may cause a intimal thickening of
   MMD vessels reported in previous study ([141]Takagi et al., 2016;
   [142]He et al., 2023b).

   NARS2 encodes asparaginyl-tRNA synthetase 2, a putative member of the
   class II family of aminoacyl-tRNA synthetases. This enzyme catalyzes
   asparagine attachment to mitochondrial tRNA and plays a critical role
   in mitochondrial protein biosynthesis ([143]Duchene et al., 2009).
   NARS2 mutations are associated with a combined oxidative
   phosphorylation deficiency 24 and autosomal recessive deafness-94
   ([144]Sissler et al., 2017). In this study, NARS2 was downregulated but
   there is a lack of research about the relationship between NARS2 and
   MMD. Due to its association with cellular energy metabolism, we
   inferred that NARS2 downregulation in MMD may lead to a deficiency in
   mitochondrial protein synthesis and affect the energy supply within the
   cell, leading to VSMC and EC dysfunction in MMD vessels.

   PTPN6 is a member of the protein tyrosine phosphatase (PTP) family, a
   group of enzymes that remove phosphate groups from tyrosine residues by
   hydrolyzing phosphoric acid monoesters. PTPs are signaling molecules
   that regulate several cellular processes, including cell growth,
   differentiation, mitotic cycle, and oncogenic transformation ([145]Amin
   et al., 2021). In this study, PTPN6 was the only upregulated among the
   four key genes. In previous studies, Clement et al. found that PTPN6
   overexpression fully restored PDGF-induced proliferation, migration,
   and signaling pathways in SMC exposed to high glucose and hypoxia
   ([146]Clément et al., 2021). As for the role of VSMC in MMD, a previous
   study reported VSMC over-proliferation in the intima and intimal
   thickening of pathological blood vessels in MMD ([147]Takekawa et al.,
   2004). This generates evidence that PTPN6 overexpression in MMD may
   promote VSMC proliferation and subsequently cause intimal thickening.
   In addition, PTPN6 expression was significantly positively correlated
   with the expression of dynamin 2 (DNM2) in this study, a GTP-binding
   protein subfamily and an important component of the cellular
   cytoskeleton ([148]Ferguson and De Camilli, 2012). Interestingly, our
   previous study also found that cytoskeletal remodeling in ECs plays an
   important role in MMD pathogenesis ([149]He et al., 2023b). This
   suggests that PTPN6 and DNM2 are involved in MMD by participating in
   regulating the cellular cytoskeleton.

   Our study has limitations. First, the sample size collected in this
   study was relatively small, which may have reduced the
   representativeness of the study and weaken the level of evidence when
   drawing conclusions about causality between identified genes and MMD
   pathogenesis. Further study needs to be conducted in a larger and
   ethnically diverse patient population to strengthen the
   generalizability of the findings. Second, the key genes and their
   functions in this study are lack of validation through in vitro
   experiments. Functional experiments are needed to validate the impact
   of the identified genes on relevant cellular processes in MMD, for
   example, the cell proliferation assays or migration assays on the gene
   overexpression or knockdown cell lines. Finally, in order to avoid the
   omissions of some potential biologically significant genes in the small
   sample size, this study used the criteria of p < 0.05 and | log2FC | >
   1 to identify DEGs, without using adjusted p-values or FDR values. This
   method may lead to a higher false positive rate, which is also a
   limitation of this study.

5 Conclusion

   In summary, this study identified four key OXPHOS-related genes
   associated with MMD: CSK, NARS2, PTPN6, and SMAD2. We systematically
   investigated the relationship between four key OXPHOS-related genes and
   the pathogenesis of MMD We also analyzed the profile of immune
   infiltration and the microenvironment of these genes. Bioinformatic
   analysis showed that the four key genes are related with angiogenesis,
   endothelial cell and VSMC proliferation, and VSMC phenotype switching,
   which may be involved in the mechanisms of MMD. Further studies will be
   needed to investigate the role of these identified genes in the
   pathogenesis of MMD in a larger and ethnically diverse patient
   population, or by in vivo and in vitro experiments. In conclusion, this
   study provides a comprehensive understanding of the important role of
   OXPHOS in MMD pathogenesis and will be conducive to the development of
   new treatment methods for MMD.

Funding Statement

   The author(s) declare that financial support was received for the
   research, authorship, and/or publication of this article. This study
   was supported by the Natural Science Foundation of China (82171887 and
   82371296 to RW).

Data availability statement

   The data presented in the study are deposited in the GEO repository
   ([150]https://www.ncbi.nlm.nih.gov/geo/info/datasets.html), accession
   number [151]GSE189993, [152]GSE157628, [153]GSE141024.

Ethics statement

   The studies involving humans were approved by the Institutional Ethics
   Committee of Beijing Tiantan Hospital, Beijing, China (KY 2020-045-02).
   The studies were conducted in accordance with the local legislation and
   institutional requirements. Written informed consent for participation
   was not required from the participants or the participants' legal
   guardians/next of kin in accordance with the national legislation and
   institutional requirements.

Author contributions

   ZH: Conceptualization, Formal Analysis, Investigation, Methodology,
   Software, Visualization, Writing–original draft, Writing–review and
   editing. JZ: Conceptualization, Formal Analysis, Investigation,
   Methodology, Software, Visualization, Writing–original draft,
   Writing–review and editing. YS: Writing–original draft. ZZ:
   Writing–original draft. YW: Writing–original draft. SX:
   Writing–original draft, Formal Analysis, Visualization. YZ:
   Supervision, Writing–original draft, Writing–review and editing. SH:
   Conceptualization, Data curation, Formal Analysis, Investigation,
   Methodology, Project administration, Resources, Software, Supervision,
   Validation, Visualization, Writing–original draft, Writing–review and
   editing. RW: Funding acquisition, Investigation, Methodology, Project
   administration, Resources, Software, Supervision, Validation,
   Visualization, Writing–review and editing.

Conflict of interest

   The authors declare that the research was conducted in the absence of
   any commercial or financial relationships that could be construed as a
   potential conflict of interest.

Publisher’s note

   All claims expressed in this article are solely those of the authors
   and do not necessarily represent those of their affiliated
   organizations, or those of the publisher, the editors and the
   reviewers. Any product that may be evaluated in this article, or claim
   that may be made by its manufacturer, is not guaranteed or endorsed by
   the publisher.

Supplementary material

   The Supplementary Material for this article can be found online at:
   [154]https://www.frontiersin.org/articles/10.3389/fgene.2024.1417329/fu
   ll#supplementary-material
   [155]Table2.XLSX^ (157.7KB, XLSX)
   [156]Table3.XLSX^ (41.9KB, XLSX)
   [157]Table1.DOCX^ (342.2KB, DOCX)
   [158]Table4.DOCX^ (38KB, DOCX)

References