Monocyte, neutrophil, and whole blood transcriptome dynamics following ischemic stroke

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

Table 1 Subject demographics and relevant clinical characteristics in each of the time windows analyzed

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

Fig. 1

Differential expression across time points. Differentially expressed genes (DEGs) across time points versus VRFC (0–24 h (TP1), 24–48 h (TP2), and >48 h (TP3)) in monocytes (A), neutrophils (B), and whole blood (C). Venn diagrams of the numbers of DEGs at each time point are shown on the left, and on the right, the corresponding bar plots are shown for the numbers of up- and downregulated DEGs found per time point

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

Fig. 2

Enriched pathways in monocytes. A Venn diagram represents all enriched pathways in monocytes at 0–24 h, 24–48 h, and >48 h (Fisher’s p-value<0.05). In total, 132 pathways overlapped all three time points. B Top 10 pathways enriched and shared at all time points. C Top 10 over-represented pathways with predicted inhibition (z≤ −2, blue shades) in every time point. D Top 10 over-represented pathways unique to each time point (Fisher’s p-value<0.05). Down arrows indicate significant predicted inhibition of pathway (z≤ −2). Blue cells indicate negative z-score, white cells indicate no direction can be predicted, and grey cells are those where no prediction can be performed. (*= Benjamini–Hochberg corrected p-value<0.05). FcγR-Phagocytosis: Fcγ receptor-mediated Phagocytosis in Macrophages and Monocytes; Ind, induction

From the pathways exclusively regulated at a single TP, suppressed oxidative phosphorylation, suppressed Wnt/Ca²⁺, 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).

Fig. 3

Enriched pathways in neutrophils. A Venn diagram represents all enriched pathways in neutrophils at 0–24 h, 24–48 h, and >48 h (Fisher’s p-value<0.05). Fourteen pathways overlapped all time points. B Top 10 common pathways enriched in all three time points. C Top over-represented pathways with predicted activation (z≥ 2, orange shades) or inhibition (z≤ −2, blue shades) at every time point. D Top 10 over-represented pathways unique for each time point (Fisher’s p-value<0.05). Down arrows indicate predicted significant inhibition (z≤ −2) and up arrows predicted significant activation (z≥ 2). Blue cells indicate negative z-score, and orange positive z-score. White cells indicate no direction can be predicted, and grey indicate that no prediction can be performed (*= Benjamini–Hochberg-corrected p-value<0.05). Reg: regulation

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

Fig. 4

Enriched pathways in whole blood. A Venn diagram represents all enriched pathways at each of the time points (0–24 h, 24–48 h, >48 h; Fisher’s p-value<0.05). There was no overlap of pathways between the time points. B Top over-represented pathways with significant z-score (z ≥ 2, predicted activation) at 0–24 h and 24–48 h. No significant z-scores were predicted for pathways at >48 h. C Top over-represented pathways unique to each TP (Fisher’s p-value<0.05). Orange cells indicate positive z-score for the pathway, and blue indicate negative z-score. White bars indicate no direction can be predicted, and grey bars indicate no prediction can be performed. Up arrows indicate predicted significant activation (z ≥ 2)

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

Fig. 5

Dynamics of gene expression based on GEDI. Tile maps based on self-organizing maps are shown for all time points for monocytes, neutrophils, and whole blood. Each tile within a cell type represents the same set of genes across all 3 time points and VRFCs. Changes in color in a tile over the time points show changes in expression of the member genes over time, while constant color over the time points shows more constant expression levels. The asterisk (*) indicates tiles where immunoglobulin genes are present, and the plus (+) indicates the presence of interleukin receptor genes. On the right, gene density maps represent the number of DEGs per tile (genes clustered based on Pearson’s correlation). Tiles behaving similarly are placed in the same mosaic neighborhood. Letters on the tiles in the gene density map depict presence of specific cell markers in the DEGs in a specific tile (M: monocytes, G: granulocyte= neutrophils, basophils and eosinophils)

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

Fig. 6

Distinctive profiles of gene expression over time after IS based on self-organizing maps (SOM). DEG profiles are shown across time points in monocytes, neutrophils, and whole blood. Dots on each profile represent 0 (VRFC), 0–24 h, 24–48 h, and >48 h time points from left to right. Rows show similar expression profiles across sample types. The top significant Gene Ontology (GO) terms for biological processes are shown. Color corresponds to the percentage of genes in the GO group present in the profiles. Abbreviations: clust.: clustering; endocardial to mesenchymal transition: endocardial cushion to mesenchymal transition; intracell.: intracellular; mol.: molecular; morphog.: morphogenesis; neg.: negative; pos.: positive; prod.: production; prolif.: proliferation; reg.: regulation; reg. of Cdc42 protein signal: regulation of Cdc42 protein signal transduction; reg. of sk. muscle contraction by Ca2+: regulation of skeletal muscle contraction by regulation of release of sequestered calcium ion; res.: response; sig.: signaling

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

Fig. 7

Differentially expressed genes (DEGs) in stroke etiologies. Venn diagrams of the numbers of DEGs at each time point versus VRFCs after cardioembolic (CE), large vessel (LV), and small vessel (SV) strokes in monocytes (A), neutrophils (B), and whole blood (C) (left panels). On the right are the corresponding bar plots for the numbers of up- and downregulated DEGs found per time point and IS cause (0–24 h (TP1), 24–48 h (TP2), and >48 h (TP3))

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

Fig. 8

Distinctive gene expression trajectories over time for different causes of IS based upon self-organizing maps (SOM). DEG trajectories across time (0 h, 0–24 h, >24 h)) are shown for cardioembolic (CE), large vessel (LV), and small vessel (SV) causes of strokes in monocytes, neutrophils, and whole blood. Dots on each profile represent every time point from left to right. Rows show similar expression profiles across stroke etiology in each sample type; empty cells indicate no matching profile was found in that sample type / etiology combination. The top enriched Gene Ontology (GO) terms for biological processes are shown for each profile. The color corresponds to the percentage of genes in the GO group present in the profiles. Abbreviations: Ag: antigen; reg: regulation; diff: differentiation; neg.: negative; pos.: positive; prom: promoter; R: receptor; recomb: recombination; RNP: ribonucleoprotein; sec.: secretion; S.: signal/signaling; TLR: Toll-like receptor; transduc: transduction; transmemb.: transmembrane

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

Fig. 9

Co-expression modules of genes significant to time after stroke, hub genes, and their functional clusters. A All WGCNA-identified gene co-expression modules in monocytes, neutrophils, and whole blood. Modules significant for time are labeled with an asterisk (*) and those additionally significant to another factor with two asterisks (**). Hub genes that varied as a function of time were clustered according to their functional organization (HG: hub gene module) in monocytes (HG-M) as shown in B (left panel), and in whole blood (HG-WB) in C (left panel). The top 3 most specific significant GO terms in each hub gene functional cluster are shown for monocytes in B (right panel) and for whole blood in C (right panel) (images modified from HumanBase). Reg: regulation

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-5_inh, IL-2_inh, and TNF_inh were predicted upstream regulators of gene expression at 0–24 h; IL-1B_inh, IL-2_inh, IL-3, IL-4_inh, IL-5_inh, IL-10_inh, IL-13, IL-15, IL-33_inh, and TNF_inh are predicted upstream regulators at 24–48 h; and IL-2_inh, IL-3, IL-4_inh, IL-13, IL-33, and TNF_inh 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-2_inh, 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 TNF_inh were predicted upstream regulators of gene expression in neutrophils at 24–48 h, and IL-2, IL-3, IL-4, IL-5_inh, IL-10 and TNF_inh 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, cAMP_act, 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.