Peripheral tissues reprogram CD8+ T cells for pathogenicity during graft-versus-host disease

Graft-versus-host disease (GVHD) is a life-threatening complication of allogeneic stem cell transplantation induced by the influx of donor-derived effector T cells (TE) into peripheral tissues. Current treatment strategies rely on targeting systemic T cells; however, the precise location and nature of instructions that program TE to become pathogenic and trigger injury are unknown. We therefore used weighted gene coexpression network analysis to construct an unbiased spatial map of TE differentiation during the evolution of GVHD and identified wide variation in effector programs in mice and humans according to location. Idiosyncrasy of effector programming in affected organs did not result from variation in T cell receptor repertoire or the selection of optimally activated TE. Instead, TE were reprogrammed by tissue-autonomous mechanisms in target organs for site-specific proinflammatory functions that were highly divergent from those primed in lymph nodes. In the skin, we combined the correlation-based network with a module-based differential expression analysis and showed that Langerhans cells provided in situ instructions for a Notch-dependent T cell gene cluster critical for triggering local injury. Thus, the principal determinant of TE pathogenicity in GVHD is the final destination, highlighting the need for target organ–specific approaches to block immunopathology while avoiding global immune suppression.

, included in a separate Excel file, lists for the 19 highly preserved F→M WGCNA gene modules, and the 2 skin-specific B6→129 WGCNA modules: the over-represented Gene Ontology categories (Supplemental Table 1A), the putative driver genes identified on the basis of high intra-modular connectivity (Supplemental Table 1B), and the transcription factors predicted to act as upstream regulators (Supplemental Table 1C). Supplemental Table 2, included in a separate Excel file, lists the clinical details for the 5 HSCT patients included in this study.
Supplemental Table 3, included in a separate Excel file, contains the results from the preprocessing performed to test for LC contaminants in LC-replete versus LC-depleted epidermal samples (Supplemental Table 3A), and lists the intestine-and skin-specific genes subtracted from the analysis (Supplemental Table 3B). Supplemental Figure 1 shows the visualization of sample relationship without gene subtraction. Supplemental Figure 2 shows the overlap between TE gene expression profiles in the F→M and B6→129 experimental models, and the similarities between F→M and B6→129 based WGCNA modules. Supplemental Figure 3 depicts the analytical pipeline followed, based upon correlation network analysis and downstream validation.
Supplemental Figure 4 shows the network representation of the B6→129 skin-specific modules MD and ME, and their overlap with the F→M skin-specific module M28. Supplemental Figure 5 shows the expression of JAG1 and DLL4 Notch ligands by LC compared to non-LC populations in the epidermis of mice developing GVHD, and the effect of in vitro Notch blockade on the capacity of male LC to stimulate IFN-g generation by activated MataHari T cells. Supplemental Figure 6 shows the kinetics of LC host-to-donor turn over and the effect host LC depletion upon TE accumulation in the epidermis in independent murine BMT models of GVHD. c) Blood -Erythrocytes were removed from whole blood samples by hypotonic lysis with distilled water. Cells were re-suspended in FACS buffer for counting and immunolabelling. d) Small intestine -The entire small intestine was flushed and rinsed with 40 ml of ice cold harvest medium (PBS, 2% FCS, 1% penicillin-streptomycin; Lonza, UK), and sectioned into ~0.5 cm pieces. The intestinal pieces were incubated with detaching medium (HBSS, 5% FCS, 1% penicillin-streptomycin, 5 mM EDTA; Lonza, UK) at 37°C with shaking at 150 rpm for 60 min.
The supernatant, containing the IEL, was passed sequentially through a 100 µm and a 40 µm cell strainer and enriched for lymphocytes by density gradient centrifugation with Ficoll Paque TM Plus (GE Healthcare, UK). The intestinal pieces were further incubated in digestion medium (RPMI, 2% FCS; Lonza, UK; 200 U/ml collagenase IV; LifeTechnologies, USA; 200 U/ml DNAse I; Sigma, UK) at 37°C with shaking at 150 rpm for 60 min, dissociated and passed sequentially through a 100 µm and a 40 µm cell strainer. The cell suspension, containing the LP cells was enriched for lymphocytes by density gradient centrifugation with Ficoll Paque TM Plus. Cells were re-suspended in FACS buffer for counting and immunolabelling. e) Skin -Epidermal and dermal immune cells were isolated from the skin in accordance with the protocol described by Henri et al. (1), with minor modifications. Briefly, the body skin was cut into ~1x1 cm pieces, after having been shaved and the subcutaneous fat removed, and the ears were s3 split in two parts (ventral and dorsal). The pieces of skin were incubated overnight in dispase medium (HBSS, 2% FCS; Lonza, UK; 2.5 mg/dl dispase II; Sigma, UK) at 4°C, and the epidermal and dermal sheets were separated and cut into ~0.5 cm fragments. The epidermal fragments were vortexed, mashed and passed sequentially through a 100 µm and a 40 µm cell strainer to disintegrate the remaining tissue and create a cell suspension. The dermal fragments further incubated in digestion medium at 37°C with shaking at 150 rpm for 60 min, dissociated and passed sequentially through a 100 µm and a 40 µm cell strainer. Cells were re-suspended in FACS buffer for counting and immunolabelling.

Flow cytometry
Cells were plated out at up to 1x10 6  s6

Immunofluorescence staining of epidermal sheets
Epidermal sheets were fixed in ice-cold acetone, blocked in a 0.25% fish gelatin, 10% goat serum solution and stained for LC with CD207 (eBioL31) (eBioscience, USA). Images were acquired with the Leica DMI4000b fluorescence microscope (Leica Microsystems, UK), using the Leica Application Suite Advanced Fluorescence v3.2 software (Leica Microsystems, UK). Raw image data were processed with Fiji (ImageJ).

Microarray analysis
A total of 36 RNA samples from the B6®129 model and 61 RNA samples from the F®M model were used for microarray analysis, representing 3 independent biological replicates for each of the tissues/compartments, treatment conditions and time points. Hybridized arrays were scanned with a GeneChip 3000 7G scanner (Affymetrix, USA) (1 batch for B6®129 samples; 3 batches for F®M samples) and the image data processed to generate .cel files. Expression Console Software, version 1.4.1 (Affymetrix, USA) was used to generate quality control statistics for each sample; samples with a high mean absolute deviation of the residual for a chip versus all chips in the data set (≥ data set average + 2 SD), a high mean absolute relative log expression (≥ data set average + 2 SD) and a low area under the curve (AUC) for a receiver operator curve (ROC) comparing the intron controls to the exon controls (≤ 0.8) were considered to have poor hybridisation quality and excluded from the analysis (1 sample in the B6®129 data set; 3 samples in de F®M data set). Raw sample expression signals were background subtracted, quantile normalized, and the probe level data were summarized using the Robust Multi-array Average algorithm (2, 3) implemented in the oligo BioConductor R package (4). The ComBat algorithm from the sva BioConductor R package (5) was employed to adjust for batch effects.
Transcripts identified through multiple probes were collapsed based on maximum expression values using the CollapseDataset module of GenePattern software (Broad Institute) (6).

Data set pre-processing and gene subtraction
Pre-processing of the data sets was performed to test for LC contamination in the microarray samples used in experiments shown in Fig. 8. We found no significant difference between LC-s7 replete and LC-depleted epidermal samples in terms of the levels of expression of LC signature genes (7) (Supplemental Table 3A). To exclude any effect of tissue-specific transcripts derived from non-T cell contamination of individual samples, transcripts identified as highly specific for the intestine and skin in the PaGenBase database [specificity measure (PMS) > 0.9] (8) were subtracted from the entire dataset prior to analysis (Supplemental Table 3B).

Samples relationship visualization
Multivariate statistical analysis methods implemented in the stats R package, in particular multidimensional scaling, were applied to perform dimensionality reduction of the datasets and visualization of the samples relationships. Similarities between groups were evaluated by hierarchical clustering, using average-linkage upon the Pearson's correlation-based dissimilarity matrix; the validity and stability of the clusters was assessed using a non-parametric bootstrap as implemented in the pvclust R package (10). Additional hierarchical clustering and heat maps were produced using the matrix visualization and analysis platform GENE-E (Broad Institute, USA).

Differential gene expression
The limma BioConductor R package was used to perform analyses of gene differential expression, using an empirical Bayes moderated t-statistic corrected for multiple hypothesis testing using Benjamini-Hochberg procedure, with a cut-off of false discovery rate (FDR) ≤ 0.05, and an absolute fold-change cut-off ≥ 2.0.  Table 4).

Weighted gene co-expression network analysis (WGCNA)
Scale-free network topology analysis of microarray expression data was performed using the WGCNA R package, as previously described (16,17). A signed hybrid weighted correlation network was constructed using a Pearson correlation matrix created from the pairwise comparison between all pairs of genes, and a soft thresholding power β=8. The topological overlap was calculated as a measure of network interconnectedness, and genes were grouped by average linkage hierarchical clustering on the basis of the topological overlap dissimilarity (1topological overlap). Module eigengenes were calculated using a dynamic tree-cutting algorithm and merging threshold function at 0.25. The modules identified were correlated to the sample traits using a binary vector representation of the tissues of origin and study groups. To validate the microarray analysis, preservation of the WGCNA modules identified in the F→M model dataset was tested against the dataset of the B6→129 model, using the R function s9 "modulePreservation" in the WGCNA R package, as previously described (18). Results were interpreted according to the following thresholds for Zsummary: if Zsummary≥10, module "strongly preserved"; if Zsummary ≥2 and <10, module "weak to moderately preserved"; if Zsummary<2, module "not preserved". Significance of the Zsummary scores was calculated by gene permutation testing.
Intramodular hub genes, which are genes that have the highest number of connections within a module, were identified on the basis of having eigengene-based connectivity (kME) > 0.8 and gene significance (GS) > 0.2. Visualisation of the modules network of gene connections was accomplished with the Cytoscape v3.5 software (19).

Prediction of upstream master regulators
The Cytoscape plugin iRegulon was used to predict the transcriptional regulatory network underlying each of the modules of co-expressed genes, as described by Janky et al. (20). Briefly, for each gene set, the regulatory sequences in 20 Kb around the transcription start site were analyzed, and motif prediction was performed using an enrichment score threshold of 2.0, a ROC threshold for AUC calculation of 0.03%, and a rank threshold for visualization of 5000. Candidate transcription factors were predicted with a maximal FDR for motif similarity of 0.05.

Overrepresentation analysis
The Web-based Gene Set Analysis Toolkit (WebGestalt), a suite of tools for functional enrichment analysis, was used to identify overrepresented KEGG pathways categories and translate gene lists into functional profiles. Statistical significance of KEGG pathways enrichment was calculated based on hypergeometric distribution statistics, adjusting the false discovery rate using the Benjamini-Hochberg procedure.

Module-based differential expression analysis (modDE)
In order to quantify potentially significant associations at the level of WGCNA-derived gene modules between tissue types and/or between conditions, we developed a gene based association test that takes advantage of the magnitude and sign of the differential expression effect size for each gene. We start by translating each gene's p value into a chi-squared statistic s10 Xi with one degree of freedom. In its simplest form, the test is designed to assess whether genes are consistently over-expressed in a set of cases compared to controls. Therefore, the test statistic T is the sum over all genes of the Xi, if Xi is positive, or 0, if Xi is negative. The distribution of this statistic under the null is obtained by assuming that each gene has a 50% chance of being over-expressed and 50% of being under-expressed. Note that this can be refined by computing the genome-wide probability of over-expressed genes, in case that the proportion differs from 50%. With either of these assumptions, the probability of observing a number K of genes overexpressed among the n genes in the module can be computed using a binomial distribution. For a given K, the statistic T is chi-squared with K degrees of freedom (because it is the sum of K one degree of freedom chi-square distributions). The probability of observing a test statistic greater than the observed value T can then be computed using a weighted sum of probabilities, summing over all possible values of K. The derivation is analogous if all genes are expected to be underexpressed. If the expectation is a mixture of over and under expressed genes, the T statistic is then summed over all genes in a manner that reflects that combination.

Module E
Module D CD8 + T cells infiltrating the epidermis in GVHD were excluded by gating. Bottom -Summary data of the difference in MFI staining between each sample and the respective isotype control (n = 2, graphs show mean ± SD; §, isotype control MFI > sample MFI). No detectable DLL1 and JAG2 expression was found on any epidermal population (data not shown). (B) Effect of in vitro LY411575 or PBS (untreated) exposure upon IFN-γ generation by concanalavin-preactivated MataHari T cells upon interaction with male or female LC. Data derived from 6 independent experiments, graphs show mean ± SD. ****p ≤ 0.0001, two-way ANOVA with Holm-Sidak correction for multiple comparisons.