Transcriptome Analysis of GVHD Reveals Aurora Kinase A as a Targetable Pathway for Disease Prevention

Sirolimus and Tac-Mtx provide graded disease protection in a NHP model of GVHD

Figure 1a shows our previously described strategy(5) for MHC-mismatched allogeneic-HCT in NHP, in addition to the GVHD prophylaxis strategies employed in the current study. As shown in the figure, this study compared untransplanted healthy control animals (‘HC’, n = 34) as well as autologously transplanted controls (‘Auto’, n = 6 for flow cytometric data and n = 3 for gene expression profiling) to MHC-mismatched allogeneic transplant recipients receiving no immunoprophylaxis (‘No Rx’, n = 4), those receiving sirolimus(20, 21) (‘Siro’, n = 4) those receiving tacrolimus-methotrexate combination immunoprophylaxis(4) (‘Tac-Mtx’, n = 3). All transplants utilized myeloablative TBI as a pre-transplant preparative regimen, transplantation of T cell-replete peripheral blood stem cell products obtained after G-CSF-mobilized leukapheresis, and supportive care with anti-bacterial, anti-fungal and anti-viral prophylaxis. Figure 1b demonstrates longitudinal analysis of clinical GVHD scores for the allo-transplant recipients, using our previously described scoring system(5). Briefly, the GVHD clinical score increases with cumulative GI-specific abnormalities (e.g. diarrhea), liver-specific abnormalities (e.g. hyperbilirubinemia) and skin-specific abnormalities (extent and character of rash). These scores demonstrate the rapid and severe clinical syndrome that accompanied unprophylaxed GVHD, the partial protection and delay in onset of severe GVHD in the recipients treated with sirolimus monotherapy and the more substantial protection from clinical disease observed in the dual therapy tac-mtx cohort. As expected, significantly shorter overall survival paralleled an increase in clinical GVHD scores with median survival times (MST) of 7.5, 17 and 49 days observed in the No Rx, Siro and Tac-Mtx cohorts, respectively, with all autologous transplant recipients surviving long-term (Figure 1c). We found a significant difference between the MST of the Tac-Mtx cohort relative to the No Rx cohort (P = 0.018, log rank test) and a difference approaching significance between the Siro and No Rx cohort (P = 0.055, log rank test). Donor engraftment (Figure 1d) occurred in all recipients, and hematopoietic reconstitution was notable for a rapid expansion of CD4+ and CD8+ T cells in the unprophylaxed recipients, that occurred faster than in either autologous controls or prophylaxed cohorts (Figure 1e).

Unsupervised analysis identifies treatment cohort as the major contributor of transcriptional variation during GVHD

To facilitate a comprehensive, unsupervised analysis of the molecules and pathways that are aberrantly regulated during GVHD and how the expression of these genes changes with immunoprophylaxis, we profiled messenger RNA (mRNA) expression from sorted CD3+ T cells isolated from healthy controls (with no history of HCT), and on day +14 or at the time of sacrifice (if sacrifice occurred prior to day +14) from transplant recipients. Principal Component Analysis (PCA), a dimensionality reducing method that is used to identify the determinants of variation in complex data(22), was then performed. This allowed us to determine the degree to which the T cell transcriptome of the three allo-HCT cohorts differed from both the autologous transplant recipients and from healthy controls, and from each other, and the degree to which transcriptional variation was shared between cohorts. As shown in Figure 2a, projections of the first and second principal components showed discrete clustering when grouped by treatment cohort. Confidence ellipses (90%) of PCA projections showed no overlap between the No Rx cohort and either healthy controls or autologous recipients and independent tests of statistical significance confirmed that the first principal component projections were different between these groups (P < 0.001 – No Rx vs HC, P < 0.001 – No Rx vs Auto; unpaired t-test corrected for multiple hypothesis testing using BH procedure). In contrast, projections from the Tac-Mtx cohort exhibited substantially more overlap with those of autologous transplant recipients and healthy controls and the first principal component of the Tac-Mtx group was not significantly different than that of the autologous controls (P > 0.05). This analysis provided strong evidence that the treatment groups could be distinguished transcriptionally from each other, and that the No Rx group could be distinguished from both Autos and HCs, and as such, provided the foundation for further analysis of the gene expression pathways and networks that define each group.

To develop a more detailed understanding of the transcriptional landscape among the treatment cohorts on a single gene level, we performed differential expression (DE) analysis. Confirming our PCA results, we found an inverse correlation between the number of DE transcripts (up- and down-regulated) compared to either HC or autologous transplant controls and the length of GVHD-free survival after allogeneic HCT (Figure 2b). Showing the complimentary comparison using a spatially-weighted, Venn diagram (Chow-Ruskey plot), Figure 2c depicts healthy controls, autologous controls and treated allogeneic HCT cohorts relative to the No Rx cohort. This analysis underscores the high degree of transcriptomic variation between unprophylaxed GVHD and all other groups, as well as the significant sharing of differentially expressed transcripts in the autologous and healthy controls (Figure 2d) when compared to the No Rx allo-HCT transcriptome. The comparative analysis of DE transcripts in Figures 2b–c demonstrated that while each treatment cohort displayed a unique transcriptional profile, there was also important sharing of DE transcripts between the cohorts.

Transcripts associated with immune effector function and cell turnover are over-represented in T cells during GVHD: Top gene analysis and Gene Set Enrichment Analysis (GSEA)

To specifically characterize the genes and pathways impacted by allo-activation, GVHD and its prophylaxis, we identified the top 100 transcripts over-represented in the No Rx group compared to Auto controls, as well as the top 100 transcripts that were comparatively under-represented in the No Rx group, (and thus enriched in the Auto controls). This analysis was performed using absolute fold change and statistical significance as ranking parameters (Figure 3 and Supplementary Table 2). The 200 identified transcripts could be largely grouped into 6 functional categories: 1) Cytokines, chemokines and their receptors; 2) Cell turnover; 3) Activation or effector function of T cells; 4) T regulatory; 5) Other; and 6) No annotation available. Figure 3 shows a heatmap of the Z-score of the expression values of these transcripts and their functional classification, with an accompanying graphic depicting the relative frequencies of the 6 categories. Supplementary Table 2 catalogues each of the 200 top transcripts with more detailed information of each transcript including fold change and statistical significance.

Working to parallel our Top Gene analysis using an approach that can overcome the bias associated with cut-off-based analyses, we performed GSEA(23–25). As depicted in the selected examples in Figure 4a–b and Supplementary Table 3, GSEA analysis identified multiple gene sets that were enriched in the leading edge (left side of figure) of the ranked gene list for No Rx vs Auto controls (transcripts enriched in T cells during GVHD-associated allo-activation). This analysis reinforced the concept that two broad categories of pathways were activated during GVHD on a system level: those involved in immune signaling (Figure 4a) and those involved in cell cycle progression and T cell proliferation (Figure 4b).

GSEA also allowed us to interrogate the degree to which T cell activation during GVHD resembled other T cell functional pathways. For this analysis, we focused on gene sets previously identified to be differentially regulated during the protective immune response to viral antigens, given the well-documented phenomenon of heterologous allo-immunity (where T cells activated by viruses have been observed to be serendipitously allo-reactive) and because of the potentially clinically significant overlap between effective downregulation of the alloimmune response, and unintended inhibition of protective immunity(26–29). The dataset for transcriptional control of the T cell response to the lymphocytic choriomeningitis virus (LCMV) is particularly robust, (30) and therefore, we performed a focused analysis of the degree to which GVHD-associated T cell activation overlapped with the immune activation triggered by LCMV. The Molecular Signatures Database (MSigDB) is a collection of annotated gene sets for use with GSEA; therefore, to perform this analysis, we probed the MSigDB database for genesets curated during both primary and secondary T cell responses in the setting of both acute and chronic LCMV infection (Gene Express Omnibus – GSE30962)(31). As shown in the representative example in Figure 4c and in Supplementary Figure 3, we found that 7 of 8 LCMV genesets curated in MSigDB were enriched during the comparison of the No Rx and auto-HCT cohorts. Pathway analysis performed on the combined leading edge genes in these enrichment plots provided key insights into the important nodes of overlap between allo- and virus-mediated T cell activation. As shown in Supplementary Table 4, this analysis demonstrated that of the nine statistically significant pathways identified (p <0.05 using Benjamini-Hochberg corrected P values), all were involved in cell proliferation. The gene lists associated with these pathways provide a new compendium of genes that are involved in both viral and allo—mediated antigen-driven T cell activation.

The gene dysregulation seen in CD3+ T cells during GVHD in NHP parallels current therapeutic and diagnostic efforts in patients

Having enumerated and classified the transcripts that were most highly over- or under-represented during GVHD compared to both HC (which identifies genes differentially expressed during both homeostatic reconstitution and allo-activation) and Auto controls (which specifically identifies those genes differentially expressed during allo-activation), we next queried our dataset individually for selected gene products that are currently, or have been recently, investigated as pharmacologic targets for GVHD prevention or therapy. Among the gene products queried were CCR5, CD2, IL12RB2 and CD38 (Figure 5a). In clinical studies, CCR5 blockade has been shown to be a potential strategy for visceral GVHD prevention(32). Consistent with this, CCR5 expression was highly enriched in the NHP GVHD transcriptome compared to HC; however, as shown in Figure 5a, CCR5 enrichment was also observed during autologous reconstitution, suggesting that its upregulation may be associated with both homeostatic and allo-mediated T cell expansion. This observation is consistent with previous studies in rhesus macaques, which have shown that CCR5 upregulation occurs during homeostatic reconstitution of CD4+T cells in animals treated with a depleting anti-CD4 antibody(33). Likewise, CD2 functions as a costimulatory molecule during antigen-specific lymphocyte activation and has been studied in a phase II trial as a potential target for the treatment of GVHD(34). Similar to what we observed for CCR5, while we found significant enrichment of CD2 in T cells during GVHD when compared to HC, this enrichment did not reach statistical significance when compared to Auto controls (Figure 5a), suggesting that while blocking CD2 may decrease T cell activation, it may not represent a specific target for allo-activated cells. In contrast, CD38, which has recently been shown to be highly correlated with GVHD in patient samples(35), was specifically enriched in NHP GVHD compared to both autologous transplants and healthy controls. Similarly, IL-12/23, which is also being investigated as a potential therapy for GVHD(36), was substantially and specifically enriched T cells during NHP GVHD. The enrichment of these receptor subunits is reflective of a high degree of Th1 skewing in the CD3 cells isolated from the No Rx cohort(37). GSEA also supported a Th1-skewed immune response, with significant enrichment for Th1 versus Th17 transcripts(38) (Figure 5b). Thus, while we did observe dysregulation of selected Th17-related transcripts in the No Rx cohort (including IL17F (but not IL17A), POU2AF1, and TSC22D3, Supplementary Figure 4), the overarching GVHD response in animals receiving no prophylaxis was Th1-polarized.

The addition of serum biomarkers to clinical GVHD staging and grading is being increasingly recognized to improve accuracy in diagnosis and prognostication. We set out to validate known, serum biomarkers of GVHD in our dataset and to provide additional information about the natural history of these markers during untreated GVHD as well as during treatment with current standards-of-care. Several of the key molecules identified from serum protein interrogation include elafin(39), ST-2(40, 41), Reg3A(42), HGF(43), IL2RA(43), TNFRSF1(43), CXCL10(44), IL-10(45), and TIM3(41), Of those proteins expressed by T cells (which excludes elafin, Reg3A and HGF) and are annotated for NHP (which excludes ST-2), we observed robust and specific upregulation of IL2RA and HAVCR2 (the gene encoding TIM3) during GVHD, with significant activation compared to both HC and Auto controls (Figure 5a). While we also observed modest upregulation of IL-10, CXCL10 and TNFRSF1B compared to healthy controls during GVHD (Figure 5a), the expression profiles of these molecules were not significantly different between the No Rx cohort and the Auto controls, suggesting that they may mark both homeostatic- and allo-proliferating T cells.

In addition to investigating genes that have been previously identified as important biomarkers or treatment targets in GVHD, this study has allowed us to begin to query for cross-over of molecules that have been previously recognized as important pathogenic biomarkers in other diseases caused by T cell activation. Consistent with the hypothesis that significant sharing of pathogenic mechanisms exists between diseases of auto- and allo-activation, as shown in Figure 5c and Supplementary Table 2, the two most highly enriched transcripts in our dataset are the ‘calprotectin’ genes, S100A8 and S100A9 (10.4- and 8.6-fold increase compared to Auto controls, respectively), which have been previously implicated as having an important pathogenic role in inflammatory bowel disease(46–48). These data provide support for the existence of common pathogenic mechanisms between these two immune-mediated disorders and point to the calprotectin gene transcripts in T cells as novel putative biomarkers of acute GHVD.

We have shown previously(5) that during NHP GVHD, T cells express high levels of Ki-67 and granzyme B. Consistent with these findings, we observed significant transcriptional upregulation of the gene encoding Ki-67 (MKI67) in the No Rx cohort compared to both HC and Auto controls, with normalization in prophylaxed animals (Figure 5d). Granzyme B gene expression was also found to be significantly upregulated in T cells during GVHD in the No Rx cohort, and was normalized by therapy. Its expression was also higher in the No Rx cohort compared to Auto controls, although this difference did not reach statistical significance. As shown in Figure 5d, flow cytometric measurements of Ki-67 and Granzyme B at the time of gene array analysis paralleled the expression data, and provided independent validation of our transcriptomic approach to defining T cell dysregulation during GVHD.

Pathophysiology of NHP aGVHD includes dysregulation of Notch, Notch targets, as well as evidence for T cell pyroptosis

Seeking to better understand the pathophysiologic drivers of primate GVHD and to provide further external validation to our systems approach, we also queried our dataset for known pathways that have been implicated in GVHD pathogenesis. These included Notch signaling and the inflammasome/pyroptosis pathways (Figure 6). Consistent with its role in T cell expansion during GVHD, the Notch pathway has been shown to be a potential target for GVHD therapy(49). In agreement with these previous studies, we observed significant and specific upregulation of NOTCH1 in CD3 T cells during acute GVHD compared to both HC and autologous controls (Figure 6a, top). Expanding this line of investigation further, we constructed an interactome, using Ingenuity Pathway Analysis (IPA), to identify molecules with which NOTCH1 has been mechanistically linked. We were able to identify a number of transcripts from the Notch interactome, (Figure 6a, bottom) (CCNE1, IL2RA, ADAM19, IFNG, VEGFA, GNB4, MCM5, TUBB4B, MCM3, PYCARD, GZMB, ENO1), that were over- or under-represented in the No Rx cohort, consistent with this pathway being a major driver of T cell allo-reactivity(50).

In addition to identifying gene expression signatures of T cell activation and proliferation during GVHD, we also identified enrichment of a striking number transcripts associated with cell death. To better understand the mechanisms of cell death that may be active during GVHD, we classified those cell death-related transcripts that were significantly over- or under-represented (adjusted P value < 0.05) in the No Rx cohort relative to both Auto controls and to HC (Figure 6b) using DeathBase(51). As shown in the figure, we observed the enrichment of many transcripts known to be pro-apoptotic as well as a relative reduction in the number of inhibitors of apoptosis in our dataset. This is consistent with other published studies of accelerated cell death during aGVHD(52). In addition, as shown in Figure 6b, significant transcriptional activation of several of the canonical pyroptosis transcripts was also observed in the No Rx cohort, including PYCARD (ASC), CASP 1, 4, and their downstream target IL-18. While the inflammasome is most commonly attributed to innate immune cells(53), recent work has documented cleaved caspase-1 in LPS-stimulated leukocytes in patients with GVHD(54), caspase-1-dependent pyroptosis in CD4+ T cells of HIV+ patients,(55, 56) and that CD4+ expression of PYCARD impacts T cell proliferation(56, 57). Together, these observations suggest that the inflammasome is also active in adaptive immune cells. Our data echo these observations and further suggest that the inflammasome may play an important role in the initiation and/or progression of GVHD.

An analysis gene expression in NHP GVHD reveals aurora kinase A as a novel potential molecular target for GVHD

To discover new putative targets for GVHD that could complement current standard of care therapy, we performed a stepwise analysis of over- or under-represented pathways using KEGG(58), Panther(59), Biocarta(60), and Reactome(61) Pathway analysis. The focus of this analysis was to determine which pathways were specific to allo-activated T cells and thus GVHD, that were not present during autologous reconstitution. We found that the pathways involved with control of cell turnover predominated this comparison (Supplementary Table 5), and further found that 4 transcripts (AURKA, RAN, KPNA2 and KIF15) associated with the formation of the mitotic cell spindle(60, 62) were significantly and specifically enriched in T cells during GVHD (Figure 7a). Both to interrogate this relationship on a systems level, and to focus our efforts on genes for which clinically relevant inhibitors are available, we constructed an AURKA interactome using IPA, which included those genes with which AURKA has been shown to interact. We then used GSEA to evaluate whether this interactome showed enrichment in the leading edge of the No Rx cohort versus autologous controls. As predicted from the DAVID Pathway analysis, the aurora kinase A interactome showed significant enrichment (Figure 7b, False Discovery Rate (FDR) = 0.006). Moreover, a closer inspection of the leading edge genes in this analysis revealed that the most enriched transcripts in the aurora kinase A interactome were also highly enriched during immunoprophylaxis with sirolimus (Figure 7c), suggesting that this may be a pathway that is left un-corrected with mTOR inhibition and thus, represents a potentially important partner with sirolimus. As shown in Supplementary Figure 5, we further found that the aurora kinase A interactome displayed little overlap with those proliferation transcripts (discussed above, and shown in Figure 4c, Supplementary Figure 3a and Supplementary Table 4) that were shared during LCMV and allo— mediated T cell activation. Taken together, these results suggest that the expression of the aurora kinase A interaction network is dysregulated during NHP GVHD to a level not seen during autologous reconstitution, and that targeting this pathway may 1) complement residual gene dysregulation that occurs during sirolimus-based GVHD prophylaxis and 2) target cell turnover pathways that may be distinct from those shared with virus-mediated T cell activation.

Inhibition of aurora kinase A mitigates GVHD and prolongs survival in a murine model of acute GVHD

To rigorously validate the Aurora Kinase A pathway, it was further explored using a multi-species analysis pipeline (Supplementary Figure 6). We first performed a proof-of-concept study in our established MHC-mismatched murine model of GVHD(63) using MLN8327, a commercially available selective, reversible small molecule inhibitor of aurora kinase A. As shown in Figure 8, inhibition of the aurora kinase A pathway with in vivo MLN-8327 administration was, indeed, able to mitigate clinical severity and lethality in the murine GVHD model, significantly extending MST from 22.5 (vehicle) to 40.5 days (MLN-8327), p <0.0001). These results predict that targeting aurora kinase A will mitigate GVHD, but that, as expected, for full protection it will likely need to be paired with adjunctive immunosuppression.

Given the results of the NHP and murine studies, we further explored whether, in human cells and patient samples, the aurora kinase A pathway represented a robust target for GVHD prevention. Two assays were performed, both of which suggest that aurora kinase A represents a highly promising pathway for further clinical investigation. In the first assay, an ex-vivo mixed lymphocyte reaction (MLR) was performed with human T cells, and the relative expression of AURKA in the allo-proliferating versus non-proliferating cells was determined. As shown in Figure 8c, allo-proliferative human CD4 and CD8 T cells were 15.3+/− 1.8-fold and 15.5+/− 2.3-fold enriched, respectively, for AURKA transcripts compared to non-proliferating CD4 and CD8 cells. Finally, we compared Aurora Kinase A gene expression in T cells purified from 11 transplant patients who developed acute GVHD and 13 control transplant patients who did not. (Figure 8d). Closely mirroring the NHP transcriptome results, this analysis confirmed enrichment of AURKA gene transcripts in patients with GVHD.

Study Design

This was a prospective cohort study in NHP designed to compare the clinical and immunologic outcomes of transplantation, and to discern the transcriptome of autologous and allogeneic transplant recipients. Experiments were performed in 4 cohorts: (1) Autologous transplants; (2) Allogeneic transplants with no GVHD prophylaxis (No Rx); (3) Allogeneic transplants using sirolimus for GVHD prophylaxis (Siro); (4) Allogeneic transplants using tacrolimus + methotrexate for GVHD prophylaxis (Tac-Mtx). Group sizes of 3–6 were chosen based on power analysis, which assumed that survival in each of the treated cohorts would be at least two-times that in the No Rx cohort. In order to comply with the IACUC ‘reduction’ mandate, based on the “three Rs” of Replacement, Refinement, and Reduction(72), statistical significance of survival data for the prophylaxis cohorts was compared to the No Rx cohort after 3 animals were enrolled, and group sizes were initially limited to 3 animals if statistical significance was obtained at that point. If other immunologic analysis required additional animals to be added to any cohort, this was performed post-hoc. The treatment of animals with immunosuppression was performed in an open-label format (no blinding). However, blinding was performed on all pathologic analysis and on the initial analysis of flow cytometry data and transcriptome sample handling and data processing. Flow cytometric immune phenotyping as well as transcriptional analysis was performed on a prospectively-designed timed sample acquisition schedule.

NHP Ethics Statement

This study used specific pathogen-free, juvenile rhesus macaques that were housed at the Yerkes National Primate Research Center and the Washington National Primate Research Center. The study was conducted in strict accordance with USDA regulations and the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. They were approved by the Emory University, University of Washington, and Seattle Children’s Research Institute’s Institutional Animal Care and Use Committees.

Rhesus HCT

We utilized our previously-described strategy for MHC-mismatched HCT in rhesus macaques(5). Briefly, apheresis was performed after G-CSF mobilization (Amgen, 50mcg/kg for 5 days), and an unmanipulated apheresis product was transplanted into MHC-mismatched recipients. For these studies, we utilized our rhesus macaque MHC typing system(23) to choose donor: recipient pairs, which were all partially MHC-matched, either as half-sibling haplo-identical pairs or as unrelated partially-matched pairs. In the case of the autologous controls, donors were harvested and apheresis products were cryopreserved for a minimum of 2 months prior to the onset of conditioning and subsequent reinfusion of the unmanipulated apheresis product. The transplanted total nucleated cell dose (TNC) and CD3+ cell doses are shown in Supplementary Table 1. The pre-HCT preparative regimen consisted of total body irradiation (TBI) of 9.6 Gy given in two fractions over two days. Irradiation was delivered with a Varian Clinac 23EX (Varian), at a dose rate of 7 cGy/min. Supportive care included antibacterial prophylaxis with polymyxin b and neomycin sulfate, (begun on Day −6 and continued until the absolute neutrophil count was greater than 500 cells/microliter) as well as enorofloxacin (Baytril®, Bayer Healthcare, 7mg/kg IM daily, begun on day −1 and continued until ANC > 500 cells/microliter). Antiviral prophylaxis (cidofovir, Gilead, 5 mg/kg IV weekly) and antifungal prophylaxis (fluconazole 5mg/kg orally, daily) were also employed. Leukoreduced (using an LRF10 leukoreduction filter, Pall Medical) and irradiated (2200 rad) platelet-rich plasma or whole blood was given for a peripheral blood platelet count of ≤ 50 × 10³ per µL or a hemoglobin < 9 g/dL, respectively, or if clinically significant hemorrhage was noted. Blood product support adhered to ABO antigen matching principles.

Four transplant cohorts were compared in this analysis: (1) allo-transplant recipients who did not receive post-transplant immunoprophylaxis, abbreviated as ‘No-Rx’ (n = 4); (2) allo-transplant recipients treated with sirolimus (LC laboratories), given daily for the length of analysis as an intramuscular formulation with doses adjusted to achieve a serum trough of 5–15 ng/dL, abbreviated as ‘Siro’ (n = 4); (3) allo-transplant recipients treated with tacrolimus (given daily for the length of analysis, with doses adjusted to achieve a serum trough of 5–15 ng/dL) and methotrexate (10mg/m2 given on days +1, +3 +8, +11 post-transplant) abbreviated as ‘Tac-Mtx’ (n = 3); (4) autologous transplant recipients who did not receive any immunoprophylaxis abbreviated as ‘Auto’ (n = 6 for flow cytometric data and n = 3 for expression profiling). GVHD clinical score was assessed weekly for allo-transplant recipients as previously described.(5) Briefly, the GVHD clinical score increases with cumulative GI-specific abnormalities (e.g. diarrhea), liver-specific abnormalities (e.g. hyperbilirubinemia) and skin-specific abnormalities (extent and character of rash). Significance of clinical scores was determined by pairwise t-test corrected for multiple hypothesis testing using Benjamini-Hochberg procedure. The Kaplan-Meier product-limit method was used to calculate survival. Differences between groups were determined using log-rank statistics.

Flow cytometry

NHP Microarray and data analysis

Following T cell purification, RNA was stabilized in T cell lysates with RLT buffer (Qiagen) and RNA was purified using RNEasy Column Kit (Qiagen). RNA was quantified using a Nanodrop Spectrophotometer (Thermo Scientific) and purity was confirmed with an RNA 6000 Nano Kit (Agilent). The purified RNA was sent to the Vanderbilt Technologies for Advanced Genomics Core and to the Oregon Health Sciences University Gene Profiling Shared Resource where RNA quantity and quality were verified, followed by cDNA/cRNA synthesis, and target hybridization to GeneChip Rhesus Macaque Genome Array (Affymetrix). The resultant fluorescent signals were processed and normalized using the Robust Multichip Averaging (RMA) Method(73). The microarrays were performed in four batches, with all batches containing samples from both healthy controls and transplanted animals. The “ComBat” algorithm was implemented to adjust for batch effects(74) (shown in Supplementary Figure 2) and probe-sets containing low signal-to-noise measurement were filtered out in order to enhance statistical testing power(75). Principal Component Analysis (PCA) was applied to summarize modes of gene array variance using the Bioconductor MADE4 package(76). Probe-sets were annotated using 1) annotation file from Dr. Robert B. Norgren Jr.(77) 2) the annotation file provided by the chip manufacturer (release 33, 10–30-12) and 3) data provided by Ingenuity Systems (Ingenuity Systems, www.ingenuity.com) for the small number of probe-sets that were not annotated by the chip manufacturer.

Analyses of gene differential expression (DE) was performed using an empirical Bayes moderated t-statistic, with a cutoff of 0.05, corrected for multiple hypothesis testing using Benjamini-Hochberg procedure and an absolute fold change cutoff > 1.4 with the Bioconductor limma package(78).

Using pathway analysis systems IPA (Ingenuity Systems, www.ingenuity.com) and DAVID(79) (encapsulating Biocarta, KEGG, and Reactome pathway analysis), we interrogated our DE gene lists for significantly enriched pathways. Significance was estimated using a right-tailed Fisher’s exact test using the Benjamini-Hochberg procedure to account for multiple testing.

Further analysis of differentiating characteristics between untreated allo-HCT and control arrays involved Gene Set Enrichment Analysis (GSEA) as described(23–25). GSEA can identify whether the members of a gene set (a collection of which are housed in the Molecular Signatures Database (MSig-DB) - Broad Institute, Boston, MA) are enriched in an independent rank-ordered profile of genes that are differentially expressed between two experimental groups. In this manner, GSEA is able to provide definitions of overrepresented biological functions without implicit bias associated with cut-off based analyses. In the current analysis, gene sets were ranked using a signal to noise ratio difference metric with 1000 permutations of gene set labels.

The gene list for cell death was obtained from DeathBase(51). The notch, aurora kinase A and sonic hedgehog interactomes were created using Ingenuity Systems. Unless otherwise noted, gene analysis was carried out using R: Bioconductor.

Murine HCT

Murine models and ethics statement: B10.BR/SgSnJ (H2k) mice were purchased from The Jackson Laboratory (Bar Harbor, ME). B6 (H2b: CD45.2) mice were purchased from the National Institutes of Health (Bethesda, MD). Mice were housed at the University of Minnesota Hospital in a specific pathogen-free facility in micro-isolator cages. Donor and recipient mice were between 8 to 10 weeks of age at the time of HCT. Murine experiments were conducted in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health, and were approved by the University of Minnesota Institutional Animal Care and Use Committee.

Transplant Regimen: Lethally irradiated (1000 cGy) B10.BR (Jackson Labs) recipients were transplanted with five different treatment regimens: (1) Bone marrow control: These animals received 1×10⁷ bone marrow (BM) cells alone from fully MHC-mismatched B6 mice (NCI, n = 15 in two independent experiments). (2) BM + 5 × 10⁶ splenocytes from B6 mice + MLN8237, a selective aurora kinase A inhibitor (Chemie Tek) (n = 18 in two independent experiments). MLN8327 was given by gavage daily at a dose of 30mg/kg from day 0→28 (80). (3) BM + splenocytes + LDE225, an inhibitor of the sonic hedgehog coreceptor Smoothened (Smo, Chemie Tek (n = 27 in three independent experiments). LDE225 was given by gavage daily at a dose of 40 mg/kg from day 0 →28 (81). (4) BM + splenocytes + the vehicle for MLN (10% 2-hydroxypropyl-b-cyclodextrin, Sigma, 0926) (n = 27 in three independent experiments). (5) BM + splenocytes + the vehicle for LDE225 (PEG300, Sigma, 91462) (n = 27 in three independent experiments). Mice were monitored for clinical signs of GVHD (with scores based on weight loss, posture, activity, fur texture and skin integrity) twice weekly as previously described(82). The Kaplan-Meier product-limit method was used to calculate survival. Differences between groups were determined using log-rank statistics.

MLR studies

PBMC were isolated from 5 healthy donors by Ficoll (GE Healthcare) fractionation. Bulk T cells were isolated using a negative immunomagnetic selection kit (Miltenyi Biotech). The negative fraction (“Responder” T cells) were labeled with Cell Trace Violet (Invitrogen) and were co-cultured with the irradiated (3500 cGy) positive immunomagnetic fraction from an allogeneic donor. MLRs were cultured for 4 days at 37°C in XVivo15 (Lonza) supplemented with 10% FCS, 2 mM glutamine, penicillin– streptomycin and gentamycin. On day 4, cells were stained for CD3, CD4, CD8 and sorted using the following gates 1) CD3+, CD4+, CTV+, 2) CD3+, CD4+, CTV-, 3) CD3+, CD8+, CTV+, 2) CD3+, CD8+, CTV−. Total RNA was obtained from cell lysates using RNAEasy Mini kit (Qiagen), RNA yields were measured using a Nanodrop Spectrophotometer (Thermo Scientific, Waltham, MA) and cDNA was synthesized using iScript kit (BioRad) normalizing the amount of input RNA across samples. AURKA, SHH, and beta-2 microglobulin (B2M) transcript levels were measured by drop digital PCR (QX200 Bio-Rad) as previously described(83) using exon-spanning primers dHsaEG5024289 for AURKA, dHsaEG5006220 for SHH, and dHsaEG5020739 for B2M (Bio-Rad, Hercules, CA). AURKA and SHH transcript levels were then normalized to B2M, and the ratio of the transcript levels in the allo-proliferating (CTV-) and non-proliferating (CTV+) cell fractions determined.

Transcriptional studies of patient samples