Introduction: This study aimed to elucidate mechanisms underlying moyamoya disease (MMD) pathogenesis and to identify potential novel biomarkers. We utilized gene co-expression networks to identify hub genes associated with the disease. Methods: Twenty-one middle cerebral artery (MCA) samples from MMD patients and 11 MCA control samples were obtained from the Gene Expression Omnibus (GEO) dataset, GSE189993. To discover functional pathways and potential biomarkers, weighted gene co-expression network analysis (WGCNA) was employed. The hub genes identified were re-assessed through differential gene expression analysis (DGEA) via DESeq2 for further reliability verification. Additional 4 samples from the superficial temporal arteries (STAs) from MMD patients were obtained from GSE141025, and a subgroup analysis stratified by arterial type (MCA vs. STA) DGEA was performed to assess if the hub genes associated with MMD are expressed significantly greater on the affected arteries compared to healthy ones in MMD. Results: WGCNA revealed a predominant module encompassing 139 hub genes, predominantly associated with the neuroactive ligand-receptor interaction (NLRI) pathway. Of those, 17 genes were validated as significantly differentially expressed. Neuromedin U receptor 1 (NMUR1) and thyrotropin-releasing hormone were 2 out of the 17 hub genes involved in the NLRI pathway (log fold change [logFC]: 1.150, p = 0.00028; logFC: 1.146, p = 0.00115, respectively). MMD-only subgroup analysis stratified by location showed that NMUR1 is significantly overexpressed in the MCA compared to the STA (logFC: 1.962; p = 0.00053) which further suggests its possible localized involvement in the progressive stenosis seen in the cerebral arteries in MMD. Conclusion: This is the first study to have performed WGCNA on samples directly affected by MMD. NMUR1 expression is well known to induce localized arterial smooth muscle constriction and, recently, type 2 inflammation which can predispose to arterial stenosis potentially advancing the symptoms and progression of MMD. Further validation and functional studies are necessary to understand the precise role of NMUR1 upregulation in MMD and its potential implications.

Moyamoya disease (MMD) is characterized by progressive stenosis of the distal internal carotid artery and often includes the proximal anterior, middle, or more rarely the posterior cerebral arteries. Although MMD can occur at any age, it typically shows a bimodal distribution with peaks at ages 5–9 years and 35–40 years. While progressive, revascularization surgery can prevent future ischemic stroke associated with MMD.

The exact etiology of MMD is unknown; however, current literature suggests that its onset might be influenced by various factors, including genetic predispositions, abnormalities in embryonic cerebrovascular development stage, autoimmune and inflammatory processes, as well as environmental factors [1]. Regarding genetic research, association studies employing computational methods have shed light on the genetic basis of MMD, which partially supports the pathological hallmarks. These studies have revealed association with alterations in genes, like MMP, VEGF, PDGFRB, and TGFB, which are substantially involved in the normal functioning of blood vessels [2].

Various computational methods, such as differential gene expression analysis (DGEA) and gene set analysis, have been utilized in the several dozen MMD association studies in the present literature [3]. Their primary goals were to identify potential biomarkers and elucidate biological pathways implicated in MMD. However, such computational methods typically focus on genes exhibiting enrichment in the diseased cohort compared to the control group. As a result, they often miss revealing the functional interconnections among the identified genes. In contrast, weighted gene co-expression analysis (WGCNA) is a powerful bioinformatic model which offers the unique ability to identify specifically hub genes. Hub genes are significant in comparison with normal genes because they tend to have a higher degree of connectivity within biological networks, indicating their central roles in regulating various cellular processes or interlinking biochemical pathways. This centrality makes hub genes influential in controlling key pathways and functions, and their dysregulation can have a more profound impact on cellular physiology and disease processes compared to less connected genes. This research represents the first application of WGCNA to samples obtained from the cerebral arteries affected by MMD. Performing WGCNA on the affected tissues themselves is important given that MMD is considered a localized disease and not systemic. Thus, studies which utilize WGCNA on MMD patients’ blood samples can and have identified strong biomarkers and, however, may not necessarily be able to identify the hub genes directly involved with the pathogenesis [4]. The primary objectives were to enhance our understanding of the underlying pathomechanism and identify novel biomarkers.

Gene Expression Datasets and Data Processing

Transcriptome-wide gene expression profiles of GSE189993 and GSE141025 were downloaded from the Gene Expression Omnibus (GEO) database. The flowchart for the present study is illustrated in Figure 1. For step 1, only the microarray dataset of GSE189993 was used which contains 21 middle cerebral artery (MCA) tissue samples from MMD patients and 11 control MCA tissue samples from non-MMD patients based on the platform of GPL16699, Agilent-039494 SurePrint G3 Human GE v2 8 × 60K Microarray 039381 (Feature Number version). The raw read count profile of GSE189993 was processed by the limma package for background correction and normalization. The “biomaRt” package was used for probe annotation. The total gene expression dataset, consisting of 16,671 genes, was further processed for WGCNA. Clinical traits were gathered from the information of GSE189993 which included disease status (state), gender, and age. All samples were divided into a control subgroup and MMD subgroup based on whether MMD was reported in the “state” group and in addition, in the same way, MMD or not in status group, female or male in gender group. The age of each patient, from 2 to 79, was recorded in age group.

Fig. 1.

Flowchart on construction of the pathogenesis-related gene of MMD occurrence. This flowchart was generated and edited with the free application online, diagrams.net (https://app.diagrams.net/). MCA, middle cerebral artery; STA, superficial temporal artery; WGCNA, weighted gene co-expression network analysis; DEGs, differential expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology; NMUR1, neuromedin U receptor 1; TRH, thyrotropin-releasing hormone.

Fig. 1.

Flowchart on construction of the pathogenesis-related gene of MMD occurrence. This flowchart was generated and edited with the free application online, diagrams.net (https://app.diagrams.net/). MCA, middle cerebral artery; STA, superficial temporal artery; WGCNA, weighted gene co-expression network analysis; DEGs, differential expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology; NMUR1, neuromedin U receptor 1; TRH, thyrotropin-releasing hormone.

Close modal

For step 2, the microarray dataset of GSE141025 provided the gene expression profile of superficial temporal artery (STA) tissue samples from 4 MMD patients based on the same previously mentioned platform of GPL16699, Agilent-039494 SurePrint G3 Human GE v2 8 × 60K Microarray 039381 (Feature Number version). These 4 samples were to be merged with 4 of the MMD MCA samples from GSE189993 to perform an MMD-only subgroup analysis stratified by artery type. This analysis was performed to further assess if certain MMD-related hub genes of interest (e.g., ones coding for receptors involved with pathways known to be associated with arterial stenosis) found in step 1 were mainly expressed in the MCA. A significantly greater expression by the MCA when compared to the STA could suggest that the hub gene might have a role in the localized pathogenesis observed at the cerebral arteries depending on the gene’s function. Differentially expressed genes (DEGs) in GSE189993 and GSE141025 were identified using the following criteria, p value <0.05 and log fold change (|logFC|)> 0.5, and are listed in supplementary Tables 3 and 4, respectively.

Patient Characteristics

The 21 patients with MMD and 11 patients with control conditions underwent neurosurgical treatment between May 2014 and September 2020 at the Nagoya University Hospital [5]. All patients with MMD were diagnosed according to the guidelines proposed by the Ministry of Health and Welfare of Japan and underwent STA-MCA anastomosis, for preventing stroke onset. The protocol was approved by the Institutional Review Board of the Nagoya University Graduate School of Medicine, and written informed consent was obtained from each participant or legal guardian by the original investigators (Mamiya et al. [5]). The control group included 6 patients with intracranial artery aneurysm and 5 patients with medically intractable seizures. Patients with ICA aneurysms underwent aneurysm trapping or proximal ligation, which required STA-MCA anastomoses for flow constructions, and patients with medically intractable seizures underwent epileptic focus resection. The control patients were enrolled in a consecutive manner.

Seven MMD patients and 2 epileptic patients would later go on to be excluded from the present study as their samples were identified as batch effects, or in other words, unwanted technical variations that are introduced in data when experiments are conducted at different times or under slightly different conditions. Thus, a total of 23 patients (14 MMD and 9 controls) were included in the study. Their individual patient characteristics are reported in online supplementary Table 1 (for all online suppl. material, see https://doi.org/10.1159/000538035), and their ongoing prescribed medications are given in online supplementary Table 2.

Intraoperative Sample Collection Procedure

The authors of the microarray dataset reported pathological characteristics of MMD samples included thinned media and thickened intima in the C7 segment of the ICA [5]; these changes have also been observed in the M4 segment of the MCA artery [6‒8]. Hence, the collected specimens of the M4 were considered to be MMD lesions. In the 14 patients with MMD and the 6 with ICA aneurysms, a small sample of the MCA in the M4 segment was collected by arteriotomy during the STA-MCA anastomosis. The vessel wall of the M4 was excised following a previously reported procedure [8]. For the 3 patients with intractable seizures, the M4 was dissected from the resected brain tissue, and in these patients, the vessel wall of the MCA was excised according to another previously reported method [6]. In patients with intractable seizures, the M4 was dissected from the resected brain tissue. Obtained samples were washed with normal saline containing 10 U/mL heparin, collected in Qiagen, and stored at −80°C until further use.

Construction of WGCNA

We used the whole expression profile of GSE189993 to perform WGCNA. A topology network of WGCNA was constructed by the R package “WGCNA” for analysis. The adjacency matrix was transformed to a topological overlap matrix; the soft-threshold power (β) was eight when 0.85 was chosen as the correlation coefficient threshold with a scale-free R2 above 0.85 and a slope near 1 [9]. The minimum number of genes in the modules was 59, and the cut height threshold was 0.25 for module detection. We constructed a module-trait relationship chart and searched for the key modules that were closely related to type (MMD or controls), sex (male or female), and age (years).

The correlation between the eigengene module and clinical traits was used to estimate the module-trait relationship and detect highly phenotypically related modules efficiently. We obtained the module significance by calculating the average absolute gene significance (GS) of all the genes involved in the module. GS is measured as the log10 transformation of the p value in the logistic regression between gene expression and clinical traits. Module member (MM) is defined as the correlation between the expression profile and each module feature value. In WGCNA, the module with the highest correlation score was defined as the key module and subjected to further analysis.

Functional Enrichment Analysis

Genes in the key module of “state” indicated a significant correlation with clinical features of MMD occurrence. Hub genes were defined as those with |GS|>0.35 and |MM|>0.75. Hub genes were selected to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. GO and KEGG pathway analyses were conducted to explore their biological functions utilizing the ShinyGO [version 0.77]. GO terms and KEGG pathways were explored with adjusted p < 0.05 and false discovery rate <0.05 set as the cutoff.

Statistical Analysis

The statistical significance of differences between two groups was analyzed using nonparametric test or t tests based on data distribution characteristics. All analyses were conducted using R Studio (v1.4.1717) software, and a p value <0.05 was considered statistically significant. The R packages DESeq2 (v1.32.0), GEOquery (v2.60.0), tidyverse (v1.3.1), CorLevelPlot (v0.99.0), limma (v3.46.0), WGCNA (v1.71), ggplot2 (v3.4.2), etc., were used in this study.

Gene Co-Expression Networks

After excluding 9 outliers (GSM5710974, GSM5710975, GSM5710967, GSM5710985, GSM5710986, GSM5710968, GSM5710969, GSM5710970, and GSM5710972) from GSE122897, the cohort consisted of 14 MCA MMD tissue samples and 9 MCA non-MMD control tissue samples as illustrated in Figure 2a. In WGCNA, we correlated each module with clinical traits in the GSE189993 dataset by calculating the module significance for each module-trait correlation. When 0.85 was used as the correlation coefficient threshold, the soft power threshold (β) was set to 18 as shown in Figure 2b. Twenty-four co-expression modules were constructed by using the average linkage hierarchical clustering algorithm as shown in Figure 2 c, d. Non-clustered genes were gathered in the gray module. The magenta module containing 217 genes was identified as having the strongest association with MMD and thus was selected as the module of interest, shown in Figure 2 d, e. For this module, 139 genes met the criteria (|MM|> 0.75 and |GS|> 0.35) to be labeled as hub genes (online suppl. Table 5). DGEA performed on the same dataset identified a total of 625 DEGs which overlapped with 17 of the 139 hub genes as shown in Figure 2f. Detailed statistics for these 17 genes, including logFC and p values, are shown in Table 1.

Fig. 2.

a Clustering of samples to detect outliers. b Scale-free topology model for finding the soft-thresholding power. c Clustering dendrogram of total genes related to MMD. d Heatmap shows the relationships between different modules and clinical traits. e Scatter plot of module membership versus GS for the magenta module to identify MEs. Scatter plot reported a Pearson correlation coefficient of −0.36 and a p value of 0.000000049. f Venn diagram shows the intersecting genes among differential expression analysis and weighted gene co-expression analysis. ME, module eigengene; p, p value; DEGs, differentially expressed genes.

Fig. 2.

a Clustering of samples to detect outliers. b Scale-free topology model for finding the soft-thresholding power. c Clustering dendrogram of total genes related to MMD. d Heatmap shows the relationships between different modules and clinical traits. e Scatter plot of module membership versus GS for the magenta module to identify MEs. Scatter plot reported a Pearson correlation coefficient of −0.36 and a p value of 0.000000049. f Venn diagram shows the intersecting genes among differential expression analysis and weighted gene co-expression analysis. ME, module eigengene; p, p value; DEGs, differentially expressed genes.

Close modal
Table 1.

List of the 17 hub genes identified by WGCNA which were confirmed as DEGs via DESeq2 analysis and expression data of 2 genes which were further analyzed in an MMD-only analysis stratified by tissue location (MCA vs. STA)

Ensembl gene IDGene symbolFull gene namelogFCp value
Hub gene expression in 12 MCA MMD versus 11 MCA non-MMD controls 
ENSG00000100105 PATZ1 POZ/BTB and AT hook containing zinc finger 1 1.430 0.00016 
ENSG00000100228 RAB36 RAB36, member RAS oncogene family 1.059 0.00046 
ENSG00000120913 PDLIM2 PDZ and LIM domain 2 1.354 0.00015 
ENSG00000136854 STXBP1 Syntaxin binding protein 1 1.105 0.00009 
ENSG00000170893 TRH Thyrotropin-releasing hormone 1.146 0.00115 
ENSG00000171596 NMUR1 Neuromedin U receptor 1 1.150 0.00028 
ENSG00000227702 LINC00111 Long intergenic nonprotein coding RNA 111 1.320 0.00075 
ENSG00000228627 NA Novel transcript 0.917 0.00113 
ENSG00000228659 LINC01201 Long intergenic nonprotein coding RNA 1201 1.142 0.00028 
ENSG00000233605 LINC01851 Long intergenic nonprotein coding RNA 1851 0.941 0.00067 
ENSG00000237728 CAMTA1-AS2 CAMTA1 antisense RNA 2 1.104 0.00118 
ENSG00000258483 LINC02251 Long intergenic nonprotein coding RNA 2251 0.831 0.00057 
ENSG00000258735 LINC00637 Long intergenic nonprotein coding RNA 637 1.213 0.00045 
ENSG00000261190 C16orf97 Chromosome 16 putative open reading frame 97 1.740 0.00045 
ENSG00000269794 NA Zinc finger protein family pseudogene 1.067 0.00096 
ENSG00000279995 NA TEC 1.306 0.00044 
ENSG00000280989 LINC00581 Long intergenic nonprotein coding RNA 581 0.976 0.00070 
Ensembl gene IDGene symbolFull gene namelogFCp value
Hub gene expression in 12 MCA MMD versus 11 MCA non-MMD controls 
ENSG00000100105 PATZ1 POZ/BTB and AT hook containing zinc finger 1 1.430 0.00016 
ENSG00000100228 RAB36 RAB36, member RAS oncogene family 1.059 0.00046 
ENSG00000120913 PDLIM2 PDZ and LIM domain 2 1.354 0.00015 
ENSG00000136854 STXBP1 Syntaxin binding protein 1 1.105 0.00009 
ENSG00000170893 TRH Thyrotropin-releasing hormone 1.146 0.00115 
ENSG00000171596 NMUR1 Neuromedin U receptor 1 1.150 0.00028 
ENSG00000227702 LINC00111 Long intergenic nonprotein coding RNA 111 1.320 0.00075 
ENSG00000228627 NA Novel transcript 0.917 0.00113 
ENSG00000228659 LINC01201 Long intergenic nonprotein coding RNA 1201 1.142 0.00028 
ENSG00000233605 LINC01851 Long intergenic nonprotein coding RNA 1851 0.941 0.00067 
ENSG00000237728 CAMTA1-AS2 CAMTA1 antisense RNA 2 1.104 0.00118 
ENSG00000258483 LINC02251 Long intergenic nonprotein coding RNA 2251 0.831 0.00057 
ENSG00000258735 LINC00637 Long intergenic nonprotein coding RNA 637 1.213 0.00045 
ENSG00000261190 C16orf97 Chromosome 16 putative open reading frame 97 1.740 0.00045 
ENSG00000269794 NA Zinc finger protein family pseudogene 1.067 0.00096 
ENSG00000279995 NA TEC 1.306 0.00044 
ENSG00000280989 LINC00581 Long intergenic nonprotein coding RNA 581 0.976 0.00070 
NLRI pathway-related hub gene expression analysis in 4 MMD MCA versus 4 MMD STA 
ENSG00000170893 TRH Thyrotropin-releasing hormone 1.287 0.82043 
ENSG00000171596 NMUR1 Neuromedin U receptor 1 1.962 0.00053 
NLRI pathway-related hub gene expression analysis in 4 MMD MCA versus 4 MMD STA 
ENSG00000170893 TRH Thyrotropin-releasing hormone 1.287 0.82043 
ENSG00000171596 NMUR1 Neuromedin U receptor 1 1.962 0.00053 

logFC, log fold change; NA, not available; NLRI, neuroactive ligand-receptor interaction.

Functional Enrichment Analysis

In our KEGG analysis which identified only one pathway (the neuroactive ligand-receptor interaction [NLRI]), we found that 6 of the 139 hub genes are involved in the NLRI pathway as shown in online supplementary Figure 1. The GO biological process functional enrichment analysis for the genes in the magenta module predominantly indicated their involvement in behavior and feeding behavior as shown in online supplementary Figure 2. However, when the 17 out of the 139 hub genes validated as DEGs were reanalyzed in the GO biological process functional enrichment, they were predominantly associated with axon target recognition, regulation of the transport of several organic substances like organic anions, organic acids, including dicarboxylic acid across the blood-brain barrier and vessels as shown in Figure 3a. These genes might also have critical roles in establishing or maintaining the transmembrane electrochemical gradient and ion transmembrane transport.

Fig. 3.

a Bubble plot of biological process (BP) of GO analysis on the 17 intersecting hub genes from magenta module. b, c Box and whisker plots on the differential expression of NMUR1 and TRH in MMD samples compared to controls. d Plot on the differential expression of NMUR1 in MMD MCA samples compared to MMD STA.

Fig. 3.

a Bubble plot of biological process (BP) of GO analysis on the 17 intersecting hub genes from magenta module. b, c Box and whisker plots on the differential expression of NMUR1 and TRH in MMD samples compared to controls. d Plot on the differential expression of NMUR1 in MMD MCA samples compared to MMD STA.

Close modal

MMD MCA versus STA-Stratified Subgroup Analysis

The NMUR1 gene underwent further examination to ascertain if its expression was significantly higher in the MCA than in other regions. From the GSE189993 dataset used in our initial analysis (step 1), 4 MMD patients (online supplementary Table 1: patient # 3, 9, 10, and 11) sampled in step 1 had both their MCA and STA sampled which were then added to GSE141025. Unlike the thyrotropin-releasing hormone (TRH) gene, which encodes a peptide hormone with systemic effects, the NMUR1 gene encodes for a receptor, implying localized effects based on expression. Hence, it was expected that the DGEA for the location-specific MMD-only subgroup would confirm no significant difference in TRH gene expression between the two arteries. Moreover, a significant difference in NMUR1 gene expression would intensify the interest in this gene and its potential role in the localized stenosis of the MCA during MMD progression. As hypothesized, our findings backed this presumption (NMUR1: p = 0.00053; TRH: p = 0.82043) which is shown in Table 1 and Figure 3.

Currently, comprehensive research on MMD at the genomic level remains sparse, and the underlying mechanisms for its onset are considered multifactorial. WGCNA while a powerful bioinformatic model generally necessitates a sample size of at least 20, making it challenging to deploy for rare diseases where sampling affected structures proves difficult. Recently, Cao et al. [10] presented findings from the inaugural bioinformatic study using WGCNA on cerebral artery samples. Notably, they utilized the dataset used in the present study (GSE189993) and supplemented it with another cohort (GSE141022) that includes 4 MMD MCAs and 4 non-MMD control MCAs. In total, 40 samples were candidates for inclusion in WGCNA. Upon closer examination, it was identified that the samples from GSE141022 correspond with those in GSE189993, displaying identical attributes such as patient age, sex, sample ID number, and sampling location. This resulted in 20% of the samples being duplicates; however, given that none of the 8 duplicates seem to have been identified as outliers during the assessment for batch effects via principal component analysis, the percentage of duplicates which were involved in WGCNA is expected to be much greater. A final dataset, GSE157628 (consisting of 11 MMD MCA, 6 IA MCA, and 3 MCA epileptic patient samples), was also used to validate the hub genes identified in WGCNA. Concerningly, all these samples also aligned with those already present in GSE189993. These overlapping datasets warrant careful reconsideration, as they indisputably influence the study conclusions.

The molecular functions of the hub genes identified in this study could offer insights into MMD diagnosis and treatment. We identified 139 hub genes in the magenta module associated with MMD onset. KEGG analysis linked these genes with the NLRI pathway, which predominantly comprises neuroreceptor genes, such as dopamine receptor and proto-oncogene. This pathway plays crucial roles in environmental information processing and signaling molecule interaction [11]. Thus, its dysregulation has been linked to various oncological pathologies including lung, renal, prostate cancer [12], and certain neurological diseases including Alzheimer’s disease [13]. To our knowledge, this study is the first to associate neurovascular disorders with the NLRI pathway.

NMUR1

The neuromedin U receptor 1 (NMUR1) predominates in the gastrointestinal tract, influencing mesenteric and portal circulation due to its ability to cause smooth muscle contraction when bound to NMU [14]. This understanding led to the first study to investigate whether NMUR1 had a role in the cardiovascular system in humans [15]. Mitchell et al. identified NMUR1 to be localized on the vascular smooth muscle cells of the saphenous vein, coronary artery, and left ventricle vascular smooth muscle cells via immunochemistry which was further backed by quantitative autoradiography showing highest receptor density in the media of large vessel. Furthermore, vasoconstriction was in agreement with the localization of NMUR1 to vascular smooth muscle cells. The results from Mitchell et al. provide valuable information to the findings in the present study regarding the overexpression of NMUR1 seen in MMD patients. The present study found that MMD patients exhibit a 2.98-fold surge in NMUR1 gene expression in their MCA compared to controls. Furthermore, our MMD-only subgroup analysis indicates a significantly heightened expression in MCA relative to STA by approximately 7.22-fold, suggesting that the MCA in MMD could experience chronic vasoconstriction, potentially leading to progressive stenosis. Interestingly, the NMU ligand, as shown by Mitchell et al. [15], mimics the potency and response of angiotensin 2, a potent vasoconstrictor. Chronic vasoconstriction due to potent hormones like angiotensin-2 is recognized to result in arterial stenosis [16]. Furthermore, histopathology of the main arteries affected by MMD when biopsied occasionally shows infiltration of T cells which contributed to the idea that MMD was also inflammatory mediated [17]. As of last year, the NMUR1 was found to be expressed by most human immune cells with higher levels in type 2 inflammatory cells including type 2 T helpers, type 2 cytotoxic T cells, group-2 innate lymphoid cells, and eosinophils and was significantly upregulated in the activated type 2 cells [18]. These authors concluded that NMU could possibly contribute to pathogenic processes of type 2 immunity-mediated diseases. It is possible that the microsamples of the MCAs from MMD patients used in this study contained a certain degree of immune cells. Thus, with the current dataset it is not feasible to conclude to what degree the increase in expression of NMUR1 was contributed by the MCA smooth muscles or lymphocytes in the samples. Regardless of whichever pathway is ultimately found to be dominant in future studies, both possible outcomes can indeed contribute to arterial stenosis.

In regard to future researchers seeking to externally validate our findings on NMUR1, it would be of great interest to rather validate via using the approach used by Mitchell et al. but on the cerebral arteries. Confirming a high density of the receptor encoded by NMUR1 via immunochemistry while confirming lower densities on other arteries unaffected by the disease could raise great interest on investigating the long-term effects of targeting the receptor of NMUR1 via administering a novel antagonist to observe if arterial stenosis in animal models has a significant effect on the rate of stenotic progression as a potential therapeutic approach to combat the high prevalence of ischemic stroke seen in MMD patients. Additionally, we suggest that future studies should include the analysis of serum for NMU ligand levels alongside NMUR1 expression. This is crucial to discern whether the overexpression of NMUR1 is primarily driven by its interaction with the NMU ligand or if it is a result of both elevated NMU ligand levels and NMUR1 overexpression in MMD patients. Understanding this distinction is key to designing experimental studies that explore how targeting this pathway might influence the progression of stenosis, especially in animal models.

Thyrotropin-Releasing Hormone

The TRH gene encodes for TRH, a tripeptide produced in the hypothalamus that regulates the secretion of thyroid-stimulating hormone from the anterior pituitary, thereby playing a crucial role in the hypothalamic-pituitary-thyroid axis and overall thyroid function. Our analysis revealed that TRH was overexpressed by approximately 2.91-fold in the MCA from MMD patients when compared to control samples. This finding aligns with previous studies that identified hyperthyroidism as a significant risk factor for MMD (OR, 10.936; p < 0.001) [19]. Thyroid hormone is known to target the vascular endothelium and exaggerate vascular reactivity including enhancing arterial vasoconstrictory response to norepinephrine [20]. Such enhanced sympathetic nervous activity might cause abnormal vessel formation and chronic vasoconstriction. Furthermore, thyrotoxicosis is linked to hyperhomocysteinemia [21], a known risk factor for arterial stenosis. Nonetheless, more functional studies are necessary to understand the precise role of TRH in MMD.

Limitations

There are several limitations worth mentioning in the present study. First, the control group of non-MMD MCA samples was not from healthy patients as 6 of the 9 harbored an internal carotid artery aneurysms and the remaining 3 were from patients with epilepsy. This limitation could have impacted the bioinformatic analysis to some extent; however, obtaining MCA samples from healthy individuals with no underlying neurological disorder is not ethically feasible. Another important factor was the age distribution. Even though 87% of the patients were adults, 3 patients were pediatrics (ages 5, 7, and 14). We know that MMD in younger pediatrics often has a more aggressive onset when compared to adults [22]. The underlying mechanism contributing to this observation has never been investigated; however, there is a possibility that differences in gene expression between the two age groups are contributing to the difference in onsets which could impact the analysis of the present study. Given the rarity of MMD, obtaining sufficient pediatric samples to stratify the analysis at present is challenging. Obtaining samples from a biobank comes with limitations as well. Patient comorbidities such as the presence of vasculitis or autoimmune diseases are important factors which can influence bioinformatic analyses conducted on vascular pathologies like MMD. Similarly, to many other bioinformatic studies on MMD and other neurovascular diseases utilizing biobanks [9, 10, 23‒25], such additional patient information was not able to be retrieved. Thus, we recommend future researchers validating our findings consider matching their study subjects with the common risk factors for atherosclerosis, such as hyperlipidemia and diabetes. While the sample size of our study, involving only 23 samples, might be seen as a limitation, it is adequate for conducting WGCNA and is comparatively large relative to other genomic studies on MMD in the literature [3]. The atrial segments of the MCA sampled in this study also raise another limitation worth mentioning. We know that the steno-occlusions seen in MMD are located more proximally (at the terminal ICA and M1 segments). The M4 segments studied in this article are more characterized by collateralization and reduced muscle layers. Thus, sampling more proximally would be ideal for future researchers attempting to improve on our work; however, this also poses the risk for ethical concerns as we rarely perform bypass surgery in the lower MCA segments. Additionally, this study lacks a control group for MCA stenosis arising from alternative etiologies, such as arteriosclerosis, indicating a gap that necessitates further exploration. In conclusion, we encourage future researchers to undertake a comprehensive external validation of the 17 potential biomarkers identified and the role of the NLRI pathway in MMD pathogenesis.

The present study identified 17 potential biomarkers for MMD and is first to link the NLRI pathway with MMD. NMUR1 expression is well known to induce vascular smooth muscle constriction and inflammation which are all factors associated with arterial stenosis. The MMD-only subgroup analysis stratified by artery type further strengthens the potential role of NMUR1 overexpression in MMD given that NMUR1 was significantly overexpressed in the MCA when compared to the STA, highlighting NMUR1 as a key gene of interest in this study. Regarding the hub gene TRH, hyperthyroidism has been identified in a recent comparative meta-analysis as a risk factor for MMD making the present study the first to support these data at the genomic level. Further validation and functional studies are essential to fully understand the specific role of NMUR1 upregulation in MMD and its potential implications.

Given that all data were obtained from a publicly available microarray dataset, this study is exempted from requiring ethics approval. Patient consent was not required as this study was based on publicly available data.

The authors have no conflicts of interest to declare.

This study was not supported by any sponsor or funder.

Conceptualization and design, statistical and bioinformatic analysis, and drafting the manuscript: Samuel D. Pettersson and Shunsuke Koga; interpretation of data: Samuel D. Pettersson, Shunsuke Koga, Shan Ali, Alejandro Enriquez-Marulanda, Philipp Taussky, and Christopher S. Ogilvy; critically revising the manuscript: Shan Ali, Alejandro Enriquez-Marulanda, Philipp Taussky, and Christopher S. Ogilvy; supervising the study: Shunsuke Koga and Christopher S. Ogilvy.

The link to the raw microarray dataset for the MCAs and STAs used can be downloaded from the links https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE189993 and https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE141025, respectively. The R script can be obtained via emailing the corresponding author, and the full gene expression and significance values for each analysis performed in the present study are reported in the online supplementary Tables 3–5.

1.
Burke
GM
,
Burke
AM
,
Sherma
AK
,
Hurley
MC
,
Batjer
HH
,
Bendok
BR
.
Moyamoya disease: a summary
.
Neurosurg Focus
.
2009
;
26
(
4
):
E11
. .
2.
Mertens
R
,
Graupera
M
,
Gerhardt
H
,
Bersano
A
,
Tournier-Lasserve
E
,
Mensah
MA
, et al
.
The genetic basis of moyamoya disease
.
Transl Stroke Res
.
2022
;
13
(
1
):
25
45
. .
3.
Guey
S
,
Tournier-Lasserve
E
,
Hervé
D
,
Kossorotoff
M
.
Moyamoya disease and syndromes: from genetics to clinical management
.
Appl Clin Genet
.
2015
;
8
:
49
68
. .
4.
Tang
Q
,
Li
W
,
Huang
J
,
Wu
Y
,
Ma
C
,
Tu
Y
, et al
.
Single-cell sequencing analysis of peripheral blood in patients with moyamoya disease
.
Orphanet J Rare Dis
.
2023
;
18
:
174
12
. .
5.
Mamiya
T
,
Kanamori
F
,
Yokoyama
K
,
Ota
A
,
Karnan
S
,
Uda
K
, et al
.
Long noncoding RNA profile of the intracranial artery in patients with moyamoya disease
.
J Neurosurg
.
2023
;
138
(
3
):
709
16
. .
6.
Takagi
Y
,
Kikuta
KI
,
Sadamasa
N
,
Nozaki
K
,
Hashimoto
N
.
Caspase-3-dependent apoptosis in middle cerebral arteries in patients with moyamoya disease
.
Neurosurgery
.
2006
;
59
(
4
):
894
901
. .
7.
Hosoda
Y
,
Ikeda
E
,
Hirose
S
.
Histopathological studies on spontaneous occlusion of the circle of Willis (cerebrovascular Moyamoya disease)
.
Clin Neurol Neurosurg
.
1997
;
99
(
Suppl 2
):
S203
8
. .
8.
Takagi
Y
,
Kikuta
KI
,
Nozaki
K
,
Fujimoto
M
,
Hayashi
J
,
Imamura
H
, et al
.
Expression of hypoxia-inducing factor-1 alpha and endoglin in intimal hyperplasia of the middle cerebral artery of patients with Moyamoya disease
.
Neurosurgery
.
2007
;
60
(
2
):
338
45
. .
9.
Zhao
C
,
Ma
Z
,
Shang
J
,
Cui
X
,
Liu
J
,
Shi
R
, et al
.
Bioinformatics analysis reveals potential biomarkers associated with the occurrence of intracranial aneurysms
.
Sci Rep
.
2022
;
12
(
1
):
13282
11
. .
10.
Cao
L
,
Ai
Y
,
Dong
Y
,
Li
D
,
Wang
H
,
Sun
K
, et al
.
Bioinformatics analysis reveals the landscape of immune cell infiltration and novel immune-related biomarkers in moyamoya disease
.
Front Genet
.
2023
;
14
:
1101612
. .
11.
Ren
C
,
Liu
X
,
Zhang
Z
,
Wang
Y
,
Duan
W
,
Li
S
, et al
.
CRISPR/Cas9-mediated efficient targeted mutagenesis in Chardonnay (Vitis vinifera L.)
.
Sci Rep
.
2016
;
6
:
32289
. .
12.
He
Z
,
Tang
F
,
Lu
Z
,
Huang
Y
,
Lei
H
,
Li
Z
, et al
.
Analysis of differentially expressed genes, clinical value and biological pathways in prostate cancer
.
Am J Transl Res
.
2018
;
10
(
5
):
1444
56
.
13.
Kelly
J
,
Moyeed
R
,
Carroll
C
,
Luo
S
,
Li
X
.
Genetic networks in Parkinson’s and Alzheimer’s disease
.
Aging
.
2020
;
12
(
6
):
5221
43
. .
14.
Hedrick
JA
,
Morse
K
,
Shan
L
,
Qiao
X
,
Pang
L
,
Wang
S
, et al
.
Identification of a human gastrointestinal tract and immune system receptor for the peptide neuromedin U
.
Mol Pharmacol
.
2000
;
58
(
4
):
870
5
. .
15.
Mitchell
JD
,
Maguire
JJ
,
Kuc
RE
,
Davenport
AP
.
Expression and vasoconstrictor function of anorexigenic peptides neuromedin U-25 and S in the human cardiovascular system
.
Cardiovasc Res
.
2009
;
81
(
2
):
353
61
. .
16.
Zieman
SJ
,
Melenovsky
V
,
Kass
DA
.
Mechanisms, pathophysiology, and therapy of arterial stiffness
.
Arterioscler Thromb Vasc Biol
.
2005
;
25
(
5
):
932
43
. .
17.
Weinberg
DG
,
Arnaout
OM
,
Rahme
RJ
,
Aoun
SG
,
Batjer
HH
,
Bendok
BR
.
Moyamoya disease: a review of histopathology, biochemistry, and genetics
.
Neurosurg Focus
.
2011
;
30
(
6
):
E20
. .
18.
Ye
Y
,
Luo
J
,
Zeng
N
,
Jiang
S
,
Chen
W
,
Hoyle
RD
, et al
.
Neuromedin U promotes human type 2 immune responses
.
Mucosal Immunol
.
2022
;
15
(
5
):
990
9
. .
19.
Ahn
JH
,
Jeon
JP
,
Kim
JE
,
Ha
EJ
,
Cho
WS
,
Park
YJ
, et al
.
Association of hyperthyroidism and thyroid autoantibodies with moyamoya disease and its stroke event: a population-based case-control study and meta-analysis
.
Neurol Med Chir
.
2018
;
58
(
3
):
116
23
. .
20.
Napoli
R
,
Biondi
B
,
Guardasole
V
,
Matarazzo
M
,
Pardo
F
,
Angelini
V
, et al
.
Impact of hyperthyroidism and its correction on vascular reactivity in humans
.
Circulation
.
2001
;
104
(
25
):
3076
80
. .
21.
Colleran
KM
,
Ratliff
DM
,
Burge
MR
.
Potential association of thyrotoxicosis with vitamin B and folate deficiencies, resulting in risk for hyperhomocysteinemia and subsequent thromboembolic events
.
Endocr Pract
.
2003
;
9
(
4
):
290
5
. .
22.
Pettersson
SD
,
Olofsson
HKL
,
Ali
S
,
Szarek
D
,
Miękisiak
G
,
Ogilvy
CS
.
Risk factors for ischemic stroke after revascularization surgery in patients with moyamoya disease: an age-stratified comparative meta-analysis
.
World Neurosurg
.
2023
;
173
:
146
57.e14
. .
23.
Jiang
Y
,
Leng
JX
,
Lin
Q
,
Zhou
F
.
Epithelial-mesenchymal transition related genes in unruptured aneurysms identified through weighted gene coexpression network analysis
.
Sci Rep
.
2022
;
12
(
1
):
225
. .
24.
Jin
F
,
Duan
C
.
Identification of immune-infiltrated hub genes as potential biomarkers of Moyamoya disease by bioinformatics analysis
.
Orphanet J Rare Dis
.
2022
;
17
(
1
):
80
. .
25.
Chen
Y
,
Tang
M
,
Li
H
,
Liu
H
,
Wang
J
,
Huang
J
.
TGFβ1 as a predictive biomarker for collateral formation within ischemic moyamoya disease
.
Front Neurol
.
2022
;
13
:
899470
. .