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

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.


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 bloodbased biomarkers.
The aim of the present IS study is to explore the blood biomarkers associated with ERS by utilizing the wholeblood 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 ERSrelated biomarkers that may serve as potential determinants for further diagnosis and treatment of IS.

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.

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 |log 2 FC | >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.

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.

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 a false discovery rate ðFDRÞ < 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.

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) 2 Computational and Mathematical Methods in Medicine   2.2   GSM922887  GSM922890  GSM922893  GSM922896  GSM922899  GSM922902  GSM922905  GSM922908  GSM922911  GSM922914  GSM922917  GSM922920  GSM922923  GSM922926  GSM922929  GSM922932  GSM922935  GSM922938  GSM922941  GSM922944  GSM922947  GSM922950  GSM922953  3100191_Stroke  3100138_Stroke  3100058_Stroke  3100039_Stroke  3100031_Stroke  3100114_Stroke  3100118_Stroke  3100042_Stroke  3100077_Stroke  3100036_Stroke  3100061_Stroke  3100126_Stroke  3100084_Stroke  212414_Control  213403_Control  213906_Control  213406_Control  212411_Control  212404_Control  212409_Control     3 Computational and Mathematical Methods in Medicine 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.

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 cutoff 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.

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 (Figures 1(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.
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 (Figures 2(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)).

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 (Figures 3(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        (Figures 4(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).

PPI Network Analysis.
To identify the interactions of the ERS-DEGs, PPI network of ERS-DEGs was established and is shown in Figure 6(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.
3.5. Network Analysis of Hub Genes. Transcriptional regulatory network analysis was performed to analyze the interaction between TFs and the hub genes (Figure 7(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.
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 (Figure 8(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.

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

13
Computational and Mathematical Methods in Medicine 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 experi-ments [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, 14 Computational and Mathematical Methods in Medicine  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. CREBrelated 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 ERSinduced 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 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. 16 Computational and Mathematical Methods in Medicine 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 experiments in 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.

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 bio-markers 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   19 Computational and Mathematical Methods in Medicine 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.

Data Availability
The datasets used and/or analyzed during the current study are available from the corresponding authors on reasonable request.