Identification of Hub Genes Associated with the Pathogenesis of Intracranial Aneurysm via Integrated Bioinformatics Analysis

Introduction

Intracranial aneurysm (IA) is a cerebrovascular disease that can be caused by several factors, such as inflammation, hemodynamic changes, genetics, smoking history, etc.^1–4 The global prevalence of IA is about 3.2%, and the main risk associated with it is that of rupture and bleeding.^5,6 The available treatment modalities for IA include surgical clipping and various endovascular therapy options (such as endovascular coiling, shunts, and intravascular devices).^6–9 However, there are still controversies about the choice of treatment modalities and the start time of treatment for IA under different conditions (ruptured or unruptured). In addition, the existing main treatment methods are all invasive procedures, which may cause various complications.¹⁰ IA and the subarachnoid hemorrhage (SAH) caused by its rupture have caused a significant burden on the global health.^5,11 Therefore, effective management of unruptured IA (UIA) and treatment of IA still remains a significant clinical challenge.¹⁰ Therefore, detailed studies on the pathogenesis of IA may potentially help to improve the treatment strategies for this devastating disease.

However, so far, no studies have clearly explained the specific mechanisms of aneurysm formation.^12–15 In recent years, with the advancement of genetic technology, several studies have performed gene expression profiling on the vascular wall tissue of IA patients in order to discover the expression of various key genes and changes in pathways involved in the formation of aneurysm. The differentially expressed genes between aneurysm as well as the normal samples were identified in a study by Yu et al¹³ and mainly involved genes regulating immune and inflammatory responses, cell growth, proliferation, differentiation and migration, and intercellular signal transduction. The validation analysis by Wang et al showed that VCAM1, MAGI2, PPP2R2, PPP2R3A genes may be related to the occurrence and development of IA, and these genes were mainly involved in encoding cell adhesion molecules (CAMs) and extracellular regulated protein kinase (ERK)/c-JunN-terminal kinase (JHK) signaling pathways.¹⁶ The team of Wei and coworkers found that CD40, CD40LG, NOS1, STRN, MCM4, MYH11 and miR-125b may play an important role in the pathogenesis of IA.¹⁴ In fact, due to the difficulty in obtaining materials from IA tissues, different data processing methods and platform conditions, the results of screening differential genes with microarray and sequencing data might vary significantly.

In order to overcome the shortcomings of limited or inconsistent results due to applications of the different technology platforms and small sample sizes, robust rank aggregation (RRA) method was used in this study to merge the different gene data sets. This method has been used in many other diseases (such as preeclampsia, acute myocardial infarction, thyroid cancer) to screen for the various differentially expressed genes (DEGs),^17–19 but it is currently used less in the studies related to IA. In this study, the gene expression datasets were downloaded from the Gene Expression Omnibus (GEO) database, and the RRA method was used to screen the significant genes between the IA tissue and the control tissue, and perform a functional enrichment analysis for the significant genes. The cytoHubba plug-in of Cytoscape software was used to identify the various hub genes. The expression of different hub genes was further verified in the sequencing dataset.

Materials and Methods

Microarray Datasets of IA

The gene expression microarray datasets of IA analyzed in this study were acquired from the Gene Expression Omnibus GEO (http://www.ncbi. nlm.nih.gov/geo). We downloaded raw data or series matrix files from GEO, including GSE15629,²⁰ GSE75436²¹ and GSE26969,²² GSE6551.²³

Identification of DEGs in IA

We used the corresponding annotation package in R software (version 4.0.3; 64-bit; https://www.R-project.org/) or annotation document to map the microarray probes to gene symbols, and to calculate the mean value when multiple probes mapped to the same gene symbol. Next, the package of “limma” in R software was used to identify the DEGs. The |log2 fold change (FC)|>1 and P-value < 0.05 were treated as the criteria to screen for the DEGs, and the DEGs were visualized with volcano plots.

RRA Analysis

The RRA method assumes that genes are randomly arranged in each dataset. When a gene ranks high in all datasets, the P value will decrease, and it is more likely to be a significant gene. The differential genes obtained in each dataset were sorted into a list of up-regulated and down-regulated genes. The “Robust Rank Aggregation” package in R software was used to rank different up-regulated and down-regulated genes and to identify the significance of the different genes. The conditions for the screening of the significant genes were set to the fold change >1 and P-value < 0.05.

Functional and Pathway Enrichment Analyses of DEGs

The “clusterProfiler” package in R was used to perform Gene Ontology (GO) functional enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of the significant genes screened out by RRA analysis.

Protein–Protein Interaction (PPI) Network Analysis

The STRING database (http://string-db.org/) was used to establish a PPI network with a confidence of 0.4 to predict the various protein–protein interactions. Thereafter, the results were imported into Cytoscape 3.8.1 to build a network model. Using cytoHubba plug-in of Cytoscape, the top 10 DEGs were selected with a high connectivity in the gene expression network as the hub genes according to the degree algorithm.

Hub Genes Expression in RNA-Seq Dataset

We further verified the hub genes by using data from an RNA-seq dataset (GSE122897,²⁴ including 44 IA, 16 controls). The RNA-seq reading count data was analyzed using DESeq2 function to identify the various DEGs. The genes with P-value <0.05 were considered to be significant.

Results

Information on Included Datasets

The microarray data GSE15629 contains 14 IA samples and 5 control samples. GSE75436 contains 15 IA samples and 15 control samples. GSE26969 has 3 IA samples and 3 control samples. GSE6551 consists of 2 IA samples and 3 control samples. The detailed information of these datasets has been shown in Table 1.

Table 1 Characteristics of the Included Microarray Datasets

Identification of DEGs

In order to ensure homogeneity, the four data sets were evaluated for quality, and the results have been shown in Supplementary Figure 1. According to the set threshold, the “limma” package in R was used to screen the DEGs. The results of the DEGs in each dataset have been represented by volcano plots in Figure 1.

Figure 1 Volcano plots of the four chip datasets. Red points: upregulated genes; Blue points: downregulated genes; Gray points: genes with no significant difference. (A) GSE15629, (B) GSE75436, (C) GSE26969, (D) GSE6551.

Results of RRA Integrated Analysis

According to the threshold conditions that were set, a total of 136 significant genes were identified by RRA integrated analysis, which included 45 up-regulated genes and 91 down-regulated genes and described in Supplementary Table 1. Then, the top ten significant genes among the up-regulated and down-regulated genes were selected to draw a heatmap, as shown in Figure 2.

Figure 2 Heatmap of the top ten up-regulated and down-regulated genes in the RRA analysis. Different colors represent the level of gene expression. Red: high expression. Blue: low expression.

Functional Enrichment Analysis of Significantly Different Genes

We performed gene enrichment analysis on the significant genes. The results showed that the degree of bio-enrichment of extracellular matrix (ECM) structural constituent (GO: 0005201, P-value = 9.02E-14) was the most significant, and the all enrichment pathways of GO are listed in Table 2. In the KEGG enrichment analysis, it was found that ECM-receptor interaction (hsa04512; P-value =0.0008) and protein digestion and absorption were significantly enriched (Table 3). Figure 3 depicts the proportion of the significant genes and the degree of enrichment in the GO and KEGG enriched pathways.

Table 2 GO Analysis of Integrated DEGs

Table 3 KEGG Analysis of Integrated DEGs

Figure 3 Bubble chart of GO and KEGG enrichment pathway. (A) GO enrichment pathway; (B) KEGG enrichment pathway.

PPI Network Analysis and Identification of Hub Genes

In order to understand the possible interactions between the proteins encoded by the significant genes in IA, a PPI network was constructed using the STRING online database, and then the Cytoscape was used to visualize the network and identify the different hub genes. Sixty-seven of 136 significantly different genes were used for the construction of the PPI network. The results of network analysis performed by Cytoscape software is shown in Figure 4A. The top ten hub genes screened by cytoHubba were VCAN, COL1A1, COL11A1, COL5A1, COL5A2, TNC, POSTN, THBS2, CHRDL1 and CDH2. The connection network between them is shown in Figure 4B.

Figure 4 PPI network and hub genes identified by cytoHubba plug-in of Cytoscape. (A) PPI network. (B) hub genes.

Analysis of the Expression Level of the Different Hub Genes

We validated the expression level of different hub genes using sequencing dataset GSE122897, and found that eight out of the 10 hub genes had significant differences in the verification dataset. Only TNC and CHRDL1 genes did not show significant differences in the validation dataset (P-value > 0.05). The difference in the expression level of hub genes between IA and the control group is shown in Figure 5.

Figure 5 The validation of the expression of the hub genes in dataset GSE122897. **P-value < 0.01; ***P-value < 0.001; ****P-value < 0.0001; ns indicates non-significant.

Discussion

Four datasets from GEO were found to be downloaded in this study, and RRA was used to perform integration analysis on them to screen for the presence of the significant genes. A total of 136 significant DEGs were screened from the four datasets, including 45 up-regulated genes and 91 down-regulated genes. The top ten hub genes screened by cytoHubba were VCAN, COL1A1, COL11A1, COL5A1, COL5A2, TNC, POSTN, THBS2, CHRDL1 and CDH2. In the sequencing dataset, 8 hub genes were verified to have differences, including VCAN, COL1A1, COL11A1, COL5A1, COL5A2, POSTN, THBS2 and CDH2.

Apoptosis and migration of the vascular smooth muscles, extracellular matrix remodeling of intracranial blood vessels, and inflammatory reactions have been widely accepted to actively participate in IA formation.^1,2,25,26

Moreover, enrichment analysis showed that extracellular matrix structural constituent and ECM-receptor interaction may be involved in the formation of IA. At the same time, it was found that the various genes encoding multiple types of collagen were involved in the structural composition of ECM. Collagen is the most important fibrin outside the cell, and it is also the skeleton of the ECM. Normal extracellular matrix components can maintain the elasticity of the arterial wall. A number of studies have found that when the synthesis of type III collagen decreases and the proportion of type I collagen increases, the blood vessel wall can undergo changes, thereby causing hemangiomas or rupture.²⁷ Yu et al found that the type I collagen gene (COL1A1) was involved in the destruction of the cellular matrix, and COL11A1 was the most significant up-regulated gene among all genes analyzed.¹³ A recent study showed that COL1A1 and COL5A2 were significantly up-regulated in the aneurysm group and could effectively regulate phosphoinositide 3-kinase (PI3K)/Akt signaling pathway.²⁸ Similar to the results of previous studies, we found that the five genes involved in encoding collagen (COL5A2, COL1A1, COL11A1, COL5A1) were all up-regulated in the intracranial aneurysm group, and may play an important role in maintaining normal extracellular matrix components.

The genes affecting ECM structural constituent also include VCAN, POSTN and THBS2. VCAN is a chondroitin sulfate proteoglycan, which has five different subtypes (V0, V1, V2, V3, V4).^29–31 VCAN can exist in a variety of the tissues and is an important component of the ECM.^29,31,32 It plays a pivotal part in regulating cell proliferation, differentiation, adhesion and migration, and can promote angiogenesis and inflammatory tumor microenvironment formation.^30,33 It has been confirmed that VCAN is over-expressed in the tumor stroma and the cancer cells of a variety of malignancies, such as breast cancer,³⁴ endometrial cancer,³⁵ testicular germ cell tumor,³⁶ bladder carcinoma,³⁷ gastric cancer.³⁸ The study by Shen et al³⁸ found that a high expression of VCAN in gastric cancer tissue was associated with a larger depth of tumor invasion and a poor prognosis. Zhang et al³⁷ also confirmed similar observations. They analyzed the lesion tissues of 417 bladder cancer patients and found that the overexpression of VCAN was related to the number of tumors, depth of invasion, lymph node metastasis, distant metastasis, and histological grade. When VCAN was silenced, it was found to effectively inhibit the growth of tumors in situ, and when combined with endostatin, it was more effective to prolong the survival time of mice.³⁹ The above studies have clearly established that an abnormal expression of VCAN in the tumors can significantly promote the tumor growth and metastasis, which is usually related to the poor prognosis of tumors. VCAN is not only related to the development and prognosis of malignant tumor diseases but also can contribute to the development of vascular diseases. VCAN affects the composition of the ECM and controls the formation of elastic fiber fibrils, which is of great significance for the remodeling of the ECM in the vascular diseases, and also plays a key role in atherosclerosis and restenosis.^40,41 ECM remodeling plays an important role in the formation of aneurysms, and the reduction of ECM is a significant feature found in IA. A previous study⁴² reported that a similar phenotype displayed by the V2 subtype of VCAN gene may help to inhibit the growth of smooth muscle and other vascular entities, thus leading to IA. Thus, VCAN gene might be an important candidate gene involved in the pathogenesis of IA. Additionally, similar to the results of previous studies, we also found that VCAN was significantly up-regulated in diseased IA and was involved in the composition of ECM.

The POSTN gene is located on the long arm of human chromosome 13 and can encode POSTN.⁴³ POSTN is a protein secreted by ECM, which is involved in the regulation of cell adhesion, differentiation, blood vessel formation and calcification, tissue fibrosis as well as other pathophysiological processes.^44–47 Some studies have found that POSTN is highly expressed in the tissues or the serum in a variety of diseases (such as pulmonary fibrosis, cartilage-like tumors, head and neck cancer), and has been reported to be associated with poor prognosis.^48–50 Therefore, POSTN may serve as a potential prognostic biomarker and treatment target for various diseases. The findings of Luo et al also established this point of view.⁵¹ They found that the serum POSTN concentration in patients with SAH increased, and it was related to the clinical severity and poor prognosis. At present, the expression of POSTN in various diseases is not fully understood, and the mechanisms involved in regulating its expression in the different diseases might also be different. Studies have found that POSTN can participate in the regulation of cell growth, tumorigenesis and drug response through activating the PI3K/Akt pathway.^50,52 Researchers found that when POSTN expression was reduced, it could significantly inhibit the formation of colorectal tumors in mice.⁴⁷ They proposed that POSTN can activate integrin-focal adhesion kinase (FAK)-Src kinase by binding to integrin, thereby activating YAK/TA2 to cause the release of interleukin-6 (IL-6), which ultimately leads to an inflammatory response. A previous study has confirmed that inflammation can serve as an important mechanism for the formation of hemangioma, and proteomics analysis has found that the deposition of POSTN in the extracellular matrix of abdominal aortic aneurysm tissue was increased.⁵³ The results reported by Yamashita et al also showed that POSTN was highly expressed in the cell infiltration and elastin degradation area of human abdominal aortic aneurysm tissue and was upregulated in the stage of abdominal aortic aneurysm expansion caused by continuous inflammation in the mice.⁵⁴ Therefore, it could be speculated that POSTN is also related to the occurrence and development of aneurysms. In our analysis, it was found that the expression of POSTN in IA tissue was significantly higher than that in the normal control group, and it was involved in the composition of ECM, which indicated that POSTN may be involved in the formation and development of IA.

THBS2 is a multifunctional stromal cell protein encoded by the THBS2 gene and belongs to the thrombin-sensitive protein family.⁵⁵ It has been reported to be involved in the regulation of angiogenesis, proliferation, apoptosis, tumor migration, autophagy and transforming growth factor (TGF) β activation.^56–59 In arteries, THBS2 is mainly secreted by the smooth muscle cells and partly by endothelial cells. A large number of studies have shown that THBS2 may be involved in the formation and development of aneurysms. For instance, a previous study has suggested that THBS2 gene polymorphism is related to thoracic aortic aneurysms, but the specific mechanism is not clear.⁶⁰ Recent studies have found that the expression of THBS2 in the aortic samples of aortic dissection (AD) patients and the concentration of THBS2 in the circulation were significantly higher than those in the control group, and the levels of THBS2 in the circulation of AD patients could positively correlate with the levels of various inflammatory factors.⁵⁶ Similarly, the circulating concentration of THBS2 in patients with abdominal aortic aneurysm was also increased, which was independently related to the risk of cardiovascular death in patients.⁶¹ Another study also confirmed that THBS2 expression was increased in abdominal aortic aneurysm tissue, and further found that THBS2 could act as a potential target of miR-195-5p, which may be involved in the formation of abdominal aortic aneurysm.⁵⁵ We found that THBS2 was highly expressed in the tissues of IA. This is similar to the results of the various previous studies, which prove that THBS2 may be involved in the formation of aneurysms, but the specific mechanism needs to be further studied.

CDH2 is a transmembrane protein that plays an important role in the cell adhesion and is encoded by the CDH2 gene.⁶² Apoptosis of the vascular smooth muscle cells can function as an important regulator of aneurysm formation. Studies have found that cadherin can regulate the behavior of the vascular smooth muscle cells through the mediation of matrix metalloproteinases.⁶³ It has been reported that the protein level of full-length nitro-cadherin in the human abdominal aortic aneurysm samples was significantly lower, but more nitro-cadherin fragments were observed.⁶³ Interestingly, they found that the combination of matrix metalloproteinase-7 deletion and overexpression of EC4-Fc (soluble N-cadherin and a smaller soluble N-cadherin truncation) significantly increased the severity of abdominal aortic aneurysms. In contrary, other study has reported that the expression of CDH2 in patients with thoracic aortic aneurysms was significantly increased.⁶⁴ After testing the cerebrospinal fluid of SAH patients with the ruptured anterior circulation aneurysm, Takase et al found that the level of soluble vascular endothelial cadherin in the cerebrospinal fluid was significantly increased, and the clinical severity of this index could positively correlate with the level of cerebrospinal fluid tumor necrosis factor-α.⁶⁵ We found that CDH2 expression was increased in the tissue samples from IA patients, which indicates that CDH2 may be involved in the formation of IA, but its involvement in specific pathways requires further analysis.

In summary, our results found that various genes differentially expressed in IA were mainly involved in the ECM structural constituent and ECM-receptor interactions. The hub genes involved in the various physiological and pathological processes, such as inflammation, angiogenesis, cell growth, proliferation, etc., are closely related to aneurysm diseases and may serve as a potential therapeutic target and new biomarker for IA. This study has certain limitations, the sample size was still small, and more experiments need to be performed to verify the expression changes of the genes and the specific mechanisms involved in the development of IA.