- Research article
- Open access
- Published:
Monocyte, neutrophil, and whole blood transcriptome dynamics following ischemic stroke
BMC Medicine volume 21, Article number: 65 (2023)
Abstract
Background
After ischemic stroke (IS), peripheral leukocytes infiltrate the damaged region and modulate the response to injury. Peripheral blood cells display distinctive gene expression signatures post-IS and these transcriptional programs reflect changes in immune responses to IS. Dissecting the temporal dynamics of gene expression after IS improves our understanding of immune and clotting responses at the molecular and cellular level that are involved in acute brain injury and may assist with time-targeted, cell-specific therapy.
Methods
The transcriptomic profiles from peripheral monocytes, neutrophils, and whole blood from 38 ischemic stroke patients and 18 controls were analyzed with RNA-seq as a function of time and etiology after stroke. Differential expression analyses were performed at 0–24 h, 24–48 h, and >48 h following stroke.
Results
Unique patterns of temporal gene expression and pathways were distinguished for monocytes, neutrophils, and whole blood with enrichment of interleukin signaling pathways for different time points and stroke etiologies. Compared to control subjects, gene expression was generally upregulated in neutrophils and generally downregulated in monocytes over all times for cardioembolic, large vessel, and small vessel strokes. Self-organizing maps identified gene clusters with similar trajectories of gene expression over time for different stroke causes and sample types. Weighted Gene Co-expression Network Analyses identified modules of co-expressed genes that significantly varied with time after stroke and included hub genes of immunoglobulin genes in whole blood.
Conclusions
Altogether, the identified genes and pathways are critical for understanding how the immune and clotting systems change over time after stroke. This study identifies potential time- and cell-specific biomarkers and treatment targets.
Background
Ischemic stroke (IS) is one of the leading causes of death and disability in the world. Brain injury follows arterial occlusions in large or small cerebral vessels. These may arise due to several different causes that ultimately deprive the tissue of necessary oxygen and glucose [1]. Effective treatments are limited to short time windows and access to stroke centers. Early diagnosis is paramount for best outcomes. Current differential diagnosis usually requires advanced brain imaging. Therefore, there is a need for tests that utilize reliable and accurate molecular biomarkers from blood.
The immune and clotting systems play critical roles in the injury and recovery from stroke. After IS, peripheral leukocytes, including monocytes and neutrophils, infiltrate the injured area, mediating the immune response that causes inflammation and subsequent resolution and repair [2, 3]. Monocytes are composed of different subsets: classical (pro-inflammatory CD14++ CD16−), intermediate (CD14++ CD16+), and non-classical (anti-inflammatory CD14+ CD16++) [4]. In addition, monocytes and neutrophils undergo polarization after IS: activated M1 monocytes or monocyte-derived macrophages and N1 neutrophils that are related to the inflammatory response can polarize to M2 or N2 phenotypes that are associated with the resolution and regenerative phase [5,6,7,8].
In models of experimental IS, neutrophils increase in the brain after 3 h and reach peak levels in the first 24 h [9, 10]. Monocytes slowly infiltrate the injury, peaking at day 3 or later after experimental stroke [2, 11]. Neutrophil levels in the brain return to near normal by a week after stroke, while monocyte increases persist for over a month [2, 9]. The prolonged monocyte/macrophage presence is likely indicative of ongoing peripheral leukocyte interaction with the injured brain associated with recovery phases. Analyzing peripheral leukocytes after stroke represents a feasible proxy to study the cellular and pathological changes that occur in response to the brain parenchyma injury. An increase of circulating neutrophils occurs promptly after stroke, and altered ratios of peripheral leukocytes (including neutrophils and monocytes) are indicators of outcome [12,13,14,15,16,17]. Successful intervention strategies in the acute and subacute phases of stroke may be improved when the pathological role of specific leukocyte types at different times is considered.
Transcriptional changes are detected promptly after IS in peripheral blood cells, showing how dynamic changes in gene expression can be revealed even in the acute phase of stroke. This results in distinct signatures depending on the cell type and stroke etiology [18,19,20,21,22]. Peripheral monocytes and neutrophils have been shown to be major cell types that display a transcriptomic response within the first 24 h after stroke [23]. In this study, the transcriptomic profiles from peripheral monocytes and neutrophils and whole blood were analyzed as a function of time and of different etiologies. Different analytical approaches (differential expression, self-organizing maps, and Weighted Gene Co-expression Network Analysis (WGCNA [24])) enabled the identification of genes that change expression following acute stroke. Though changes of gene expression in whole blood have been described within 0 to 24 h following ischemic stroke using microarrays [23, 25], this is the first study to analyze the transcriptional profiles of monocytes, neutrophils, and whole blood with RNA-seq at times ranging from 0 to >48 h. The focus on the response over time in different cell types is crucial for the eventual development of diagnostic biomarkers and cell- and time-tailored treatments.
Methods
Subjects
Thirty-eight ischemic stroke (IS) patients and 18 vascular risk factor control (VRFC) subjects were recruited at the University of California at Davis Medical Center under a study protocol reviewed and approved by the Institutional Review Board (IRB ID 248994-41). The study adheres to federal and state regulations for protection of human research subjects, The Common Rule, Belmont Report, and Institutional policies and procedures. Written informed consent was obtained from all participants or a legally authorized representative.
The criteria for recruitment are detailed in our previous study [18]. Briefly, IS diagnoses (cardioembolic (CE), large vessel (LV), and small vessel/ lacunar (SV)) were confirmed by two independent neurologists based on history, exam, brain CT or MRI, and other testing. The exclusion criteria were as follows: anticoagulation therapy (using coumadin, heparin, or any NOACs), immunosuppressive therapy, current or recent (2 weeks) infection, and hematological malignancies. Vascular risk factor control (VRFC) subjects had no history of stroke, myocardial infarction, or peripheral vascular disease, and they were recruited based on the presence of vascular risk factors including hypertension, hypercholesterolemia, and/or type 2 diabetes.
Whole blood for RNA analysis was drawn directly into PAXgene RNA stabilizing tubes for subsequent batch isolation. Blood for immune cell populations was collected in citrate tubes for immunomagnetic isolation by RoboSep (StemCell Technologies, Inc.). Cell isolation was performed as described in Carmona-Mora et al. [18]. Monocytes were positively selected using antibodies to CD14 to a purity of >93%, and neutrophils were enriched by negative selection to a purity of >99% as previously validated by flow cytometry.
RNA sequencing and differential gene expression analyses
RNA isolation and cDNA library preparation were performed as previously described [18]. In summary, total RNA was extracted from isolated monocytes and neutrophils using the Zymo Direct-zol RNA mini-prep kit (Zymo Research) according to the manufacturer’s protocol. This was followed by treatment with DNase (QIAgen). Total RNA from whole blood samples was extracted using QIAcube with PAXgene Blood miRNA Kit (QIAgen). Ribosomal RNA and globin transcripts were depleted using InDA-C (aka AnyDeplete) during library preparation by the NuGEN Ovation Universal RNA-Seq system (Tecan Genomics, Inc.). RNA sequencing yielded an average of 200 M ± 10 M 2×150 bp reads per sample. Raw data were processed to generate counts as previously described [18]. Briefly, raw reads were processed using expHTS [26] to trim low-quality sequences, for adapter contamination, and to remove PCR duplicates. Trimmed reads were aligned to the human GRCh38 primary assembly genome (GENCODE v25 annotation (http://www.gencodegenes.org/), using STAR v. 2.5.2b aligner [27]. Raw counts by-gene were generated using featureCounts of the Subread software v.1.6.0 [28], and normalized (transcripts per million, TPM) on Partek Flow software (Partek Inc.). Partek Genomics Suite was used for differential expression with an analysis of covariance (ANCOVA) model with REML variance estimate using the model: Y ijklmn = μ + Diabetes i + Diagnosis j + Hypercholesterolemia k + Hypertension l + Time (h) + Diagnosis*Time Point (TP) jm+ ε ijklmn, where Y ijklmn represents the nth observation on the ith Diabetes, jth Diagnosis, kth Hypercholesterolemia, lth Hypertension, mth Time Point (TP), μ is the common effect for the whole experiment, and ε ijklmn represents the random error component [29].
To identify differentially expressed genes, subjects were split into time points (TPs) from stroke onset (TP1= 0–24 h; TP2= 24–48 h; and TP3 = > 48 h; mean and SD for time (h) in every time point are available in Table 1). Vascular risk factor control (VRFC) subjects were assigned time point zero (TP0). Contrasts included time after stroke, interaction of diagnosis x TP, and major risk factor categories (diabetes, hypercholesterolemia, and hypertension) with cutoffs for significance set to p < 0.02 and fold-change > |1.2| to create lists of genes. Fisher's least significant difference (LSD) was used for individual contrasts [30].
Gene clustering
Gene Expression Dynamics Inspector (GEDI) v2.1 [31] was used to create mosaic grids (tiles) of self-organizing maps for visualization of differentially expressed genes over time (https://github.com/midas-wyss/gedi). Two phases of training iteration were used (40 and 100) with linear initialization. Grid sizes were chosen depending on the total number of differentially expressed genes to analyze per sample type to keep a similar number of genes per tile in all mosaics (5×7, 7×8, and 4×6, for monocytes, neutrophils, and whole blood samples respectively). Tiles corresponding to gene clusters of like-behaving genes were formed based on Pearson’s correlation. Tiles are composed of the same genes across time points, and mosaics for monocytes, neutrophils, and whole blood have different tile composition.
Self-organizing maps (SOM) [32] were implemented in Partek Genomics Suite (alpha value set to 0.1, with random initialization, exponential decay function, Gaussian neighborhood, and rectangular topology). This was done to examine trajectories of gene expression over time. The input for SOM consisted of differentially expressed genes that are present in at least two of the studied time points. In total, 500,000 training iterations were performed. Map height and width were set to 4×4 (16 profiles) in monocytes, neutrophils, and whole blood in the analyses irrespective of IS cause. For SOM in DEGs analyzed per IS etiology, map height and width were set as follows: 2×2 (CE), 2×3 (LV), and 2×2 (SV) in monocytes; 2×2 (CE), 3×3 (LV), and 2×3 (SV) in neutrophils; 2×3 (CE), 2×2 (LV), and 2×2 (SV) in whole blood. Differentially expressed gene expression was standardized by shifting to a mean of 0 (standard deviation (SD) of 1). Profiles were summarized and represented with ± 1SD. Profiles with similar dynamics were merged (based on similar directionality in every time point) for gene ontology analyses and visualization.
Cell-specific markers
The presence or enrichment (p-value < 0.05 for significant enrichment) of gene lists with blood cell type-specific genes was assessed by comparing to previously described blood cell type-specific genes [33, 34]. Enrichment analyses were performed using hypergeometric probability testing (R function phyper).
Pathway and gene ontology analyses
Pathway enrichment analyses were performed using Ingenuity Pathway Analysis (IPA, Ingenuity Systems®, QIAgen). For input, differentially expressed genes and their fold-changes from every time point and sample type, with a p < 0.05 and fold-change > |1.2|, were used. Pathways and predicted upstream regulators (Ingenuity Upstream Regulator analysis in IPA, white paper, Ingenuity Systems®, QIAgen) with Fisher’s exact test p < 0.05 were considered statistically over-represented, and those that also have a Benjamini–Hochberg False Discovery Rate (FDR) correction for multiple comparisons are indicated in the figures. IPA also computes significant pathway activation or suppression (z ≥ 2.0 and z ≤ −2.0, respectively), by using the z-score—which is based on comparisons between input data and pathway patterns, causal relationships, and curated literature. Gene ontology (GO) enrichment was explored as implemented in Partek Genomics Suite in the Gene Set Analysis, using Fisher’s exact test and FDR correction for multiple comparisons, with significance set at p<0.05.
Weighted gene co-expression network construction and analysis
Separate weighted gene co-expression networks were generated for isolated monocyte (MON network) and neutrophil (NEU network) data, as well as for whole blood (WB network). VRFC samples were excluded, and genes below a minimum of 40 counts in every sample were filtered out for MON and NEU, and below a total of 80 counts for WB. The MON network was generated using the 14,955 detected genes after filtering across 35 IS samples; the NEU network was generated using 13,921 genes across 31 IS samples; the WB network was generated using 15,360 genes across 37 IS samples. Data were imported into R and checked for missing or zero-variance counts using the function goodSamplesGenes.
Networks were generated with the Weighted Gene Co-Expression Network Analysis (WGCNA) package using a Pearson correlation to measure co-expression [24]. An approximate scale-free topology was depicted by the data. Soft thresholding powers (β) of 14, 8, and 16 were selected for the MON, NEU, and WB networks, respectively, to maximize strong correlations between genes while minimizing weak correlations [35]. A signed network was used to consider both positive and negative correlations [36]. The cutreeDynamic function (method = tree; deepsplit = 1; minimum module size = 50) was used to form modules due to its adaptability to complex dendrograms and ability to identify nested modules [36]. Hub genes were defined as the top 5% by interconnectivity and may represent genes with important regulatory or molecular signaling roles [37, 38].
Module association with clinical parameters
Module-parameter associations were determined using Kruskal-Wallis and Spearman Ranked Correlation tests in Partek Genomics Suite for categorical and continuous variables, respectively. Ranked statistical tests were utilized to minimize the impact of outliers. Parameters were associated with the module eigengene, or first principal component of expression of genes within a module. Continuous time since event and all clinical parameters (age, diabetes, hypercholesterolemia, hypertension, race, and sex) were examined on the complete datasets. A p-value < 0.05 was considered significant.
Hub gene analyses
Functional modules were detected with HumanBase [39] in a tissue-specific manner, using all hub genes per sample type as input (monocytes and whole blood). Briefly, this tool is based on shared k-nearest-neighbors and the Louvain community-finding algorithm. Gene ontology terms in the results are considered significant based on a corrected P-value < 0.05 (one-sided Fisher’s exact tests with a Benjamini–Hochberg FDR correction for multiple comparison).
Results
Cohort demographics
Individuals were binned into time points (TPs) from stroke onset (TP1= 0–24 h, average ~15 h; TP2= 24–48 h, average ~32 h; and TP3= >48 h, average ~56 h), and vascular risk factor controls (VRFC) were assigned TP0 (Table 1). The cohort demographics and clinical characteristics per time point are presented in Table 1, as well as the total number of subjects per subgroup. The statistical analysis of variables showed that age, diabetes, and hypertension were not significantly different (p < 0.05) between TPs and controls in any of the sample types (monocytes, neutrophils, or whole blood). NIHSS at admission was significantly higher in subjects for samples at TP3 when comparing to VRFC in monocytes. Hyperlipidemia was significantly different for TP1 in all sample types when comparing to VRFC, while sex is significantly different in whole blood at TP3. In this cohort, 5 out of 38 IS patients received thrombolytic therapy (recombinant tissue plasminogen activator, rtPA) within 4.5h of their stroke. None of the IS cases developed hemorrhagic transformation.
Gene expression dynamics in the acute and subacute post-stroke phase
The dynamic changes in differential gene expression across different time windows were assessed in monocytes (M, Additional file 1: Tables S1A-C), neutrophils (N, Additional file 1: Tables S1 D-F), and whole blood (WB, Additional file 1: Tables S1G-I) using a significance cutoff P<0.02 (Fig. 1). In monocytes, there were 290, 645, and 352 DEGs at 0–24 h (TP1), 24–48 h (TP2), and >48 h (TP3) compared to controls, respectively (Fig. 1). In neutrophils, 508 (TP1), 745 (TP2), and 547 (TP3) genes were differentially expressed compared to controls. A total of 610 (TP1), 260 (TP2), and 227 (TP3) DEGs were identified in whole blood compared to controls (Fig. 1 and Additional file 1: Table S1). Across the time points, the number of up- or downregulated DEGs in each sample type showed distinct and dynamic trends. For example, in monocytes, most DEGs are downregulated, with the most at TP2 (Fig. 1A, right panel). Conversely, in neutrophils most genes are upregulated, peaking at TP2 (Fig. 1B, right panel). In whole blood, the highest upregulation occurs within the first 24 h (Fig. 1C, right panel). Biotypes of DEGs were relatively similar in monocytes and neutrophils, while in whole blood more lncRNAs and T-cell receptor genes were detected at >48 h post-stroke (Additional file 2: Fig. S1).
Pathway enrichment analyses also show distinctive shifts in molecular function and cellular roles of the DEGs over time per sample type. Most of the over-represented pathways in monocytes are shared across at least two TPs (Fig. 2A) (average of 88.8% shared pathways) and most of these pathways are predicted to be significantly suppressed (Z ≤ −2.0), including chemokine, IL-8 signaling, and production of NO and ROS in monocytes/macrophages (Additional file 1: Tables S2A-C and Additional file 2: Fig. S2B and 2C). The only enriched pathways with predicted activation are PTEN signaling (at 24–48 h and >48 h) and RhoGDI signaling (at 24–48 h) (Additional file 1: Tables S2B-C). Calcium signaling, CCR3, chemokine, CREB, and CXCR4 signaling were the top enriched pathways in monocytes at 0–24h (Fig. 2C). Actin, B cell receptor, ephrin, Erb, and ERK/MAK signaling were the top pathways over-represented in monocytes at 24–48 h (Fig. 2C). Similar pathways were regulated in monocytes at times >48 h in addition to estrogen, FLT3, GM-CSF, GNRH, and HIF signaling (Fig. 2C).
From the pathways exclusively regulated at a single TP, suppressed oxidative phosphorylation, suppressed Wnt/Ca2+, heparan sulfate, and nNOS pathways top the list in the first 24 h (Fig. 2D, left panel). Ephrin B, Tec kinase, reelin, TGF-β, and IL-6 signaling are amongst the suppressed pathways at TP2 (Fig. 2D, middle panel), while tryptophan degradation and IL-15 signaling are suppressed at times over 48h (Fig. 2D, right panel).
In neutrophils, fewer enriched pathways were shared across at least two TPs (Fig. 3A, Additional file 1: Tables S2D-F) (average of 62.5% of shared pathways in TPs pairwise comparison). There were increased numbers of pathways at TP2 and TP3 compared to TP1 (Fig. 3A, C). EIF2 signaling and leukocyte extravasation were over-represented at all TPs (Fig. 3B). CD40 and HIF pathways are over-represented at TP2 and TP3. Several interleukin pathways are enriched in the later TPs (Fig. 3C and Additional file 1: Tables S2D-F). IL22 was expressed at TP2 and TP3 (Fig. 3C). IL-1, 4, 7, and 8 were exclusively over-represented in TP2, while IL-2, 10, and 15 pathways were only enriched after 48 h (Fig. 3D).
In contrast to the isolated cell samples, whole blood showed no shared pathways across all TPs (Fig. 4A). Most pathways were enriched at only one TP in whole blood (Fig. 4C and Additional file 1: Tables S2 G-I), and many of these were specific for 0–24 h including p53, AMPK and ATM signaling, and FXR/RXR activation (Fig. 4C).
Upstream regulators were computed to identify drivers of gene expression at specific TPs (Additional file 1: Table S3). Similar to what was seen with canonical pathway enrichment, there were more regulators shared across TPs for monocytes (Additional file 2: Fig. S2, Additional file 1: Tables S3A-C) and neutrophils (Additional file 2: Fig. S3; Additional file 1: Tables S3 D-F) compared to whole blood (Additional file 2: Fig. S4; Additional file 1: Tables S3 G-I). Upstream regulators (both activators and suppressors) for monocytes included many molecules common to all time points: 15-LOX, ACKR3, BTNL2, CSF2, ERG, ERK, filgrastim, IKZF1, lipopolysaccharide, Msx3, OGA, PCYT2, PGR, PLA2R1, Progesterone, Rhox5, TNF, U18666A, and WAC (Additional file 2: Fig. S2; Additional file 1: Tables S3A-C). Upstream regulators for neutrophils predicted at all time points included the following: calcitriol, CIP2A, ciprofloxacin, CSF3, dexamethasone, ESR1, filgrastim, FOXO3, FOXO4, interferon beta-1a, LARP1, MRTFB, MYCN, NUPR1, OSM, RPS15, RRP1B, torin1, TP53, and YAP1 (Additional file 2: Fig. S3; Additional file 1: Tables S3 D-F).
In contrast, the vast majority of upstream regulators are unique for every TP in whole blood (Additional file 2: Fig. S4; Additional file 1: Tables S3 G-I). In addition, there are more regulators in the first 24 h in whole blood (Additional file 2: Fig. S4; Additional file 1: Tables S3 G-I), consistent with the proportion of whole blood DEGs seen at this time point.
Identification of time-dependent DEG clustering profiles
Tile mosaics based on self-organizing maps were constructed for DEGs with GEDI to allow an overall view of gene clusters across time by correlated expression. These mosaics reveal the distinctive dynamics of gene expression changes across every time point in monocytes (Additional file 1: Table S4A), neutrophils (Additional file 1: Table S4B), and whole blood (Additional file 1: Table S4C) (Fig. 5). Every tile clusters a specific group of DEGs, and the tile mosaics for every time point show how that group of genes changes over time (Fig. 5 and Additional file 1: Table S4). Tracked at each coordinate, most tiles in the mosaic maps show evidence of dynamic changes at each time point (Fig. 5, 0 to >48 h).
To understand the trajectory of the DEG clusters over time following stroke and to identify key genes per time window, we analyzed expression profiles using self-organizing maps (SOM) (Additional file 1: Table S5), where every cluster was plotted separately to obtain a profile. For each cell type including monocytes (Additional file 1: Table S5A), neutrophils (Additional file 1: Table S5B), and whole blood (Additional file 1: Table S5C), there were different patterns of gene expression that included increases over time, decreases over time, peaks at 24–48 h, and valleys at 24–48 h (Fig. 6). Though comparable expression patterns between monocytes, neutrophils, and whole blood are seen over time, these were associated with unique functions (GO) in each cell type (Fig. 6, Additional file 1: Tables S6A-L). Examples include neutrophil degranulation which increases over time in whole blood, decreases of IL-1 secretion over time in monocytes, leukocyte migration peaking in neutrophils at 24 h, and a valley for regulation of platelet activation at 24–48 h in monocytes (Fig. 6).
Time-dependent gene expression changes per stroke etiology
To identify expression changes associated with IS causes, the same model and criteria for DEGs were used, but time points were considered as 0–24 h and over 24 h (>24 h) post-stroke in order to assure enough IS cases per etiology and time point. The cohort characteristics and analyses for the clinical variables after this re-stratification can be found in Additional file 1: Table S7. The cohort was split according to the main IS etiologies: cardioembolic (CE, Additional file 1: Tables S8 A, B, G, H, M, N), large vessel (LV, Additional file 1: Tables S8 C, D, I, J, O, P), and small vessel (SV, Additional file 1: Tables S8 E, F, K, L, Q, R).
More DEGs were identified in all IS causes vs. VRFC in the >24 h period than in the first 24 h for monocytes and neutrophils (Fig. 7, Additional file 1: Table S8). In contrast, in whole blood the first day post-IS showed higher DEGs in CE and LV stroke, pointing to other critical cell types in the response within this time window (Fig. 7 and Additional file 1: Table S8). Gene expression tended to be suppressed in monocytes following CE, LV, and SV strokes at all time points, whereas gene expression was generally upregulated in neutrophils in all stroke etiologies at all time points (Fig. 7). Gene expression was increased in whole blood at 0–24 h in CV and LV strokes (Fig. 7). tPA treatment on 4–5 patients (depending on sample type) did not affect the gene expression changes identified. On PCA (principal component analysis), tPA subjects did not cluster specifically when assessing the DEGs in both 0–24h and >24 h (Additional file 2: Fig. S5).
Key DEG clusters per IS cause were identified using SOM for gene expression trajectories (Additional file 1: Table S9). Comparable expression trajectories over time are seen for CE, LV, and SV strokes in monocytes (Additional file 1: Tables S9 A, B, C), neutrophils (Additional file 1: Tables S9 D, E, F), and whole blood (Additional file 1: Tables S9 G, H, I) (Fig. 8), although not all trajectory types are consistently present in every IS cause and in every sample type (Fig. 8).
GO enrichment revealed the associated functions for similar profile trends in every sample type, which tended to be unique for each IS etiology and cell type (Additional file 1: Tables S10 A-W; Fig. 8). In neutrophils, for example, the GO terms neutrophil migration and positive regulation of defense response are enriched in the profile of genes that decrease expression throughout TPs in LV and SV strokes, while it is not present in CE (Fig. 8). Interestingly, regulation of neutrophil migration is enriched in the profile of DEGs that consistently increase expression in LV stroke in the whole blood samples (Additional file 1: Table S10).
Time-associated gene expression networks/modules in monocytes, neutrophils, and whole blood
To analyze the time-dependent changes that occur after IS from a genome-wide perspective, we used Weighted Gene Co-expression Network Analysis (WGCNA). The gene expression counts for monocytes, neutrophils, and whole blood samples from IS patients were used to generate three separate networks (Additional file 2: Fig. S6), and significant modules of co-expressed genes were identified in relation to time (as continuous variable in hours post-IS) (Fig. 9A, Additional file 1: Table S11).
For monocytes, 36 modules of co-expressed genes and a module of unassigned genes (MON-Grey) were identified within this cell type network. Sixteen monocyte modules were found to be significant for time: MON-DarkRed, MON-MidnightBlue, MON-GreenYellow, MON-Sienna3, MON-Red, MON-Pink, MON-LightCyan, MON-Violet, MON-Yellow, MON-DarkMagenta, MON-SteelBlue, MON-Blue, MON-Grey60, MON-Tan, MON-PaleTurquoise, and MON-Magenta (Fig. 9A). The gene lists for each of these 16 modules are provided in Additional file 1: Table S11A. Of these, MON-GreenYellow and MON-Blue were also significantly related to diabetes. In neutrophils, 48 modules of co-expressed genes and a module of unassigned genes (NEU-Grey) were identified within the network (Fig. 9A). One neutrophil module (NEU-DarkSlateBlue) was significant for its relationship to time. The gene list associated with this one neutrophil module is provided in Additional file 1: Table S11B.
Within the whole blood network, 51 modules of co-expressed genes and a module (WB-Grey) of unassigned genes were identified. Ten of these whole blood modules were significantly related to time: WB-Brown, WB-Grey60, WB-Tan, WB-Turquoise, WB-MidnightBlue, WB-SteelBlue, WB-DarkRed, WB-SaddleBrown, WB-Pink, and WB-MediumPurple3. Of these, WB-Tan and WB-DarkRed were also related to sex, while WB-SteelBlue and WB-MediumPurple3 were also related to age. The gene lists for each of these 10 whole blood modules are provided in Additional file 1: Table S11C.
Genes that are highly interconnected in each time-associated module were also identified for all sample types. These represent “hub” genes with high potential for functional or regulatory significance (Additional file 1: Table S12). Examples of hub genes from the different modules in monocytes include the following: SCAF11, CREBZF, SETX, JAK1, EIF3F, HNRNPK, UBC, DICER1, CAPN2, RAB10, SF3B1, DAZAP2, UTY, SPARC, PPP4C, and RAB11FIP1. Examples of hub genes from the different modules in whole blood include the following: TNFRSF18, PCSK9, IGSF9, MTHFSD, DCAF12, BCL2L1, MAU2, TLL2, FCRL6, PARVG, and VASP.
A tissue-specific network-based functional characterization of the hub genes was performed [39] to gain perspective on the hub genes identified in monocytes (Fig. 9B, Additional file 1: Table S13A) and whole blood (Fig. 9C, Additional file 1: Table S13C). Neutrophils were not analyzed since only two hub genes were found: RP11-35J10.7 and AP000593.6 (Additional file 1: Table S12B).
In monocytes, eight functional HumanBase modules were identified and included the following terms: response to virus and immune effector process (HG-M1, Hub Genes in Monocytes cluster 1), RNA splicing (HG-M2), RNA metabolism (HG-M3), DNA repair (HG-M4), cell morphogenesis (HG-M5), regulation of cell cycle (HG-M3 and HG-M6), cell adhesion and secretion (HG-M7), and protein transport (HG-M8) (Fig. 9B and Additional file 1: Table S13). The terms for HG-M2-4 are unique to monocytes (lists in Additional file 1: Table S13B). The network-based functional interpretation of the whole blood hub genes in time-related WGCNA modules identified five functional HumanBase modules that included the following terms: response to virus and immune effector process (HG-WB1, Hub Genes in Whole Blood cluster 1), leukocyte activation (HG-WB2), filopodium assembly (HG-WB3), regulation of cellular response to transforming growth factor beta receptor (HG-WB4), gene silencing (HG-WB4), and cytoskeleton organization (HG-WB5) (Fig. 9C and Additional file 1: Table S13D). The first two terms in HG-WB1 were shared between whole blood and monocytes HG-M1. The rest were only identified in whole blood (overlaps not shown, all lists available in Additional file 1: Table S13).
There was a significant correlation between the expression of some monocyte and whole blood hub genes with NIHSS on admission (Additional file 1: Tables S14A-B). Examples of hub genes in monocytes that correlate with NIHSS included MIER1, DICER1, FBXW5, SELPLG, TIA1, and BCKDK (p < 0.01, Additional file 1: Table S14A). There is no overlap between the NIHSS significantly correlated hub genes from monocytes and whole blood (Additional file 1: Table S14). The hub genes identified in neutrophils included RP11-35J10.7 which is a novel transcript that is sense intronic to OVCH2 (ovochymase-2) (Additional file 1: Table S12B). The second was AP000593.6, a U2 small nuclear RNA auxiliary factor 1 (U2AF1), currently annotated as a pseudogene (Additional file 1: Table S12B). Since there were only two neutrophil hub genes, no additional functional analyses were performed for neutrophils.
Gene markers delineating M1/M2 monocyte and N1/N2 neutrophil subsets were identified in the time-related modules. Classical or inflammatory monocyte markers like CD14, CCR2, CSF1R, and non-classical or intermediate marker FCGR3A [40, 41] were found in three time-associated modules in monocytes (Additional file 1: Table S11A) along with M2/N2 polarization-related markers CCL2 and STAT3 [8, 42,43,44]. The presence of CCR2 in a time-related module, the main chemoattractant for monocytes to the injury site in IS [45, 46], may also suggest an active recruitment of this leukocyte type. In neutrophils, none of these markers were present in the modules, while in whole blood CCR2 and TNF (M2 [44, 47] and M1/N1 markers [43, 48,49,50], respectively) were present in two significant time modules (Additional file 1: Table S11C).
Discussion
Ischemic stroke elicits specific responses in peripheral blood, which can be seen at the transcriptional level in monocytes and neutrophils [18] as well as other cell types. This is the first study to analyze those differences based on time trajectories in the early time window after stroke using RNA-seq. A summary of the flow of approaches utilized is presented in Additional file 2: Fig. S7. Hundreds of genes change expression at different time points following stroke, with most genes in neutrophils being upregulated and most genes in monocytes being downregulated over time. The genes that change expression at each time point and in each cell type are associated with specific pathways that are usually unique for each cell type and time. Self-organizing maps identified specific trajectories of gene expression for monocytes, neutrophils and whole blood that were similar but were associated with differing signaling pathways. These pathways also differed as a function of cardioembolic, large vessel, and small vessel causes of stroke. We also show modules of genes that are co-expressed and change with time following stroke using WGCNA in neutrophils, monocytes, and whole blood. Specific hub genes (potential master regulators) were identified for modules that were associated with time after stroke. These analyses help understand genes and functions changing across the early post-ischemic stroke period and groups of genes that behave in a concerted fashion. Both aspects are important in the search for better understanding of the complex molecular and immune interactions after ischemic stroke to identify optimal treatment targets and their optimal time windows.
Gene expression dynamics after ischemic stroke
Differential gene expression changes in ischemic stroke (IS) versus controls were assessed by splitting subject samples into different time windows. Monocytes displayed downregulation of most of their DEGs, while neutrophils were generally upregulated. This is in accordance with the up- and downregulation patterns seen in our previous study [18] that did not consider time. For both monocytes and neutrophils, the largest number of differentially expressed genes was observed at the 24–48 h time point, whereas the greatest number of regulated genes in whole blood was observed at the 0–24 h time point. This suggests that additional cell types in whole blood may drive the early peripheral immune response in the first 24 h.
Canonical pathways associated with monocytes over time were mostly suppressed, similar to our previous findings [18]. Many of the over-represented functional pathways are shared in monocytes across time points, with very few shared pathways over time in neutrophils, and none in whole blood. These findings suggest more dynamic changes of genes and pathways in neutrophils and leukocytes, with much less dynamic change in monocytes at least over the first few days after a stroke. When analyzing enrichment of canonical pathways and prediction of upstream regulators, a z-score is calculated with Ingenuity software, where z ≥ 2 is activated (denoted ahead with subscript act) and z ≤ −2 is inhibited (denoted ahead with subscript inh).
Cytokine expression changes, and specifically those of interleukins, were prominent and some changed with time in monocytes. Since changes of many interleukins have been directly quantified in blood of IS patients over time [51, 52], it was not surprising that we saw many changes of cytokine and interleukin genes representing their receptors. The IL-4 signaling pathway that is essential for M2 polarization [53] was enriched across all time points with IL-4 itself predicted as inhibited based on the observed overall gene expression profile in monocytes at 24–48 h and >48 h in this study. The IL-8 signaling pathway inh was also enriched across all time points in monocytes. Pro-inflammatory cytokine IL-8 has a role in monocytic recruitment, by promoting monocyte adhesion to the vascular endothelium [54]. TGF-ꞵ (produced by M2 macrophages) has a protective role in stroke [55] and is inhibited in monocytes only at 24–48 h. IL-6 signaling was enriched and significantly suppressed only at 24–48 h, and IL-15 signaling is suppressed at all time points. IL-6 (pro- and anti- inflammatory cytokine) is a marker of infarct size and functional outcome and is part of M2 polarization signaling promoting monocytic differentiation into macrophages [56,57,58,59,60,61,62]. IL-15 is mostly produced by monocytes and its blockade reduces brain injury after ischemia [63, 64]. In monocytes, IL-5inh, IL-2inh, and TNFinh were predicted upstream regulators of gene expression at 0–24 h; IL-1Binh, IL-2inh, IL-3, IL-4inh, IL-5inh, IL-10inh, IL-13, IL-15, IL-33inh, and TNFinh are predicted upstream regulators at 24–48 h; and IL-2inh, IL-3, IL-4inh, IL-13, IL-33, and TNFinh are upstream regulators at >48 h.
The release of granules from neutrophils results from IL-8 stimulation [65, 66], a pathway that was over-represented only at 24–48 h and predicted as inhibited. Other pathways uniquely enriched between 24 and 48 h post-IS included those for N2 markers TGF-β and IL-4. IL-4 activates neutrophils [67], and along with pro-inflammatory IL-1 and IL-7, are key for T-cell proliferation and homeostasis [68].
The IL-15 pathway inh is enriched at 24–48 h and >48 h, and IL-15 is identified as an upstream regulator itself at 24–48 h. IL-10 pathway and IL-2 pathway inh were only present in neutrophils at times >48 h. IL-10 is an anti-inflammatory cytokine and N2 marker which decreases inflammation and apoptosis [62]. IL-15 induces changes in neutrophils and increases their phagocytic activity [69]. IL-2inh, a T-cell growth factor and activator, is predicted as inhibited in this study which corresponds with the temporal profile for this cytokine in IS [51, 70]. Additionally, the N1 marker GM-CSF signaling pathway, which is protective in experimental stroke [71], was enriched at > 48 h after stroke. IL-2, IL-3, IL-4, IL-5, IL-6, IL-10, IL-11RA, IL-15, IL-16, IL-23A, IL-32, and TNFinh were predicted upstream regulators of gene expression in neutrophils at 24–48 h, and IL-2, IL-3, IL-4, IL-5inh, IL-10 and TNFinh are predicted upstream regulators of gene expression profiles seen in neutrophils at >48 h.
Whole blood over-represented pathways were most robust at 0–24 h after stroke, and rather than being related to cytokines or chemokines, they appear to relate to specific biological responses to ischemia including p53, ATM, cAMPact, and AMP-activated protein kinase signaling [72,73,74,75,76,77]. Ephrin receptor signaling act that is key in angiogenesis after stroke is over-represented in whole blood at 24–48 h after stroke [78]. In general, fewer canonical pathways were over-represented in whole blood, despite having as many or more DEGs as the monocyte and neutrophil datasets. This could be due to a wide range of biological processes that might be regulated in opposite directions in a heterogeneous sample of multiple cell types from peripheral blood. This highlights the need to further delineate temporal changes of DEGs in each immune cell type following stroke. Though no cytokines were predicted as upstream regulators of gene expression in whole blood, there were a significant number of microRNAs at the three times that were identified as predicted upstream regulators.
The molecules identified as upstream regulators of gene expression in the different cell types might be particularly important therapeutic targets since they have potential to affect so many downstream genes with concerted expression changes after stroke [79]. For example, in neutrophils calcitriol is identified as an upstream regulator at all three time points, but with a trend towards suppression between 0 to 48 h, followed by activation at >48 h after IS. Though calcitriol has shown neuroprotective effects in experimental models of IS [80,81,82], our human data might point to loss of efficacy outside of an early time window. Pterostilbene, an antioxidant, anti-inflammatory, and anti-apoptotic compound with a beneficial effect after stroke [83,84,85,86], was predicted as an upstream regulator of the observed changes in neutrophil gene expression but only at times >48h after IS.
Time-associated gene expression modules and networks in monocytes, neutrophils, and whole blood
By analyzing correlation of co-expression patterns as a continuous variable over time, different modules and their associated hub genes were identified. Hub genes potentially drive the gene expression changes for the different modules of co-expressed genes. The most time-associated modules were associated with monocytes followed by whole blood and then neutrophils.
Monocyte polarization markers are present in four time-associated modules, including the MON-LightCyan module where CD14 (Classical or inflammatory monocyte marker) and STAT3 (M2 marker) genes are present. In this module, hub genes are also significantly enriched for monocyte-specific genes, suggesting that the expression dynamics of these genes could be critical for the monocyte response after stroke. Furthermore, the vast majority of the MON-LightCyan module hub genes display higher expression levels in classical monocytes (The Human Protein Atlas [41]). Other key markers present in time-related modules are CCR2 (Classical monocyte marker), CSF1R (M2 polarization), and CCL2 (M2 marker). Together, these results may indicate transformation of monocytes to the restorative M2 type over time, even within the ~72 h time period of this study. Nonetheless, the correlation of CCR2 with time is consistent with active recruitment of monocytes from the bone marrow, which is in line with experimental data where monocytes begin to accumulate around day 3 or later after ischemic stroke [87].
In neutrophils, the genes in the only time-significant module did not overlap with cell type or polarization markers. In contrast, the gene markers TNF-α (M1 and N1) and CCR2 (M2) were present in whole blood modules associated with time. As expected, different polarization types are present in whole blood, the same is true for cell type markers, which were significantly enriched in several significant to time modules (monocyte, granulocyte, erythroblast, natural killer, and T cell markers). CCR2 expression in neutrophils is not expected in a resting state, but has been linked to altered neutrophil programming in inflammatory states [88] and promoting chemotactic attraction to the injury site [88].
The functional and biological roles of the highly interconnected time-associated hub genes identified through WGCNA point to unique processes and drivers of gene expression in monocytes and whole blood. Hub genes in monocytes are enriched significantly for RNA splicing and RNA metabolism functions. This may reflect active formation of specifically spliced gene transcripts in the response to injury and timing of polarization which likely changes as monocytes move from inflammatory to restorative subtypes (M1 and M2 monocytes, respectively).
Hub gene functional modules in whole blood represent a composite of the gene expression changes in individual immune cell types. Enrichment of a wide host of functions including leukocyte activation, filopodium assembly, cytoskeleton organization, regulation of cellular response to transforming growth factor beta receptor, and gene silencing point to the breadth of cell types in blood undergoing cell proliferation, activation, and migration in the first 72 h after ischemic stroke. Through different approaches, Li et al. [89] identified several stroke-related hub genes, where one of them is common to our analysis in whole blood hub genes (TSPAN14). This gene is involved in the regulation of Notch signaling pathway [90] and may be a promising target in CE stroke.
Several time-associated hub genes were positively correlated with NIHSS (NIH stroke scale, a measure of stroke severity) at admission. In monocytes, a positive significant correlation with NIHSS was found for TIA1 (T-cell-restricted intracellular antigen-1 aka cytotoxic granule associated RNA binding protein), an anti-inflammatory protein in peripheral tissues and a repressor of TNF-α expression. It is a M1 marker and a key regulator of the innate immune response of the CNS during stress [91,92,93]. MIER1, a transcriptional repressor [94], and RMI1, a DNA repair protein [95], were also positively correlated with stroke severity. Variants in these genes are associated with monocyte counts [96] and myeloid white cell counts [96], respectively; and expression positively correlates with the infiltration of monocytes, macrophages, and other immune cells in gastric carcinoma [97]. Other NIHSS-time-associated hub genes (positive correlation) included the following: SELPLG, high affinity receptor for cell adhesion molecules in leukocytes [98]; GNAI3, associated with the response to intracerebral hemorrhage [99] and involved in VEGF-induced angiogenesis [100]; and FBW5, an E3 ubiquitin ligase and negative regulator of MAP 3K7/TAK1 signaling in the interleukin-1B (IL1B) signaling pathway [101]. IL1B is a key player in the pathogenesis of brain damage after ischemia [102,103,104] .
In neutrophils, time-associated hub genes did not show a significant correlation with NIHSS at admission. Only one gene in the time-significant module, MCTS1, negatively correlated with stroke severity. MCTS1 is a translation enhancer [105] and is a target of let7i, which is involved in leukocyte attachment and recruitment to the endothelium in the brain [106, 107].
In whole blood, immunoglobulin constant region and variable region genes are hub genes in time modules. IGLV1-40, IGLV3-27, IGKV1-12, IGHV3-30, and IGLC3 showed positive correlation with NIHSS at admission. This highlights changes in the humoral immune response across early times after stroke and could relate to stroke outcomes [108]. However, given that most of the genes code for variable region chains, this response may also be variable across patients. Further studies are needed to examine these immune humoral profiles and their relationship to evolving stroke injury and repair.
Refining key genes and responses after IS by analyzing gene clusters
GEDI
To generate DEG clusters based on self-organizing maps, tiled mosaics were constructed for each cell type over the three time windows. The GEDI [31] maps allow visualizing coordinated gene expression changes across time for smaller groups of DEGs (1–77 genes per tile/cluster, averaging 11.9 (M), 9.9 (N), 8.4 (WB)), and to visually compare distinct organization of the mosaics between monocytes, neutrophils, and whole blood. In the mosaics, tiles of interest are those that have marked differences in expression at a specific time point, those that change consistently through time, and those that maintain constant expression levels across time. From the DEGs grouped in different tiles, monocyte and neutrophil-specific markers were identified. Most tiles containing these markers were neighboring each other, likely because they were correlated by expression through the dimension reduction used to construct the GEDI maps. Nonetheless, in whole blood, the cell-specific marker genes present in some tiles (monocytes, granulocytes, erythroblast, B cells, and megakaryocytes) do not group in the same neighborhoods. Analyzing distant cell marker tiles can be critical to refine different functional relevance for the DEGs in those tiles.
In the monocyte GEDI maps, the lower left tile is one case of a cell marker-containing tile that is unrelated to other tiles in the mosaic with other cell-specific markers. This tile/group of DEGs has sustained low expression (namely 1,7) and contains the monocyte marker CACNA2D4 (Calcium voltage-gated channel auxiliary subunit alpha2delta 4). This gene displays higher expression in classical monocytes (http://www.proteinatlas.org) [41, 109] and belongs to the TCR signaling pathway. CACNA2D4 gene expression dynamics, as visible in its GEDI tile, could indicate a switch to non-classical monocytes. Closer examination of the same tile shows other genes including SURF1 (Cytochrome C Oxidase Assembly Factor) and TSPAN14, which are associated with monocyte counts in genetic studies [96]. Other genes in this cluster include DNPH1 (2'-deoxynucleoside 5'-phosphate N-hydrolase 1), ZBTB5 (Zinc finger and BTB domain containing 5), MTBP (MDM2 binding protein), and six other uncharacterized transcripts, which may share similar roles to the other genes in the tile. These expression correlations still do not fully define shifts towards a cell subtype, like classical, non-classical, and intermediate monocytes. For example, TSPAN14 is highly expressed in non-classical monocytes more than in other subtypes, while DNPH1 is predominant in intermediate monocytes. Altogether, looking at specific tiles can implicate key genes and cellular processes that change with time after stroke.
In monocytes, “neighborhoods” of tiles with monocyte markers can be seen on the upper and lower parts of the mosaic. These tiles showed different trends across time, from decreasing expression (tile 1,1) to unchanged gene expression (i.e., 1,7; 5,1; and 5,2). The grouping of cell-specific genes is also seen and even more accentuated in neutrophils, where granulocyte markers cluster in four opposite corners of the mosaic, displaying increasing or consistent high expression. Also in these neutrophil mosaics, two upper tiles in opposite corners (namely 2,1—left and 7,1—right) rapidly change expression across TPs 1 and 2, and then decrease to reflect levels more like TP0/controls. These tiles also have interleukin receptor genes (Fig. 5, depicted by +). Tile 2,8, which contains an immunoglobulin gene (depicted with *), shows a pattern of increasing expression across time. DEGs of interest present in tile 2,8 of the neutrophil mosaic include TLR5, an activator of the innate immune response [110], and platelet aggregation gene PDK1, key for cell division in hypoxic conditions [111]. This tile also includes relevant genes that have robust expression in healthy granulocytes or neutrophils, but in stroke are more suppressed in TP3, like KREMEN1, a negative regulator of Wnt/β catenin pathway [112] and PPP4R2, a modulator of neuronal differentiation and survival [113].
Additionally, we consistently see dynamic profiles of immunoglobulin expression in our samples including the whole blood mosaic (left lowermost tiles). These clusters of DEGs peak in the first 24 h and decrease thereafter. Several immunoglobulin genes are expressed as a function of time across samples. B cell activation and increased immunoglobulin production have been demonstrated after stroke [114,115,116]. Further work could elucidate whether the immunoglobulin genes detected in this study are related to the expected immune response after stroke, or to the production of autoantibodies seen in ischemic stroke patients [117,118,119,120]. The latter will be interesting to explore, given the cellular response to brain antigens is associated with infarct size and outcome [121, 122].
Similar expression dynamics by SOM
Self-organizing map (SOM) profiles enable grouping of the DEGs into clusters based on trajectory/directionality of expression changes and their functional associations across the time following stroke. After analyzing GEDI tiles, this is also useful because ontology analyses are more precise in larger gene groups and SOM profiles draw from larger and/or more balanced clusters. These profiles are crucial for understanding which molecular pathways and cell types show progressive activation or suppression over time, or whether they have “peaks” or “valleys” of expression over time.
In monocytes and neutrophils, the profiles of DEGs that peak only in the first 24 h after stroke are enriched for myeloid leukocyte activation and leukocyte migration, respectively. Interleukin-1 beta secretion in monocytes, and positive regulation of reactive oxygen species metabolism in neutrophils, are over-represented in profiles that decrease expression over time. Genes that are suppressed between 0 and 24 h and then recover to levels seen in VRFCs are enriched for terms related to JNK and MAPK cascades, platelet activation, and macrophage chemotaxis in monocytes. This pattern also has over-representation of chemokine production and leukocyte aggregation genes in neutrophils. Opposite trends can also be seen between monocytes and neutrophils: Cdc42-associated signaling (critical in cell growth and differentiation [123]) is enriched in DEGs that decrease expression over time in monocytes, while in neutrophils Cdc42-associated signaling is associated with genes that increase expression over time.
Furthermore, the trajectory of the SOM profiles detected in the DEGs at different time points could be used to prioritize diagnostic biomarker candidates. A diagnostic panel that could be used to diagnose stroke in the first 3 days should primarily include DEGs that either consistently increase or consistently decrease over the 3 days after IS. Such panels could indicate not only that an ischemic stroke had occurred, but could indicate roughly how much time had elapsed following the stroke—something that cannot be estimated with current methodology.
Refining genes for cardioembolic, large vessel, and small vessel strokes
There were large numbers of DEGs expressed in monocytes and neutrophils for each cause of stroke (CE, LV, SV) with somewhat more DEGs 24 h after stroke. In contrast, there were many more DEGs expressed in whole blood at 0–24 h in CE, LV, and SV stroke. This suggests that other cells (e.g., B and T cells) in whole blood (in addition to neutrophils and monocytes) were contributing to the whole blood responses to stroke at 0–24 h.
The trajectories of SOM profiles in strokes of all causes combined were similar to profiles for each stroke cause including CE, LV, and SV causes of strokes. Though the temporal profiles were similar, the DEGs and enriched GO terms were generally different for CE, LV, and SV causes of stroke in monocytes, neutrophils, and whole blood. There were exceptions, such as the “Regulation of Golgi to plasma membrane protein transport” pathway, which progressively decreased in expression in monocytes in CE, LV, and SV stroke, a pathway shown to affect outcomes in experimental stroke [124]. Another example was “NADH oxidation” which peaked at 24 h in neutrophils in CE, LV, and SV strokes and is known to play a role in experimental stroke [125]. There were few shared GO terms enriched in DEGs in whole blood CE, LV, and SV strokes, likely related to the cellular heterogeneity of whole blood.
SOM profiles identified important attributes of DEGs from opposite expression trajectories and between IS etiologies. In neutrophils, the toll-like receptor signaling pathway is associated with a profile that decreases at less than 24 h and increases after 24 h in CE stroke. In contrast, for LV, the toll-like receptor pathway is enriched in a profile that decreases at all times compared to controls. Toll-like receptor signaling modulates critical immunomodulatory NFkB signaling and is a promising target for treating cardiovascular disease [126,127,128] (clinicaltrials.gov ID# NCT04734548). TLR signaling impacts downstream pro- or anti-inflammatory molecules, including TNF-α, interleukins, interferons, and TGF-β [129].
There are also differences between SOM profiles in whole blood for different causes of stroke. “Positive regulation of IL-1β secretion” genes peak in the first 24 h in CE and LV strokes, whereas in SV stroke IL-1 genes decrease at 24 h and then recover at times after 24 h. An example of the complexity of the profiles is observed for whole blood in SV strokes where both a peaking of positive regulation of cytokines is noted at 24 h with decreases thereafter, as well as consistent increase over time of pathways associated with negative regulation of cytokines. Thus, either different cell types in whole blood are responding differently, and/or different cytokines are being regulated differently. This serves to emphasize how complex the temporal and cellular responses are to stroke, and to warn against time-agnostic approaches to treatment at single times of strokes of different causes. Although the analyses subgrouping by IS cause and time point could be considered as exploratory due to smaller sample size, they highlight an important consideration for future studies to determine biomarker and treatment targets.
Limitations
Though it was possible to identify markers for monocyte and neutrophil polarization in either direction (inflammatory or anti-inflammatory types), it was not possible to define a unique shift towards a specific subtype since we expect that both M1/M2 and N1/N2 phenotypes are present in the samples of pooled cells in the ~3 day time window studied here. It is possible that studying longer times after stroke might allow detection of these shifts in gene expression responsible for the evolution into polarized M2 or N2 phenotypes. Moreover, future single-cell RNA-seq studies should enable a better identification of the different cell subpopulations present at each time point.
This study employed single samples from single patients at different times. Thus, the changes of gene expression represent the average of multiple patients over the times stipulated. We have shown that individual genetic variation plays a key role in the specific expression response after stroke [130]. Our approach presented here based on the availability of samples would include inter-individual variability of expression as compared to differences measured in a single subject over time. Thus, future studies should consider a serial sample approach in every subject over longer periods of time to investigate gene expression dynamics in every subject and see how these vary between subjects as a function of time and its relation to infarct volume, evolving NIHSS, and other clinical variables.
A strength of the study was the multiple analytical approaches used, including differential expression, GEDI, SOM, and WGCNA, with each providing insight to the complexity and dynamics of gene expression changes after stroke. SOM in particular was able to demonstrate how pathways changed over time for each cell type and cause of stroke. WGCNA identified modules of co-expressed genes and associated hub genes that changed over time. However, even WGCNA had some weaknesses in that few hub genes and only one time-associated module were identified in neutrophils. Both these results likely arose from differences in size of input list of genes or specific parameters used in the models, since we indeed observed slightly more time-regulated DEGs in neutrophils compared to monocytes when studied by time point. Since WGCNA identified modules that correlated with the continuous time parameter, it may have missed modules with complex dynamic profiles (like the identified SOM profiles). Pathway and functional analyses helped to interpret aggregate shifts of groups of genes, while reducing the impact of potential false positives of individual genes. Though these results are likely to be more reliable, they are hampered in the current study by the small sample size for many of the time points particularly for those considering different causes of stroke in the three sample types examined (isolated monocytes, isolated neutrophils, whole peripheral blood). For example, subgrouping of TP3 resulted in a group of female only subjects, which may impact the genes discovered for that time point in particular, as sex influences gene expression and immune response after stroke [25]. Because of demographics and group size, gene identification for the first 48 h is more robust than over 48 h. Limitations such as this could be addressed in future studies with much larger sample sizes of blood samples collected over multiple times in a true longitudinal fashion from each subject. Still, the use of a parallel analysis approach like co-expression networks, which accounted for change of time as a continuous variable, also supported gene discovery for later times. Overall, however, each of these analytical approaches provided unique insights into the pathophysiology of stroke based upon cells in blood which modulate the peripheral clotting and immune responses in stroke.
Conclusions
-
We identified key genes associated with time in the response from leukocytes after ischemic stroke. Some of these correlated with stroke severity.
-
Differentially expressed genes have distinctive trajectories for all IS etiologies analyzed, which allows refinement of critical genes and functions at specific time points after stroke.
-
Altogether, the identified changes in gene expression and pathways over time are critical for understanding how the immune and clotting systems change dynamically after stroke, and point to the complexities of identifying biomarkers and treatment targets for stroke.
Availability of data and materials
Gene expression data and clinical parameters of the subjects studied are available at https://doi.org/10.5281/zenodo.6977557 [131]. The code used for analyses is available at https://github.com/paulicm-UCD/gene_networks_stroke. Other data used and/or analyzed during the current study are available from the corresponding author on reasonable request.
Abbreviations
- ATM:
-
Ataxia-telangiectasia mutated kinase
- CE:
-
Cardioembolic
- DE:
-
Differentially expressed
- DEGs:
-
Differentially expressed genes
- GEDI:
-
Gene Expression Dynamics Inspector
- GO:
-
Gene ontology
- IL:
-
Interleukin
- IS:
-
Ischemic stroke
- LV:
-
Large vessel
- NIHSS:
-
National Institutes of Health Stroke Scale
- SOM:
-
Self-organizing maps
- SV:
-
Small vessel
- TP:
-
Time point
- VRFC:
-
Vascular risk factor control
- WGCNA:
-
Weighted Gene Co-expression Network Analysis
References
BCV C, De Silva DA, Macleod MR, et al. Ischaemic stroke. Nat Rev Dis Primers. 2019;5:70.
Qiu Y-M, Zhang C-L, Chen A-Q, et al. Immune Cells in the BBB Disruption After Acute Ischemic Stroke: Targets for Immune Therapy? Front Immunol. 2021;12:678744.
Zera KA, Buckwalter MS. The Local and Peripheral Immune Responses to Stroke: Implications for Therapeutic Development. Neurotherapeutics. 2020;17:414–35.
ElAli A, Jean LBN. The Role of Monocytes in Ischemic Stroke Pathobiology: New Avenues to Explore. Front Aging Neurosci. 2016;8:29.
Park J, Chang JY, Kim JY, et al. Monocyte Transmodulation: The Next Novel Therapeutic Approach in Overcoming Ischemic Stroke? Front Neurol. 2020;11:578003.
Jian Z, Liu R, Zhu X, et al. The Involvement and Therapy Target of Immune Cells After Ischemic Stroke. Front Immunol. 2019;10:2167.
Cai W, Liu S, Hu M, et al. Functional Dynamics of Neutrophils After Ischemic Stroke. Transl Stroke Res. 2020;11:108–21.
Cuartero MI, Ballesteros I, Moraga A, et al. N2 neutrophils, novel players in brain inflammation after stroke: modulation by the PPARγ agonist rosiglitazone. Stroke. 2013;44:3498–508.
Gelderblom M, Leypoldt F, Steinbach K, et al. Temporal and spatial dynamics of cerebral immune cell accumulation in stroke. Stroke. 2009;40:1849–57.
Garcia JH, Liu KF, Yoshida Y, et al. Influx of leukocytes and platelets in an evolving brain infarct (Wistar rat). Am J Pathol. 1994;144:188–99.
Wattananit S, Tornero D, Graubardt N, et al. Monocyte-Derived Macrophages Contribute to Spontaneous Long-Term Functional Recovery after Stroke in Mice. J Neurosci. 2016;36:4182–95.
Gill D, Sivakumaran P, Aravind A, et al. Temporal Trends in the Levels of Peripherally Circulating Leukocyte Subtypes in the Hours after Ischemic Stroke. J Stroke Cerebrovasc Dis. 2018;27:198–202.
Pektezel MY, Yilmaz E, Arsava EM, et al. Neutrophil-to-Lymphocyte Ratio and Response to Intravenous Thrombolysis in Patients with Acute Ischemic Stroke. J Stroke Cerebrovasc Dis. 2019;28:1853–9.
Bi Y, Shen J, Chen S-C, et al. Prognostic value of neutrophil to lymphocyte ratio in acute ischemic stroke after reperfusion therapy. Sci Rep. 2021;11:6177.
Wang A, Quan K, Tian X, et al. Leukocyte subtypes and adverse clinical outcomes in patients with acute ischemic cerebrovascular events. Ann Transl Med. 2021;9:748.
Chen Y, Ren J, Yang N, et al. Eosinophil-to-Monocyte Ratio is a Potential Predictor of Prognosis in Acute Ischemic Stroke Patients After Intravenous Thrombolysis. Clin Interv Aging. 2021;16:853–62.
Dong X, Nao J, Gao Y. Peripheral monocyte count predicts outcomes in patients with acute ischemic stroke treated with rtPA thrombolysis. Neurotox Res. 2020;37:469–77.
Carmona-Mora P, Ander BP, Jickling GC, et al. Distinct peripheral blood monocyte and neutrophil transcriptional programs following intracerebral hemorrhage and different etiologies of ischemic stroke. J Cereb Blood Flow Metab. 2021;41:1398–416.
Stamova B, Ander BP, Jickling G, et al. The intracerebral hemorrhage blood transcriptome in humans differs from the ischemic stroke and vascular risk factor control blood transcriptomes. J Cereb Blood Flow Metab. 2019;39:1818–35.
Xu H, Tang Y, Liu D-Z, et al. Gene expression in peripheral blood differs after cardioembolic compared with large-vessel atherosclerotic stroke: biomarkers for the etiology of ischemic stroke. J Cereb Blood Flow Metab. 2008;28:1320–8.
Jickling GC, Stamova B, Ander BP, et al. Profiles of lacunar and nonlacunar stroke. Ann Neurol. 2011;70:477–85.
Dykstra-Aiello C, Jickling GC, Ander BP, et al. Intracerebral Hemorrhage and Ischemic Stroke of Different Etiologies Have Distinct Alternatively Spliced mRNA Profiles in the Blood: a Pilot RNA-seq Study. Transl Stroke Res. 2015;6:284–9.
Tang Y, Xu H, Du X, et al. Gene expression in blood changes rapidly in neutrophils and monocytes after ischemic stroke in humans: a microarray study. J Cereb Blood Flow Metab. 2006;26:1089–102.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. 2008;9:559.
Stamova B, Jickling GC, Ander BP, et al. Gene expression in peripheral immune cells following cardioembolic stroke is sexually dimorphic. PLoS One. 2014;9:e102550.
Streett DA, Petersen KR, Gerritsen AT, et al. expHTS: analysis of high throughput sequence data in an experimental framework. In: Proceedings of the 6th ACM Conference on Bioinformatics, Computational Biology and Health Informatics. New York: Association for Computing Machinery; 2015. p. 523–4.
Dobin A, Davis CA, Schlesinger F, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.
Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30:923–30.
Thompson WA. The Problem of Negative Estimates of Variance Components. Ann Math Stat. 1962;33:273–89.
Tamhane A, Dunlop D. Statistics and data analysis: From elementary to intermediate, https://www.scholars.northwestern.edu/en/publications/statistics-and-data-analysis-from-elementary-to-intermediate-2(2000, Accessed 22 Sept 2021).
Eichler GS, Huang S, Ingber DE. Gene Expression Dynamics Inspector (GEDI): for integrative analysis of expression profiles. Bioinformatics. 2003;19:2321–2.
Kohonen T. Self-Organizing Maps. Springer Science & Business Media; 2001.
Chtanova T, Newton R, Liu SM, et al. Identification of T cell-restricted genes, and signatures for different T cell responses, using a comprehensive collection of microarray datasets. J Immunol. 2005;175:7837–47.
Watkins NA, Gusnanto A, de Bono B, et al. A HaemAtlas: characterizing gene expression in differentiated human blood cells. Blood. 2009;113:e1–9.
Zhang B, Horvath S. A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol. 2005;4:Article17.
Langfelder P. Signed vs. unsigned topological overlap matrix: technical report. 2013. https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/TechnicalReports/signedTOM.pdf. Accessed March 2022.
Yang Y, Han L, Yuan Y, et al. Gene co-expression network analysis reveals common system-level properties of prognostic genes across cancer types. Nat Commun. 2014;5:3231.
Langfelder P, Mischel PS, Horvath S. When is hub gene selection better than standard meta-analysis? PLoS One. 2013;8:e61505.
Krishnan A, Zhang R, Yao V, et al. Genome-wide prediction and functional characterization of the genetic basis of autism spectrum disorder. Nat Neurosci. 2016;19:1454–62.
Gordon S, Taylor PR. Monocyte and macrophage heterogeneity. Nat Rev Immunol. 2005;5:953–64.
Uhlen M, Oksvold P, Fagerberg L, et al. Towards a knowledge-based Human Protein Atlas. Nat Biotechnol. 2010;28:1248–50.
Piccard H, Muschel RJ, Opdenakker G. On the dual roles and polarized phenotypes of neutrophils in tumor development and progression. Crit Rev Oncol Hematol. 2012;82:296–309.
Wang N, Liang H, Zen K. Molecular mechanisms that influence the macrophage M1–M2 polarization balance. Front Immunol. 2014;5:614.
Sierra-Filardi E, Nieto C, Domínguez-Soto A, et al. CCL2 shapes macrophage polarization by GM-CSF and M-CSF: identification of CCL2/CCR2-dependent gene expression profile. J Immunol. 2014;192:3858–67.
Greter M, Lelios I, Croxford AL. Microglia versus myeloid cell nomenclature during brain inflammation. Front Immunol. 2015;6:249.
Yan Y-P, Sailor KA, Lang BT, et al. Monocyte chemoattractant protein-1 plays a critical role in neuroblast migration after focal cerebral ischemia. J Cereb Blood Flow Metab. 2007;27:1213–24.
Goerdt S, Politz O, Schledzewski K, et al. Alternative versus classical activation of macrophages. Pathobiology. 1999;67:222–6.
Fridlender ZG, Sun J, Kim S, et al. Polarization of tumor-associated neutrophil phenotype by TGF-beta: ‘N1’ versus ‘N2’ TAN. Cancer Cell. 2009;16:183–94.
Ohms M, Möller S, Laskay T. An attempt to polarize human neutrophils toward n1 and n2 phenotypes in vitro. Front Immunol. 2020;11:532.
Lin S-H, Chuang H-Y, Ho J-C, et al. Treatment with TNF-α inhibitor rectifies M1 macrophage polarization from blood CD14+ monocytes in patients with psoriasis independent of STAT1 and IRF-1 activation. J Dermatol Sci. 2018;91:276–84.
Nayak AR, Kashyap RS, Kabra D, et al. Time course of inflammatory cytokines in acute ischemic stroke patients and their relation to inter-alfa trypsin inhibitor heavy chain 4 and outcome. Ann Indian Acad Neurol. 2012;15:181–5.
Perini F, Morra M, Alecci M, et al. Temporal profile of serum anti-inflammatory and pro-inflammatory interleukins in acute ischemic stroke patients. Neurol Sci. 2001;22:289–96.
Liu X, Liu J, Zhao S, et al. Interleukin-4 is essential for Microglia/Macrophage M2 polarization and long-term recovery after cerebral ischemia. Stroke. 2016;47:498–504.
Gerszten RE, Garcia-Zepeda EA, Lim YC, et al. MCP-1 and IL-8 trigger firm adhesion of monocytes to vascular endothelium under flow conditions. Nature. 1999;398:718–23.
Pantoni L, Sarti C, Inzitari D. Cytokines and cell adhesion molecules in cerebral ischemia: experimental bases and therapeutic perspectives. Arterioscler Thromb Vasc Biol. 1998;18:503–13.
Hotter B, Hoffmann S, Ulm L, et al. IL-6 Plasma Levels Correlate With Cerebral Perfusion Deficits and Infarct Sizes in Stroke Patients Without Associated Infections. Front Neurol. 2019;10:83.
Bustamante A, Sobrino T, Giralt D, et al. Prognostic value of blood interleukin-6 in the prediction of functional outcome after stroke: a systematic review and meta-analysis. J Neuroimmunol. 2014;274:215–24.
Chomarat P, Banchereau J, Davoust J, et al. IL-6 switches the differentiation of monocytes from dendritic cells to macrophages. Nat Immunol. 2000;1:510–4.
Mitani H, Katayama N, Araki H, et al. Activity of interleukin 6 in the differentiation of monocytes to macrophages and dendritic cells. Br J Haematol. 2000;109:288–95.
Casella G, Garzetti L, Gatta AT, et al. IL4 induces IL6-producing M2 macrophages associated to inhibition of neuroinflammation in vitro and in vivo. J Neuroinflammation. 2016;13:139.
Braune J, Weyer U, Hobusch C, et al. IL-6 Regulates M2 polarization and local proliferation of adipose tissue macrophages in obesity. J Immunol. 2017;198:2927–34.
Dugue R, Nath M, Dugue A, et al. Roles of Pro- and Anti-inflammatory Cytokines in Traumatic Brain Injury and Acute Ischemic Stroke. In: Abreu GEA, editor. Mechanisms of Neuroinflammation. Rijeka: IntechOpen; 2017.
Lee GA, Lin T-N, Chen C-Y, et al. Interleukin 15 blockade protects the brain from cerebral ischemia-reperfusion injury. Brain Behav Immun. 2018;73:562–70.
Neely GG, Robbins SM, Amankwah EK, et al. Lipopolysaccharide-stimulated or granulocyte-macrophage colony-stimulating factor-stimulated monocytes rapidly express biologically active IL-15 on their cell surface independent of new protein synthesis. J Immunol. 2001;167:5011–7.
Detmers PA, Lo SK, Olsen-Egbert E, et al. Neutrophil-activating protein 1/interleukin 8 stimulates the binding activity of the leukocyte adhesion receptor CD11b/CD18 on human neutrophils. J Exp Med. 1990;171:1155–62.
Detmers PA, Powell DE, Walz A, et al. Differential effects of neutrophil-activating peptide 1/IL-8 and its homologues on leukocyte adhesion and phagocytosis. J Immunol. 1991;147:4211–7.
Boey H, Rosenbaum R, Castracane J, et al. Interleukin-4 is a neutrophil activator. J Allergy Clin Immunol. 1989;83:978–84.
ElKassar N, Gress RE. An overview of IL-7 biology and its use in immunotherapy. J Immunotoxicol. 2010;7:1–7.
Girard D, Paquet ME, Paquin R, et al. Differential effects of interleukin-15 (IL-15) and IL-2 on human neutrophils: modulation of phagocytosis, cytoskeleton rearrangement, gene expression, and apoptosis by IL-15. Blood. 1996;88:3176–84.
Shevach EM. Application of IL-2 therapy to target T regulatory cell function. Trends Immunol. 2012;33:626–32.
Lanfranconi S, Locatelli F, Corti S, et al. Growth factors in ischemic stroke. J Cell Mol Med. 2011;15:1645–87.
Xie G-H, Dai H-J, Liu F, et al. A dual role of atm in ischemic preconditioning and ischemic injury. Cell Mol Neurobiol. 2020;40:785–99.
Yin K-J, Chen S-D, Lee J-M, et al. ATM gene regulates oxygen-glucose deprivation-induced nuclear factor-kappaB DNA-binding activity and downstream apoptotic cascade in mouse cerebrovascular endothelial cells. Stroke. 2002;33:2471–7.
Sermeus A, Michiels C. Reciprocal influence of the p53 and the hypoxic pathways. Cell Death Dis. 2011;2:e164.
Filichia E, Shen H, Zhou X, et al. Forebrain neuronal specific ablation of p53 gene provides protection in a cortical ischemic stroke model. Neuroscience. 2015;295:1–10.
Xin M, Feng J, Hao Y, et al. Cyclic adenosine monophosphate in acute ischemic stroke: some to update, more to explore. J Neurol Sci. 2020;413:116775.
Jiang S, Li T, Ji T, et al. AMPK: Potential Therapeutic Target for Ischemic Stroke. Theranostics. 2018;8:4535–51.
Elgebaly MM. Ephrin-Eph signaling as a novel neuroprotection path in ischemic stroke. J Mol Neurosci. 2020;70:2001–6.
Krämer A, Green J, Pollard J Jr, et al. Causal analysis approaches in ingenuity pathway analysis. Bioinformatics. 2014;30:523–30.
Evans MA, Kim HA, Ling YH, et al. Vitamin D3 Supplementation Reduces Subsequent Brain Injury and Inflammation Associated with Ischemic Stroke. Neuromolecular Med. 2018;20:147–59.
Fu J, Xue R, Gu J, et al. Neuroprotective effect of calcitriol on ischemic/reperfusion injury through the NR3A/CREB pathways in the rat hippocampus. Mol Med Rep. 2013;8:1708–14.
Sadeghian N, Shadman J, Moradi A, et al. Calcitriol protects the Blood-Brain Barrier integrity against ischemic stroke and reduces vasogenic brain edema via antioxidant and antiapoptotic actions in rats. Brain Res Bull. 2019;150:281–9.
Liu H, Wu X, Luo J, et al. Pterostilbene attenuates astrocytic inflammation and neuronal oxidative injury after ischemia-Reperfusion by Inhibiting NF-κB Phosphorylation. Front Immunol. 2019;10:2408.
McCormack D, McFadden D. A review of pterostilbene antioxidant activity and disease modification. Oxid Med Cell Longev. 2013;2013:575482.
Liu J, Xu J, Mi Y, et al. Pterostilbene alleviates cerebral ischemia and reperfusion injury in rats by modulating microglial activation. Food Funct. 2020;11:5432–45.
Tu Q, Le D, Wang C, et al. Pterostilbene attenuates ischemic stroke by modulating miR-21-5p/PDCD4 axis in vivo and in vitro. J Funct Foods. 2020;75:104275.
Kovacs-Litman A, Vonderwalde I. Monocyte-derived macrophages modulate inflammation and promote long-term functional recovery in a mouse model of ischemia. J Neurosci. 2016;36:9757–9.
Speyer CL, Gao H, Rancilio NJ, et al. Novel chemokine responsiveness and mobilization of neutrophils during sepsis. Am J Pathol. 2004;165:2187–96.
Li Q, Gao X, Luo X, et al. Identification of Hub Genes Associated with Immune Infiltration in Cardioembolic Stroke by Whole Blood Transcriptome Analysis. Dis Markers. 2022;2022:8086991.
Dornier E, Coumailleau F, Ottavi J-F, et al. TspanC8 tetraspanins regulate ADAM10/Kuzbanian trafficking and promote Notch activation in flies and mammals. J Cell Biol. 2012;199:481–96.
LeBlang CJ, Medalla M, Nicoletti NW, et al. Reduction of the RNA Binding Protein TIA1 Exacerbates Neuroinflammation in Tauopathy. Front Neurosci; 14. Epub ahead of print 2020. https://doi.org/10.3389/fnins.2020.00285.
Piecyk M, Wax S, Beck AR, et al. TIA-1 is a translational silencer that selectively regulates the expression of TNF-alpha. EMBO J. 2000;19:4154–63.
Gueydan C, Droogmans L, Chalon P, et al. Identification of TIAR as a protein binding to the translational regulatory AU-rich element of tumor necrosis factor α mRNA. J Biol Chem. 1999;274:2322–6.
Ding Z, Gillespie LL, Paterno GD. Human MI-ER1 alpha and beta function as transcriptional repressors by recruitment of histone deacetylase 1 to their conserved ELM2 domain. Mol Cell Biol. 2003;23:250–8.
Fang L, Sun X, Wang Y, et al. RMI1 contributes to DNA repair and to the tolerance to camptothecin. FASEB J. 2019;33:5561–70.
Chen M-H, Raffield LM, Mousas A, et al. Trans-ethnic and ancestry-specific blood-cell genetics in 746,667 individuals from 5 global populations. Cell. 2020;182:1198–1213.e14.
Chen S, Liu W, Huang Y. Identification and external validation of a prognostic signature associated with DNA repair genes in gastric cancer. Sci Rep. 2021;11:7141.
Xu T, Liu W, Yang C, et al. Lipid raft-associated β-adducin is required for PSGL-1-mediated neutrophil rolling on P-selectin. J Leukoc Biol. 2015;97:297–306.
Durocher M, Ander BP, Jickling G, et al. Inflammatory, regulatory, and autophagy co-expression modules and hub genes underlie the peripheral immune response to human intracerebral hemorrhage. J Neuroinflammation. 2019;16:56.
Sun J, Huang W, Yang S-F, et al. Gαi1 and Gαi3mediate VEGF-induced VEGFR2 endocytosis, signaling and angiogenesis. Theranostics. 2018;8:4695–709.
Minoda Y, Sakurai H, Kobayashi T, et al. An F-box protein, FBXW5, negatively regulates TAK1 MAP3K in the IL-1beta signaling pathway. Biochem Biophys Res Commun. 2009;381:412–7.
Boutin H, LeFeuvre RA, Horai R, et al. Role of IL-1alpha and IL-1beta in ischemic brain damage. J Neurosci. 2001;21:5528–34.
Touzani O, Boutin H, Chuquet J, et al. Potential mechanisms of interleukin-1 involvement in cerebral ischaemia. J Neuroimmunol. 1999;100:203–15.
Murray KN, Parry-Jones AR, Allan SM. Interleukin-1 and acute brain injury. Front Cell Neurosci. 2015;9:18.
Ahmed YL, Schleich S, Bohlen J, et al. DENR-MCTS1 heterodimerization and tRNA recruitment are required for translation reinitiation. PLoS Biol. 2018;16:e2005160.
Jickling GC, Ander BP, Shroff N, et al. Leukocyte response is regulated by microRNA let7i in patients with acute ischemic stroke. Neurology. 2016;87:2198–205.
Bernstein DL, Jiang X, Rom S. let-7 microRNAs: Their Role in Cerebral and Cardiovascular Diseases, Inflammation, Cancer, and Their Regulation. Biomedicines; 9. Epub ahead of print 26 May 2021.https://doi.org/10.3390/biomedicines9060606.
Engler-Chiurazzi EB, Monaghan KL, Wan ECK, et al. Role of B cells and the aging brain in stroke recovery and treatment. Geroscience. 2020;42:1199–216.
Uhlén M, Fagerberg L, Hallström BM, et al. Tissue-based map of the human proteome. Science,https://science.sciencemag.org/content/347/6220/1260419?ijkey=0c35588049707564cc6251a2021287a5f153dffd&keytype2=tf_ipsecsha(2015, Accessed 22 Sept 2021).
Hayashi F, Smith KD, Ozinsky A, et al. The innate immune response to bacterial flagellin is mediated by Toll-like receptor 5. Nature. 2001;410:1099–103.
Münzer P, Walker-Allgaier B, Geue S, et al. PDK1 Determines Collagen-Dependent Platelet Ca2+ Signaling and Is Critical to Development of Ischemic Stroke In Vivo. Arterioscler Thromb Vasc Biol. 2016;36:1507–16.
Mao B, Wu W, Davidson G, et al. Kremen proteins are Dickkopf receptors that regulate Wnt/beta-catenin signalling. Nature. 2002;417:664–7.
Bosio Y, Berto G, Camera P, et al. PPP4R2 regulates neuronal cell differentiation and survival, functionally cooperating with SMN. Eur J Cell Biol. 2012;91:662–74.
Selvaraj UM, Poinsatte K, Torres V, et al. Heterogeneity of B cell functions in stroke-related risk, prevention, injury, and repair. Neurotherapeutics. 2016;13:729–47.
Ortega SB, Torres VO, Latchney SE, et al. B cells migrate into remote brain areas and support neurogenesis and functional recovery after focal stroke in mice. Proc Natl Acad Sci U S A. 2020;117:4983–93.
Zbesko JC, Frye JB, Becktel DA, et al. IgA natural antibodies are produced following T-cell independent B-cell activation following stroke. Brain Behav Immun. 2021;91:578–86.
Bornstein NM, Aronovich B, Korczyn AD, et al. Antibodies to brain antigens following stroke. Neurology. 2001;56:529–30.
Prüss H, Iggena D, Baldinger T, et al. Evidence of intrathecal immunoglobulin synthesis in stroke: a cohort study. Arch Neurol. 2012;69:714–7.
Iadecola C, Anrather J. The immunology of stroke: from mechanisms to translation. Nat Med. 2011;17:796–808.
Javidi E, Magnus T. Autoimmunity after ischemic stroke and brain injury. Front Immunol. 2019;10:686.
Kalev-Zylinska ML, Symes W, Little KCE, et al. Stroke patients develop antibodies that react with components of N-methyl-D-aspartate receptor subunit 1 in proportion to lesion size. Stroke. 2013;44:2212–9.
Becker KJ, Kalil AJ, Tanzi P, et al. Autoimmune responses to the brain after stroke are associated with worse outcome. Stroke. 2011;42:2763–9.
Feng Q, Cerione RA. Chapter 218 - Cdc42 and Its Cellular Functions. In: Bradshaw RA, Dennis EA, editors. Handbook of Cell Signaling. 2nd ed. San Diego: Academic Press; 2010. p. 1785–94.
Yuan D, Hu K, Loke CM, et al. Interruption of endolysosomal trafficking leads to stroke brain injury. Exp Neurol. 2021;345:113827.
Wang X-X, Wang F, Mao G-H, et al. NADPH is superior to NADH or edaravone in ameliorating metabolic disturbance and brain injury in ischemic stroke. Acta Pharmacol Sin. Epub ahead of print 24 June 2021.https://doi.org/10.1038/s41401-021-00705-5.
Nalamolu KR, Challa SR, Fornal CA, et al. Attenuation of the Induction of TLRs 2 and 4 Mitigates Inflammation and Promotes Neurological Recovery After Focal Cerebral Ischemia. Transl Stroke Res. 2021;12:923–36.
Gao W, Xiong Y, Li Q, et al. Inhibition of Toll-Like receptor signaling as a promising therapy for inflammatory diseases: a journey from molecular to nano therapeutics. Front Physiol. 2017;8:508.
Ashayeri Ahmadabad R, Mirzaasgari Z, Gorji A, et al. Toll-like receptor signaling pathways: novel therapeutic targets for cerebrovascular disorders. Int J Mol Sci; 22. Epub ahead of print 7 2021. https://doi.org/10.3390/ijms22116153.
Gesuete R, Kohama SG, Stenzel-Poore MP. Toll-like receptors and ischemic brain injury. J Neuropathol Exp Neurol. 2014;73:378–86.
Amini H, Shroff N, Stamova B, et al. Genetic variation contributes to gene expression response in ischemic stroke: an eQTL study. Ann Clin Transl Neurol. 2020;7:1648–60.
Carmona-Mora P, Knepp B, Jickling GC, Zhan X, Hakoupian M, Hull H, et al. Gene expression in monocytes, neutrophils and whole blood after stroke. Zenodo. 2022. https://doi.org/10.5281/zenodo.6977557.
Acknowledgements
We thank the patients and families that participated in the study. The sequencing was carried out by the DNA Technologies and Expression Analysis Core at the UC Davis Genome Center, supported by NIH Shared Instrumentation grant 1S10OD010786-01. Monocyte and neutrophil clip art from Servier https://smart.servier.com/ and licensed under CC-BY 3.0 Unported https://creativecommons.org/licenses/by/3.0/. We thank the Biostatistics group from the Clinical and Translational Science Center at UC Davis for their support on statistical considerations and analyses (supported by the National Center for Advancing Translational Sciences, NIH grant number UL1 TR001860).
Funding
These studies were supported by grants from the National Institutes of Health (NS097000, NS101718, NS075035, NS079153, and NS106950 to FRS, BSS, BPA, and GCJ). The funding agencies had no influence on the study design, analyses, interpretation of data, or in writing the manuscript.
Author information
Authors and Affiliations
Contributions
Conceived and designed the experiments, and edited the manuscript: PC-M, BPA, BS, FRS. Performed the experiments and/or reviewed data: PC-M, BK, MH, HH, XZ, NA, HA. Clinical characterization: GCJ, FRS. Wrote the first draft: PC-M. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
This study was reviewed and approved by the University of California, Davis Institutional Review Board (IRB ID# 248994-41). All participants or a legally authorized representative provided written informed consent to participate in this study.
Consent for publication
Not applicable
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Additional file 1: Table S1.
Differentially expressed genes (DEGs) significant to time point (TP) p-value < 0.02 on the contrasts IS TP vs. VRFC (TP0). Table S2. Enriched canonical pathways overrepresented in all time points (TP) (Fisher's p-value < 0.05). Table S3. Predicted upstream regulators in all time points (TP) (Fishers' p-value < 0.05). Table S4. Differentially expressed genes (DEGs) distributed per GEDI map tile based on Pearson's correlation of gene expression values. Table S5. Differentially expressed genes (DEGs) distributed per SOM profile based on Euclidean distance. Table S6. GO terms for biological processes enriched in the SOM profiles from Fig. 6. Table S7. Subject demographics and relevant clinical characteristics in the time bins analyzed after regrouping per IS etiology. Table S8. Differentially expressed genes (DEGs) significant to time point (TP) p-value < 0.02 on the contrasts IS TP-IS etiology vs. VRFC (TP0). Table S9. DEGs distributed per SOM profile in all IS etiologies, based on Euclidean distance. Table S10. GO terms for biological processes enriched in the SOM profiles from Fig. 9, for all IS etiologies. Table S11. Genes in WGCNA modules associated time (h). Table S12. Hub genes in the modules significant to time (h). Table S13. HumanBase functional clustering of hub genes from time-associated modules and their corresponding gene ontology terms. Table S14. Correlation between NIHSS at admission and expression of time-associated hub genes.
Additional file 2: Supplemental Figure 1.
Percentage of DEGs classified per biotype. For every time point and comparison between time points, the percentage of biotype categories for DEGs are shown in the bar plots and corresponding table for (A) monocytes, (B) neutrophils and (C) whole blood. Ig: immunoglobulin; TcR: T-cell receptor; TEC: To be Experimentally Confirmed. Supplemental Figure 2. Predicted upstream regulators for monocytes. (A) Venn diagram represents all upstream regulators predicted at 0-24 h, 24-48 h, and >48 h (Fisher’s p-value<0.05). (B) Top over-represented regulators common to all time points (Fisher’s p-value<0.05). (C) Top over-represented regulators that are specific for each time point (Fisher’s p-value<0.05). Up arrows indicate predicted significant activation (z ≥ 2) or down arrow - significant inhibition (z ≤ −2). White cells indicate no direction can be predicted. (*= Benjamini-Hochberg corrected p-value<0.05). R: receptor; N.R: nuclear receptor. Supplemental Figure 3. Predicted upstream regulators for neutrophils. (A) Venn diagram represents all upstream regulators predicted at 0-24 h, 24-48 h, and >48 h (Fisher’s p-value<0.05). (B) Top over-represented regulators shared at all time points (Fisher’s p-value<0.05). (C) Top over-represented regulators that are specific for each time point (Fisher’s p-value<0.05). Up arrows indicate predicted significant activation (z ≥ 2) or down arrow - significant inhibition (z ≤ −2). White cells indicate no direction can be predicted. (*= Benjamini-Hochberg corrected p-value<0.05). R: receptor; endo: endogenous; reg: regulator. Supplemental Figure 4. Predicted upstream regulators for whole blood. (A) Venn diagram represents all upstream regulators predicted at 0-24 h, 24-48 h, and >48 h (Fisher’s p-value<0.05). (B) Only 2 predicted regulators (Fisher’s p-value<0.05) are shared for the 3 time points. (C) Top over-represented regulators that are specific for each time point (Fisher’s p-value<0.05). Up arrows indicate predicted significant activation (z ≥ 2) or down arrow - significant inhibition (z ≤ −2). White cells indicate no direction can be predicted. R: receptor; reg: regulator; endo: endogenous. Supplemental Figure 5. tPA administration is not associated with differential expression. A) Characteristics of the five cases where tPA was administered. Principal Component Analyses (PCA) using the DEGs from the time points 0-24 and >24 h in LV and SV strokes in B) monocytes, C) neutrophils, and D) whole blood samples. Cases where tPA was administered are colored in orange. Supplemental Figure 6. Co-expression network construction. WGCNA Soft-thresholding power scale free topology fit (left panel) and mean connectivity (right panel) plots for all genes analyzed in (A) monocytes, (B) neutrophils and (C) whole blood. Reference lines are at 0.8 (red) and 0.9 (green) on the left panels and at 100 (red) and 200 (green) on the right panels. Supplemental Figure 7. Workflow of the study. Monocytes and neutrophils were isolated by flow cytometry (scatter plots adapted in part from Carmona-Mora et al. [18]. RNA-seq from the isolated cell samples and whole blood was used for two parallel approaches: differential expression and co-expression networks. For the first one, patients were split into time point groups and self-organizing maps were constructed with the differentially expressed genes. For the construction of co-expression networks, time was considered as a continuous variable.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Carmona-Mora, P., Knepp, B., Jickling, G.C. et al. Monocyte, neutrophil, and whole blood transcriptome dynamics following ischemic stroke. BMC Med 21, 65 (2023). https://doi.org/10.1186/s12916-023-02766-1
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12916-023-02766-1