# Integrative Analyses of Biomarkers Associated with Endoplasmic Reticulum Stress in Ischemic Stroke

Authors: Xiaoting Zhang; Xi Li; Jinyan Gu; Jingpei Guo; Jiayao Chen; Ke Zhang; Junfeng Liu; Jiani Liu; Chao Peng; Hanwei Liu; Bin Zhou
Journal: Computational and Mathematical Methods in Medicine
 (2022)

---

## Abstract

Background. Neuronal apoptosis, which is the primary pathological transform of cerebral injury following ischemic stroke (IS), is considered to be induced by endoplasmic reticulum stress (ERS) by numerous reports. However, ERS biomarkers in IS have not been fully identified yet. Consequently, the present study is aimed at exploring potential blood biomarkers by investigating the molecular mechanisms of ERS promoting neuronal apoptosis following IS development. Methods. A comprehensive analysis was performed with two free-accessible whole-blood datasets (GSE16561 and GSE37587) from the Gene Expression Omnibus database. Genetic information from 107 IS and 24 healthy controls was employed to analyze the differentially expressed genes (DEGs). Genes related to ERS (ERS-DEGs) were identified from the analysis. Enrichment analyses were performed to explore the biofunction and correlated signal pathways of ERS-DEGs. Protein-protein interaction (PPI) network and immune correlation analyses were performed to identify the hub genes along with their correspondent expressions and functions, all of which contributed to incremental diagnostic values. Results. A total of 60 IS-related DEGs were identified, of which 27 genes were confirmed as ERS-DEGs. GO and KEGG enrichment analysis corroborated that upregulated ERS-DEGs were principally enriched in pathways related to immunity, including neutrophil activation and Th17 cell differentiation. Moreover, the GSEA and GSVA indicated that T cell-related signal pathways were the most considerably immune pathways for ERS-DEG enrichment. A total of 10 hub genes were filtered out via the PPI network analysis. Immune correlation analysis confirmed that the expression of hub genes is associated with immune cell infiltration. Conclusions. By integrating and analyzing the two gene expression data profiles, it can be inferred that ERS may be involved in the development of neuronal apoptosis following IS via immune homeostasis. The identified hub genes, which are associated with immune cell infiltration, may serve as potential biomarkers for relative diagnosis and therapy.


---

## Body

## 1. Introduction

Stroke is one of the leading causes of mortality and disability worldwide. Ischemic stroke (IS) attributes to over 85% morbidity of total stroke cases [1]. Injury ensuing consequent to the interruption in cerebral blood flow, known as “ischemia and brain hypoxia,” can ultimately result in cell degeneration, apoptosis, and neurological dysfunction. Pathology of the transversion entangles a variety of cellular and biochemical molecular mechanisms, including apoptosis, neuroinflammation, oxidative stress, glutamate excitotoxicity, and energy depletion, in which neuronal apoptosis was considered to be the primary pathological metamorphosis following IS [2]. Early and effective neuroprotective therapies are needed to decrease the disability and mortality rates of IS [3]. However, the complexity of its pathophysiological mechanisms poses a major challenge to clinical treatments, and effective neuroprotectors have not been established yet [4]. Therefore, it is of great demand to explore the molecular mechanisms that regulate neuronal apoptosis and identify potential biomarkers for improvement of the current diagnostic and therapeutic procedures.Numerous reports have confirmed that cerebral ischemia-induced ERS is a critical precipitating factor of neuronal apoptosis [5]. Endoplasmic reticulum (ER) conducts the proper folding and assembly of secretory and membranous proteins. However, when ER homeostasis is disturbed, unfolded and misfolded proteins accumulate in the ER lumen, resulting in ER stress (ERS). ERS can trigger neuronal apoptosis by multiple mechanisms, including inducing calcium overload, suppressing protein synthesis, and activating the inflammatory signal pathways [6]. Previous reports have shown that ERS may be a potential target for the diagnosis and treatment of IS; however, the key biomarkers of ERS in IS injury have not been identified yet [7]. Among all kinds of biomarkers, blood-based biomarkers contain a notable preponderance, such as minimal invasiveness, cost-effectiveness, high patient acceptability, and availability in different clinical settings [8]. The combination of microarray techniques and bioinformatics analysis was classified as a high-throughput, effective, and comprehensive method to detect the objective biomarkers [9]. Thus, we used bioinformatics analysis for the exploration of ERS blood-based biomarkers.The aim of the present IS study is to explore the blood biomarkers associated with ERS by utilizing the whole-blood samples from the Gene Expression Omnibus (GEO) database. Differentially expressed genes (DEGs) between the patients with IS and healthy controls were identified, and the DEGs that correlated to ERS (ERS-DEGs) were further screened. Moreover, the biological functions of ERS-DEGs were analyzed via four enrichment analyses. Hub genes were screened out from ERS-DEGs via PPI analysis. Furthermore, correlation analysis and immune infiltration were used to assess the diagnostic value of the hub genes. In conclusion, the present study identified a few ERS-related biomarkers that may serve as potential determinants for further diagnosis and treatment of IS.


## 2. Materials and Methods

### 2.1. Data Download and Preprocess

GEO (http://www.ncbi.nlm.nih.gov/geo) [10] is a free-accessible database that provides reliable profiles of IS, from where two gene expression profiles, GSE16561 [11] and GSE37587 [12], were searched with the keyword “ischemic stroke” and downloaded originally via GEOquery package [13] of R software V.3.6.5 (http://www.r-project.org/). Data type description is termed “expression profiling by array,” and the species is Homo sapiens. The platform for both databases is “GPL6883 Illumina HumanRef-8 v3.0 expression beadchip,” which comprises a total of 131 human whole-blood samples. Specifically, GSE16561 consists of the whole-blood samples of 39 patients with IS and 24 healthy controls, and GSE37587 contained 68 whole-blood samples of patients with IS. The raw datasets of GSE16561 and GSE37587 were combined. The sva package [14] was employed for background correction and data normalization, and the obtained result was demonstrated using the boxplot diagrams. Furthermore, diagrams of principal component analysis (PCA) were drawn using the ggplot2 package.


### 2.2. Identification of Differentially Expressed Genes

DEGs in the combined dataset were screened out via limma package [15] according to the inclusion criteria of adjusted p value < 0.05 and ∣log2FC∣>0.1. Visualization of DEGs was plotted using the ggplot2 package along with volcano plot and heatmap. The list of ERS-related genes (ERSGs) was procured from the GeneCards database (https://www.genecards.org/) [16]. Venn diagram was applied to illustrate the intersection of corresponding genes, which are related to DEGs and ERS and thus, contributing to the filtrate, ERS-DEGs.


### 2.3. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Analysis

GO analysis is a conventional and pragmatic method for interpreting the characteristic biological functions of genes [17], including annotation of biological processes (BP), molecular functions (MF), and cellular components (CC). KEGG analysis provides all known biochemical pathways documented in a comprehensive biological database [18]. GO and KEGG pathway enrichment analyses for the ERS-DEGs were implemented via clusterProfiler package [19]. An adjusted p value < 0.05 would be considered statistically significant. Furthermore, we selected the GO/KEGG pathway terms with the highest enrichment degree to draw the network maps, respectively.


### 2.4. Gene Set Enrichment Analysis (GSEA) and Gene Set Variation Analysis (GSVA)

To further evaluate the significant alterations in the functions and pathways of the gene expression matrix of ERS-DEGs, GSEA was performed using the clusterProfiler R package. Reference gene sets were selected as “c2.cp.kegg.v7.0.symbols.gmt.” Results with afalsediscoveryrateFDR<0.25 and adjusted p<0.05 were set as the threshold for determining the significance of enrichment. Based on the results obtained from the GSEA, we explored the signal pathways of ERS-DEGs in IS via the GSVA R package [20]. Enrichment was considered significant for adjusted p<0.05.


### 2.5. Protein–Protein Interaction Network Analysis and Hub Gene Identification

The protein-protein interaction (PPI) network of ERS-DEGs was constructed by the Retrieval of Interacting Genes (STRING) database (http://string-db.org) [21], which is an online tool to identify functional interactions. Furthermore, “cytoHubba,” a plugin in Cytoscape software, was employed to screen hub genes within the PPI network.


### 2.6. Network Analysis of Hub Genes

The online visual analytics platform NetworkAnalyst (https://www.networkanalyst.ca/) [22] was used to predict latent transcription factors (TFs) and miRNAs associated with the hub genes. The interaction of ERS-DEGs and potential TFs was evaluated using the obtained results of Gene Regulatory Networks (GRN) and gene-targeted record of ENCODE ChIP-seq data. The TF-miRNA coregulatory network and RegNetwork databases were used to explore the interactions between ERS-DEGs and potential miRNAs. The interaction data of ERS-DEGs and potential drugs was obtained from the DrugBank database and the analysis of protein-drug interactions in diseases, drugs, and chemicals.


### 2.7. Correlation Analysis of Hub Genes and Immune Infiltration

CIBERSORT is a tool used for deconvolution of the expression matrix of transcriptome, based on the principle of linear support vector regression. The tool can also be used to estimate the compositions and abundances of immune cells in a mixed cell population [23]. Thus, the gene expression data was uploaded to CIBERSORT, and the cut-off criteria for statistical significance was set at a p value less than 0.05. Consequently, the matrix data of infiltrating immune cells was obtained. The distribution of 22 types of infiltrating immune cells in each sample was presented on the heatmap plot using the heatmap R package. The correlation of 22 types of infiltrating immune cells was also visualized in the heatmaps, drawn by the corrplot package in R software.


## 2.1. Data Download and Preprocess

GEO (http://www.ncbi.nlm.nih.gov/geo) [10] is a free-accessible database that provides reliable profiles of IS, from where two gene expression profiles, GSE16561 [11] and GSE37587 [12], were searched with the keyword “ischemic stroke” and downloaded originally via GEOquery package [13] of R software V.3.6.5 (http://www.r-project.org/). Data type description is termed “expression profiling by array,” and the species is Homo sapiens. The platform for both databases is “GPL6883 Illumina HumanRef-8 v3.0 expression beadchip,” which comprises a total of 131 human whole-blood samples. Specifically, GSE16561 consists of the whole-blood samples of 39 patients with IS and 24 healthy controls, and GSE37587 contained 68 whole-blood samples of patients with IS. The raw datasets of GSE16561 and GSE37587 were combined. The sva package [14] was employed for background correction and data normalization, and the obtained result was demonstrated using the boxplot diagrams. Furthermore, diagrams of principal component analysis (PCA) were drawn using the ggplot2 package.


## 2.2. Identification of Differentially Expressed Genes

DEGs in the combined dataset were screened out via limma package [15] according to the inclusion criteria of adjusted p value < 0.05 and ∣log2FC∣>0.1. Visualization of DEGs was plotted using the ggplot2 package along with volcano plot and heatmap. The list of ERS-related genes (ERSGs) was procured from the GeneCards database (https://www.genecards.org/) [16]. Venn diagram was applied to illustrate the intersection of corresponding genes, which are related to DEGs and ERS and thus, contributing to the filtrate, ERS-DEGs.


## 2.3. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Analysis

GO analysis is a conventional and pragmatic method for interpreting the characteristic biological functions of genes [17], including annotation of biological processes (BP), molecular functions (MF), and cellular components (CC). KEGG analysis provides all known biochemical pathways documented in a comprehensive biological database [18]. GO and KEGG pathway enrichment analyses for the ERS-DEGs were implemented via clusterProfiler package [19]. An adjusted p value < 0.05 would be considered statistically significant. Furthermore, we selected the GO/KEGG pathway terms with the highest enrichment degree to draw the network maps, respectively.


## 2.4. Gene Set Enrichment Analysis (GSEA) and Gene Set Variation Analysis (GSVA)

To further evaluate the significant alterations in the functions and pathways of the gene expression matrix of ERS-DEGs, GSEA was performed using the clusterProfiler R package. Reference gene sets were selected as “c2.cp.kegg.v7.0.symbols.gmt.” Results with afalsediscoveryrateFDR<0.25 and adjusted p<0.05 were set as the threshold for determining the significance of enrichment. Based on the results obtained from the GSEA, we explored the signal pathways of ERS-DEGs in IS via the GSVA R package [20]. Enrichment was considered significant for adjusted p<0.05.


## 2.5. Protein–Protein Interaction Network Analysis and Hub Gene Identification

The protein-protein interaction (PPI) network of ERS-DEGs was constructed by the Retrieval of Interacting Genes (STRING) database (http://string-db.org) [21], which is an online tool to identify functional interactions. Furthermore, “cytoHubba,” a plugin in Cytoscape software, was employed to screen hub genes within the PPI network.


## 2.6. Network Analysis of Hub Genes

The online visual analytics platform NetworkAnalyst (https://www.networkanalyst.ca/) [22] was used to predict latent transcription factors (TFs) and miRNAs associated with the hub genes. The interaction of ERS-DEGs and potential TFs was evaluated using the obtained results of Gene Regulatory Networks (GRN) and gene-targeted record of ENCODE ChIP-seq data. The TF-miRNA coregulatory network and RegNetwork databases were used to explore the interactions between ERS-DEGs and potential miRNAs. The interaction data of ERS-DEGs and potential drugs was obtained from the DrugBank database and the analysis of protein-drug interactions in diseases, drugs, and chemicals.


## 2.7. Correlation Analysis of Hub Genes and Immune Infiltration

CIBERSORT is a tool used for deconvolution of the expression matrix of transcriptome, based on the principle of linear support vector regression. The tool can also be used to estimate the compositions and abundances of immune cells in a mixed cell population [23]. Thus, the gene expression data was uploaded to CIBERSORT, and the cut-off criteria for statistical significance was set at a p value less than 0.05. Consequently, the matrix data of infiltrating immune cells was obtained. The distribution of 22 types of infiltrating immune cells in each sample was presented on the heatmap plot using the heatmap R package. The correlation of 22 types of infiltrating immune cells was also visualized in the heatmaps, drawn by the corrplot package in R software.


## 3. Results

### 3.1. Preprocessing of Datasets and Identification of DEGs

The two gene expression datasets (GSE16561 and GSE37587) of blood samples were first merged. The resulting single dataset was normalized and standardized prior to analysis. Illustration of datasets before and after calibration is shown in the boxplots (Figures1(a) and 1(b)). We proceeded to evaluate the variations between healthy controls and the IS group via PCA and sample clustering analysis. As shown in Figures 1(c) and 1(d), the clustering of each group was more marked following the data preprocessing. Hence, the sample data were considered as a reliable source.Figure 1
Illustration of GSE16561 and GSE37587 datasets before and after calibration. (a, b) The boxplot diagrams of the datasets with the interbatch difference were removed before and after to calibration. (c, d) The PCA diagrams of the datasets with the interbatch difference were removed before and after calibration.
(a)(b)(c)(d)Following data preprocessing, a total of 60 DEGs; including 38 downregulated genes and 22 upregulated genes, were extracted from the gene expression matrix. This is shown in the heatmap and volcano plots (Figures2(a) and 2(b)). Subsequently, a total of 7092 ERSGs were sorted out from the GeneCards database, and Venn diagram was applied to visualize the overlap and differences between DEGs and ERSGs. The two datasets showed an overlap of 27 ERS-DEGs (Figure 2(c)).Figure 2
ERS-DEGs were obtained from the intersection of DEGs and ERS-related genes. (a) Heatmap of DEGs; (b) volcano plot of the distributions of DEGs. Red and blue dots represent upregulated and downregulated genes, respectively. No significantly changed genes were marked as grey dots; (c) Venn diagram of ERS-DEGs. Venn diagram analysis executed an intersection of DEGs and ERSGs to reveal the ERS-DEGs.
(a)(b)(c)


### 3.2. KEGG and GO Enrichment Analysis of ERS-DEGs

GO analysis results indicate that the upregulated ERS-DEGs were significantly enriched in BP pathways, including immune response-activating cell surface receptor signaling pathway, immune response-activating signal transduction, neutrophil degranulation, and neutrophil activation; all of which are involved in immune homeostasis. As for MF and CC analysis, the ERS-DEGs were mainly implicated with the external side of plasma membrane and MHC class II receptor activity, respectively (Figures3(a)–3(d)). The results are detailed in Table 1. KEGG pathway analysis demonstrated that the upregulated ERS-DEGs were predominantly enriched in immune-related pathways, including leishmaniasis, staphylococcus aureus infection, Th17 cell differentiation, systemic lupus erythematosus, phagosome, and tuberculosis (Figures 3(e) and 3(f)). These results are shown in detail in Table 2. The pathways that were most significantly enriched from ERS-DEGs in GO and KEGG analysis are shown in diagrams (Figures 3(g) and 3(h)).Figure 3
GO/KEGG enrichment analyses of ERS-DEGs. (a) Advanced bubble chart shows GO functional enrichment significance items of ERS-DEGs in three functional groups: BP, MF, and CC. Thex-axis and y-axis labels represent the gene ratios enriched in the GO analysis group and GO terms, respectively. The color of dots represents the adjusted p value: the redder the color, the lesser the adj. p value; the bluer the color, the greater the adj. p value. The size of dots represents gene counts. (b–d) Chord diagrams show the distribution of ERS-DEGs in the three different GO-enriched functions. The color of dots represents the log FC of the genes. (e) Advanced bubble chart shows enrichment of ERS-DEGs in signal pathways. (f) The chord diagram shows the distribution of ERS-DEGs in different KEGG pathways. (g) The pathway diagram shows the immune response-activating cell surface receptor signaling pathway with the highest enrichment degree in GO analysis. (h) The pathway diagram shows that the leishmaniasis pathway has the highest enrichment degree in KEGG analysis.
(a)(b)(c)(d)(e)(f)(g)(h)Table 1
GO enrichment analysis of ERS-DEGs.
OntologyIDDescriptionGeneRatioBgRatiop valuep.adjustq valueBPGO:0043312Neutrophil degranulation12/49485/186703.18e-091.67e-061.30e-06BPGO:0002283Neutrophil activation involved in immune response12/49488/186703.41e-091.67e-061.30e-06BPGO:0042119Neutrophil activation12/49498/186704.28e-091.67e-061.30e-06BPGO:0002446Neutrophil mediated immunity12/49499/186704.38e-091.67e-061.30e-06BPGO:0030098Lymphocyte differentiation9/49353/186702.94e-078.99e-057.00e-05CCGO:0009897External side of plasma membrane11/50393/197173.19e-094.09e-072.76e-07CCGO:0042581Specific granule7/50160/197171.51e-079.69e-066.53e-06CCGO:0070820Tertiary granule6/50164/197173.55e-061.51e-041.02e-04CCGO:1904724Tertiary granule lumen4/5055/197171.13e-053.54e-042.38e-04CCGO:0034774Secretory granule lumen7/50321/197171.56e-053.54e-042.38=-04MFGO:0019864IgG binding2/4911/176974.07e-040.0730.058Table 2
KEGG enrichment analysis of ERS-DEGs.
OntologyIDDescriptionGeneRatioBgRatiop valuep.adjustq valueKEGGhsa04640Hematopoietic cell lineage5/3299/80763.87e-050.0040.004KEGGhsa05140Leishmaniasis4/3277/80762.24e-040.0130.010KEGGhsa05340Primary immunodeficiency3/3238/80764.34e-040.0150.012KEGGhsa05150Staphylococcus aureus infection4/3296/80765.22e-040.0150.012KEGGhsa04659Th17 cell differentiation4/32107/80767.87e-040.0180.014


### 3.3. GSEA and GSVA

GSEA was performed to identify the biological pathways with significant alterations towards IS pathology. A total of 10 pathways were enriched by this method. Interestingly, five of them were closely associated with T cells, including T cytotoxic, T helper, IL17, NO2IL12, and CTL pathways (Figures4(a)–4(j)). The results are detailed in Table 3. Furthermore, GSVA was conducted to further explore the signal pathways in IS, including GNF2_CD33, GSE 3688 STAT5 AB knockin vs WT T cell IL 2 T-related 17h DN, and GSE29618_monocyte_vs_PDC_up (Figure 5).Figure 4
The top 10 significantly enriched pathways based on GSEA. (a) WP purine metabolism; (b) BioCarta T cytotoxic pathway; (c) BioCarta T helper pathway; (d) Reactome mRNA decay by 3′ to 5′ exoribonuclease; (e) BioCarta IL17 pathway; (f) Reactome transport of nucleoside-free purine and pyrimidine bases across the plasma membrane; (g) KEGG primary immunodeficiency; (h) BioCarta NO2IL12 pathway; (i) BioCarta CTL pathway; (j) Reactome mitochondrial calcium ion transport.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)Table 3
GSEA in IS.
IDDescriptionSet sizeEnrichment scoreNESp valuep.adjustq valuesRankGO:0031982Vesicle4850.1727146973.4112344380.0014326650.0041335450.001924525700GO:0005576Extracellular region4100.1749941963.2667012790.0014492750.0041335450.001924525435GO:0005739Mitochondrion2430.1422813612.2130932510.0014641290.0041335450.001924525983GO:0005615Extracellular space3290.190210293.2741988360.0014749260.0041335450.001924525252GO:0043230Extracellular organelle2630.220613933.5410205510.0014771050.0041335450.001924525621GO:0070062Extracellular exosome2610.2206988693.5249952460.0014771050.0041335450.001924525621GO:1903561Extracellular vesicle2630.220613933.5410205510.0014771050.0041335450.001924525621GO:0030054Cell junction2010.165794812.4084827940.0014814810.0041335450.001924525189GO:0031410Cytoplasmic vesicle3220.1457290152.487680470.0014814810.0041335450.001924525589GO:0097708Intracellular vesicle3220.1457290152.487680470.0014814810.0041335450.001924525589GO:0031967Organelle envelope2040.1868419672.7303292620.0014858840.0041335450.001924525788GO:0031975Envelope2040.1868419672.7303292620.0014858840.0041335450.001924525788GO:0030141Secretory granule1500.2257188682.8554146860.0014970060.0041335450.001924525404GO:0098796Membrane protein complex1500.1621314822.0510142490.0014970060.0041335450.001924525656GO:0099503Secretory vesicle1710.2220231932.9948720730.0015105740.0041335450.001924525404GO:0005740Mitochondrial envelope1370.1900353122.3069467260.0015290520.0041335450.001924525616GO:1990904Ribonucleoprotein complex1290.1796543652.1313439140.0015313940.0041335450.001924525902GO:0031966Mitochondrial membrane1300.1779346722.1169076690.0015337420.0041335450.001924525616GO:0019866Organelle inner membrane1000.244613892.6241194150.0015723270.0041335450.001924525609GO:0070161Anchoring junction1040.2515477142.7352353030.0015898250.0041335450.001924525149GO:0012506Vesicle membrane1210.1743482661.9974747970.003110420.0073979230.0034443771110GO:0030659Cytoplasmic vesicle membrane1180.1773064262.008090750.003129890.0073979230.0034443771110GO:0031090Organelle membrane4740.1024336962.0166298440.0043415340.0098156420.004570036656GO:0000228Nuclear chromosome129-0.149453948-1.9909431440.0057306590.0124164280.005780928647GO:1902494Catalytic complex2040.1251160271.8283255960.0133729570.027815750.0129506531192GO:0045202Synapse1220.1474493791.6933664670.0203442880.0406885760.018944074167GO:0000790Nuclear chromatin107-0.137337172-1.6669836440.0244565220.0471014490.021929825621Figure 5
Heatmap of 10 signal pathways of ERS-DEGs based on GSVA.


### 3.4. PPI Network Analysis

To identify the interactions of the ERS-DEGs, PPI network of ERS-DEGs was established and is shown in Figure6(a). The identified hub genes, predicted via cytoHubba plugin is also shown in the aforementioned figure. The top 10 hub genes in an increasing order of interaction degree are as follows: HIF1A, CREBBP, EP300, ARNT, TP53, PPARG, VHL, EPAS1, HIF1AN, and CREB1 (Figure 6(b)). The hub genes may be key ERS-related biomarkers of IS, which require further clinical trials for validation.Figure 6
Protein–protein interaction (PPI) networks. (a) PPI network analysis of ERS-DEGs. Red and blue colors represent high and low expression levels of genes, respectively. (b) The top 10 hub genes from the PPI network were predicted and shown on the subnetwork. Node color reflects the degree of accuracy in prediction with a gradient shift from yellow to red. Illustratively, the redder the color, the higher the accuracy of prediction.
(a)(b)


### 3.5. Network Analysis of Hub Genes

Transcriptional regulatory network analysis was performed to analyze the interaction between TFs and the hub genes (Figure7(a)). The top 3 hub genes associated with TFs were TP53, EPAS1, and VHL. Moreover, a network of miRNAs and hub genes was constructed (Figure 7(b)). The top 3 hub genes associated with miRNAs were TP53, CREB1, and HIF1A. As potential therapeutic targets, relative drugs and molecular compounds of the hub genes were also explored by the drug-gene interaction networks. As shown in Figure 7(c), a total of 7 drugs or molecular compounds were found to interact with hub genes of HIF 1A, CREB1, and HIF1AN. They are carvedilol, 2-methoxyestradiol, N-[(1-chloro-4-hydroxyisoquinolin-3-YL)carbonyl] glycine, naloxone, adenosine monophosphate, D-tartaric acid, and N-(carboxyvarbonyl)-D-phenylalanine.Figure 7
Network analyses of the hub genes. (a) Network of hub gene-TFs. Dots in blue to green represent the hub genes, and dots in red to yellow represent the TFs. (b) Network of hub gene-miRNA. Red and blue nodes are hub genes and miRNAs, respectively. (c) Networks of hub gene-related drug. Red and blue circle nodes are hub genes and drugs/molecular compounds, respectively.
(a)(b)(c)


### 3.6. Immune-Related Gene Identification and Functional Correlation Analysis

The result of the CIBERSORT analysis revealed that the T and NK cell subtypes were predominant among the 22 types of infiltrating immune cells (Figure8(a)). The correlation between 22 immune cell subtypes was also performed and presented via heatmap. The result revealed that M0 macrophages and resting NK cells had the strongest and weakest correlation with other immune cells, respectively (Figure 8(b)). Furthermore, the hub genes, ARNT, CREB1, CREBBP, EP300, EPAS1, HIF 1A, HIF 1AN, and TP53, were identified, and they correlated with 11 different types of immune cells (Figure 8(c)). We found a significant correlation between the number of neutrophil and T cell subtypes and the expression level of hub genes. The results suggested that neutrophils and T cells play a critical role in the development of ERS following IS.Figure 8
The nose cone is detachable upon impact. Evaluation and visualization of infiltrating immune cells. (a) Bar diagram displayed the proportion of infiltrating immune cells based on CIBERSORT algorithm. (b) Heatmap of the correlation between 22 subtypes of immune cells. Blue and red areas indicate positive and negative correlations, respectively. The darker the color, the stronger the correlation. (c) The correlations between hub genes and immune cells.
(a)(b)(c)


## 3.1. Preprocessing of Datasets and Identification of DEGs

The two gene expression datasets (GSE16561 and GSE37587) of blood samples were first merged. The resulting single dataset was normalized and standardized prior to analysis. Illustration of datasets before and after calibration is shown in the boxplots (Figures1(a) and 1(b)). We proceeded to evaluate the variations between healthy controls and the IS group via PCA and sample clustering analysis. As shown in Figures 1(c) and 1(d), the clustering of each group was more marked following the data preprocessing. Hence, the sample data were considered as a reliable source.Figure 1
Illustration of GSE16561 and GSE37587 datasets before and after calibration. (a, b) The boxplot diagrams of the datasets with the interbatch difference were removed before and after to calibration. (c, d) The PCA diagrams of the datasets with the interbatch difference were removed before and after calibration.
(a)(b)(c)(d)Following data preprocessing, a total of 60 DEGs; including 38 downregulated genes and 22 upregulated genes, were extracted from the gene expression matrix. This is shown in the heatmap and volcano plots (Figures2(a) and 2(b)). Subsequently, a total of 7092 ERSGs were sorted out from the GeneCards database, and Venn diagram was applied to visualize the overlap and differences between DEGs and ERSGs. The two datasets showed an overlap of 27 ERS-DEGs (Figure 2(c)).Figure 2
ERS-DEGs were obtained from the intersection of DEGs and ERS-related genes. (a) Heatmap of DEGs; (b) volcano plot of the distributions of DEGs. Red and blue dots represent upregulated and downregulated genes, respectively. No significantly changed genes were marked as grey dots; (c) Venn diagram of ERS-DEGs. Venn diagram analysis executed an intersection of DEGs and ERSGs to reveal the ERS-DEGs.
(a)(b)(c)


## 3.2. KEGG and GO Enrichment Analysis of ERS-DEGs

GO analysis results indicate that the upregulated ERS-DEGs were significantly enriched in BP pathways, including immune response-activating cell surface receptor signaling pathway, immune response-activating signal transduction, neutrophil degranulation, and neutrophil activation; all of which are involved in immune homeostasis. As for MF and CC analysis, the ERS-DEGs were mainly implicated with the external side of plasma membrane and MHC class II receptor activity, respectively (Figures3(a)–3(d)). The results are detailed in Table 1. KEGG pathway analysis demonstrated that the upregulated ERS-DEGs were predominantly enriched in immune-related pathways, including leishmaniasis, staphylococcus aureus infection, Th17 cell differentiation, systemic lupus erythematosus, phagosome, and tuberculosis (Figures 3(e) and 3(f)). These results are shown in detail in Table 2. The pathways that were most significantly enriched from ERS-DEGs in GO and KEGG analysis are shown in diagrams (Figures 3(g) and 3(h)).Figure 3
GO/KEGG enrichment analyses of ERS-DEGs. (a) Advanced bubble chart shows GO functional enrichment significance items of ERS-DEGs in three functional groups: BP, MF, and CC. Thex-axis and y-axis labels represent the gene ratios enriched in the GO analysis group and GO terms, respectively. The color of dots represents the adjusted p value: the redder the color, the lesser the adj. p value; the bluer the color, the greater the adj. p value. The size of dots represents gene counts. (b–d) Chord diagrams show the distribution of ERS-DEGs in the three different GO-enriched functions. The color of dots represents the log FC of the genes. (e) Advanced bubble chart shows enrichment of ERS-DEGs in signal pathways. (f) The chord diagram shows the distribution of ERS-DEGs in different KEGG pathways. (g) The pathway diagram shows the immune response-activating cell surface receptor signaling pathway with the highest enrichment degree in GO analysis. (h) The pathway diagram shows that the leishmaniasis pathway has the highest enrichment degree in KEGG analysis.
(a)(b)(c)(d)(e)(f)(g)(h)Table 1
GO enrichment analysis of ERS-DEGs.
OntologyIDDescriptionGeneRatioBgRatiop valuep.adjustq valueBPGO:0043312Neutrophil degranulation12/49485/186703.18e-091.67e-061.30e-06BPGO:0002283Neutrophil activation involved in immune response12/49488/186703.41e-091.67e-061.30e-06BPGO:0042119Neutrophil activation12/49498/186704.28e-091.67e-061.30e-06BPGO:0002446Neutrophil mediated immunity12/49499/186704.38e-091.67e-061.30e-06BPGO:0030098Lymphocyte differentiation9/49353/186702.94e-078.99e-057.00e-05CCGO:0009897External side of plasma membrane11/50393/197173.19e-094.09e-072.76e-07CCGO:0042581Specific granule7/50160/197171.51e-079.69e-066.53e-06CCGO:0070820Tertiary granule6/50164/197173.55e-061.51e-041.02e-04CCGO:1904724Tertiary granule lumen4/5055/197171.13e-053.54e-042.38e-04CCGO:0034774Secretory granule lumen7/50321/197171.56e-053.54e-042.38=-04MFGO:0019864IgG binding2/4911/176974.07e-040.0730.058Table 2
KEGG enrichment analysis of ERS-DEGs.
OntologyIDDescriptionGeneRatioBgRatiop valuep.adjustq valueKEGGhsa04640Hematopoietic cell lineage5/3299/80763.87e-050.0040.004KEGGhsa05140Leishmaniasis4/3277/80762.24e-040.0130.010KEGGhsa05340Primary immunodeficiency3/3238/80764.34e-040.0150.012KEGGhsa05150Staphylococcus aureus infection4/3296/80765.22e-040.0150.012KEGGhsa04659Th17 cell differentiation4/32107/80767.87e-040.0180.014


## 3.3. GSEA and GSVA

GSEA was performed to identify the biological pathways with significant alterations towards IS pathology. A total of 10 pathways were enriched by this method. Interestingly, five of them were closely associated with T cells, including T cytotoxic, T helper, IL17, NO2IL12, and CTL pathways (Figures4(a)–4(j)). The results are detailed in Table 3. Furthermore, GSVA was conducted to further explore the signal pathways in IS, including GNF2_CD33, GSE 3688 STAT5 AB knockin vs WT T cell IL 2 T-related 17h DN, and GSE29618_monocyte_vs_PDC_up (Figure 5).Figure 4
The top 10 significantly enriched pathways based on GSEA. (a) WP purine metabolism; (b) BioCarta T cytotoxic pathway; (c) BioCarta T helper pathway; (d) Reactome mRNA decay by 3′ to 5′ exoribonuclease; (e) BioCarta IL17 pathway; (f) Reactome transport of nucleoside-free purine and pyrimidine bases across the plasma membrane; (g) KEGG primary immunodeficiency; (h) BioCarta NO2IL12 pathway; (i) BioCarta CTL pathway; (j) Reactome mitochondrial calcium ion transport.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)Table 3
GSEA in IS.
IDDescriptionSet sizeEnrichment scoreNESp valuep.adjustq valuesRankGO:0031982Vesicle4850.1727146973.4112344380.0014326650.0041335450.001924525700GO:0005576Extracellular region4100.1749941963.2667012790.0014492750.0041335450.001924525435GO:0005739Mitochondrion2430.1422813612.2130932510.0014641290.0041335450.001924525983GO:0005615Extracellular space3290.190210293.2741988360.0014749260.0041335450.001924525252GO:0043230Extracellular organelle2630.220613933.5410205510.0014771050.0041335450.001924525621GO:0070062Extracellular exosome2610.2206988693.5249952460.0014771050.0041335450.001924525621GO:1903561Extracellular vesicle2630.220613933.5410205510.0014771050.0041335450.001924525621GO:0030054Cell junction2010.165794812.4084827940.0014814810.0041335450.001924525189GO:0031410Cytoplasmic vesicle3220.1457290152.487680470.0014814810.0041335450.001924525589GO:0097708Intracellular vesicle3220.1457290152.487680470.0014814810.0041335450.001924525589GO:0031967Organelle envelope2040.1868419672.7303292620.0014858840.0041335450.001924525788GO:0031975Envelope2040.1868419672.7303292620.0014858840.0041335450.001924525788GO:0030141Secretory granule1500.2257188682.8554146860.0014970060.0041335450.001924525404GO:0098796Membrane protein complex1500.1621314822.0510142490.0014970060.0041335450.001924525656GO:0099503Secretory vesicle1710.2220231932.9948720730.0015105740.0041335450.001924525404GO:0005740Mitochondrial envelope1370.1900353122.3069467260.0015290520.0041335450.001924525616GO:1990904Ribonucleoprotein complex1290.1796543652.1313439140.0015313940.0041335450.001924525902GO:0031966Mitochondrial membrane1300.1779346722.1169076690.0015337420.0041335450.001924525616GO:0019866Organelle inner membrane1000.244613892.6241194150.0015723270.0041335450.001924525609GO:0070161Anchoring junction1040.2515477142.7352353030.0015898250.0041335450.001924525149GO:0012506Vesicle membrane1210.1743482661.9974747970.003110420.0073979230.0034443771110GO:0030659Cytoplasmic vesicle membrane1180.1773064262.008090750.003129890.0073979230.0034443771110GO:0031090Organelle membrane4740.1024336962.0166298440.0043415340.0098156420.004570036656GO:0000228Nuclear chromosome129-0.149453948-1.9909431440.0057306590.0124164280.005780928647GO:1902494Catalytic complex2040.1251160271.8283255960.0133729570.027815750.0129506531192GO:0045202Synapse1220.1474493791.6933664670.0203442880.0406885760.018944074167GO:0000790Nuclear chromatin107-0.137337172-1.6669836440.0244565220.0471014490.021929825621Figure 5
Heatmap of 10 signal pathways of ERS-DEGs based on GSVA.


## 3.4. PPI Network Analysis

To identify the interactions of the ERS-DEGs, PPI network of ERS-DEGs was established and is shown in Figure6(a). The identified hub genes, predicted via cytoHubba plugin is also shown in the aforementioned figure. The top 10 hub genes in an increasing order of interaction degree are as follows: HIF1A, CREBBP, EP300, ARNT, TP53, PPARG, VHL, EPAS1, HIF1AN, and CREB1 (Figure 6(b)). The hub genes may be key ERS-related biomarkers of IS, which require further clinical trials for validation.Figure 6
Protein–protein interaction (PPI) networks. (a) PPI network analysis of ERS-DEGs. Red and blue colors represent high and low expression levels of genes, respectively. (b) The top 10 hub genes from the PPI network were predicted and shown on the subnetwork. Node color reflects the degree of accuracy in prediction with a gradient shift from yellow to red. Illustratively, the redder the color, the higher the accuracy of prediction.
(a)(b)


## 3.5. Network Analysis of Hub Genes

Transcriptional regulatory network analysis was performed to analyze the interaction between TFs and the hub genes (Figure7(a)). The top 3 hub genes associated with TFs were TP53, EPAS1, and VHL. Moreover, a network of miRNAs and hub genes was constructed (Figure 7(b)). The top 3 hub genes associated with miRNAs were TP53, CREB1, and HIF1A. As potential therapeutic targets, relative drugs and molecular compounds of the hub genes were also explored by the drug-gene interaction networks. As shown in Figure 7(c), a total of 7 drugs or molecular compounds were found to interact with hub genes of HIF 1A, CREB1, and HIF1AN. They are carvedilol, 2-methoxyestradiol, N-[(1-chloro-4-hydroxyisoquinolin-3-YL)carbonyl] glycine, naloxone, adenosine monophosphate, D-tartaric acid, and N-(carboxyvarbonyl)-D-phenylalanine.Figure 7
Network analyses of the hub genes. (a) Network of hub gene-TFs. Dots in blue to green represent the hub genes, and dots in red to yellow represent the TFs. (b) Network of hub gene-miRNA. Red and blue nodes are hub genes and miRNAs, respectively. (c) Networks of hub gene-related drug. Red and blue circle nodes are hub genes and drugs/molecular compounds, respectively.
(a)(b)(c)


## 3.6. Immune-Related Gene Identification and Functional Correlation Analysis

The result of the CIBERSORT analysis revealed that the T and NK cell subtypes were predominant among the 22 types of infiltrating immune cells (Figure8(a)). The correlation between 22 immune cell subtypes was also performed and presented via heatmap. The result revealed that M0 macrophages and resting NK cells had the strongest and weakest correlation with other immune cells, respectively (Figure 8(b)). Furthermore, the hub genes, ARNT, CREB1, CREBBP, EP300, EPAS1, HIF 1A, HIF 1AN, and TP53, were identified, and they correlated with 11 different types of immune cells (Figure 8(c)). We found a significant correlation between the number of neutrophil and T cell subtypes and the expression level of hub genes. The results suggested that neutrophils and T cells play a critical role in the development of ERS following IS.Figure 8
The nose cone is detachable upon impact. Evaluation and visualization of infiltrating immune cells. (a) Bar diagram displayed the proportion of infiltrating immune cells based on CIBERSORT algorithm. (b) Heatmap of the correlation between 22 subtypes of immune cells. Blue and red areas indicate positive and negative correlations, respectively. The darker the color, the stronger the correlation. (c) The correlations between hub genes and immune cells.
(a)(b)(c)


## 4. Discussion

IS is the primary cause of mortality and disability among Chinese residents, and its morbidity has increased over the last decade, thus, posing a huge burden to national medical and economic systems [24]. Neuronal apoptosis is the key etiological factor of neuronal damage following cerebral ischemia. Activated ERS has been reported to play a critical role in inducing neuronal cell apoptosis in IS [25]. Therefore, exploring potential biomarkers related to ERS is indispensable for the elucidation of neuronal damage following IS. However, ERS biomarkers in IS have not yet been fully identified. In this study, we comprehensively analyzed two mRNA microarray datasets (GSE16561 and GSE37587), including 107 and 24 whole-blood samples from reliable profiles of IS and healthy controls, respectively. A total of 60 DEGs were distinguished from the two datasets, including 38 and 22 downregulated and upregulated genes, respectively. Subsequently, the intersection analysis of DEGs and ERSGs was performed to obtain 27 ERS-DEGs. Enrichment analysis revealed that the modules and pathways ERS-DEGs enriched in were closely related to immunity. The data suggest that local inflammation may be involved in the development of IS. Moreover, we constructed a PPI network, and 10 hub genes with the highest scores in ERS-DEGs were identified. Our results suggested that the identified hub genes may interact with each other in the network and are also associated with immune cell infiltration.ERS is a pathological condition related to hypoxia, starvation, calcium imbalance, and free radical excess. Contrarily, ERS can trigger the unfolded protein response, which restores homeostasis by decreasing protein translation and upregulation of ER chaperone gene expression [26]. Interestingly, unregulated ERS can initiate neuronal apoptosis by activating apoptosis signal pathways directly and also further exacerbates neuronal damage through activating inflammatory signal pathways. Recent studies have reported that ERS is a major cause of neuronal apoptosis [27]. A previous study has revealed that the promotion of ERS via knockdown Hes1 induced vast neuronal apoptosis in a mouse model with middle cerebral artery occlusion (MCAO) [25]. In addition, Li et al. reported that the inhibition of ERS signal pathway by γ-glutamylcysteine attenuated neuronal apoptosis in the ischemic brain of mice [28]. The aforementioned results support the idea that ERS may be a potential target in IS injury. Therefore, exploring biomarkers of ERS in blood samples of patients with IS is necessary for early diagnosis and treatment. Furthermore, it may aid the identification of novel therapeutic targets.To investigate the involvement of ERS-DEGs in the biological mechanism of IS, GO and KEGG enrichment analyses were performed, and the results revealed that ERS-DEGs are mainly enriched in immune-related pathways. In the BP annotations, neutrophil degranulation and neutrophil activation involved in immune response were found to be significantly related to the development of IS. Neutrophils are the first inflammatory responders to be recruited to the site of cerebral infarction after IS [29]. Enhanced ERS was suggested to be associated with severe neutrophil inflammation [30]. Notably, activated neutrophils release neutrophil argininase-1, which induces cell apoptosis through the ERS pathway [31]. Clinical researches have established that neutrophil overexpression is a positive indicator of stroke progression [32]. Further research has shown that inhibition of neutrophil infiltration into ischemic lesions decreases infarct size and mitigates stroke pathology [33]. These results are similar to our predictions.Similarly, GSEA and GSVA enrichment analyses also revealed that ERS-DEGs were significantly enriched in T cell-related pathways, including T cytotoxic, T helper, IL17, NO2IL12, and CTL pathways. The results above indicate that T cells play an important role in IS. Few researchers have observed significant increases in CD4+ and CD8+ T cells in the peri-infarct area one month following IS experiments [34]. Persistent neuronal ERS has been reported to promote CD4+ and CD8+ T cell priming by inducing activation of STING signal cascade [35]. Several studies have reported that CD8+ T cells could induce neuronal damage directly by the cytotoxic function, while CD4+ T cells could aggravate the local inflammation by activating the effector T cell. Researches supported that T cell-targeted therapy could effectively reduce the infarct size of brain tissues by depleting T cells [34, 36]. In summary, our results further confirmed that T cell responses associated with ERS play an integral role in post-IS neuroinflammation.In addition, among the 20 KEGG-enriched signaling pathways, we observed that the pathway with high enrichment was Th17 cell differentiation. In previous studies, Th17 cells were observed to be differentiated from CD4+ T cell in response to TGF-β and IL-6 and can as well produce proinflammatory cytokines. However, the unstable and plastic Th17 cells can also be transdifferentiated into Treg cells, which can attenuate inflammation [37]. Hence, T cells have a dual role in the inflammatory response after stroke. ERS was reported one of the mechanisms that deregulate Th17/Treg cells [38]. Therefore, suppressing ERS may be the effective therapeutic approach that regulated the balance between proinflammatory and anti-inflammatory mechanisms. In summary, ERS is indispensable for inflammatory responses, and our results are consistent with these findings.PPI molecular interaction network was employed to analyze the protein interaction of ERS-DEGs. Hypoxia inducible factor 1 subunit alpha (HIF1A) was identified as one of the top hub genes in this study, which is an important transcription factor in the hypoxic or ischemic brain and aids cells in the adaptation to hypoxic conditions [39]. During hypoxia, the upregulation of HIF1A can induce glycolysis and inhibit mitochondrial respiration to decrease oxygen consumption, thus, improving cell proliferation. Researches support the important role of HIF1A as a neuroprotector in ischemic brain injury [40]. Furthermore, other members in the hypoxia-inducible factor family, including EPAS1 (HIF2A) and ARNT (HIF1β) genes, were also identified as hub genes. Notably, cAMP response element-binding protein 1 (CREB1), was also reported as one of the hub genes that are related to neuronal survival after ischemia. CREB-related therapeutics have extensive application prospects in cerebral protection after IS [41]. Although the CREB1 inhibitor, naloxone, is conventionally used for the treatment of opioid addiction, its neurorestorative effects have been observed in patients with IS and MCAO mouse models [42]. Hence, our results indicate that the hypoxia-inducible factor family and CREB family genes are potential therapeutic targets and biomarkers of ERS in IS.Our results indicate that CREB binding protein (CREBBP) is the top hub gene in the PPI network analysis. It has been reported to have an increased expression after IS and to be involved in neurovascular reconstruction and neuron protection [43]. Moreover, persistent ERS can promote the stability of CREBBP [44]. These findings are consistent with our results, which suggest that CREBBP may be a potential biomarker of ERS in IS. An analogue of CREBBP, EP300, was also identified as a hub gene in the PPI network. It participates in mediating the ubiquitination and degradation of CHOP protein and prevents ERS-induced apoptosis under hypoxic conditions [45]. This hub gene has been previously reported as a potential biomarker of IS, which is consistent with our report. However, the function of EP300 in IS pathology has not yet been elucidated [46]. Therefore, as a potential biomarker of IS, future investigation of EP300 is worthwhile.In the correlation analysis between hub genes and TF or hub genes and miRNA, TP53 was found to have rich networks. Previous reports have shown that ERS can induce significant upregulation of TP53 [47]. Rapid and substantial accumulation of TP53 was observed in the ischemic penumbra, which resulted in neuronal apoptosis [48]. Therefore, TP53 may be a potential therapeutic target of IS. Interestingly, it has been reported that the inhibition of TP53 can considerably decrease infarct size [49]. Our results are consistent with the aforementioned reports, indicating that TP53 may be a potential biomarker of IS. Further investigation of the molecular mechanisms of hub genes in IS may provide novel innovations for the diagnosis and treatment of IS.However, there are several limitations to this study. First, the results were not backed up with experimentsin vitro and in vivo, to further verify the biological functions of these potential ERSGs after IS. Secondly, the research was not able to evaluate the correlation between hub genes and severity of pathology in patients with IS, due to a lack of corresponding clinical studies. We also observed that the batch differences cannot be avoided and removed by the analysis, due to high-content datasets. Thus, a larger sample size and clinical study setting are required for further studies, in order to strengthen the statistical power and obtain more reliable results. Related experimental evidence is also required to fully elucidate the role of the hub genes at the molecular cellular level in IS. In addition, considering that protein biomarkers are more reliable and applicable than gene biomarkers in clinical practices, proteomics data is required for further validation of our results.


## 5. Conclusions

In conclusion, we performed a comprehensive analysis of the differentially expressed genes associated with ERS following IS, by combining two datasets (GSE16561 and GSE37587) to explore the molecular mechanism and biomarkers of ERS, which are involved in neuronal damage of IS. ERS-DEGs and hub genes were identified, which may be potential therapeutic targets and biomarkers of IS. Enrichment analysis indicated that the immune-related signal pathway is an assignable mechanism in the progression of IS. Network and immune correlation analyses further confirmed that the identified hub genes are closely related to immune infiltration. However, further clinical studies should be conducted in order to validate our findings.



---

*Source: 4212180-2022-08-25.xml*