Differential blood transcriptome modules predict response to corticosteroid therapy in alcoholic hepatitis.

BACKGROUND & AIMS: In patients with severe alcoholic hepatitis (SAH), little is known about the profile of peripheral blood mononuclear cells (PBMCs) at baseline and during corticosteroid therapy, among those who can be treated successfully with steroids (steroid-responders [R] and those who cannot (steroid-non-responders [NR]); 2 groups with different outcomes.
METHODS: We performed RNA-seq analysis in PBMCs from 32 patients with definite SAH, at baseline and after 7 days of corticosteroids. The data were sorted into R and NR (n = 16, each group) using the Lille model and 346 blood transcription modules (BTMs) were identified. BTMs are predefined modules of highly co-expressed PBMC genes, which can determine specific immune cell types and cellular functions. The activity of each BTM was taken as the mean value of its member genes.
RESULTS: At baseline, 345 BTMs had higher activity (i.e. were upregulated) in NR relative to R. The 100 most upregulated BTMs in NR, included several modules related to lymphoid lineage (T, B, and natural killer [NK] cells), modules for cell division and mitochondrial respiratory electron transport chain (ETC, relating to energy production), but only a few modules of myeloid cells. Correlation studies of BTM activities found features of significantly greater activation/proliferation and differentiation for T and B cells in NR relative to R. After 7 days of corticosteroids, NR had no significant changes in BTM activities relative to baseline, whereas R had downregulation of BTMs related to innate and adaptive immunity.
CONCLUSIONS: At baseline and during corticosteroid therapy, increased activity in the PBMCs of gene modules related to activation/proliferation and differentiation of T and B cells, NK cells, and mitochondrial ETC, is a hallmark of SAH patients who are steroid-non-responders. LAY
SUMMARY: Patients with severe alcoholic hepatitis receive steroid therapy as the main line of treatment; however, this treatment is ineffective in some patients. This only becomes apparent after 7 days of steroid therapy. We have developed an approach where it can be estimated if a patient is going to respond or not to steroid therapy using the gene expression information of blood cells. This method will allow clinicians to assess the response of patients to steroids earlier, and will help them in adopting alternate strategies if the treatment is found to be ineffective in a particular patient.

Introduction

The mainstay of treatment in severe alcoholic hepatitis (SAH) is corticosteroid therapy, which is targeted towards suppression of inflammation caused by alcohol overdose. The response to corticosteroids is assessed after 7 days of therapy according to the Lille model, which is based on characteristics which can be assessed at the patient’s bedside. However, ~40% of patients with SAH do not respond to corticosteroids (hereafter, called non-responders [NR]) and have worse outcomes than patients who do respond to this therapy. However, in responders to corticosteroids (R), the treatment reduces the risk of death within 28 days of treatment, but not in the following 6 months. In addition, in patients with SAH, corticosteroid therapy increases the risk of bacterial infection and infection-related death., Together, these findings indicate that there is an urgent need to develop new therapies for patients with SAH. A prerequisite to finding new therapies could be the identification of molecular and cellular differences in the components of the immune and inflammatory responses between R and NR, at baseline and during corticosteroid therapy. We have shown previously that NR have a baseline underexpression of the glucocorticoid receptor, which is attributable to a regulator protein BAG1 (BCL2 associated athanogene), suggesting a molecular basis for resistance to endogenous and exogenous corticosteroids. In the same study, the baseline total leucocyte count was greater in NR than R, whilst there was no between-group difference in the differential neutrophil count, suggesting a contribution of mononuclear cells, that is monocytes or lymphocytes or both, to the higher baseline total leucocyte count in NR. Of note, it has been shown that, at baseline, peripheral blood lymphocytes from NR had a greater propensity to proliferate than lymphocytes from R.

High-throughput methods such as transcriptomics have been used to capture the variability of the immune system in various liver diseases., Recently, using high-throughput transcriptomics of whole blood in patients with various autoimmune or infectious diseases, or of peripheral blood mononuclear cells (PBMCs) in healthy subjects receiving vaccination,, has led to the development of gene modules (made of co-expressed genes across various immunological conditions) that were more accurate than the classical identification of differentially expressed genes (DEGs) in determining blood cell subsets which are primary players in the response to immunological stimuli. Such co-expressed gene groups (with fixed gene membership) have been identified as blood transcriptome modules (BTMs). In previous studies, the activity of each BTM module has been taken as the mean value of its member genes.,

In the current study, results of PBMC transcriptomics obtained by RNA sequencing (RNA-seq) were used to determine if there were broad immunological differences discernible at the baseline, to identify patients who would eventually respond or not respond to steroid therapy. In addition, the PBMC transcriptome was assessed in the same patients after 7 days of corticosteroid therapy to evaluate the impact of therapy on immune cells and functions.

Patients and methods

Patient selection

In a cross-sectional study, 102 consecutive patients with clinical and analytical features of SAH were shortlisted at the Institute of Liver and Biliary Sciences, New Delhi, between October 2013 and March, 2015. Patients having gastrointestinal haemorrhage in the past 15 days, moderate to severe ascites, hepatocellular carcinoma or other forms of malignancy, serological marker positivity for HBV, HCV, HIV, or peptic ulcer before/at presentation were excluded. All patients were screened for infection and sepsis by blood cultures, urine cultures, procalcitonin levels, and chest X-ray before enrolment. Transjugular liver biopsy was performed to confirm histological evidence of SAH. Fifty-four patients were characterised as having alcoholic hepatitis histologically and having a Maddrey’s discriminant function (MDF) score of 32 or more, recent onset of jaundice, history of long-standing heavy alcohol consumption that continued until the past 30 days before the onset of jaundice, and significantly raised aspartate aminotransferase and alanine aminotransferase. Baseline demographic profiles were recorded and blood samples were collected (before the start of steroid therapy, prednisolone 40 mg/day) and after therapy at Day 7. The laboratory staff processing the samples were blinded regarding the clinical details. All patients were managed according to the standard of care including intensive care monitoring, high-calorie diet (35–40 cal/kg/day), intravenous albumin, broad-spectrum antibiotics, and corticosteroid therapy for 7 days. The severity of liver disease was assessed using the MDF, Child-Pugh (CP), and model for end-stage liver disease (MELD) scores at enrolment and at follow-up. The response to corticosteroid therapy was assessed at Day 7, by calculating the Lille score. The institutional ethics review board approved the study (protocol number NCT01820208) and informed consent was obtained from all patients enrolled in the study. The patients in the R and NR groups did not show significant differences in clinical characteristics except for total leucocyte count, platelet counts, Lille score, and mortality, as shown in Table S1 in our previous publication.

Of the 54 patients with liver biopsy who underwent corticosteroid therapy, 38 were R) and 16 were NR to this treatment. All 16 NR were included in the transcriptomic study, whereas for logistical and cost reasons only 16 R of the total of 38 were randomly included in the RNA sequencing (RNA-seq) study. There was no difference in baseline characteristics between the 16 R included in the study and the remaining 22 who were not included. From the whole blood, PBMCs were isolated following the Ficoll Hyper gradient protocol. From the isolated PBMCs, total RNA was isolated using the RNAEasy kit (Cat#52304, Qiagen, USA) method.

RNA preparation from PBMCs

PBMCs were obtained at baseline and Day 7 post-therapy in all patients. Briefly, total RNA was extracted from PBMCs using the Trizol method (Invitrogen, Carlsbad, CA, USA). RNA was assessed with the Agilent BioAnalyzer for quantity and quality. For library preparation, we used TruSeq mRNA library preparation kit version 3.0 (Illumina, San Diego, CA, USA) as per manufacturer’s instructions and sequencing was performed on a HiSeq (Illumina) platform. Reads were aligned to the Human Reference Sequence GRCh37.p13 build using Bowtie21. The gene expression data was normalised for within and between samples using the transcripts per million and trimmed mean of M values methods. More details are provided in the Supplementary material.

Use of RNA-seq data

The RNA-seq data obtained in our cohort of 32 patients with SAH (16 R and 16 NR) were used for different types of analyses described below.

Analyses of differential expression among patients

DESeq2 software, version 1.24.0,21 was used for differential expression analysis. The between-group differential expression was calculated using custom in-house script on normalised expression values of every gene including fold-change and log fold-change values (1.5-fold or more), whilst estimating their significance (p <0.05). A list of DEGs was established in each of the following pairwise comparisons: 1 was cross-sectional (NR vs. R, at baseline) and 2 were longitudinal (R at Day 7 vs. R at baseline; and NR at Day 7 vs. NR at baseline). Analysis of functional groups of genes, as defined using Reactome (with Enrichr, https://maayanlab.cloud/Enrichr/), was used to identify enriched pathways among the genes of each of the 3 lists of DEGs.

BTM analysis

We used the 346 BTMs developed by Li et al. that were publicly available (Table S1). Each BTM contained a variable number of co-expressed genes and had received an identification number. A large number of BTMs were related to specific immune-cell subsets; for example, the BTMs ‘M.37.1: enriched in neutrophils (I)’; ‘M11.0: enriched in monocytes (II)’; ‘M7.0: enriched in T cells (I)’; ‘M7.2: enriched in natural killer [NK] cells (I)’; and ‘M47.0: enriched in B cells (I)’. Other modules were related to immune signalling; for example, ‘M16: TLR and inflammatory signalling’; M37.0: immune activation – generic cluster’ (Table S1). For each patient, the activity of each BTM module was taken as the mean value of its member genes. BTM activities were used for cross-sectional comparisons (NR vs. R, both at baseline) and longitudinal comparisons (R at Day 7 vs. R at Day 0; NR at Day 7 vs. NR at Day 0). Significantly correlated BTMs were identified using Morpheus (https://software.broadinstitute.org/morpheus/). To capture the highly correlated BTMs, unsupervised hierarchical clustering analysis was performed using the R package ‘linkcomm’.

PBMC RNA-seq deconvolution analysis

The ABIS software (https://github.com/giannimonaco/ABIS) was used to deconvolute our RNA-seq data. Indeed, ABIS through RNA-seq signatures normalised by mRNA abundance allows absolute deconvolution of human immune cell types present in PBMCs. CIBERSORTx was used for determining the correlation between the blood cell counts and those inferred from the RNA-seq analysis. Inferred CIBERSORTx lymphocyte counts were calculated as the sum of the counts of naive B cells, memory B cells, CD8 T cells, naive CD4 T cells, resting memory CD4 T cells, and activated memory CD4 T cells. The counts of monocytes, M0 macrophages, M1 macrophages, and M2 macrophages were summed to infer CIBERSORTx monocyte counts.

Flow cytometry

Frozen PBMCs obtained at baseline and after 7 days of corticosteroid therapy in 20 additional patients with SAH (10 R and 10 NR) were used for flow cytometry studies. Because results obtained with RNA-seq data indicated several differences in the lymphoid lineage between NR and R, flow cytometry focused on this lineage. Samples were washed twice with RPMI media and 1× PBS and further stained with fluorochrome-labelled monoclonal antibodies for B and T lymphocyte-specific subtypes. The lymphocytes population was gated using FSC and SSC plots. The populations were further subdivided and expressed as the percentage frequency of cells: total B cell population (CD19+), transitional B cells (CD24+highCD38+high), plasmablasts (CD20-CD38+), activated B cells (CD19+CD23+), B10 cells (CD19+CD23+IL-10+), B regulatory cells (CD19+CD23+CD25+), naive B cells (CD19+CD27-), memory B cells (CD19+CD27+), T helper (Th) cells (CD3+CD4+), and CD8 T cells (CD3+CD8+). All the cell populations were acquired on BD FACS (flow-assisted cytometry) verse platform and data were analysed using FlowJo software. The personnel handling the PBMCs were kept blinded to the source of the sample. Results obtained with flow cytometry were used for cross-sectional (NR vs. R at baseline) and longitudinal comparisons (R at Day 7 vs. R at baseline; and NR at Day 7 vs. NR at baseline).

Statistical analysis

R version 3.5.3 (2019) and RStudio (version 1.2.1335, 2019; R Foundation for Statistical Computing) were used for statistical analysis. An unpaired or paired Student t test was used, where appropriate, to estimate the significance between the study groups. False discovery rates, wherever mentioned, were estimated using the default method available in the applications used for the analysis. For measuring the distance between the BTMs for clustering, BTM activity scores were standardised for each BTM, by subtracting the BTM activity score mean value and dividing it by standard deviation. Regression was performed using Pearson correlation. To determine the accuracy of BTMs to distinguish the potential steroid response at baseline, the area under the receiver-operating-characteristic curve was assessed for all the 346 BTMs in comparison with the Lille score for estimating the deterministic accuracy of identifying the R and NR at baseline.

Results

Baseline clinical blood counts

First, we examined the results of routine clinical blood counts obtained in our patients. At baseline, the total blood leucocyte count and differential neutrophil count were abnormally increased in both eventual NR and R indicating leucocytosis and neutrophilia in the 2 groups. However, baseline leucocytosis was greater in eventual NR than R, whereas there was no between-group difference in the differential blood count of not only neutrophils (as expected), but also monocytes and lymphocytes (Table 1). Nevertheless, the greater total leucocyte count in NR relative to R may result from the cumulative effects of slightly greater absolute counts (frequency) of neutrophils, monocytes, and/or lymphocytes. Moreover, peripheral mononuclear cells include different subsets which are not assessed by routine techniques, and can differ between groups. Together, these findings indicate that other approaches are required to address the question of potential differences in peripheral mononuclear cells between NR and R.

Table 1

Total leucocyte count and differential blood count of leucocytes at baseline and after 7 days of corticosteroids in patients with severe alcoholic hepatitis who are responders (R) or non-responders (NR) to corticosteroids∗.

Variable	Baseline		After 7 days of corticosteroids		p value
	R (1)	NR (2)	R (3)	NR (4)	(1 vs. 2)	(3 vs. 4)	(1 vs. 3)	(2 vs. 4)
Total leucocyte count (×103 cells/mm3)	12.9 (4.0–31.9)	15.2 (7.9–33.0)	15.4 (4.3–19.4)	17.3 (5.0–34.9)	0.013	0.046	0.805	0.741
Differential blood count (%)
Neutrophils	78 (46–90)	81 (67–90)	82 (48–96)	83 (69–87)	0.713	0.62	0.976	0.42
Monocytes	8 (4–20)	7 (4–13)	8 (4–18)	8 (5–13)	0.67	0.84	0.132	0.089
Lymphocytes	13 (6–30)	13.5 (5–31)	14 (5–27)	15 (5–24)	0.92	0.71	0.657	0.703

Part of these results have been published in Sharma et al.

We considered that analysis of PBMC RNA-seq data may be a relevant approach to address this question. Indeed, RNA-seq (CIBERSORTx)–inferred white cell counts were also positively correlated with clinical laboratory measurements of complete blood counts (Fig. S5), which suggested that PBMC RNA-seq provided information that correlated with criterion-standard clinical measurements of blood counts.

Analysis of DEGs

To gain insights into eventual baseline differences in circulating mononuclear cells between NR and R, we first analysed PBMC RNA-seq data to determine features that may segregate R and NR at baseline. We first performed an unsupervised clustering analysis of the samples using the 18,650 expressed PBMC genes, which showed significant overexpression in NR relative to R (Fig. S1A). Enrichment analysis of the overexpressed gene set identifies multiple biological processes/molecular functions/cellular components, thus losing specificity (Table S2). As can be observed, the enrichment analysis does not highlight specific processes which can be used as markers for segregating SAH patients at baseline. Also, although there were DEGs between baseline and Day 7 in each group (11,773 genes for R [Fig. S1B] and 10,244 genes for NR [Fig. S1C]), their enrichment analysis could not categorically identify the immune profile associated with steroid R and NR (Tables S3 and S4). Moreover, association with an immune cell type could not be drawn from this analysis. Thus, we used the approach of analysing BTMs, which can reveal immunological features that are obscured in the classical enrichment analysis.

Searching between-group differences in immune cell subsets at baseline

Global increases in the baseline activities of blood BTMs in NR relative to R

As identified by Li et al.,, using RNA-seq data we were able to estimate the BTM activities for all of the 346 BTMs at baseline. We performed principal components analysis on BTM activities and found 2 distinct groups (R and NR) clearly distinguishable by the first principal component (PC1, 77.2% of the total variation; Fig. 1A). Unsupervised hierarchical cluster analysis showed that the baseline activity of the 346 BTMs differed between NR and R (Fig. 1B); the differences in BTM activities were significant (Table S5).

Fig. 1

Blood transcription modules (BTMs) are able to differentiate responders (R) and non-responders (NR) to steroid therapy at baseline.

(A) Principal component analysis of 346 BTMs using their activity scores showed that BTMs can segregate the NR and R, at baseline. The colour of the dots is only to identify the R or NR group. (B) Unsupervised clustering also confirmed the ability of the 346 BTMs to segregate NR and R, at baseline, with the exception of 1 sample (NR). The brown colour shows an increase and blue colour shows a decrease in gene expression in different samples (in columns). (C) Heatmap showing unsupervised clustering of the top 207 differentially expressed (>2-fold and p <0.05) BTMs between NR and R, which were able to segregate robustly NR and R, at baseline. The red colour shows an increase and green colour shows a decrease in gene expression in different samples (in columns).

BTMs related to lymphoid lineage have greater baseline activities in NR relative to R

Of the 346 BTMs, 207 had an upregulation of their activity by 2-fold or more in NR relative to R (Fig. 1C), with 38 modules directly related to mononuclear immune cells (Table S5, shaded grey). Of the 38 BTMs, 7 (18.4%) were related to myeloid cells and 31 (81.6%) to cells of the lymphoid lineage (B, T, and NK cells), a BTM related to NK cells (M61.1) showing the highest between-group difference (5.28 increase; Table S5). Among the top 100 most upregulated BTMs in NR relative to R, 11 were related to the lymphoid lineage and only 2 specifically to myeloid cells indicating that the upregulation of gene modules of lymphoid lineage was important in differentiating NR from R.

NR have evidence of increased T-cell populations at baseline

As the BTMs correspond to diverse cell type and cellular process or pathway type, we determined if the baseline BTM upregulation (increased activities by 2-fold or more, in NR relative to R), could identify the processes associated with specific blood cell subsets. The BTMs with highly similar activities across samples are likely to be coherent and perform together. Thus, an unsupervised hierarchical cluster analysis sorted the 100 upregulated BTMs into 15 clusters (Fig. 2). Among these, T cell-related BTMs including gene modules for ‘T cell surface signature’, ‘enriched in T cells (I)’, ‘CD4 T cell surface signature Th2-stimulated’, and ‘T cell differentiation (Th2)’ clustered together. BTMs related to cell proliferation including gene modules for ‘cell division in stimulated CD4 T cells’, clustered with those for cell cycle and cell division. Of note, clusters enriched in gene modules related to cell division and E2F transcription factor network were closely associated (Fig. 2). Interestingly, gene module of ‘CD4 T cell surface signature Th1-stimulated’ clustered with the module of ‘enriched in calcium signalling proteins’ (known to be involved in T cell activation). Together these findings suggest that activation and proliferation of CD4 T cells, and their differentiation into effector subtypes (here, with a mixed Th profile, Th1/Th2), were more marked in NR than in R.

Fig. 2

Correlation of cell type and function of blood transcription molecules (BTMs).

Dendrogram of unsupervised clustered top 100 BTMs to determine the correlations between the BTMs. The BTMs aggregated into 15 clusters, with BTMs belonging to B- and T-cell types clustering together. The colour of the edges in the dendrogram is to identify the clusters and do not refer to expression values. TBA, to be announced, i.e. those BTMs whose exact cellular/functional characterisation is pending; ZZ, break in dendrogram branches to make it fit in one page.

Next, using the ABIS software package, a means of deconvoluting RNA-seq data, we found independent evidence of significantly higher baseline CD4 T cell populations (including naïve and memory CD4 T cells; Fig. S2) in NR relative to R. In addition, the RNA-seq-inferred baseline naive CD8 T cells populations were found to be significantly higher at baseline in NR relative to R (Fig. S2). At baseline, among unconventional T cells, the RNA-seq-inferred MAIT (mucosal-associated invariant T cell) population was significantly higher in NR, while γδ T cells were not different between the 2 groups (Fig. S2). The findings related to CD4 and CD8 T cells obtained with RNA-seq analysis were strongly supported by results of our FACS analysis in isolated PBMCs from 20 new patients with SAH (Table S6) at baseline (Table S3; 10 of these patients were subsequently found to be R and 10 NR, showing a greater baseline frequency of both CD4 and CD8 T cells (Fig. 3).

Fig. 3

Baseline peripheral blood mononuclear cell flow cytometry confirms cell types identified by blood transcription modules.

Comparative estimates of cell frequency (%) for B, T, and NK cells types at baseline. Significant difference between transitional, activated and regulatory B cells, and CD4+ and CD8+ T-cells and NK cells were observed (Mann-Whitney U test). NR, non-responder; R, responder.

NR have evidence of increased B-cell populations at baseline

Fig. 2 shows that the largest cluster centred around B cell-related BTMs including gene modules for ‘enriched in B cells (VI)’, ‘B cell surface signature’, ‘memory B cell surface signature’, ‘naive B cell surface signature’, and module M217 ‘TBA (source: B cells)’ in association with BTMs related to cell proliferation including those relating to cell cycle, c-MYC translational network (a known amplifier of lymphocyte genes), E2F1 targets, amino acid metabolism, and transport. Moreover, cluster for gene modules related to naïve and memory B cells clustered closely. Together, these findings suggest a greater baseline activation and proliferation of B cells in NR relative to R.

Consistent with these results, the RNA-seq (ABIS)-inferred B-cell populations (including naïve B cells, memory B cells, and plasmablasts) were higher at baseline in NR relative to R (Fig. S2); a finding suggesting B-cell activation in NR. Our results of FACS analysis at baseline in the new patients’ cohort showing decreased frequency of transitional B cells and increased frequencies of activated B cells and regulatory B cells in NR relative to R were also consistent with an increased B-cell activation in NR. Of note, the frequency of naïve and memory B cells assessed by FACS did not differ between NR and R, a finding contrasting with the increases of these two B-cell subsets using RNA-seq data (see earlier and Fig. 3B and S2). A decrease in transitional B cells associated with increased activated B cells is seen during the early steps of B cell activation, whereas plasmablast production occurs at a later stage. Therefore, it is possible that patients enrolled for FACS and those enrolled for RNA-seq were not investigated exactly at the same stage of B-cell activation.

NR have evidence of increases in NK cells at baseline

Fig. 2 shows that a cluster was composed of a single gene module for ‘enriched in NK cells, KIR cluster’, suggesting that the increased activity of this module in NR reflected an increased frequency in NK cells. Consistent with this hypothesis, the RNA-seq-inferred NK populations at baseline were higher in NR than R (Fig. S2) and FACS analysis showed that baseline NK-cell frequency was significantly higher in NR relative to R (Fig. 3B).

NR have evidence of increased myeloid cells at baseline

In Fig. 2, BTMs related to myeloid cell formed a separate cluster, that is the gene module of ‘activated dendritic cells’ along with modules for ‘nucleotide metabolism’, ‘purine nucleotide biosynthesis’, ‘nuclear pore, transport; mRNA splicing, processing’, and ‘nuclear pore complex’, suggesting a close association of dendritic cell activation with increased nucleotide synthesis and utilisation in NR. Thus, activation of dendritic cells, is an anabolic process requiring nucleotides for nucleic acid synthesis, in particular mRNA encoding for a broad variety of cytokines, chemokines, among others., Consistent with these results, RNA-seq-inferred subpopulations of myeloid and plasmacytoid dendritic cells were higher in NR than R (Fig. S2). Of note, the RNA-seq-inferred monocyte population was increased in NR whereas the population of conventional and intermediate monocytes did not differ between the 2 groups. These findings suggest that, at baseline, a subpopulation of monocytes may be electively increased in NR vs. R, a hypothesis that should be addressed in future studies.

Increased baseline BTM activities related to the mitochondrial ETC in NR

Fig. 2 shows a cluster including BTM modules M231, M219, and M216, which are all gene modules of ‘respiratory electron transport chain (ETC) mitochondrion’. Because the mitochondrial ETC plays a crucial role in oxidative phosphorylation (OxPhos) and therefore energy production, our results suggest greater energy production in PBMCs from NR relative to PBMCs from R. Above-mentioned clusters were not associated with specific cell subsets, suggesting gene modules related to ETC may be expressed by different cell subsets.

Corticosteroid therapy associates with PBMC gene expression reprogramming in R but not in NR

Results of the routine laboratory values for leucocytes during corticosteroid therapy did not show differences between NR and R. To understand the variations in BTM activities arising with therapy, within the R and NR groups, we performed pairwise longitudinal comparison of the activity of each BTM at Day 7 of corticosteroid therapy relative to the corresponding value at baseline; the comparison being hereafter referred to as Day 7/0. In R, at Day 7/0, there were 64 BTMs with 2-fold change in their activities (3 upregulated, 61 downregulated). Among the upregulated BTMs at Day 7/0 in R, the gene module of ‘enriched in NK cells (KIR cluster)’ was the top one with a 22-fold increase in activity (Table S7A). An increase in NK cells was observed in R at Day 7 compared to baseline (Fig. 4C), although the increase was not as much as it was with the related BTM. Among downregulated BTMs at Day 7/0 in R, there were several gene modules related to innate immunity including those relating to mitogen-activated protein kinase signalling and AP-1 transcription factors, inflammasome, lysosome, monocytes, and antigen presentation. Interestingly, even though transcriptomics was measured on PBMCs, the gene module of ‘enriched in neutrophils (II)’ was also downregulated. There was also a downregulation of gene modules related to adaptive immunity including those relating to T cell signalling and co-stimulation, mitosis in CD4 T cells, B cell development, B cell-receptor signalling, cytoskeleton, and cell migration (Table S7A). Of note, there was also a downregulation of gene modules for cell migration. Together, these findings indicate that corticosteroid therapy is associated with extensive reprogramming of gene expression in innate and adaptive immune cells that compose PBMCs from R, giving a molecular basis to the clinical phenotype of ‘responder to steroids’.

Fig. 4

Flow cytometry-based assessment confirms variations in blood transcription module (BTM) identified cells associated with 7-days of steroid therapy.

(A) Comparison of B and T cell type frequencies between baseline and day 7 shows a significant reduction in naive, transitional, activated and regulatory B cells in responders. A similar trend is observed in non-responders, with the exception of activated B-cells, which were not significantly different post-therapy from baseline. (B) An assessment of T-cells showed an increase in CD4+ and CD8+ cells post-therapy in R but not in NR, suggesting defective cell types in NR, which could be the reason for patients not responding to therapy. (C) Baseline cell frequency of NK cells corroborates with significant increase in BTM 61.1 (NK cell related KIR cluster activity; Table S4) among NR. (Mann-Whitney U test.)

In contrast to R, NR did not exhibit significant changes in BTMs at Day 7/0 (Table S7B). It is noteworthy that an unsupervised cluster analysis of BTM activities measured at Day 7 revealed that these activities were globally greater in NR than R (Fig. 5A). Interestingly, principal components analysis on BTM activities at Day 7 still found 2 distinct groups (R and NR) clearly distinguishable by the first principal component (PC1, 88.2% of the total variation; Fig. 5B). Together, these findings indicate that during corticosteroid therapy, unlike PBMCs from R, PBMCs from NR did not exhibit gene expression reprogramming the activity of every BTM (reflecting their richness in expressed genes), giving a molecular basis to the clinical phenotype of ‘non-responders to steroid’.

Fig. 5

Variations between responders and non-responders in blood transcription module (BTM) scores after 7 days of steroid therapy.

(A) Heatmap for unsupervised clustering of BTM activity scores of NR and R, at day 7 post-steroid therapy. Note the persistent increase in BTM activity among NR and reduction among R post-therapy. The red colour shows an increase and green colour shows a decrease in gene expression in different samples (in columns). (B) Principal components analysis of BTMs post-therapy confirmed their ability to robustly segregate NR and R, at day 7 also. The colour of the dots is only to identify the R or NR group. NR, non-responder; R, responder.

Next, we sought to determine whether corticosteroid therapy was associated with changes in RNA-seq-inferred immune cell populations, by performing longitudinal comparisons (Day 7/0), in each group, R and NR. Surprisingly, the RNA-seq (ABIS)-inferred profile in cell-subtype changes was very similar in NR and R. Indeed, at Day 7/0, T-cell subpopulations were not significantly affected in the 2 groups, whereas the populations of monocytes and naive B cells were significantly decreased in NR and R (Fig. S3 and S4, respectively). The results of FACS analysis in PBMCs obtained from the patients of the additional cohort after 7 days of corticosteroids, enabled to refine changes in B cells that occurred in each group between Day 7 and baseline. Again, changes were similar in the 2 groups and consisted in increases in the transitional B cell population and decreases in the frequency of naïve B cells, activated B cells, and regulatory B cells (Fig. 4A). Because of the similar pattern of changes in several B cell subset frequencies in R and NR, corticosteroids did not seem to exert a control at the level cell subset frequency (i.e. homeostasis or differentiation). However, the frequencies of CD4 and CD8 T cell increased significantly in R but not NR (Fig. 4B), suggesting that corticosteroid-related mechanisms may positively regulate the T cell frequency in R. The fact that gene modules related to T cells were downregulated in R (see earlier), suggests that in T cells, corticosteroids may control gene expression level and cell subset frequency through different mechanisms. Correlations between the blood cell counts by FACS and the RNA-seq data (using CIBERSORTx) demonstrated positive correlations (Fig. S5). The correlation coefficients were moderate, probably owing to the low sample size.

Discriminating accuracy of baseline BTM activities in predicting the response to corticosteroid therapy

Distinction provided by BTM activities in segregating R and NR at baseline prompted us to investigate the discriminating accuracy of the top 100 BTMs (at baseline) with that of the Lille score (7 days post-therapy). The Lille score of course showed an AUC of 1, as it is the parameter on which the patients were segregated. All the BTMs were considered for this analysis. All the 100 BTMs had a significant AUC for non-response, of which most were of more than 0.8 (Table S8). The best performers had an AUC of 1 and included gene modules for ‘cell division – E2F transcription network’, ‘Rho GTPase cycle’, ‘cell division’, ‘chemokine cluster II’, ‘enriched in B cells VI’, ‘TBA (M41.2)’, ‘enriched in T cells II’, ‘TBA (M41.0)’, and ‘lipid metabolism, endoplasmic reticulum’, with an AUC of 1. It is important to note that the baseline clinical parameters such as the CTP score, MDF, and MELD score were not significantly different in R and NR.

Discussion

This is the first study using results of PBMC RNA-seq to investigate whether the composition of PBMCs may differ between NR and R, at baseline and during corticosteroid therapy. The use of 346 BTMs and their activities (according to Li et al.),, distinct from the conventional use of classical analysis of DEGs, allowed us to clearly distinguish NR from R, at baseline and after 7 days of corticosteroids. Although it was surprising to find a large number of BTMs to be overexpressed in the NR compared with R, it should be noted that the BTMs constitute a certain set of genes. Given that there were large numbers of genes that were overexpressed in NR at baseline (Fig. S1), it is not surprising that the BTMs derived from such gene sets would also reflect similar patterns. Also, these BTMs achieved a significance of p <0.05.

The most important results obtained at baseline were as follows: First, analysis of BTM activities was consistent with greater activation/proliferation and differentiation of T cells and B cells, in NR relative to R; these findings being supported by 2 independent approaches, including deconvolution of RNA-seq data using the ABIS software and FACS analysis in a new cohort of patients with SAH. Second, NR had greater activity in a gene module for NK cells that FACS results suggested to be related to an increased frequency of NK cells. Third, NR had greater activities of gene modules related to ETC and therefore OxPhos and energy production.

Activation of immune cells consists in the stimulation of anabolic pathways (see earlier) which are energetically expensive.[20], [21], [22] However, the mechanisms for energy production depend on the subset of immune cells. In activated myeloid cells, the production of energy relies on aerobic glycolysis (Warburg effect) whereas in T and B cells, activation relies on aerobic glycolysis but also OxPhos, which is crucial., Therefore, our BTM and FACS findings of greater baseline activation/proliferation in T and B cells from NR relative to R likely explains the upregulation of BTMs related to ETC in PBMCs from NR.

After 7 days of corticosteroids, R had downregulation of several BTMs related to both innate and adaptive immunity, indicating that a profound gene expression reprogramming forms the molecular basis of the appropriate PBMC response to steroids. In contrast, NR did not exhibit any features of gene expression reprogramming because they had no significant changes in their BTM activities relative to baseline. In other words, it is likely that lymphocyte activation/proliferation observed at baseline in NR, persisted in these patients after 7 days of corticosteroids.

Future studies should address several questions which might be inter-related. What are the molecular mechanisms for the baseline lymphoid activation in NR? Are there commonalities between these baseline molecular mechanisms and those which are involved in the ‘PBMC resistance’ to corticosteroid therapy? Does the lymphocyte activation, which exists at baseline and persists despite corticosteroid administration, play a role in the poor outcome of NR? Because lymphocyte activation is energetically expensive, its persistence over time in NR may have a metabolic cost resulting in reallocation of nutrients to fuel immune cell metabolism at the expense of the homeostasis of several non-haematopoietic tissues and leading to their dysfunction. Addressing these questions is important because the answer might produce clues to novel therapeutic approaches, for example which might interfere with immune cell metabolism.

The robustness of BTMs in identifying the response to therapy is evident from their ability to accurately discriminate the patients as effectively as the current standard Lille score, but only with an advantage of doing so at baseline itself. It is interesting to note that out of the significantly associated 100 BTMs, the majority demonstrated a discriminating accuracy of more than 0.8, which was not observed with other clinical parameters such as MELD, CTP, or MDF. In our efforts to determine a marker for early detection of response to steroids, the BTMs have proven to be most effective when compared with hepatic or metabolic markers. Our study highlights the importance of BTMs in determining the immunological phenotypes and serves as a starting point for further investigations using BTMs in liver diseases. It also helps to identify corticosteroid NR before initiation of the therapy.

Clinically actionable biomarker(s) predicting the response to corticosteroids in patients with SAH are still an unmet medical need. Our results suggest that either a gene signature composed of a limited number of genes which are easy to measure in clinical laboratories or lymphocyte subpopulation(s) identified by FACS, could serve as clinically actionable biomarkers of the response to corticosteroids in these patients. Potential biomarkers should be investigated in future studies.

Financial support

The current study was supported by Grant no. 4903 of the .

Authors’ contributions

Performed the experiments, generated, and analysed the data: S.S., J.S.M., S.B., S.M.S. Prepared the manuscript: S.S., J.S.M., S.B., R.M., S.K.S. Conceptualised and supervised the study, and finalised the manuscript: R.M., S.K.S.

Data availability statement

Data are available upon request to the corresponding author.

Conflicts of interest

The authors disclose no conflicts.

Please refer to the accompanying ICMJE disclosure forms for further details.