Stress responses, M2 macrophages, and a distinct microbial signature in fatal intestinal acute graft-versus-host disease

Patients with severe gastrointestinal (GI) acute graft-versus-host disease (aGVHD) who do not respond to steroids have a dismal prognosis. For decades, such patients have been treated with escalated immunosuppression based upon the clinical assumption that the donor T cells responsible for aGVHD are persistently damaging tissues despite high-dose steroids. Unfortunately, few patients with steroid-refractory aGVHD respond to escalated T cell–targeted immunosuppression and most ultimately succumb to unresolved organ damage or infection (1, 2).

Little data derived from human GVHD target tissue samples exist to guide novel therapies in steroid-refractory aGVHD. Detailed histologic evaluation of GI aGVHD has demonstrated crypt loss, mucosal atrophy, and mucosal hemorrhage as potential hallmarks of severe aGVHD (3–6). Inflammatory findings are often nonspecific, providing little information as to the functional changes occurring in the tissue (6). This paucity of target tissue biological data may have contributed to the relative lack of progress in therapy for steroid-refractory aGVHD for many years.

To uncover novel molecular insights driving fatal human steroid-refractory aGVHD, we performed RNA-Seq on archived FFPE rectosigmoid biopsies obtained during the routine clinical care of patients with a diagnosis of lower GI stage 3–4 aGVHD at onset (n = 22) compared with their repeat biopsies obtained for evaluation of steroid-refractory aGVHD. In doing so, we could identify transcriptome changes in the severely affected intestinal mucosa when the patients’ clinical status deteriorated because of apparent inadequate steroid response. We also compared aGVHD samples with 10 normal rectosigmoid biopsies to gain additional insight into the pathophysiology of life-threatening aGVHD itself, validating our findings using NanoString. Both microbial and human RNA are sequenced with RNA-Seq, and simultaneous analysis of host–microbe interactions is increasingly feasible albeit largely limited to in vitro experiments to date (7, 8). We analyzed microbial gene expression in our intestinal mucosa biopsies and noted a distinct profile in rapidly fatal aGVHD.

Patient descriptions and outcomes. Patient demographics are detailed in Table 1. This was an adult population undergoing hematopoietic cell transplantation (HCT) for hematologic malignancy. The majority of patients were male (73%) and underwent reduced intensity conditioning (68%). Peripheral blood stem cells were the graft source for 50% of the recipients. The second rectosigmoid biopsy performed for evaluation of steroid-refractory aGVHD occurred at a median of 19 days (IQR 12–49 days) after the first biopsy. Seventeen patients received antithymocyte globulin as second-line therapy in the cohort; the remaining five patients received etanercept as their second-line therapy. Eighteen of twenty-two patients (82%) in this steroid-refractory GI aGVHD cohort died at a median of 82.5 days (IQR 50–160 days) from the first rectosigmoid biopsy, with the majority of patients dying from GVHD (59%).

Table 1

Patient demographics

Differential expression of single genes. We show the overall structure of the data set as demonstrated by principal component analysis of RNA-Seq and the nCounter Human Immunology Panel (NanoString) in Figure 1, A and B, respectively. Overall, the normal samples tended to cluster, with no clear pattern in the data structure when comparing onset versus steroid-refractory GI GVHD biopsies.

Figure 1

Principal component analysis. In both the RNA-Seq dataset (A) and NanoString dataset (B), there was no clear separation of sample types.

Several immune- and damage-relevant genes differed between aGVHD onset and normal samples in single-gene analyses detailed in Table 2. The most highly differentially expressed gene in aGVHD compared with normal was chitinase 3-like 1 (CHI3L1; fold change 3.4; adjusted P = 0.003), a recently described antigenic target in Crohn’s disease (9). The most significantly decreased gene in aGVHD compared with normal was aquaporin-8 (AQP8; fold change –5.0; adjusted P < 0.001), a water channel with decreased expression in inflammatory bowel disease (10). Comparing onset aGVHD versus steroid-refractory aGVHD in paired analyses, several RNA-related and metallothionein (MT) genes were increased, while chemokine (C-C motif) ligand 18 (CCL18) and intelectin-1 (ITLN1) were decreased (Table 3). Descriptions of the functions of the top 25 up- and downregulated genes with annotations summarized from RefSeq (11) are detailed in Table 2 and Table 3.

Table 2

Top 50 differentially expressed genes along with their principal protein functions from RefSeq (11) in onset GVHD versus normal colon

Table 3

Top 50 differentially expressed genes along with their principal protein functions from RefSeq (11) between onset and steroid-refractory GVHD

In supervised analyses of immune checkpoint genes that can inhibit aGVHD, we found a significant increase in indoleamine 2,3-dioxygenase (IDO1) and programmed cell death ligand 1 (PDL1) at the onset of aGVHD that was subsequently lost in the steroid-refractory setting (Figure 2), both validated with NanoString. We attempted to but could not identify PDL1 expression at the protein level in either aGVHD or normal GI tissue using the antibody in clinical use for solid tumors (Supplemental Figure 1; supplemental material available online with this article; https://doi.org/10.1172/jci.insight.129762DS1). We found a similar pattern of increased gene expression of TIGIT at the onset of aGVHD that significantly decreased in refractory aGVHD (only accessible for analysis using NanoString) (P < 0.05).

Figure 2

Hypothesis-driven single gene analyses involving immunologic tolerance and gene expression stimulated by microbes. RNA-Seq (A) and NanoString (B) single-gene analyses of programmed death ligand 1 (PDL1), aryl hydrocarbon receptor (AhR), indoleamine 2,3-dioxygenase (IDO1), amphiregulin (AREG), and T cell immunoreceptor with Ig and immunoreceptor tyrosine-based inhibition motif (ITIM) domains (TIGIT). Box plots show differences between groups (onset n = 22, refractory n = 22, and normal n = 10), with medians compared by Wilcoxon’s rank-sum test. The box plots depict the minimum and maximum values (whiskers), the upper and lower quartiles, and the median. The length of the box represents the IQR.

As for potential interaction of host mucosa with microbiota, we found a significant increase in AhR gene expression at aGVHD onset that was also elevated in the steroid-refractory setting, validated by NanoString (Figure 2). Similarly, we found that AREG expression was significantly higher in onset and steroid-refractory aGVHD compared with normal rectosigmoid biopsies, although there was no difference between the GVHD groups (Figure 2, AREG only accessible for analysis by RNA-Seq).

Gene set enrichment analysis. We compared gene sets using the hallmark, curated, and immunologic signature gene sets from the MSigDB collections to determine whether previously described gene sets could distinguish between the groups of biopsies. Figure 3A shows a heatmap of top differentially expressed genes within these data sets. Mitochondrial gene HMGCS2 was the most significantly downregulated whereas SETD7 was the most highly upregulated in aGVHD onset compared with normal samples. The top gene set increased in aGVHD was GSE8835, indicating changes in gene expression of T cells (Figure 3B), consistent with known aGVHD pathophysiology. However, gene set enrichment analysis (GSEA) showed no pathways that were statistically significantly different between onset versus steroid-refractory GI aGVHD, with the heatmap of the top differentially expressed genes showing less distinct boundaries (Figure 4C).

Figure 3

Comparison of acute graft-versus-host disease (aGVHD) versus normal control gene expression profiles. (A) Heatmap of differentially expressed genes in onset acute graft-versus-host disease (aGVHD) versus controls. (B) Enrichment plot of top expressed pathways in aGVHD versus controls.

Figure 4

Comparison of acute graft-versus-host disease (aGVHD) onset versus steroid-refractory aGVHD gene expression profiles. GO biological processes as determined by top 25 increased (A) and top 25 decreased (B) genes in onset versus steroid-refractory aGVHD. (C) Heatmap of differentially expressed genes in onset versus steroid-refractory aGVHD.

Enrichr analysis. The Enrichr analysis of the top 25 increased genes in aGVHD compared with normal is located at http://amp.pharm.mssm.edu/Enrichr/enrich?dataset=48mwo, and analysis of the top 25 decreased genes is located at http://amp.pharm.mssm.edu/Enrichr/enrich?dataset=48mwq The most highly enriched Gene Ontology (GO) biological process in aGVHD samples was inflammatory response (GO:0006954, P < 0.001), consistent with the known pathophysiology of aGVHD. The most enriched biological process based upon decreased gene expression in aGVHD was positive regulation of respiratory burst (GO:0060267, P < 0.001).

Biological processes differing between onset and steroid-refractory aGVHD based upon the top 25 increased and decreased genes are located at http://amp.pharm.mssm.edu/Enrichr/enrich?dataset=48mwc and http://amp.pharm.mssm.edu/Enrichr/enrich?dataset=48mwk, respectively, and summarized in Figure 4. We observed that multiple biological processes involving MTs were significantly increased, whereas the only process that was significantly decreased in steroid-refractory GI aGVHD was regulation of T cell apoptosis (GO:0070234), related to decreased expression of IDO1 and PDL1 in the steroid-refractory state (Table 4 and Table 5).

Table 4

Top GO biological processes differentiating onset aGVHD versus steroid-refractory aGVHD

Table 5

Top GO biological processes differentiating onset aGVHD versus steroid-refractory aGVHD

Cell deconvolution analysis. Using CIBERSORT to estimate cell composition, we observed more CD4⁺ activated memory T cells and M0 macrophages in onset GI aGVHD compared with normal rectosigmoid biopsies (Figure 5), consistent with known aGVHD pathophysiology. Both onset and steroid-refractory aGVHD had increased M1 macrophages compared with normal. However, M2 macrophages were the highest in steroid-refractory GI aGVHD compared with both aGVHD onset and normal biopsies. Plasma cells were also enriched in steroid-refractory aGVHD compared with aGVHD onset. We observed no differences in activated T cell subsets between onset and steroid-refractory aGVHD. Resting CD4⁺ memory T cells and memory B cells were lower in onset and steroid-refractory aGVHD compared with controls. Resting mast cells were lower in steroid-refractory aGVHD compared with normal. Other comparisons did not reach statistical significance.

Figure 5

Immune cell deconvolution results using CIBERSORT. Contour plots show differences between groups (onset n = 22, refractory n = 22, and normal n = 10) with medians compared by Wilcoxon’s rank-sum test.

Association of gene expression with early mortality. Survival time after diagnosis of aGVHD is the most important clinical outcome. We sought to identify a gene expression pattern suggestive of early mortality by comparing differential gene expression in patients’ survival above and below the median of 82.5 days. We found those patients with an early death had higher expression of 4 genes: glycerol kinase (GK); phosphoglycerate kinase 1 pseudogene 2 (PGK1P2), histone cluster 1 H1 family member A (HIST1H1A), and cytochrome C oxidase subunit 7B pseudogene 2 (COX7BP2; Figure 6). Based upon cell signaling analysis by Enrichr, these genes indicate a greater degree of DNA damage and stress in patients with poor survival (Supplemental Figure 2). Patients with survival above the median generally had higher expression of genes (Figure 6) associated with activation of Smoothened (Supplemental Figure 3; P = 0.005; Q = 0.04), a Hedgehog pathway gene shown to constrain damage due to colitis in a preclinical model (12).

Figure 6

Unsupervised hierarchical cluster of top 10 genes differentiating early (<82.5 days) versus later (>82.5 days) death at aGVHD onset. The top 10 genes comprise the rows, while individual patients and their survival time in days is indicated in the columns.

We also explored microbial sequences using PathSeq (13) to determine whether there was a microbial pattern at aGVHD onset that was associated with early mortality. This approach subtracts potential microbial genes from the human sequences for estimation of candidate microbes within whole-genome sequencing data. We identified the top 10 microbes with the greatest differences between the 2 groups based on survival time. In unsupervised hierarchical clustering of these top 10 candidate microbes, we observed that the majority of patients with poor survival had low presence of candidate microbes compared with those with longer survival (Figure 7). This suggests a loss of mucosal microbes in general, as opposed to a pattern of differing microbial species, may be relevant to survival in aGVHD.

Figure 7

Unsupervised hierarchical cluster of predicted microbes differentiating early (<82.5 days) versus later (>82.5 days) death at aGVHD onset. The top 10 microbes comprise the rows, while individual patients and their survival time in days is indicated in the columns.

Multiple T cell–independent mechanisms contribute to the pathophysiology of steroid-refractory aGVHD. Our results suggest that steroid-refractory aGVHD may be characterized by loss of host tissue tolerance signals, increased stress-induced MT expression, and accumulation of M2 macrophages and plasma cells. Early death was associated with higher expression of genes related to DNA damage at the onset of aGVHD, whereas prolonged survival may be associated with modulation of the Hedgehog pathway via Smoothened. These potentially novel findings broaden the horizon for clinical investigation into a fatal disease. Notably, we did not identify significant differences in T cell content in onset versus steroid-refractory aGVHD. Our data are consistent with a recently described mouse model of steroid-refractory aGVHD that did not show differences in donor T cell–driven inflammation (14). Despite the lower quality input RNA from FFPE, our RNA-Seq results could be validated by NanoString, a platform well suited for robust gene expression analysis in FFPE clinical samples (15).

We made several potentially novel observations comparing aGVHD onset with normal biopsies in this analysis. For example, we observed a significant elevation of CHI3L1, also known as YKL-40, a proinflammatory factor with known pathogenic roles in neovascularization, macrophage recruitment, and bacterial adhesion, in aGVHD (16–19). The potential role of CHI3L1 in clinical allogeneic HCT has otherwise not been described aside from a report identifying elevated plasma CHI3L1 in transplant recipients with a very high HCT comorbidity index of 5+ (20). Increased CHI3L1 may suggest a longer duration of inflammation in aGVHD. In mice, CHI3L1 expression is increased the most during the chronic phase of colitis, in which it functions to promote intestinal epithelial cell proliferation and survival (21). On the other hand, CHI3L1 may worsen outcomes of aGVHD by facilitating bacterial adhesion to colonic epithelial cells and possibly even serving as a neutrophilic antigenic target, as recently described in Crohn’s disease (9, 22). Fecal CHI3L1 can be quantified, and its concentration is closely associated with mucosal damage in inflammatory bowel disease (23). Recently, Li et al. described that knockout of CHI3L1 in donor T cells increased aGVHD lethality in a murine model (24). Our results may not be comparable because our study was not limited to isolating T cell effects. Nonetheless, future GVHD models should further delineate the role of CHI3L1 expression in host intestinal epithelial cells, and its potential role as an alloantigen in aGVHD should be explored.

Our analysis also identified 5-fold decreased gene expression of AQP8 at the onset of aGVHD. AQP8 is a water channel important for hydrogen peroxide transport and mitochondrial function known to be decreased in inflammatory bowel disease (IBD), suggesting additional molecular parallels between the inflammatory conditions of aGVHD and IBD (10, 25–27). In our GSEA, an additional mitochondrial gene, HMGCS2, was significantly lower in aGVHD compared with controls. In rats, expression of HCGCS2 is dependent upon microbiota-derived butyrate (28). Although this relationship is not yet established in humans, the role of loss of butyrate as the major microbiota-derived energy source for enterocytes in aGVHD pathogenesis is of intense interest (29). These findings regarding AQP8 and HMGCS2, along with our findings of lower microbial RNA in patients with early mortality, suggest that fatal aGVHD may be associated with mitochondriopathy related at least in part to dysbiosis. Indeed, a critical role of butyrate improving mitochondrial respiration and preventing autophagy in colonocytes has been identified (30). A recent large-scale RNA-Seq analysis of over 600 patients with ulcerative colitis also indicated reduced mitochondrial energy function along with shifts in microbial community structures (31). Further functional studies will be necessary to confirm this potential microbiota–mitochondrial crosstalk in aGVHD.

A major strength of this study is paired analyses in patients at diagnosis of aGVHD compared with when they became steroid refractory. We observed that processes involved in MT-associated stress responses are increased in steroid-refractory aGVHD compared with aGVHD onset. MTs are widely expressed stress response proteins that bind heavy divalent metals and mitigate oxidative stress (32). MT-positive cells can be observed infiltrating the lamina propria in patients with inflammatory and infectious colitis, with the number of MT-positive cells correlated with the severity of colitis (33). Anti-MT antibodies have been developed for use in mouse models of IBD and represent a potential novel therapeutic target for patients (34). It cannot be concluded from our analysis that an increase in MT gene expression is solely linked to stress because glucocorticoid exposure can also induce MT gene expression, thus possibly indicating a treatment effect in these biopsies (35).

Beyond inflammation, steroid-refractory versus onset aGVHD appears to be associated with a tissue repair response with increased AhR, AREG, and SETD7 expression. AhR activation is protective in colitis by regulating intestinal crypt stem cell differentiation (36). AhR signaling can also drive immune responses in favor of either regulatory T cell or Th17 differentiation, suggesting ligand- and cell type–specific effects that could either ameliorate or worsen aGVHD (37, 38). AREG is similarly protective in models of colitis and associated with Treg- and T9-mediated tissue repair (39–43). AREG protein expression in the GI tract during clinical aGVHD can be variable, and the prognosis is largely driven by whether AREG is found in the circulation (unfavorable) or not (favorable) (44, 45). SETD7, the most highly differentially expressed gene comparing onset aGVHD with normal biopsies in our GSEA, is required for Wnt-dependent intestinal regeneration in chemical- or irradiation-induced damage (46). Overall, an increase in expression of AhR, AREG, and SETD7 likely represents molecular processes supporting mucosal restitution after damage, although an increase in these genes did not lead to improved clinical outcomes in these patients.

One possible explanation for a failure of mucosal healing could be decreased expression of immune checkpoint genes, including IDO1, PDL1, and TIGIT. However, it is also possible that this observation is due to low interferon-γ related to low T cell content, but this cannot be conclusively determined through RNA-Seq. Previous studies have shown an increase in IDO1 expression in the setting of aGVHD (47, 48). Our results are consistent with these reports and extend these findings by demonstrating that IDO1 expression decreases when patients transition to the steroid-refractory state. Similarly, PDL1 expression can be increased in colonic epithelium in the setting of inflammation, and its loss on host tissues is associated with poorer outcomes due to increased gut T cell homing, impaired T cell apoptosis, and impaired IL-22 production from CD11c⁺CD11b⁺ lamina propria cells (49, 50). TIGIT is involved in regulating experimental aGVHD, with TIGIT-Fc attenuating aGVHD severity and T cell TIGIT deficiency accelerating it (51, 52). Identifying therapies that can sustain expression of inhibitory checkpoints may prove protective in life-threatening intestinal aGVHD.

Our analysis found that M2 macrophages but not T cells accumulated during steroid-refractory aGVHD. M2 macrophages have also been shown to accumulate with chronicity in IBD and impair enterocyte differentiation through activation of Wnt signaling, suggesting their presence may indicate incomplete tissue repair (53). Future functional and single-cell studies of the mucosal immune cells will be required to clarify the role of M2 macrophages in steroid-refractory aGVHD.

Our study is limited in that our findings cannot be validated using immunohistochemistry because the samples were consumed during RNA extraction. Additionally, our analysis was on heterogeneous populations of cells, and subtle shifts within smaller immune cell populations may not be detectable. A large, independent data set with clinical samples adequate for both single-cell sequencing and histology will be necessary to validate these findings. Such a study may not be possible within single institutions because of small sample size and most patients not undergoing repeat diagnostic endoscopy in the setting of steroid-refractory aGVHD. Nonetheless, focused efforts designed to deconstruct the molecular mechanisms underlying steroid refractoriness in aGVHD should continue.

In summary, our analysis of the human mucosal aGVHD transcriptome expands the paradigm of T cell–centric therapies for steroid-refractory GI aGVHD. We have highlighted multiple potential mechanisms for scientific investigation and therapeutic exploration for this fatal condition. Ultimately, multifaceted and personalized approaches that support tolerance, enhance tissue repair, restore commensal microbiota, and improve mitochondrial function will be needed to reduce the substantial risk of mortality due to severe intestinal aGVHD.

All patients in this data set underwent allogeneic HCT at the University of Minnesota between 2008 and 2016. We performed RNA-Seq on FFPE rectosigmoid biopsies from 22 adult patients with clinical stage 3–4 GI aGVHD at the onset of GI aGVHD, repeated at the diagnosis of steroid-refractory GI aGVHD. Rectosigmoid biopsies have 70% sensitivity, 87% specificity, and 97% positive predictive value for the diagnosis of aGVHD (54), therefore representing an excellent, consistent site of biopsy for RNA-Seq analysis. All 22 patients had histologic evidence of GI aGVHD and did not have confounding infectious diagnoses (e.g., cytomegalovirus or concurrent C. difficile colitis). These patients were the only adult HCT recipients at our institution with stage 3–4 severity at onset who underwent repeat biopsy and had FFPE blocks available for analysis at the time of the study. We compared these 22 paired GI aGVHD biopsies with FFPE rectosigmoid biopsies obtained from 10 healthy controls who had histologic analyses previously published (44).

We extracted RNA from the FFPE biopsies using PureLink FFPE RNA Isolation Kit (Thermo Fisher Scientific). The entire biopsy was consumed in the process of RNA extraction. FFPE samples yielded an adequate quantity of RNA for sequencing (median 11,707 ng, range 290–33,220 ng), although as expected RNA integrity numbers (RINs) were low (median RIN 2.2, range 1.1–3.3). After extraction, we prepared the libraries for analysis with the SMARTer Stranded Total RNA Kit (Takara Bio) and performed sequencing using the Illumina HiSeq 2500 platform, high-output mode, paired-end reads, 50 cycles. For validation of key immunologic findings, we analyzed residual RNA using the nCounter Human Immunology Panel (NanoString).

We aligned and assembled sequences using TopHat and Cufflinks, respectively; performed GSEA using MySigDB gene sets in GSEA v3.0; analyzed pathways using Enrichr; and estimated immune cell composition using CIBERSORT. For estimation of microbial RNA content, we analyzed reads that were not mapped to the human genome. Specifically, we used hisat–un-conc-gz to output unmapped pair ends as FASTQ files; then we aligned unmapped pair ends using the hisat2 algorithm. We then used PathSeq as a mapping and taxonomic classification algorithm to estimate abundance of candidate microbes using read counts assigned to each organism (13). This approach has previously been used to assess mucosal dysbiosis in Crohn’s disease (55). The raw data for both RNA-Seq and PathSeq are available at the Gene Expression Omnibus (accession GSE134662). To assess PDL1 staining in colon biopsies, we performed staining in a single test case (rectosigmoid biopsy of aGVHD from a patient not in this cohort) and control case (normal colon biopsy from a patient not in this cohort) using Food and Drug Administration–approved antibody Ventana SP263 (Roche Diagnostics) in the University of Minnesota Medical Center’s clinical pathology laboratory according to the manufacturer’s instructions.

Statistics. We assessed the overall structure of the data using principal component analysis. We compared gene expression across groups, adjusting for multiple comparisons with Bonferroni’s correction within the Limma package. Our single-gene expression data set was large; many dependent variables could have led to significantly low P values by chance (47). Limma provided Bonferroni’s adjusted P values in its contrast matrix, accounting for sample size and prior test statistics used in the appropriate design matrix. In hypothesis-driven analyses, we made comparisons between groups using Wilcoxon’s rank-sum test without correction for multiple comparisons. All P values were 2 sided.

Study approval. This study of previously archived biological specimens and data was approved by the University of Minnesota Institutional Review Board. We received a waiver of consent for accessing archived samples and data, as the majority of patients were deceased.