Exploration of the Shared Gene Signatures and Molecular Mechanisms Between Ischemic Stroke and Atherosclerosis

Introduction

Atherosclerosis (AS) is an important pathogenetic basis for the high incidence of cardiovascular and cerebrovascular diseases in the middle-aged and elderly populations, and researchers have done a great deal of work on it. It is a common and frequent disease that directly affects human health, and it is the first cause of death in the global population;¹ Ischemic stroke (IS) refers to the occurrence of brain tissue infarction caused by the blockage of cerebral arteries. This condition is characterized by the impairment of neurons, astrocytes, and oligodendrocytes. IS is considered the primary vascular event within the central nervous system that significantly contributes to mortality and disability in contemporary society. Furthermore, it is recognized as a major global cause of both disability and death.² Atherosclerosis of the carotid arteries leads to carotid artery stenosis, which in turn is one of the most important triggers of IS.

This research presents a collaborative analysis of two diseases, AS and IS, where the related genes of each disease were screened independently to identify the common related genes between the two disorders. The WGCNA technique was employed to immediately identify the modules that exhibit substantial associations with the diseases. Subsequently, by establishing a network and pinpointing the crucial genes inside the network, we ultimately discovered the relevant genes that are associated to the two diseases.

The objective of our work is to examine the inherent connection between AS and IS, with the intention of discovering specific biological clock markers that are linked to both disorders. While we conducted clinical PCR validation, it is important to note that this does not guarantee the reliability of our study. To establish greater confidence, a longer duration and a bigger sample size for follow-up studies are necessary.

Methods

GEO Dataset Download and Process

We searched the NCBI GEO database for AS and IS gene expression profiles using “atherosclerosis” or “ischemic stroke”.³ The following criteria filter the dataset: First, gene expression profiling must include cases and controls. Firstly, it is imperative that the gene expression profiling incorporates both cases and controls. Furthermore, the selection of the organization should be based on blood-related criteria. Additionally, it is imperative to maintain a minimum sample size of 10 in each group to uphold the precision of the Weighted Gene Co-expression Network Analysis (WGCNA). Finally, the GEO dataset numbered GSE9874, GSE22255, GSE23746 and GSE16561,⁴ GSE22255⁵ GSE23746, and GSE16561^6–8 were selected. The Series Matrix Files, furnished by contributors, were received and the probes were aligned with their gene symbols following the annotation document of the respective platforms. The details of the four datasets, including GSE number, detection platforms, samples, and RNA source types, are consolidated in Table 1. We designated GSE9874 and GSE22255 as a discovery cohort, and GSE23746 and GSE16561 as a validation cohort for the crucial DEG analysis.

Table 1 Summary of Those Four GEO Datasets Involving as and is Patients

Weighted Gene Co-Expression Network Analysis

The weighted gene co-expression network analysis (WGCNA) is a potent algorithm capable of identifying co-expressed gene modules of high biological significance and investigating the correlation between gene networks and diseases.⁹ We used WGCNA package Version 1.61¹⁰ in R.4.0.3 software (https://cran.r-project.org/web/packages/WGCNA/index.html) to obtain the AS and IS associated modules in GSE9874 and GSE22255, respectively. Initially, the hierarchical clustering analysis was conducted using the Hclust Function. Subsequently, the suitable soft powers β (ranging from1 to 30) were chosen using the “pickSoftThreshold” function, adhering to the scale-free network standard. The soft power value β and gene correlations matrix, calculated by Pearson analysis for all gene pairs, were then used to construct the adjacency matrix. The topological overlap matrix (TOM) and the corresponding dissimilarity (1-TOM) were derived from the adjacency matrix. A hierarchical clustering dendrogram was then constructed, and similar gene expressions were categorized into different modules (Parameter: minModuleSize = 100, CutHeight = 0.25). Finally, the module eigengene (ME) was used to summarize the expression profiles of each module, and the association between the ME and clinical characteristics was computed. As a result, modules with a high correlation coefficient with clinical features (|Cor|>0.3 and p<0.05) were prioritized, and the genes in these modules were selected for future investigation. The shared genes associated with both AS and IS were overlapped using Jvenn.¹⁰

Identification of DEGs in as and is

We utilized the Limma package (https://cran.r-project.org/web/packages/WGCNA/index.html) Version 1.61¹¹ in R.4.0.3 software to screen the differentially expressed genes (DEGs) between case samples and normal CTRL samples in GSE9874 and GSE22255, respectively. The FDR < 0.05 and |log2fold change (FC)| > 0.5 were chosen as the cut-off criteria. Then we kept the intersection of DEGs in the two comparison groups as the DEGs which are related with both of AS and IS.

Construction of PPI Network and Identification of Hub DEGs

We evaluated information on protein-protein interactions (PPIs) using the STRING database, which is accessible at https://cn.string-db.org/¹² (Confidence score > 0.4 was set as cut-off criteria). In addition, to explore the relationship between DEGs, we converted the results visually by using Cytoscape (http://www.cytoscape.org/) Version 3.9.0.¹³ Cytoscape plugin ClueGO Version 2.5.9 categorised non-redundant GO terms and depicted them as a functionally grouped network to analyze these AS and IS shared genes’ probable roles in the PPI network.¹⁴ We used the Cytoscape plugin CytoHubba Version 0.1¹⁵ to discover hub DEGs utilizing MCC, MNC, DEGREE, and EPC. After rating the top 20 genes in the network using four methods, we utilized venn analysis to find the hub genes.

Finally, we compared the AS and IS-related genes from the WGCNA result and the hub genes, they are the most important hub genes related to both diseases.

Construction of as and is-Related miRNA-mRNA Network

We further explore whether some microRNAs are regulating risk genes in AS and IS. The HumanmicroRNA Disease Database v3.2 (HMDD) (http://www.cuilab.cn/hmdd) is a database that collects evidence for the relationships between human miRNAs and diseases. We used this database to acquire the miRNAs linked with AS and IS.¹⁶ To explore miRNAs’ function, we used the mirPath v3.0 software (https://dianalab.e-ce.uth.gr/html/mirpathv3/) in the DIANA tool¹⁷ to perform the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis. The KEGG pathways with p-values <0.05 were considered significant. Then, we got the target genes regulated by related miRNAs from the miRWalk 3.0 database¹⁸ (http://129.206.7.150/). The previous step’s intersection of the relevant miRNAs’ target genes and the most important hub genes related to both diseases was used to establish the miRNAs-mRNAs controlled connection. Following that, we performed KEGG analysis of the regulated hub genes using the KEGG Orthology-Based Annotation System (KOBAS) (http://kobas.cbi.pku.edu.cn/index.php). Finally, the KEGG pathways of AS and IS-related miRNAs and regulated hub genes were used to construct miRNA-mRNA regulation networks.

Validation of Hub Signatures in Independent Cohort

In order to verify the consistency of hub gene expression level distribution, We performed the expression distribution display in both discovery cohort datasets (GSE9874 and GSE22255) and additional validation datasets (GSE23746 and GSE16561). A flow diagram of the study is shown in Figure 1.

Figure 1 The workflow chart of the study.

Quantitative Real-Time Reverse Transcriptions PCR

The venous blood samples were obtained at Liaocheng People’s Hospital. The samples were collected between 6–8 a.m. on the day after a 12 p.m. cutoff time for food intake. The blood samples were promptly transported to the hospital’s laboratory for centrifugation following their collection.

Peripheral blood samples were collected from a cohort of 30 individuals, and total RNA was extracted using Total RNA Extraction Reagent (19201ES60, Yeasen, China). Real-time quantitative PCR was subsequently conducted using the One-Step qRT-PCR Kit (D7268S, Beyotime, China) and specific primers designed for this purpose (Table 2). The gene used as a reference in this study was glyceraldehyde-3-phosphate dehydrogenase (GAPDH). The statistical analyses were conducted and visualized using GraphPad Prism 8 software. A significance level of less than 0.05 was deemed as indicative of a significant difference.

Table 2 Primers Used for RT-PCR in This Study

Results

Cohort of Discovery: AS and is Co-Expression Modules Determined by GSE9874 and GSE22255 Expression Profiles

The expression profiles of GSE9874 and GSE22255 yielded over 20,000 genes per dataset. However, the majority of these genes did not exhibit differential expression across samples. Consequently, we focused on the top quartile of genes that displayed significant variation, which consisted of approximately 5000 genes. In this study, the soft threshold β was 4 in the WGCNA analysis of AS and 12 in IS. Through the application of WGCNA, we discerned a total of eight modules in both GSE9874 and GSE22255, each module represented by a distinct color. Subsequently, we constructed a heatmap to illustrate the relationships between modules and traits. This was based on the Spearman correlation coefficient, which was used to assess the connection between each module and the disease (Figures 2A and B). Three modules: turquoise, brown, yellow and three modules: yellow, turquoise, black have high positive association with IS and AS, respectively; At the same time, three modules: black, blue, green and two modules: blue, red have high negative association with IS and AS, respectively.

Figure 2 (A). Left: The cluster dendrogram of co-expression genes in GSE22255 (IS). Right: Module–trait relationships in IS. Each cell contains the corresponding correlation and p-value. (B) Left: The cluster dendrogram of co-expression genes in GSE9874 (AS). Right: Module–trait relationships in IS. Each cell contains the corresponding correlation and p-value.

We then compared positive and negative correlation module genes. Figure 3 shows that 236 genes overlapped in AS and IS positivity modules and 29 genes in negativity modules. They may be closely linked to AS and IS pathogenesis.

Figure 3 The shared genes between Positive (A) and Negative (B) modules in AS and IS.

Identification of DEGs in as and is

A gene signature was identified by comparing mRNA expression levels in AS and IS versus CTRL in the GSE9874 and GSE22255 datasets. AS versus CTRL gene expression analysis revealed 608 DEGs, 254 upregulated and 354 downregulated in GSE9874 (Figure 4A) and 552 DEGs, 201 upregulated and 351 downregulated in GSE22255 (Figure 4B). Then we compared the two DEG sets and kept the intersection them, as shown in Figure 5, we finally shared 30 downregulated DEGs and 43 upregulated DEGs in both AS and IS.

Figure 4 The DEGs of GSE9874 (A) and GSE22255 (B). The blue and red dots represent DEGs filtered based on the cutoff criteria of adjusted |log2 (fold change)| >0.5 and FDR < 0.05, While the grey dots represent genes that do not satisfy the cutoff criteria.

Figure 5 Comparison Venn of DEGs in AS and IS. The graph on the left shows up-regulated genes and the graph on the right shows down-regulated genes.

Construction of PPI Network and Identification of Hub DEGs

PPI networks of the shared DEGs which contained 70 nodes and 291 edges were constructed using STRING and visualized in Cytoscape (Figure 6).

Figure 6 PPI network of shared DEGs.Green and pink nodes represent downregulated and upregulated shared DEGs.

We next used GlueGo to enrich GO (Figure 7). CCR1 chemokine receptor binding, T-helper cell differentiation regulation, and arsenic response were the top three enriched GO keywords. This function may be relevant in AS and IS, as CCR1 chemokine receptor binding accounted for 35% of total GO keywords and was related with 17 genes. Next, we used the CytoHubba models in Cytoscape to rank the top 20 network genes, including Degree, MCC, MNC, and EPC (Figure 8). Notably, the top 20 genes from topological analysis algorithms included 19 intersection hub genes: IRF4, TRIB1, CCL7, KRAS, CCL4, IER3, PHLDA1, CSF1, KLF10, JUNB, PPP1R15A, TIMP, DUSP2, ATF3, TNFAIP3, CCL3, CTSL, ZC3H12A, IER2, which may involve in the development of AS and IS. We compared the AS and IS-related genes from the WGCNA result (236 positive-related genes, 29 negative-related genes) and the 19 hub genes. Finally, we got 6 overlapped genes: ATF3, CCL3, CCL4, JUNB, KRAS, ZC3H12A. They are the most important hub genes related to both AS and IS based on our results.

Figure 7 The interaction network of GO terms generated by the Cytoscape plug-in ClueGO and the significant term of each group is highlighted; Proportion of each GO terms group in the total.

Figure 8 Venn display of top 20 genes based on Degree, MCC (Maximal Clique Centrality), MNC (Maximum Neighborhood Component), EPC (Edge Percolated Component).

Construction of as and is-Related miRNA-mRNA Network

Based on the findings from the HMDD database, a total of 116 miRNAs were identified to have associations with AS, whereas 25 miRNAs were found to be linked with IS. Then we searched miRNAs that can regulate 6 hub genes in miRWalk 3.0, we got 27 paired miRNA-mRNA relationships, including 26 miRNAs, four of which were overlapped with AS and IS-related miRNAs: hsa-miR-330-5p, hsa-miR-143-3p, hsa-miR-16-5p, hsa-miR-152-3p. Next, we performed KEGG enrichment analysis for the miRNAs and the 6 hub genes, we got 31 KEGG pathways related to the miRNAs from mirPathv3, and 58 KEGG pathways related to the 6 hub genes through KOBAS. We compared the KEGG pathways, there are 9 overlapped terms: ErbB signaling pathway, Choline metabolism in cancer, Neurotrophin signaling pathway, Thyroid hormone signaling pathway, FoxO signaling pathway, Insulin signaling pathway, Proteoglycans in cancer, MAPK signaling pathway, as shown in Figure 9.

Figure 9 miRNA-6 hub genes-KEGG pathway network. Red, yellow nodes represent 6 DEGs and related miRNAs, blue and green nodes represent 6 genes related to KEGG and overlapped 9 KEGG pathways; the Green line represents miRNA-KEGG, the grey line represents miRNA-DEGs, the red line represents KRAS- overlapped 9 KEGG pathways.

Hsa-miR-330-5p, hsa-miR-143-3p, hsa-miR-16-5p, and hsa-miR-152-3p are related to AS or IS, and they all regulate a common hub gene: KRAS. What’s more, the four miRNAs and KRAS are all related to 8 KEGG pathways: ErbB signaling pathway, Choline metabolism in cancer, Neurotrophin signaling pathway, Thyroid hormone signaling pathway, FoxO signaling pathway, Insulin signaling pathway, Proteoglycans in cancer, MAPK signaling pathway, shown in Figure 10.

Figure 10 miRNA-KRAS-KEGG pathway network. Red, yellow nodes represent KRAS and AS or IS-related miRNAs, green nodes represent overlapped KEGG pathways; the Green line represents miRNA-KEGG, the grey line represents miRNA-DEGs, the red line represents KRAS- overlapped 9 KEGG pathways.

Validation of Hub Signatures in Independent Cohort and Clinical Samples

To validate our results, to verify the consistency of hub gene expression level distribution, we first performed the expression distribution display in GSE9874 and GSE22255, shown in Figure 11, all six genes are upregulated in AS and IS samples. In addition, we also investigated the expression level of the six hub genes in additional validation datasets (GSE23746 and GSE16561), all of the expression levels were consistent with the distribution in the discovery dataset. In IS, the up-regulation trend of all the genes was significant, while in AS, the up-regulation trend of ATF3, CCL3, and KRAS was significant, although the expression levels of CCL4, JUNB, and ZC3H12A were not significant between AS and CTRL samples, the direction of differences were consistent with the discovery data set. Their expression trends were consistent with the results of peripheral blood qPCR validation (Figure 12), Table 3 lists all sample details.

Table 3 Clinical Features of in vivo Validation

Figure 11 Six hub genes expression levels in IS-related datasets (GSE22255 and GSE16561) (A) and AS-related datasets (GSE9874 and GSE23746) (B). * p<0.05, ** p<0.01, ***p<0.001.

Figure 12 In vivo verification Relative Expression. Levels of Hub Genes in Patients Clinically Diagnosed with IS (A) and AS (B). * p<0.05, ** p<0.01, ***p<0.001.

Discussion

Ischemic stroke and carotid atherosclerosis, both of which have a high incidence, are high-risk factors for neurologic deterioration (ND) and often occur simultaneously. The pathophysiology of ischemic stroke and carotid atherosclerosis share many common features, and both have similar causative factors and influences, such as hypertension, hyperlipidemia, and hyperglycemia. Studies have shown that chronic inflammatory processes involving the vessel wall may induce atherosclerosis and are strongly associated with the development and prognosis of patients with ischemic stroke.¹⁹ Inflammation, oxidative stress, and vascular aging are risk factors for atherosclerosis formation and acute stroke. Platelet-related biological markers, such as PAg-ADP, PAg-COL, and FIB index levels, are strongly associated with atherogenesis and acute stroke;²⁰ Positive correlation between calcitoninogen levels and carotid intima-media thickness formation in patients with acute ischemic stroke;²¹ Elevated risk of hemorrhagic transformation in patients with acute stroke with atherosclerosis.²² The above studies have contributed to our understanding of the relationship between ischemic stroke and carotid atherosclerosis, and have helped to reveal common molecular mechanisms between the two. However, there seem to be few studies exploring the pathogenesis between ischemic stroke and carotid atherosclerosis at the genetic level, and to our knowledge, this is the first time that data from public databases on ischemic stroke and carotid atherosclerosis have been integrated to explore common mechanisms between the two.

The WGCNA algorithm and DEG analysis revealed that IS and carotid atherosclerosis may share chemokine receptor binding and T helper cell differentiation. Chemokines, low-molecular-weight proteins grouped into four subfamilies—C, CC, CXC, and CX3C—impact leukocyte trafficking during inflammation. The CCR1 gene encodes a beta chemokine receptor, a G-protein-coupled receptor. Several brain studies have found CCR1 on neurons, microglia, and astrocytes. CCR1, a chemokine membrane receptor, is abundantly expressed in the brain. CCR1 is activated early in acute ischemic stroke and higher in cerebral hemorrhage and Alzheimer’s disease patients.^23,24 Expressed by macrophages and lymphocytes, increased exogenous CCR1 expression attenuates plaque formation;²⁵ Blockade of the CCR1 in vivo chemokine pathway reduces the progression of atherosclerosis in a mouse model of hypercholesterolemia.²⁶ CCR1 increases the activation and inflammatory response of spinal microglia and enhances pain sensitivity in rats, causing neuropathic pain.²⁷ The above studies are consistent with our study, showing that CCR1 promotes neuroinflammation and facilitates the deterioration of neurological disorders, as well as increasing the formation of atherosclerotic plaques, regulation of T-helper cell differentiation, a type of T-cells on the surface of which an antigen-recognition receptor, the TCR, interacts with antigen-presenting cells, the APCs, for type II MHC molecules, which can transmit the immune signals. CD4 helper T cells are also important in driving inflammation in atherosclerotic plaques, and helper T cells are recruited to and infiltrate plaques, leading to an increase in proinflammatory cytokines.^28,29 T-helper affects ischemic stroke by participating in the immune/inflammatory response.²⁸

T-helper cells act by coordinating and modulating the immune response, and acute ischemic stroke promotes an acute inflammatory cascade during which helper T cells are activated, leading to post-ischemic brain injury.^29–31

As mentioned above, we found that chemokine receptor binding and regulation of T-helper cell differentiation were enriched in both atherosclerosis and ischemic stroke, and by integrating the multivariate databases, we ended up with 6 hub genes (ATF3, CCL3, CCL4, JUNB, KRAS, and ZC3H12A), which were significantly up-regulated in both atherosclerosis and ischemic stroke groups as compared to the control group. The ATF3 gene is a gene that encodes a modifiable transcription factor that is involved in the regulation of a variety of physiological processes such as cell growth, apoptosis, and stress responses. ATF3 is not expressed in healthy and intact neurons but is expressed when axons are damaged.³² ATF3 Deficiency Inhibits Cell Viability and Induces Apoptosis Can Reduce the Risk of Ischemic Stroke.^31,33 ATF3 is a key regulator of AS progression, directly affecting atherosclerotic lesions not only through the regulation of vascular homeostasis but also through systemic glycolipid metabolism and inflammatory responses. ATF3 regulates endothelial inflammation, foam cell production, and vascular remodeling, all of which have a direct impact on atherosclerotic lesions. Extravascular ATF3 expression is also involved in AS by regulating glycolipid metabolism and inflammatory responses. CCL3 and CCL4, members of the CC subfamily of chemokines, elevate intravascular atherosclerotic adhesion, promote atherosclerotic plaque formation, and are biological markers of carotid atherosclerotic plaque.^34–36 In ischemic stroke, levels of CCl3 and CCL4 are elevated and can recruit microglia, neutrophils, and monocytes for their chemotactic function.^37,38 JUNB can right cerebral vascular smooth muscle migration and atherosclerosis,^39,40 elevated expression after ischemic stroke.^41,42 ZC3H12A is induced by cholesterol and is involved in cholesterol-induced HUVEC DNA damage that exacerbates atherosclerotic plaque formation;^43,44 Whereas studies have shown that upregulation of ZC3H12A expression has a mitigating effect on ischemic stroke,⁴⁵ however, the mechanism of ZC3H12A’s role in atherosclerosis and ischemic stroke needs further investigation. Previous Bioinformatics Analysis Shows Elevated KRAS in Ischemic Stroke,⁴⁶ and can promote atherosclerotic plaque formation.⁴⁷ The mechanism of KRAS in ischemic stroke and atherosclerosis is not clear. By constructing a miRNA regulatory network, we obtained four miRNAs known to be associated with both IS and AS: hsa-miR-330-5p, hsa-miR-143-3p, hsa-miR-16-5p, hsa-miR-152-3p, and all regulate KRAS, suggesting that KRAS may be the most critical gene acting in atherosclerosis and ischemic stroke, possibly through ErbB signaling pathway, Choline metabolism in cancer, Neurotrophin signaling pathway, Thyroid hormone signaling pathway, FoxO signaling pathway, Insulin signaling pathway, Proteoglycans in cancer, MAPK signaling pathway.

We collected clinical samples to validate our identified pivotal genes, and the results showed that the expression levels of ATF3, CCL3, CCL4, JUNB, KRAS, and ZC3H12A were high in the ischemic stroke and atherosclerosis groups compared to the control group, which was in general agreement with our bioinformatic analysis described above.

We integrated ischemic stroke and carotid atherosclerosis public database data for the first time. We found that ATF3, CCL3, CCL4, JUNB, KRAS, and ZC3H12A may be common mechanisms affecting both, and these metrics are expected to be novel candidate genes for biomarkers or potential therapeutic targets for both. We validated the identified pivotal genes by collecting clinical samples and also obtained expression trends consistent with our analyses. This study has limitations. Due to the single-center nature of the study, the inclusion of clinical patients was limited. Therefore, researchers will need to carry out more multi-center and large-sample investigations in order to validate these findings in the future. Furthermore, there is a dearth of research on the genetic mechanisms of ATF3, CCL3, CCL4, JUNB, KRAS, and ZC3H12A. Therefore, it is imperative that we do further cellular and animal experiments to investigate these genes.