Potential Role of the Chemotaxis System in Formation and Progression of Intracranial Aneurysms Through Weighted Gene Co-Expression Network Analysis

Introduction

Cardiovascular disease is still the leading cause of death worldwide. Among cardiovascular diseases, aneurysm is not as widely known and has been studied less well than coronary heart disease, but the serious consequences of aneurysm cannot be underestimated. An aneurysm is an abnormal bulge or swelling in an arterial blood vessel. The underlying disease progression involves weakening of the arterial wall, leading to the risk of progressive expansion or even rupture of the aneurysm.^1,2 Intracranial aneurysm (IA) develops in the cerebral artery and is characterized by local dilation of the cerebral artery.³ Saccular aneurysm is the most common, and is the main cause of spontaneous subarachnoid hemorrhage (SAH).⁴ Multifactorial pathophysiological changes after SAH usually lead to irreversible brain damage, and cause a high incidence of medical complications and mortality.^5–7 Individuals affected by SAH are younger than those affected by other forms of stroke, with an average age of 50 years,^8,9 and survivors often experience dysfunction and/or cognitive impairment for a long time.^10,11

Despite numerous studies over the past several decades, the pathophysiological mechanisms underlying the formation, development, and rupture of aneurysms have yet to be fully elucidated, hindering optimal treatment.^12–14 The pathological features of IA are considered to be due to the destruction of extracellular matrix (ECM), leading to damage to the structural integrity of the arteries.^15,16 Inflammation is considered to be a core response throughout the occurrence and development of disease.^17,18 The complex inflammatory response involves a series of pathological processes, including matrix metalloproteinase (MMP) activation, vascular smooth muscle cell (VSMC) apoptosis, macrophage infiltration, and activation of intracellular oxidative stress.^19–22 Macrophages infiltrate the arterial wall due to the actions of inflammatory chemokines.²³ Macrophages express and release MMPs, destroying the ECM of the arterial wall, which in turn leads to the recruitment of other inflammatory cells, aggravating the degeneration and weakening of the arterial wall, and ultimately leading to the formation and growth of aneurysms.^17,23,24 These data indicate the critical role of immune and inflammatory processes in the formation and development of IA.

With advances in microarrays, high-throughput sequencing, and bioinformatics analysis, the transcriptome changes and molecular mechanisms of many diseases can be readily determined. Through in-depth analysis of differentially expressed genes (DEGs) between different samples, such as using GO (gene ontology) and KEGG (Kyoto Gene and Genomic Encyclopedia) to analyze the functions or pathways enriched by DEGs, using protein–protein interactions (PPI) to predict the interaction between proteins, using CIBERSORT to study the degree of immune cell infiltration, and its like, which enable us to a more comprehensive and profound understanding of the disease.^25–27 By constructing a scale-free weighted network, weighted gene co-expression network analysis (WGCNA) can be used to study biologically significant gene sets related to sample characteristics, and explore the core genes of internal modules that are closely related in co-expression modules.^28,29 In this study, the microarray data of IA patients and normal arteries from public databases were examined by WGCNA to identify important modules closely related to IA, which may provide a better understanding of the molecular mechanisms underlying IA.

Methods

Data Collection and Preprocessing

This study was performed using transcriptome data of IA patients available from a public database (Gene Expression Omnibus, [GEO], https://www.ncbi.nlm.nih.gov/geo/). Datasets related to IA—GSE13353, GSE54083, and GSE75436—were included in the analysis and downloaded from the GEO database. The three datasets contained a total of 72 samples, including 25 normal samples and 47 IA samples. The raw data were processed by the RMA (Robust Multi-Array Average) method in the affy (version 1.70.0) R package.³⁰ The missing values were filled using the impute (version 1.66.0) R package. After merging the three datasets, the batch effects were removed by limma (version 3.48.1) and sva (version 3.40.0) R packages.^31,32 Additionally, approval for this study was obtained from the Ethics Committee of the First Affiliated Hospital of Nanchang University. The workflow of the present study was in Figure 1.

Figure 1 Work flow chart of this study.

Differentially Expressed Gene Analysis

After the above processing of the datasets, we used the limma (version 3.48.1) R package to filter the DEGs between the IA and normal samples.³¹ The threshold values were set as adjusted, p < 0.05 and |log₂FC| > 1.5. Thermogram analysis of DEGs was conducted using the pheatmap (version 1.0.12) R package.³³

Functional and Pathway Enrichment Analyses and Principal Component Analysis

GO annotation contains three categories, ie, biological process (BP), cell component (CC), and molecular function (MF), which can be used to identify the biological attributes of genes and gene sets for all organisms.³⁴ The KEGG is a comprehensive database integrating information on genome, chemistry, and system function.³⁵ The KEGG database associates the gene catalogs obtained from completely sequenced genomes with the system functions at the cell, species, and ecosystem levels. To explore the biological functions of the candidate gene modules, the online tool, WebGestalt (version 6.8, http://www.webgestalt.org/), was used for GO and KEGG enrichment analyses, for the genome protein-coding was used as reference set.³⁶ In the analyses, p < 0.05 was taken to indicate significant enrichment of DEGs. The results of the analyses were downloaded from the web and the figures of enrichment analyses were drawn using the ggplot2 (version 3.3.5) and dplyr (version 1.0.7) R packages.³⁷ Principal component analysis (PCA) is a multivariate statistical method to examine the correlations between multiple variables.³⁸ By deriving a few principal components from the original variables, they retained as much information about the original variables as possible and were not related to each other. PCA was conducted on the all genes set and DEGs subset to identify differences in the transcriptome patterns between the IA and normal samples.

Co-Expression Network Construction by Weighted Gene Co-Expression Network Analysis

The WGCNA (version 1.70–3) R package was used to construct the weighted gene co-expression network.²⁹ Genes with the top 25% variance were selected for the subsequent analysis step. The soft power threshold β was determined by the “pickSoftThreshold” function and the screening criterion of β was an R² of scale-free fit > 0.85. Topological overlap measure (TOM) matrices were obtained from adjacency matrices to detect gene modules of the co-expression networks. Based on the TOM matrices, the average linkage hierarchical clustering method was used to construct clustering dendrograms with a minimum module size of 50. Similar gene modules were merged, with a threshold of 0.25. The expression profiles of each module were summarized using the eigengene (ME) module and the correlation between the module and clinical status was calculated. The Pearson’s correlation and Student’s asymptotic p-value were calculated by the relative functions in the WGCNA (version 1.70–3) R package.²⁹ Key modules will be detected and further analyzed for GO, KEGG and protein–protein interactions (PPI).

Hub Gene Identification and Receiver Operating Characteristic Curve Analysis

The online tool, String (version 11.0, https://string-db.org/), which is designed for predicting protein–protein interactions (PPI) based on the evidence of experiments and database information, was used to construct a PPI network of all the DEGs and DEGs overlapped with the WGCNA modules.³⁹ Cytoscape software (version 3.8.2) was used for further analysis of PPI networks and hub gene selection.⁴⁰ The top 10 hub genes in DEGs and each module were selected with the maximal clique centrality (MCC) method using the cytoHubba plugin in Cytoscape.⁴¹ The genes with the top 10 MCC values were considered hub genes and visualized together with their connected nodes. Hub genes were used to conduct receiver operating characteristic (ROC) curve analysis with GraphPad software (version 8.0) to show the potential diagnostic value of these genes. Hub genes were validated using the GSE122897 dataset.

Evaluation of Aneurysm-Infiltrating Immune Cells

CIBERSORT is a deconvolution algorithm based on gene expression, which can accurately quantify the relative levels of different immune cell types in a complex gene expression mixture through a set of barcode gene expression values.⁴² The analysis of infiltrating immune cells on the genes set was carried out using R (version 4.1.0) and the signature genes of 22 types of infiltrating immune cells (LM22 (22 immune cell types)) and CIBERSORT software source code in R were downloaded from CIBERSORT (https://cibersort.stanford.edu/). A filter to remove samples with p > 0.05 from the data was deployed to ensure accurate prediction of the samples. After data processing, three normal and 24 IA samples were included in the follow-up analysis.

Statistical Analysis

Data preprocessing, DEG screening, WGCNA, and immune cell infiltration analysis were performed in R (version 4.1.0). Functional annotation analysis was performed by WebGestalt (version 6.8, http://www.webgestalt.org/). PPI network analyses were conducted using String (version 11.0, https://string-db.org/). Hub genes were mined using MCC with CytoHubba in Cytoscape (version 3.8.2). The ROC curves were drawn with GraphPad (version 8.0). In all analyses, p < 0.05 was taken to indicate statistical significance.

Results

DEGs and Related Pathways

According to the adjusted p < 0.05 and |log₂FC| > 1.5, we screened 266 DEGs, of which 162 were upregulated and 104 were downregulated in the IA group compared to the normal group (Figure 2A and B). The top 10 DEGs based on fold change are shown in Table 1 (information for all DEGs is listed in Table S1). To investigate the biological functions of these DEGs, we performed enrichment analysis on GO-BP (Figure 2C) and KEGG (Figure 2D) gene sets (GO-BP and KEGG IDs and further detailed information are provided in Table S2, GO-CC and GO-MF are provided in Figure S4A and B). The results of GO-BP enrichment analysis indicated that inflammatory response, defense response, immune response, and leukocyte activation were significantly enriched, suggesting that immune processes may play essential roles in IA formation and development. Correspondingly, KEGG pathway enrichment analysis suggested that the chemokine signaling pathway, Toll-like receptor signaling pathway, and complement and coagulation cascades may participate in the IA process.

Table 1 DEGs with Top-10 Fold-Change Between Normal and IA Samples

Figure 2 Identification and analysis of differentially expressed genes (DEGs). (A) Heatmap for DEGs. (B) Volcano plot for DEGs. (C) Gene set enrichment analysis of Gene Ontology biological process (GO-BP) among all the DEGs. (D) Gene set enrichment analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) among all the DEGs.

Principal Component Analysis

PCA was performed on the overall genes and DEGs to determine the expression patterns in IA compared to normal samples. According to PCA on all genes, PC1 and PC2 accounted for 52.1% and 7.7% of gene expression value variation, respectively (Figure S1A, Table S3). The variation proportions of PC1 and PC2 were 44.9% and 20.2%, respectively, in PCA of the DEGs (Figure S1B, Table S3). In addition, the PCA visualization graphs illustrated that there was indeed a difference between the IA and normal groups (Figures 3A, B, S1C and D).

Figure 3 Principal Component Analysis (PCA) on the DEGs. (A) 2D-PCA plot of DEGs. (B) 3D-PCA graph of DEGs.

Weighted Gene Co-Expression Network Construction

WGCNA was performed on 25 normal samples and 47 IA samples. The sample clustering dendrograms are shown in Figure 4A and all samples were located in the clusters and passed the cutoff thresholds. A scale-free network was constructed with a soft threshold of 7 and a correlation coefficient threshold of 0.85 (Figure 4B). The gene modules were detected based on the TOM, and five modules were finally detected (Figure 4C and D). Three modules (red, pink, and blue) were markedly correlated with IA (Figure 4E) (detailed information of genes in each module are provided in Table S4). The genes within blue and pink modules were significantly correlated with IA (Figure 4F–H). Therefore, the blue (1789 genes) and pink (486 genes) modules were defined as the key modules that were highly correlated with IA.

Figure 4 Construction of weighted gene co-expression network on IA and normal conditions. (A) Detection of outliers with sample clustering. All samples are located in the clusters and pass the cutoff thresholds. (B) Soft-threshold power analysis implemented to obtain the scale-free fit index of the network topology. (C) Identification of gene modules. (D) Heatmap of the topological overlap matrix (TOM) of genes selected for WGCNA. (E) Correlation heat map of gene modules and phenotypes. (F) Scatter plot of gene and module correlation within red module. (G) Scatter plot of gene and module correlation within blue module. (H) Scatter plot of gene and module correlation within pink module.

Enrichment Analysis of Key Modules

To explore the physiological and pathological pathways related to IA, the DEGs were intersected with the genes in the blue and pink modules, respectively. GO-BP and KEGG pathway enrichment analysis were performed. The results showed enrichment of immune response, inflammatory response, leukocyte activation, the chemokine signaling pathway, and Toll-like receptor signaling pathway in the blue module (Figure 5A and B, detailed information is shown in Table S5, GO-CC and GO-MF are provided in Figure S4C and D). Enrichment analysis in the pink module suggested that inflammatory response, defense response, myeloid leukocyte cytokine production, cytokine production involved in immune response, regulation of mast cell cytokine production, regulation of inflammatory response, positive regulation of immune system process, leukocyte-mediated immunity, mast cell cytokine production, regulation of interleukin-10 production, and complement and coagulation cascades may be related to IA (Figure 5C and D, detailed information is shown in Table S6, GO-CC and GO-MF are provided in Figure S4E and F).

Figure 5 Gene set enrichment analysis of GO-BP and KEGG ontologies and pathways among DEGs in blue and pink modules. (A and B) GO-BP and KEGG of DEGs in the blue module. (C and D) GO-BP and KEGG of DEGs in the pink module.

Identification of the Hub Genes

All DEGs and DEGs in the blue and pink modules were input into the String database to construct PPI networks (Figure S2). The results were downloaded and input into Cytoscape to identify the hub genes using the MCC algorithm in CytoHubba, and the genes with the top 10 MCC values were considered hub genes and visualized together with their connected nodes. Among the hub genes of DEGs and the blue module, CCR5 and CCL20 overlapped (Figure 6A and B), and FPR3 and CCL4 overlapped in the DEGs and the pink module (Figure 6A and C). Therefore, CCR5, CCL20, CCL4, and FPR3 were selected as key genes.

Figure 6 Hub genes identification by the CytoHubba algorithm. (A) The top 10 hub genes ranked by the MCC algorithm and their neighbors among all the DEGs. (B) The top 10 hub genes ranked by the MCC algorithm and their neighbors in the blue module. (C) The top 10 hub genes ranked by the MCC algorithm and their neighbors in the pink module. Red-orange-yellow represents key genes, and the darker the color, the stronger the association with other genes.

Different Immune Infiltrative Patterns Between IA and Normal Conditions

Gene set enrichment analysis of GO-BP and KEGG ontology and pathways indicated that pathways related to the immune process were significantly enriched. Therefore, CIBERSORT was used to evaluate the patterns of immune infiltration between IA and normal conditions. As mentioned in the Methods section, three normal samples and 24 IA samples were included in the CIBERSORT analysis, and the relative proportions of 22 immune cells are shown in Figure 7A and B. The proportion of macrophages was higher in normal samples, while M0 and M2 macrophages showed the opposite trend. The proportion of neutrophils was higher in normal samples and that of activated mast cells was higher in IA samples (Figure 7C).

Figure 7 Immune infiltrative patterns in IA and normal conditions. (A and B) Histogram of 22 types of immune cells in each IA and normal tissue. (C) Violin chart showing the differences in infiltrating immune cells between the two groups.

Validation of Key Genes and ROC Analysis

Key genes were validated using the GSE122897 dataset, and the results indicated significant differences in expression levels of all key genes between the normal and IA samples (Figure 8A–D). Finally, ROC analysis was conducted to explore potential biomarkers. The ROC curves for the key genes are shown in Figure 8E–H. The top 10 hub genes of all DEGs were also analyzed by ROC analysis, and the areas under the curves (AUCs) for all of these genes were > 0.7 with p < 0.05. CCL21, NPY1R, and ADORA3 had the highest AUCs (AUCs = 0.907, 0.898, and 0.866, respectively) (Figure S3).

Figure 8 ROC analysis and validation of key genes. (A–D) Validation of CCR5, CCL4, CCL20 and FPR3. (E–H) ROC analysis of CCR5, CCL4, CCL20 and FPR3. *P<0.05; **P<0.01; ***P<0.001.

Discussion

In this study, 266 DEGs were identified, including 162 that were upregulated and 104 that were downregulated in the IA samples. Consistent with previous studies, enrichment analysis of DEGs on GO-BP and KEGG suggested that inflammation and immune-related pathways may be most closely related to IA. To identify the basic modules closely related to IA, we performed WGCNA on IA and normal samples together and obtained two modules related to IA. In addition, functional enrichment analysis of intersecting genes between modules and DEGs also suggested that inflammation and immune processes may play important roles in the development and progression of IA. The hub genes in DEGs and the two modules were identified separately. Among the hub genes, CCR5, CCL4, CCL20, and FPR3 overlapped between the DEGs and the modules, and were selected as key genes. As the key genes were closely related to inflammation and immune processes, we performed immune infiltration analysis in IA and normal samples, and found that the proportions of M2 macrophages and activated mast cells were significantly higher in IA samples. Finally, all of these key genes were validated by another GEO dataset, and ROC analysis indicated that the top ten DEGs have potential diagnostic utility and are potential biomarkers.

A number of studies in human and animal models have confirmed immune cell infiltration dominated by macrophages in IA tissues, which is consistent with the obvious infiltration of macrophages in IA samples observed in the present study.

Chemokines are proinflammatory cytokines, which are the main mediators of macrophage chemotaxis and play a pivotal role in inflammation by activating cell chemotaxis.⁴³ Chemokines can chemoattract inflammatory cells, including macrophages, neutrophils, lymphocytes, fibroblasts, and other inflammatory cells to infiltrate sites of inflammation and exert an inflammatory effect.⁴⁴ In addition, several chemokines may play different roles in the formation of aneurysms. Abdominal aortic aneurysm (AAA) and IA have similar risk factors, including age, smoking, hypertension, and familial tendencies, and are similar in many of their pathogenic mechanisms and pathological features, including inflammation, oxidative stress, endothelial dysfunction, vascular smooth muscle cell apoptosis, leukocyte infiltration, and MMP-mediated degradation of the ECM leading to weakening of the arterial wall.^9,17 CCL2, an effective chemokine for monocytes and macrophages, shows significantly elevated expression in IA, and plays an important role in the formation and development of IA.^45–47 Iida et al reported that MKEY, a peptide CXCL4-CCL5 heterodimer inhibitor, suppressed the formation and progression of AAA in mice.⁴⁸ CCR5, a chemokine G protein-coupled receptor (GPCR) responsible for some of the functions of the chemokines CCL3, CCL4, and CCL5, has mainly been detected in macrophages, T lymphocytes, and dendritic cells.^49,50 MacTaggart et al reported that deficiency of CCR5 had no significant effect on the formation or development of AAA,⁵¹ while Ishida et al reported that knockout of CCR5 increased the formation of aneurysms in mice, suggesting a preventive effect of CCR5 on aneurysms.⁵² Another study showed that inhibiting the CCL5/CCR5 axis could reduce macrophage infiltration while suppressing AAA formation and progression, suggesting that CCR5 may have a positive effect on aneurysms.⁴⁸ The specific role of CCR5 in aneurysms has not yet been definitively elucidated. As the upstream molecule of CCR5, CCL4 is highly expressed in AAA, but further studies are required to determine whether it exerts a role through CRR5 and whether the CCL4/CCR5 axis plays a role in inhibiting or promoting aneurysms.⁵³ CCL20 is expressed in various human tissues and immune cells, mainly observed in the lymph nodes, lungs, and liver, and is produced by cells related to inflammation and autoimmune reactions, such as macrophages, neutrophils, and natural killer (NK) cells, Th17 cells, and B cells.^54,55 CCL20 and its specific receptor, CCR6, have been shown to play important roles in autoimmune diseases (rheumatoid arthritis, psoriasis, etc) and tumor diseases (hepatocellular carcinoma, pancreatic cancer, etc).^56–58 Co-culture with M2 macrophages was shown to enhance the invasiveness of pancreatic cancer cells, while CCR6 deficiency significantly reduced the invasiveness of the cells, and addition of CCL20 enhanced the epithelial–mesenchymal transition (EMT) of these cells.⁵⁸ CCL20 expression is also increased in AAA, but further studies are required to determine its role in disease progression.⁵⁹

Formyl peptide receptor (FPR) belongs to the GPCR family, and three FPRs have been identified, to date, in humans: FPR1–FPR3.⁶⁰ FPRs are mainly expressed in several types of innate immune cells, including neutrophils and monocytes/macrophages.⁶¹ The activation of FPR1 and FPR2 by chemotactic agonists was shown to trigger a series of signaling events, leading to myeloid cell migration, mediator release, increased phagocytosis, and gene transcription, all of which mediate inflammation and immunity.^62,63 Although the overall function of FPR3 remains unclear, recent studies have shown that activation of FPR3 can cause the chemotactic migration of neutrophils and improve the survival rate of murine sepsis models, suggesting that FPR3 may play an important role in the inflammatory response.⁶⁴ However, there have been few studies on the relationship between FPR3 and aneurysms. Further studies are therefore required to determine whether FPR3 is involved in the formation and progression of aneurysms.

Mast cells are the main participants in the inflammatory response and, together with other inflammatory cells, form a specific tissue response to infection and injury. By releasing various mediators, including tryptase, chymotrypsin, cathepsin, and interleukins, activated mast cells can affect vascular inflammation and remodeling.^65–67 Zhang et al reported that mast cell tryptase deficiency could attenuate mouse AAA formation.⁶⁸ Mesenchymal stem cells can suppress the progression of intracranial aneurysms by inhibiting the activation of mast cells.^69,70 Furukawa et al suggested that mast cells play a key role in promoting aneurysm rupture but not formation.⁷¹ Although mast cells play an important role in the progression of aneurysms, further studies are required to determine the specific underlying molecular mechanisms and whether chemokines are involved in the regulation of mast cells.

This study had several limitations. First, due to the nature of the secondary analysis, it was difficult to assess the reliability of the original samples. Although we included the data of three data sets that were normalized to remove batch effects, the origins of the samples and the processing that had been done previously could not be fully standardized, which would affect subsequent data analysis. For easier filtering of key genes, we only performed follow-up analysis between DEGs in different modules. Although we attempted to avoid this limitation by defining a smaller threshold, to include more candidate genes, we may have lost some information contained in genes that were not significantly differentially expressed. After estimating infiltrating immune cells of aneurysm by CIBERSORT, only three normal and 24 IA samples were included in the follow-up analysis, and this small sample size may have influenced the results of the analysis. Finally, all analyses were based on bioinformatics algorithms, and further experimental and clinical studies are necessary to verify the results. In summary, further studies are needed to confirm the molecular mechanism of IA and its relationship with immune infiltration.

Conclusion

In this study, we introduced an effective method to evaluate the co-expression network of IA and normal conditions. Several DEGs and two modules significantly related to the disease were identified. GO and KEGG functional enrichment analyses showed enrichment of pathways related to inflammatory and immune processes. Hub genes related to chemotaxis may play pivotal roles in the formation and progression of IA. Further analyses suggested that macrophages and mast cells may be the driving factors of IA, which should be examined in further studies.