The overall workflow of this study is shown in Fig. 1. Briefly, the GSE23163 dataset was collected from the GEO database (http://www.ncbi.nlm.nih.gov/geo) and converted into a suitable format for analysis, and differentially expressed genes (DEGs) were identified. Subsequently, the dataset was subjected to WGCNA to delineate the modular gene nest most associated with CIR, followed by the immune infiltration analysis. Then the pivotal genes were identified by intersecting the prepared immune gene profiles, DEGs, and modular genes. The differential expression of pivotal genes was verified after successfully constructing a CIR animal model. Afterward, a protein interaction network (modeling of microRNA-DEGs and transcription factor (TF)-DEGs) was constructed. Finally, a molecular docking simulation of potential small drug molecules to treat CIR was performed.

Fig. 1.

Study workflow. This study includes seven key stages. In phase I, an appropriate dataset was collected and converted into a format suitable for analysis. In phase II, DEGs were delineated in the target dataset. In phase III, the modular gene nest most associated with CIR was targeted by WGCNA. In Phase IV, the dataset was subjected to immune infiltration analysis. Subsequently, the prepared immune gene profiles, DEGs, and modular genes were intersected to obtain pivotal genes. In Phase V, the differential expression of pivotal genes was verified after successfully constructing a CIR animal model. In Phase VI, a protein interaction network (modelling of miRNA-DEGs and TF-DEGs) was constructed. Finally, in Phase VII, the molecular docking simulation of potential small drug molecules to treat CIR was performed. DEGs, differentially expressed genes; miRNA, microRNA; GEO, Gene Expression Omnibus; GO, Gene Ontology; PPI, protein–protein interaction; cMap, Connectivity Map; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, Gene Set Enrichment Analysis; WGCNA, Weight gene co-expression network analysis; qPCR, quantitative real-time PCR; IHC, immunohistochemistry; IFC, immunofluorescence staining; WB, Western blotting.

The gene expression profiling dataset GSE23163 (annotation platform GPL6885) contained 48 CIR and 16 sham surgery samples. The downloaded raw data were converted into gene names with the assistance of probes. Then, the sequencing data were log transformed and normalized using the “sva” package into expression datasets that could be directly analyzed.

The DEGs between the CIR and sham samples were determined using the “limma” package. The p-values lower than 0.05 and logFC absolute values $>$1 were defined as significant, and the magnitude of logFC values and p-values were visualized hierarchically as volcano plots and heat maps. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed on CIR and normal samples using the GSEA 4.2.3 software (Cambridge, MA, USA) http://www.gsea-msigdb.org/gsea/index.jsp. The random combination was set at 1000 times. To investigate the CIR-related risk signature’s molecular characteristics and pathways, we conducted GSEA to display differences in 159 gene sets and 186 KEGG gene sets. Gene sets that had enrichment scores with a normalized p-value 0.05 and a false discovery rate of 0.25 were considered significantly enriched.

PPIs are a target of cell biology research and a prerequisite for systems biology studies. PPI networks are defined as graphs where nodes and edges represent proteins and their interactions, respectively [25, 26]. To investigate the potential molecular mechanisms of CIR from the perspective of protein interactions, a PPI network of DEGs was constructed using the STRING database (https://string-db.org/) and visualized using Cytoscape (https://cytoscape.org/) [27]. The “cytohubba” (sorting by “degree” size) and “MCODE” (basic cut-off, criterion degree cutoff = 2, node score cutoff = 0.2, k-core = 2, maximum depth = 100) plugins were used to further validate the plausibility of the hub genes [28]. TF-miRNA co-regulatory interactions were collected from the RegNetwork repository (http://www.regnetworkweb.org) to detect miRNAs and TFs regulating the DEG of interest at the post-transcriptional and transcriptional levels. The DEG-miRNA network was validated using the experimentally validated miRNA-target interactions databases TarBase and miRTarBase, based on the selected genes, and checked on the platform. The DEG-miRNA target genes were validated using the JASPAR (https://jaspar.genereg.net/) and GEPIA (http://gepia.cancer-pku.cn/) databases.

The relevant network was built using the “WGCNA” package to identify the key molecules of CIR. The main steps were: (i) A power of $\beta$ = 7 was selected according to the scale-free topology criterion; (ii) Hierarchical clustering trees were constructed based on the weighted correlation coefficients of the genes, and the related modules were merged; (iii) Identification of hub gene clusters based on module-phenotype correlations; (iv) GO and KEGG enrichment analyses according to the modular gene ensemble.

CIBERSORT is a gene expression profiling-based method used for quantifying immune cell infiltration in tissues. The CIBERSORT analysis is based on LM22, the annotated gene signature matrix comprising 22 functionally defined immune cell subtypes. The 22 immune cell types included macrophage subsets (M0, M1, and M2), T cells (CD8+ T cells, naïve CD4+ T cells, memory resting CD4+ T cells, memory activated CD4+ T cells, Tfh cells, regulatory T cells, and $\gamma$$\delta$ T cells), natural killer cells (resting and activated), mast cells (resting and activated), B cells (naive and memory), dendritic cells (resting and activated), monocytes, plasma cells, neutrophils, and eosinophils. In the present study, we downloaded the LM22 gene matrix from the CIBERSORT portal (http://cibersort.stanford.edu/). To improve the accuracy of the inverse fold product algorithm, the p-value and root mean square error of CIBERSORT were counted for each sample file.

The sequences of potential hub genes corresponding to TF-binding sites were obtained from the UCSC database (https://genome.ucsc.edu/). The sequence information was then imported into JASPAR to retrieve possible TFs. Candidates that did not meet the inclusion criteria (p 0.000001) were eliminated based on their transcriptional direction and binding power. Finally, expression correlation was verified using the GEPIA database.

The molecular docking simulations were performed using the following steps. First, the up- and downregulated DEGs were introduced into the clue.io platform to screen the therapeutic molecules. Then, the two- and three-dimensional structures of the compounds were extracted from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). Afterward, the structure of the protein receptor encoded by the target gene was elucidated by screening the molecules in the UniProt (https://www.uniprot.org/) and PDB databases (https://www.rcsb.org/). Finally, CHEM 3D (PerkinElmer, Waltham, MA, USA) and AutoDock Vina-1.5.7 (The Scripps Research Institute, SanDiego, CA, USA) were used for the spatial processing of ligands and receptors, and PYMOL (http://www.pymol.org/pymol) was used for docking and annotating the drug and protein molecules.

All animal experiments were performed in accordance with protocols approved by the Ethics Committee of Nanjing Medical University and complied with the Guidelines of Laboratory Animals for Biomedical Research published by the National Institutes of Health (NIH publication, revised in 2011).

Eight-week-old male C57BL/6J mice (n = 12) weighing 20–25 g were obtained from Beijing Weitong Lihua Company (Beijing, China). The mice were housed under a 12-hour light/dark cycle with free access to food and water. The MCAO/R model was established as described previously, with minor modifications [29]. Briefly, mice were anesthetized with 2% isoflurane and placed on a heating pad to maintain body temperature at 37 °C. The left common carotid and external carotid arteries were exposed and ligated permanently to block blood flow. A small incision was made in the common carotid artery, from which a monofilament nylon suture (final tip diameter 0.20 $\pm$ 0.01 mm, Cat.1620A4, Beijing Cinontech Co. Ltd, Beijing, China) was inserted and advanced towards the internal carotid artery. The internal carotid artery was occluded at the level of the middle cerebral artery branches for 60 min. The suture was then carefully removed to allow reperfusion for 24 h. Mice with excessive bleeding or that died within 24 h after reperfusion were excluded from the study data.

The neurological deficit score of ischemic mice was evaluated 24 h after reperfusion based on the five-point scale model established previously [30]. Score 0: no neurological deficit. Score 1: failure to extend the contralateral forelimb. Score 2: circling the contralateral side. Score 3: falling to the contralateral side. Score 4: no spontaneous motor activity or death.

The mice were euthanized 24 h after reperfusion, and their brains were quickly removed. Five 2 mm-thick consecutive coronal slices were made and incubated in 2% 2,3,5-Triphenyl tetrazolium chloride (Cat.T8877, Sigma-Aldrich, St. Louis, MO, USA) solution at 37 °C for 30 min in the dark. Afterward, the slices were imaged using a digital camera, and the infarct area of each slice was measured using ImageJ software.1.8.0 (National Institutes of Health, Bethesda, MD, USA). The percentage of each cerebral infarct volume was calculated according to the following formula: total cerebral infarct volume/total brain volume $\times$ 100.

Total RNA from the cortex tissues of the removed brains was isolated using TRIzol reagent (Cat.9208, Takara, Takara Bio, Inc., Otsu, Shiga, Japan), following the manufacturer’s protocol. Total RNA (500 ng) was reversely transcribed into cDNA using PrimeScriptTM RT Master Mix (Cat. RR036A; Takara, Japan). Real-time PCR was performed on the system using SYBR Green dye (Cat.11199ES03, Yeasen, Shanghai, China). The following primer pairs for murine were used: C-C motif chemokine ligand 3 (Ccl3; Mus): 5${}^{\prime }$-CCAGCCAGGTGTCATTTT-3${}^{\prime }$, 5${}^{\prime }$-GCATTCAGTTCCAGGTCAG-3${}^{\prime }$; Ccl4 (Mus):5${}^{\prime }$-AACCTAACCCCGAGCAA-3${}^{\prime }$, 5${}^{\prime }$-AAACAGCAGGAAGTGGGA-3${}^{\prime }$; Ccl7 (Mus):5${}^{\prime }$-CAGAAGTGGGTCGAGGAG-3${}^{\prime }$, 5${}^{\prime }$-AAGAACAGCGGTGAGGAA-3${}^{\prime }$; CXC motif chemokine ligand 1 (Cxcl1;(Mus):5${}^{\prime }$-AACCGAAGTCATAGCCACA-3${}^{\prime }$, 5${}^{\prime }$-GGGGACACCTTTTAGCATC-3${}^{\prime }$; Cd14 (Mus):5${}^{\prime }$-CTGTCTTTCACTGGGCTGA-3${}^{\prime }$, 5${}^{\prime }$-GCTCATCTGGGCTAGGG-3${}^{\prime }$; Suppressor of cytokine signaling (Socs3; Mus):5${}^{\prime }$-GTCGGTCGCCTCTCCTC-3${}^{\prime }$, 5${}^{\prime }$-ATCTGCCCCTGGCTCAC-3${}^{\prime }$; Hspa1a (Mus):5${}^{\prime }$-GCTTGTGTCGGGTCCTTC-3${}^{\prime }$, 5${}^{\prime }$-TGCTGTCACTTCACCTCCA-3${}^{\prime }$; Gapdh (Mus):5${}^{\prime }$-AGGTCGGTGTGAACGGATTTG-3${}^{\prime }$, 5${}^{\prime }$-TGTAGACCATGTAGTTGAGGTCA-3${}^{\prime }$.

Total protein was extracted from tissues according to the instructions for Radioimmunoprecipitation Assay Lysis Buffer (Beyotime, Shanghai, China). Next, total protein concentration was measured using the BCA assay (Beyotime, China). Total proteins were separated by Sodium Dodecyl Sulfate Polyacrylamide Gel Electrophoresis on 12.5% polyacrylamide gels and electrophoretically transferred to PVDF membranes, which were subsequently sealed with 5% BSA. The membranes were incubated overnight at 4 ${}^{\circ }$C with diluted primary antibodies for Ccl3 (1:1000 dilution, Affinity Biosciences, Cincinnati, OH, USA), Ccl4 (1:1,000 dilution, Affinity Biosciences, USA), Ccl7 (1:1,000 dilution, Affinity Biosciences, USA) Cxcl1 (1:1000 dilution, Affinity Biosciences, USA), Hapa1a (1:1000 dilution, Wanleibio, Shenyang, China), Cd14 (1:1000 dilution, Wanleibio, China), and Socs3 (1:1000 dilution, Wanleibio, China). The membranes were then further incubated with secondary antibody (ZSGB-BIO, Beijing, China) for 1 hour at room temperature. Finally, the protein bands were analyzed and quantified using ImageJ software (National Institutes of Health, Bethesda, MD, USA).

The Human Protein Atlas Database (HPAD; https://www.proteinatlas.org/) was used to verify the immunohistochemistry of brain tissues and fluorescence staining of cells. This database is helpful for systematically studying the transcription and translation levels of coding genes in different tissue types. The localization of target genes in tissues and cell protein expression was further analyzed.

A total of 25 DEGs associated with CIR were identified in the GSE23163 dataset, of which 23 were upregulated, and 2 were downregulated. Genes with a $|$logFC$|$ of $>$2 were additionally labeled in the volcano plot, namely CCl4 and Timp1. The difference in expression between the CIR and sham groups was also observed in the heatmap (Fig. 2).

Fig. 2.

Identification of DEGs in CIR. Graphical visualization of the results of the variance analysis. (A,B) Gene expression profiles of GSE23163 are visualized in (A) volcano plots and (B) heat maps. DEGs, differentially expressed genes; CIR, cerebral ischemia-reperfusion injury; Con, Control group; Rep, Ischemia-reperfusion group.

Based on WGCNA, the gene nests were divided into six modules, each with a specific color. The correlation coefficients of the modules were 0.69 (MEblue, p 0.001), 0.76 (MEbrown, p 0.001), –0.69 (MEturquoise, p 0.001), –0.63 (MEgreen, p 0.001), –0.65 (MEyellow, p 0.001), and 0.8 (MEgrey, p 0.001). The MEbrown module comprising 47 genes contained the most relevant gene set for CIR (Fig. 3).

Fig. 3.

Gene co-expression modules. Screening for the gene module that best matches the trait. (A) Scale-free index and mean connectivity analyses for various soft-threshold powers. (B) Hierarchical cluster dendrogram of CIR-related genes based on one dissimilarity measure. The color band shows the results obtained from the automatic single-block analysis. (C) Module–trait relationships in CIR; each cell contains the corresponding correlation and p-value. The value between –1 and 1 represents the correlation between the module and clinical features.

Box line plot demonstrated significant infiltration in 11 of the 22 immune cells in CIR mice compared to that in normal control mice. However, medium B cell memory, T cell CD8, T cell CD4 naïve, T cell CD4 memory resting, T cell regulatory Tregs, monocytes, macrophages M2, dendritic cells resting, activated dendritic cells, eosinophils, and neutrophils were not infiltrated (Fig. 4). The intersection of the immune gene dataset (containing 1793 known immune-related genes), DEGs (n = 25), and MEbrown module genes (n = 47 phenotypically related genes) identified seven pivotal genes, including Ccl3, Ccl4, Ccl7, Cxcl1, Hspa1a, Cd14, and Socs3.

Fig. 4.

Immune infiltration analysis and intersection dataset. To accurately evaluate the composition of immune cells and the establishment of key molecules in the microenvironment after ischemia reperfusion. (A) Venn diagram of the intersection of the brown module (MEbrown), immune genes, and DEGs. (B) The proportion of 22 types of immune cells corresponding to immune infiltration analysis. *p 0.05, **p 0.01, ****p 0.0001.

GO analysis (Fig. 5) demonstrated that the genes in the MEbrown module were involved in biological processes such as “cellular response to biotic stimulus”, “regulation of ERK1 and ERK2 cascade”, “regulation of inflammatory response”, “response to chemokines”, “cellular response to chemokines”, and “cell junction assembly”. The cellular component GO terms of these genes were enriched in “membrane rafts”, “membrane microdomains”, “receptor complexes”, “stress fibers”, and “contractile actin filament bundles”. The molecular function terms were mainly enriched in “chemokine activity”, “protein tyrosine/threonine phosphatase activity”, “protein folding chaperone”, and “transmembrane receptor protein tyrosine kinase activity”.

Fig. 5.

Enrichment analyses of gene modules. Enrichment analysis of individual modules distinguished by WGCNA, including blue, brown, green, turquoise and yellow. Enrichment analyses of gene modules for (A) biological process; (B) cell component; (C) molecular function; and (D) KEGG pathway. KEGG, Kyoto Encyclopedia of Genes and Genomes.

KEGG analysis identified the enrichment of the pathways involved in “stress fibers”, “membrane rafts”, “membrane microdomains”, “membrane regions”, and “receptor complexes”. The enrichment scores for CIR mice vs control mice were obtained by screening the genes based on their $|$normalized enrichment scores (NES)$|$ $>$0.7, p 0.05, and q 0.25). The GSEA (Supplementary Fig. 1) identified enrichment of biological processes terms related to “regulation of acute inflammatory response”, “regulation of eosinophils chemotaxis”, “regulation of eosinophils migration”, “regulation of lymphocyte chemotaxis”, and “positive regulation of interleukin 1 $\beta$ production”. The identified cellular components terms were “clathrin-coated vesicle”, “cytochrome complex”, “NADH dehydrogenase complex”, and “phagophore assembly site”. The molecular function terms were “ccr chemokine receptor binding”, “chemokine activity”, “chemokine receptor binding”, and “NADH+ nucleosidase activity”. The KEGG-enriched pathways were related to “type 1 diabetes mellitus”, “systemic lupus erythematosus”, “the nod-like receptor signaling pathway”, and “Leishmania infection”. These results indicated a considerable similarity between the results of the enrichment analysis of the MEbrown module and gene set, supporting the reliability of the results (Table 1).

Mining the binding sites for TFs using the UCSC database, followed by the screening of TFs using the JASPAR database, identified seven significant TFs. These TFs included JUNB, BACH2, TBP, ZNF384, RELA, PATZ1, ZNF460, and ZNF24. Of these, JUNB (R = 0.51, p 0.001), TBP (R = 0.54, p 0.001), ZNF384 (R = 0.62, p 0.001), and ZNF460 (R = 0.5, p 0.001) were strongly correlated with gene expression in the myocardial tissue (Fig. 6).

Fig. 6.

Binding sites of transcription factors and expression correlation verification. After initial prediction of potential regulatory molecules, expression correlation of transcription factors and regulated genes is verified in public databases. (A) Binding site of JUNB. (B) Binding site of TBP. (C) Binding site of ZNF384. (D) Binding site of ZNF460. (E–H) Correlations of the expression levels of (E) JunB and Ccl4; (F) Tbp and Ccl7; (G) Znf384 and Cd14; (H) and Znf460 and Socs3 genes in myocardial tissue.

The constructed PPI network contained 21 nodes and 61 connection lines. We identified the central genes according to the “degree” method and found that Cxcl1 had the highest linkage value, followed by Ccl3, Ccl4, Ccl7, Cd14, Socs3, and Timp1. These results are concordant with the results obtained by the intersection of the modules and immune infiltration data. Moreover, the set of hub genes identified by the cytohubba and MCODE plugins also contained these seven genes (Fig. 7).

Fig. 7.

PPI and TF-gene-miRNA networks. Constructing and visualizing regulatory networks through Cytoscape. (A) PPI network constructed using identified DEGs. Nodes represent DEGs, and edges represent the connection between DEGs. The depth of the color and the size of the circle reflects the degree of connectivity. (B) The pale-yellow nodes represent genes, the orange nodes represent TFs, the blue nodes represent miRNA, and the connecting lines represent the regulatory relationship between them. PPI, protein-protein interaction; TF, transcription factor; miRNA, microRNA; DEGs, differentially expressed genes.

Additionally, we determined the network of action between TFs, miRNAs, and hub genes. The node color and size were based on the strength of the relationship in the TF database. The TF-miRNA co-regulatory network comprised 133 nodes and 167 edges, where the nodes comprised seven hub genes, 72 miRNAs, and 54 TFs (Fig. 7).

The triggering and antagonistic drug molecules were identified based on the similarity between the expression profiles of DEG and small molecules. The closer the score value is to 100, the more similar the DEG is to the small-molecule treatment record, and the closer it is to –100, the more dissimilar the DEG is to the small-molecule treatment record. Chaetocin, auranofin, menadione, calyculin, rilmenidine, guanabenz, escitalopram, and nifedipine were the highest-scoring compounds (Supplementary Table 1). Nifedipine, the most likely CIR antagonist, successfully molecularly docked with proteins encoded by Ccl3, Ccl4, Ccl7, Cd14, Cxcl1, Hspa1a, and Socs3, with molecular binding energies of –5.4, –7.7, –4.2, –5.0, –6.0, –4.7, and –5.3 kcal/mol, respectively (Fig. 8).

Fig. 8.

Molecular docking simulations with nifedipine. Molecular docking of the screened small molecules for possible therapeutic use with the proteins encoded by each target gene, and if a stable binding site can be formed, it indicates a possible therapeutic effect. (A) CCL3. (B) CCL4. (C) CCL7. (D) CD14. (E) CXCL1. (F) HSPA1A. (G) SOCS3.

The representative 2,3,5-triphenyl tetrazolium chloride staining images, infarct volume percentage, and neurological deficit score 24 h after MCAO/R are shown in Fig. 9. The mean percentage of cerebral infarct volume was 33.12%, and the average neurological deficit score was 2.33 in the MCAO/R group. The real-time PCR revealed remarkably higher expression of all DEGs (including Ccl3, Ccl4, Ccl7, Cxcl1, Cd14, Socs3, and Hspa1a) in the cortex tissues of the MCAO/R group than that in the cortex tissues of the sham group. Each DEG showed a significant statistical difference with p 0.0001. In vivo validation results were consistent with those of our previous bioinformatic analysis (Fig. 9 and Supplementary Fig. 2). At the translational level, Western blot (WB) was used to detect the changed proteins of target protein expression. Compared with the normal group, the expression levels of Ccl3, Ccl4, Ccl7, Hspa1a, and Socs3 were significantly changed in the CIR group; however, there were no significant differences in the expressions of Cxcl1 and Cd14 (Fig. 10). In addition, the immunohistochemistry results of CD14, SOCS3, and HSPA1A in human brain tissues and their localization in brain cells in the database clearly distinguished their expression levels and distribution in the nucleoplasm (Supplementary Fig. 3).

Fig. 9.

Evaluation of animal models and real-time PCR results of hub gene expression validation. Tissues from successfully constructed ischemia-reperfused mice were subjected to RNA extraction, reverse transcription, and PCR to verify the differences in expression of individual genes. (A) The schematic representation of the model used in this study. (B) 2,3,5-Triphenyltetrazolium chloride staining images of brain sections. (C) proportion of infarct volume (n = 6), and (D) neurological deficit scores 24 h after middle cerebral artery occlusion and reperfusion (n = 6). (E) Expression of hub genes estimated using real-time PCR (n = 6, The data are shown as the means $\pm$ SEM, ****p 0.0001, unpaired t-test).

Fig. 10.

Western blotting results of hub gene expression validation. Genes that reflect differential expression in PCR validation are validated at the translational level and presented quantitatively. (A) The expression levels of CCL3, CCL7, HSPA1A, CCL4, SOCS3, CD14 and CXCL1 were assessed by western blotting. (B) The protein levels were quantified by normalizing to the GAPDH levels (n = 6, The data are shown as the means $\pm$ SEM, **p 0.01, ***p 0.001, ****p 0.0001, ns, no significance, unpaired t-test).