Spatially resolved immune exhaustion within the alloreactive microenvironment predicts liver transplant rejection

Allograft rejection is a frequent complication following solid organ transplantation, but defining specific immune subsets mediating alloimmunity has been elusive due to the scarcity of tissue in clinical biopsy specimens. Single cell techniques have emerged as valuable tools for studying mechanisms of disease in complex tissue microenvironments. Here, we developed a highly multiplexed imaging mass cytometry panel, single cell analysis pipeline, and semi-supervised immune cell clustering algorithm to study archival biopsy specimens from 79 liver transplant (LT) recipients with histopathological diagnoses of either no rejection (NR), acute T-cell mediated rejection (TCMR), or chronic rejection (CR). This approach generated a spatially resolved proteomic atlas of 461,816 cells derived from 98 pathologist-selected regions of interest relevant to clinical diagnosis of rejection. We identified 41 distinct cell populations (32 immune and 9 parenchymal cell phenotypes) that defined key elements of the alloimmune microenvironment (AME), identified significant cell-cell interactions, and established higher order cellular neighborhoods. Our analysis revealed that both regulatory (HLA-DR+ Treg) and exhausted T-cell phenotypes (PD1+CD4+ and PD1+CD8+ T-cells), combined with variations in M2 macrophage polarization, were a unique signature of TCMR. TCMR was further characterized by alterations in cell-to-cell interactions among both exhausted immune subsets and inflammatory populations, with expansion of a CD8 enriched cellular neighborhood comprised of Treg, exhausted T-cell subsets, proliferating CD8+ T-cells, and cytotoxic T-cells. These data enabled creation of a predictive model of clinical outcomes using a subset of cell types to differentiate TCMR from NR (AUC = 0.96 ± 0.04) and TCMR from CR (AUC = 0.96 ± 0.06) with high sensitivity and specificity. Collectively, these data provide mechanistic insights into the AME in clinical LT, including a substantial role for immune exhaustion in TCMR with identification of novel targets for more focused immunotherapy in allograft rejection. Our study also offers a conceptual framework for applying spatial proteomics to study immunological diseases in archival clinical specimens.

INTRODUCTION T-cell mediated rejection (TCMR) remains the most frequent complication after liver transplantation (LT), occurring within the rst six months in up to 35% of adult LT recipients [1][2][3] . While TCMR is generally responsive to treatment with pulse corticosteroids, adjustment of maintenance immunosuppression regimens is key for preventing future TCMR episodes 4 . Ultimately, up to 10% of patients will develop steroid resistance and have recurrent episodes of TCMR. The diagnosis of TCMR hinges upon histological examination of a core biopsy stained with hematoxylin & eosin by a clinical pathologist using Rejection Activity Index (RAI), a composite score ranging from 0-9 based on severity of portal in ammation, bile duct in ammation, and venous endotheliitis 5,6 . After its inception following a Banff consensus conference in 1995, the RAI has become the gold standard to establish the diagnosis of TCMR and guide treatment strategies in clinical LT. There have been minimal changes in the RAI since it was rst introduced, with additional criteria for antibody mediated rejection (AMR), a rare entity in LT, in 2016 6 .
In parallel, options for both induction and maintenance immunosuppression in LT have not changed substantially since the 1990s and rely on therapeutics that cause non-speci c suppression of entire leukocyte populations. For instance, the two mainstay treatments broadly suppress the T-cell compartment (calcineurin inhibitors), or they function by globally inhibiting both macrophages and T-cells (corticosteroids) 7 . Thus, the absence of speci c targeting for TCMR-associated immune subpopulations in LT results in both suboptimal prevention and treatment of TCMR episodes, as well as a variety of unintended, and often severe, adverse medication side effects.
Improving our understanding of the complex alloimmune microenvironment (AME) in clinical LT would enable development of novel, focused, and personalized immunotherapies. Donor-derived antigen presenting cells (APCs) expressing allograft antigen on both MHC I and II can activate host CD8 + and CD4 + T-cells via the direct pathway, ultimately leading to tissue damage via Fas-FasL or granzyme/perforin production and secretion of pro-in ammatory cytokines 8 . The indirect pathway, which has been implicated in late TCMR, is mediated by recipient APCs in ltrating the allograft over time 8 . However, deeper characterization of graft-in ltrating leukocytes driving TCMR in clinical LT is still needed.
When compared to experimental heart and kidney transplantation, small and large animal models of LT are more technically demanding while offering a lower threshold for tolerance induction and thus less opportunity to recapitulate alloimmunity in clinical TCMR 9 . Examination of clinical samples has been limited by the tiny amount of tissue available from a core needle biopsy specimen. The INTERLIVER study examined over 200 clinical LT biopsies using genome microarrays and archetypal analysis, and differentially expressed (DE) genes involving both effector T-cell and injury-related pathways were identi ed in the small subset of biopsies with TCMR. Supervised and unsupervised molecular classi ers based on the top 30 DE genes had only a modest ability to predict histological TCMR (AUC 0.57 and 0.70, respectively) 10 . A more recent histologic study analyzing post-LT biopsies demonstrated that CD8 + T-cells form an immune synapse with APCs, with an association between segregation of CD3 and CD45 molecules on CD8 + T-cells, immunosuppression weaning failure, and development of TCMR 11 . Thus, key features driving the intrahepatic alloimmune response, including composition and phenotype of alloreactive T-cell subpopulations and interactions between innate and adaptive cells, remain elusive.
CD4 + CD25 + FoxP3 + regulatory T-cells (Tregs) have been a central focus in both experimental and clinical LT 12 . Despite substantial evidence that Tregs are central mediators of rejection and immune tolerance, clinical trials designed to expand Treg either via therapeutic intervention or cellular therapies have not yet resulted in positive clinical outcomes 13 . The programmed death 1 (PD1) pathway has also emerged an important physiologic immune checkpoint to maintain peripheral T-cell tolerance and regulate adaptive immune responses particularly during chronic antigen exposure 14 . PD1 can be expressed on both B and T-cell populations, including Tregs upon activation, with constant high expression levels following sustained antigen exposure. The PD1 pathway antagonizes T-cell receptor (TCR) engagement and CD28 co-stimulation signals, attenuating downstream cytokine production, proliferation, cell metabolism and survival; thus, ultimately moderating T-cell activity 15,16 . The role of PD1 signaling in transplantation is not well de ned, with preliminary studies on heart and kidney allografts implicating PD1-related signals as markers of allograft rejection [17][18][19][20][21][22] . However, in a recent study in clinical LT, ow cytometric analysis failed to demonstrate a difference in PD1 expression in allograft-in ltrating T-cells isolated from liver explant (n = 5), rejection (n = 7), and no-rejection liver biopsies (n = 7) 23 . Thus, detailed study of the relationship between different regulatory and in ammatory immune cell populations is critical for de ning the important aspects of the AME, optimizing identi cation of predictive biomarkers of TCMR, and identifying focused targets for immunotherapy.
Here, using multiplexed proteomics-based Imaging Mass Cytometry (IMC) analysis, we developed an immune cell phenotyping and analysis pipeline that enabled granular, single-cell characterization of over 30 discrete immune cell types, with spatial assessment of the AME in a large population of post-LT patients with no rejection (NR), TCMR, and progression to CR. We de ned signi cant cell-to-cell interactions and identi ed spatial motifs as well as predicted single cell phenotypes associated with TCMR. This approach revealed that within the AME the evolution of the immune response in TCMR was associated with intragraft presence of speci c T-cell subpopulations exhibiting an exhaustion phenotype.
However, these PD1 + T-cells were lost in CR, suggesting an important role as mediators and potential biomarkers of TCMR. Furthermore, we showed that lymphocytes and macrophages are spatially organized into aggregates in which strong interactions among PD1 + and effector T-cells exist as well as between CD8 + T-cells and speci c macrophage subpopulations. Collectively, our data offer a detailed and spatially conscious atlas of immune in ltrates in the liver AME during TCMR episodes that represent putative in situ biomarkers of rejection. Our data provide a framework for histologic assessment of complex immune microenvironments at single cell resolution in archival clinical samples, which can inform development of novel clinical assays improving treatment speci city and support development of novel targets for immunotherapy.

METHODS
This study was approved by the Health Science Campus Institutional Review board of the University of Southern California (HS-18-00708).

Patients
LT recipients were retrospectively identi ed using our institutional transplant database. Patients > 18 years at the time of transplant who underwent biopsy of their liver allograft to rule out suspected TCMR or patients with CR undergoing re-transplantation between 1/2000 and 12/2021 met inclusion criteria. Patients were excluded if the histologic diagnosis was associated with reactivation or concurrent viral infection (i.e., Hepatitis C or cytomegalovirus), anatomic causes of graft dysfunction (i.e., vascular stenoses and/or biliary strictures), or advanced brosis (bridging brosis based on Trichrome staining). Pathology reports were reviewed by a pathologist with expertise in LT to prioritize selection of patients with RAI ≥ 4 for the TCMR group (n = 41 patients, 58 ROIs, median RAI of 5 (Interquartile range (IQR) 5-6)). LT recipients who did not have evidence of rejection on their biopsy (RAI = 0) were selected for the NR group (n = 24). The CR patients (n = 14) were identi ed at the time of re-transplant for CR with histologic con rmation of CR in the explant.

Clinical Data and Demographics
Demographics and clinical parameters were obtained via comprehensive chart review and included age, sex, ethnicity, race, age at transplant, serum biochemistries, immunosuppression regimens, and all biopsy data including indication for biopsy, timing of biopsy in relation to LT, and pathology reports. Relevant demographic variables are summarized for the cohort in Extended Data Table 1. For consistency, RAI score and detailed breakdown of sub-scores were independently reviewed by a liver specialized pathologist. This review showed close agreement with the pathologic evaluation performed at the time of biopsy.

Imaging Mass Cytometry
Formalin-xed para n embedded (FFPE) tissue sections of liver biopsy specimens or explants (4µm) were selected by the pathologist to identify 1mm 2 regions of interest (ROI) for IMC acquisition focusing on representative periportal regions of the biopsies used in the clinical assessment of RAI. The SC2 Core Facility at Children's Hospital-Los Angeles performed all staining and image acquisition for this study. Slides were stained using a custom 22-marker antibody panel. Structural markers included two nuclear intercalator dyes, collagen, CD31 (vascular endothelium), and CK7 (bile ducts). Immune lineage markers included CD3, CD4, CD8, CD20, CD68, CD11b and functional or phenotypic markers included CD279 (PD1), FoxP3, Ki67, and Granzyme B among others (Extended Data Table 2). IMC staining was performed using techniques described previously 24 . Ninety-six ROIs (average of 1.2 ROI/patient) were ablated using the Hyperion Imaging System (Standard Biotools) at a power range of 3.5-4.5 with a laser frequency at 200Hz. Data were supplied as .txt and .mcd les for use in segmentation and downstream analyses.

Image Pre-Processing and Segmentation
Pre-processing steps were completed using the MATLAB package MAUI (MBI Analysis User Interface) 25 .
CD68 was used as the basis for channel spillover correction and noise removal and channel aggregate removal steps were implemented individually on each channel. Pre-processing was conducted in three batches by clinical outcome (NR, TCMR, CR), to account for staining background and noise differences between disease states. Cell segmentation was performed using Mesmer (DeepCell) and following the Bodenmiller Steinbock pipeline 26 . Image pre-processing was performed in MATLAB (version R2022b) and Python (version 3.10.8).

Single Cell Phenotyping
Cell segmentation outputs were loaded into R (version 4.2.2) to perform downstream analysis. Patient ID and clinical group identi ers were added to the Single Cell Experiment object 27 . Data were arcsin transformed using a cofactor of 5 and standardized by channel to account for differences in signal intensities.
Metaclusters were identi ed using a supervised clustering approach outlined in Extended Data Table 3.
Labelling accuracy was veri ed by reviewing concurrent metacluster label and channel expression on tissue sections. Masks were used to visualize cell labels ("cytomapper::plotCells") 28 . TIFF images were scaled, and channel signals were normalized and visualized individually ("cytomapper::plotPixels"). For each patient, metacluster proportions were calculated using the overall cell count as baseline and statistically compared across clinical groups. Subclustering was performed on the ve most relevant immune metaclusters (CD4 + T-cells, CD8 + T-cells, B cells, macrophages, and monocytes) using a semisupervised approach. A total of 30 subclusters were identi ed, leading to a nal 32 immune clusters and 9 non-immune clusters in the overall dataset. Dimensionality reduction was performed using t-Distributed Stochastic Neighbor Embedding (t-SNE) to visualize meta-and subclusters by clinical outcome 29 . T-SNE was also used to visualize possible batch effects between patients. Batch correction was performed using the mutual nearest neighbors (MNN) method, but ultimately not used for downstream analysis to avoid possibly also eliminating biological differences present in the data 30 .

Trajectory Inference
To investigate whether cell phenotypes identi ed via IMC represented a pseudotemporal evolution of the alloimmune microenvironment in LT rejection, we performed trajectory inference on each metacluster. Metacluster-speci c dimensionality reduction was performed by Uniform Manifold Approximation and Projection (UMAP) 31 . Trajectories were identi ed by Slingshot, an algorithm that can model branching lineages in single-cell data 32 . To ensure proper orientation of each trajectory, a coarse clustering was performed using k-means (k = 2, except for CD8 + T-cells where k = 5) and the cluster with the highest proportion of cells from NR samples was set as the initial cluster.

Spatial Analysis
The k-nearest neighbor approach (k = 10) was used to create the cell-cell interaction graph, which was visualized on tissue using the "imcRtools::plotSpatial" function 26 . Neighborhood analysis ("imcRtools::testInteractions") was implemented on each clinical subset to analyze pairwise interactions between metaclusters and between subclusters and to compare differences across clinical outcomes. Cell-cell interactions were calculated using permutation testing (1,000 permutations, α = 0.01) to determine whether cell types interact more (attraction) or less (avoidance) frequently than random permutations. Graph network analysis using igraph was used to visualize subcluster interactions 33 .
Cellular neighborhood analysis was implemented ("imcRtools::aggregateNeighbors") using the constructed k-nearest neighbor spatial graph (k = 10). Cells were re-clustered based on the cell types in their direct spatial neighborhood to obtain 9 cellular neighborhoods (CN). Cell type abundance of each CN was visualized on a heatmap to aid CN annotation. For each patient, CN proportions were calculated, visualized, and statistically tested to detect any differences across clinical group. CNs were also visualized on the tissue to detect any visual differences in spatial composition across clinical group.
Due to the spatial relevance of TCMR in ltrates to the vascular endothelium (based off clinical Banff criteria RAI scoring), the median distance of each cell type (meta-and subclusters) to endothelial cells was calculated and compared across clinical groups.

Predictive Modeling
Logistic LASSO (least absolute shrinkage and selection operator) regression was used to build predictive models of NR vs TCMR, TCMR vs CR, and NR vs CR. LASSO is a shrinkage method that aids in feature selection and avoids over tting. LASSO adds an L1 regularization term (sum of absolute values of the coe cients) so that the selected coe cients minimize the loss function , where y is the vector of the binary clinical outcome, X is the feature matrix, is the vector of coe cients, and is the regularization coe cient. A 5-fold cross-validation (CV) technique was used to nd the optimal value. For each comparison, model building was done using those cell types found to be statistically signi cant in the pairwise comparisons as input. Bootstrapping, a sampling with replacement technique, was implemented to rank the importance of all features (5,000 iterations). In each iteration, logistic LASSO regression was implemented on a subset of the data and non-zero coe cients were stored. Variable frequency was determined, and variables with ≥ 50% frequency were selected for the nal model. To evaluate model performance, data was cross validated by randomly splitting into training and validation sets at a 75/25 ratio with 1,000 iterations. In each iteration, the model was trained on the training set using the features identi ed during bootstrapping. Clinical outcome was then predicted on the validation set and stored alongside performance metrics (sensitivity, speci city, accuracy, and area under the curve [AUC]). Final model coe cients were obtained by averaging all coe cients. Final model performance was calculated using the evaluation metrics obtained from all iterations (mean ± SD). Receiver operating characteristic (ROC) curve was calculated using the median prediction of each patient. Correlation between the median prediction and actual clinical outcome was calculated using Spearman correlation and signi cance was tested using the Wilcoxon Rank-Sum test.
nCounter Transcriptomic and TCR Expression Analysis A subset of representative tissue samples (4 NR and 4 TCMR) was selected, and at least ve 5µm FFPE sections per block were combined for RNA extraction using the Rneasy Kit (Qiagen). Extracted RNA was quanti ed using the NanoDrop system (Thermo Fisher Scienti c), and 200 ng of total RNA was used for gene expression analysis. Samples were processed using the nCounter Nanostring platform and the PanCancer Immune Pro ling and T-cell repertoire panels according to the manufacturer's guidelines (NanoString Technologies). Raw counts were normalized using internal positive standards and housekeeping genes with the nSolver Analysis 4.0 and Advanced Analysis 2.0 software (NanoString Technologies). Expression of scaled log2 gene counts were visualized using heatmaps to determine expression differences between NR and TCMR samples. Raw counts of TRAV, TRBV, and TRGV and TRDV genes were expressed as a proportion among total TRAV, TRBV, or TRGV and TRDV gene counts, respectively, for each patient and normalized using the median value from the healthy control group. Publicly available data from a study of 6 patients who underwent IL-2 therapy and subsequently had rejection episodes within 6 months post treatment was used to validate some of these ndings 34 . The fold-change in mean gene expression between NR and TCMR as well as baseline and 4-weeks post treatment was compared to show similarity in gene upregulation.

Statistical Analysis
The Shapiro-Wilk test was used to test for normality. One-way ANOVA and two-sample t-test were used to analyze parametric data. Kruskal-Wallis and Wilcoxon Rank-Sum tests were used to analyze nonparametric data. P-values were corrected for multiple testing using the Holm method. Number of cells for each meta-and subcluster were reported as median [Q1, Q3]. A 0.05 p-value cut-off was used throughout the analysis to determine statistical signi cance. All statistical tests were carried out in R (version 4.2.2).

Major cell types and proportions in liver allografts with and without rejection
We applied IMC to 24 NR liver core biopsies, 41 biopsies with proven TCMR, and 14 CR samples using our customized analysis pipeline (Fig. 1a). By segmenting the acquired 96 multiplexed images, we generated a single cell atlas of the AME containing a total of 461,816 cells (average 4,811 ± 2,291 cell/ROI) which were classi ed into 10 main cell populations or 'metaclusters'. We evaluated raw image signals, postsegmentation dimensionality reduction (t-SNE) of individual markers, immune metaclusters by patient, and difference in mean fold change expression of all markers among the three clinical groups. (Extended Data Fig. 1). We rst projected metaclusters onto tissue sections, separating out non-immune metaclusters (hepatocytes, cholangiocytes, endothelial cells) and immune metaclusters (CD4 + T-cells, CD8 + T-cells, macrophages, monocytes, neutrophils, B cells and plasma cells, Fig. 1b). We then quanti ed the number of cells within each metacluster and evaluated scaled marker expression of lineage markers (Fig. 1c). Hepatocytes were the most common non-immune cell type, representing 62.6% of all cells identi ed, while macrophages were the most common immune cell type, representing 9.9% of all cells identi ed. Rare populations were also identi ed, including cholangiocytes (1.8% of all cells) and B cells (1% of all cells). Next, t-SNE was used to visualize differences in cell metaclusters between clinical groups (Fig. 1d) Fig. 2i-j). While MHC II molecules are constitutively expressed on human cholangiocytes, the in ammatory state of several diseases including primary biliary cirrhosis, primary sclerosing cholangitis, graft versus host disease, and even liver TCMR has been associated with MHC II overexpression on cholangiocytes, which may function as APCs in the liver 35,36 . Within immune metaclusters, there was an increase in CD8 + T-cells between NR and TCMR as well as NR and CR (p < 0.01) with a subtle increase in monocytes from NR to TCMR (p < 0.01, Fig. 1e). Despite macrophages being the most common immune metacluster, which is consistent with their pivotal role in regulating liver immune function, there were no differences in abundance between clinical groups (Fig. 1e, Extended Data Data Fig. 3b-c). Compared to NR, TCMR had a greater proportion of CD3 high CD4 + T-cells, naïve CD4 + Tcells, and activated CD4 + T-cells, which is consistent with acute alloreactivity (Fig. 2f). While their overall frequency was rare, we observed a concomitant increase in regulatory cell types, including HLADR + Tregs and PD1 + CD4 + T-cells, in the TCMR group when compared to NR, suggesting that their expansion counters effector alloreactive T-cell activity (Fig. 2f). We also determined that the CD3 high CD4 + T-cell subset represented most of the CD4 + T-cells in CR, with a signi cant decrease in resident memory T-cells and higher proportion of activated T-cells when compared to NR (Fig. 2f). Unlike TCMR, there was no expansion of the regulatory HLADR + Treg or PD1 + CD4 + T-cell populations in CR. To understand the lineage maturation trajectory of CD4 + T-cells in the alloimmune microenvironment, pseudotime reconstruction was performed (Fig. 2g) 32 . This provides further evidence that NR is primarily associated with CD4 + resident memory T-cells and suggests that CD4 + T-cell subpopulations increased during TCMR and CR originate and proliferate from circulating CD4 + T-cells (Fig. 2g). These data also suggest that the expanded Treg and PD1 + CD4 + T-cells observed in TCMR represent late-stage effector cells unique to this phase of alloimmunity.
Macrophages: Among the immune meta-clusters, macrophages were the most abundant cell type (Extended Data Fig. 1d) in all clinical groups, which highlights their key role in liver homeostasis, disease, and injury processes 37,38 . Indeed, macrophages can participate in robust in ltration of the AME during severe rejection episodes; however, their role has rarely been investigated in TCMR and CR in clinical LT 24,41 . We have previously shown that CR is characterized by a discrete macrophage phenotype absent in NR 24 . Thus, to obtain a detailed representation of the macrophages complexity and heterogenous activity in LT, we rst divided macrophages M1 and M2 based on their expression of CD163 (Extended Data Fig. 4a) 42 . The overall distribution of M1 and M2 did not differ among NR, TCMR, and CR, nor did the ratio of M2:M1 macrophages (Extended Data Fig. 4b- (Ki67 + CD68 + CD163 Hi ) and 'HLADR + M2' (HLADR + CD68 + CD163 Hi ), Fig. 4a-c, Extended Data Fig. 4e). Consistent with the activation of an in ammatory process, a greater percentage of proliferating M1 macrophages was observed in TCMR compared to NR and CR (Fig. 4f). We found one M1 and one M2 macrophage subset each expressing CD16. Both NR and TCMR exhibited a greater percentage of CD16 + M1 macrophages when compared with CR. The CD16 + M2 macrophage subcluster was most abundant in NR and appeared to become progressively depleted from TCMR to CR. These cells might represent a population of regulatory and anti-in ammatory macrophages (M2b), capable of IL-10 secretion ( Fig. 4f) 43,44 . A subpopulation of HLADR + M2 macrophages showed the opposing pattern to CD16 + M2 cells and was more abundant in both TCMR and CR than NR (Fig. 4f). These HLADR + M2 macrophages might represent a different activation state compared with the generic M2 macrophage subpopulations or suggest a unique specialization of those cells such as for antigen presentation.
Classical monocytes represented the most abundant subset across all clinical groups (Extended Data Fig. 5e-f), and the comparison of percentage across the three allograft states showed that intermediate and non-classical monocytes comprised a greater proportion of the monocyte metacluster in NR compared with TCMR and CR (Extended Data Fig. 5g).
B-cells and Plasma cells: B cells represented the smallest metacluster in the overall dataset (4,881 or 1% of all cells identi ed, Fig. 1c). Comparison of the three B cell subpopulations identi ed did not highlight any difference across different alloimmune states (Extended data Fig. 6). Finally, the small fraction of plasma cells identi ed, approximately 1.2% of all cells contained in the dataset, showed a higher proportion of those cells in TCMR than CR (Fig. 1e).
Spatial interaction and multicellular functional motifs de ne liver allograft rejection pathology Next, we examined the spatial data layer from our single cell proteomic IMC atlas to assess pairwise relationships between immune subpopulations within each clinical group by applying neighborhood and correlation analysis to characterize the statistical probabilities of cell-cell interactions (Fig. 5). Overall, a greater number of interactions, either via avoidance or attraction, were observed in TCMR when compared to NR and CR ( Fig. 5a-b, Extended Data Fig. 7a-b). The AME in TCMR was characterized by CD3 + CD8 + Tcells showing attraction to APCs including proliferating and M1 macrophages, classical monocytes, HLADR + M2 macrophages and B-cells, as well as CD3 + CD4 + T-cells supporting the idea of complex multicellular interactions characterizing this pro-in ammatory state (Fig. 5a). The neighborhood analysis revealed the presence of exhausted T-cells (PD1 + CD4 + , PD1 + CD28 + CD8 + , and PD1 + CD8 + T-cells) and Tregs in the vicinity of effector T-cells, which established a greater number of positive interactions when compared to NR, suggesting that a close crosstalk between those two ends of the spectrum T-cell phenotypes occurs in TCMR ( Fig. 5a-b, Extended Data Fig. 7a). Conversely, resident memory CD4 + T-cells showed no contact or avoidance with exhausted phenotypes in NR and TCMR respectively ( Fig. 5a-b, Extended Data Fig. 7a). In CR, HLADR + M2 macrophages surrounded HLADR + hepatocytes and M1 macrophages, while reciprocal strong interaction between CD16 + M2 macrophages and M2 macrophages were observed, likely representing a niche in which further differentiation of M2 macrophages occurs (Extended Data Fig. 7b).
Due to the relationship of the RAI score for TCMR to endothelial in ammation, we evaluated the distributions of distance to endothelial cells and each immune subpopulation across clinical groups (Extended Data Fig. 7c). Most pro-in ammatory subpopulations including CD3 + CD4 + and CD3 + CD8 + T-cells, proliferating and cytotoxic CD8 + T-cells, and classical monocytes resided near CD31 + endothelial cells, while resident memory CD4 + T-cells, CD16 + M1 macrophages, and CD16 + M2 macrophages were distributed throughout the tissue (Extended Data Fig. 7c).
To evaluate higher order spatial motifs as potential functional units associated with liver allograft pathology, we performed cellular neighborhood analysis, where patterns of higher order cellular structures (with 10 nearest neighbors) were clustered into nine novel cellular neighborhoods (CNs). These were labeled according to the cell types in each cluster as shown in the heatmap in Fig. 5c: hepatocytes, vasculature, granulocyte enriched, activated macrophages, CD8 enriched, CD16 + T-helper enriched, Thelper enriched, B-cell and monocyte enriched, and bile ducts (Fig. 5c-d). We then visualized and compared the proportions of these CNs across clinical groups (Fig. 5d-e, Extended Data Fig. 7d). For the non-immune predominant CNs (hepatocyte, vascular, and bile duct) there were few differences between clinical groups except for a slightly smaller proportion of the hepatocyte CN in TCMR (likely a direct consequence of the increase of immune cell enriched CNs in TCMR, Extended Data Fig. 7d). The CD8 enriched CN (which included HLADR + CD4 + Treg and both PD1 + CD8 + T-cell subpopulations) was expanded in TCMR when compared with NR and CR (p < 0.001 and p < 0.01 respectively). Additionally, the 'B cell and monocyte' CN was most abundant in TCMR and was also expanded in CR compared to NR (p < 0.001). Conversely, the CD16 + T-helper enriched CN was more abundant in NR when compared with TCMR and CR (p < 0.001), suggesting that this CN could be a marker of allograft tissue homeostasis (Extended Data Fig. 7d).

Evaluation of exhaustion markers in liver biopsy tissue via RNA Sequencing
Our single cell atlas across the spectrum of rejection in LT identi ed several unique cell types increased during TCMR, including Tregs, PD1 + CD4 + T-cells, PD1 + CD8 + T-cells, and HLADR + M2 macrophages.
To obtain further evidence for the functional status of these populations, including whether the identi ed PD1 + T-cells represent a terminally differentiated, activated CD4 + or CD8 + T-cell versus an exhausted effector T-cell population, we performed bulk transcriptomic analysis using the nCounter platform 45 . Liver core biopsies from 8 samples, comprising the most representative 4 TCMR and 4 NR cases were used for analysis. We selected 23 genes de ning T-cell phenotypes including helper function, exhaustion, and cytotoxic activity (Fig. 6a). By comparing DE genes between NR and TCMR, we identi ed an overall upregulation of genes typically associated with cytotoxic activity and as well as upregulation of PDCD1 (Programmed cell death 1 or PD1) gene expression in TCMR samples, which is consistent with the higher percentage of PD1 + CD4 + and CD8 + T-cell identi ed in our IMC dataset (Fig. 6a). We also observed increased DE of PDCD1LG1 (Programmed cell death-ligand 1 or PD-L1), PDCD1LG2 (Programmed cell death-ligand 2 or PD-L2), CTLA4 (Cytotoxic T-Lymphocyte Associated Protein 4), Lag3 (Lymphocyte activating 3), and CD160 (CD160 antigen) genes in TCMR compared to NR, con rming upregulation of both ligands for PD1 and T-cell exhaustion markers (Fig. 6a). We selected 18 genes which differentiate the diverse macrophage polarization in M1, M2a, M2b and M2c (Fig. 6b). Genes for both pro-in ammatory cytokines such as CXCL9 (C-X-C motif chemokine ligand 9) and CXCL10 (C-X-C motif chemokine ligand 10), mainly expressed by M1 macrophages, and anti-in ammatory cytokines including IL-10 (Interleukin 10), CCL22 (C-C motif Chemokine Ligand 22), CCL24 (C-C motif Chemokine Ligand 24) mostly associated with M2 macrophage polarization, were upregulated in TCMR when compared to NR.
We have not yet identi ed a reliable NK marker for IMC in liver tissue, so we also evaluated markers of NK cells including IL21R (Interleukin 21 receptor), XCL1 (X-C motif Chemokine Ligand 1) and NCR1 (natural cytotoxicity triggering receptor 1), which were upregulated in TCMR samples (Fig. 6b). We also identi ed an overall upregulation of several genes associated with neutrophils, B cells, mast cells, and dendritic cells (Extended data Fig. 8a). We validated these ndings using a publicly available database including of 6 cases of biopsy-proven TCMR in clinical LT, demonstrating upregulation of the same DE genes from our analysis (Extended data Fig. 8b-c) 34 .
Finally, we evaluated the intragraft T-cell receptor (TCR) repertoire to determine whether alloreactive Tcells during TCMR exhibited oligoclonal expansion via bulk analysis of TCR diversity using the 129 gene TCR Diversity nCounter assay. By using the average expression of TCR genes in NR samples as baseline, we analyzed which TCR variable region genes were upregulated or downregulated in TCMR compared to NR. We found an upregulation in the majority of TCR alpha (TRAV), TCR beta (TRBV), and TCR delta (TRDV) and TCR gamma (TRGV) variable region genes (Fig. 6c). These data suggest that TCMR is associated with an increase in the number of T-cells clones rather than expansion of few pathogenic alloreactive clones, which is similar to prior reports 46 .

Cell composition is highly predictive of rejection in liver allografts
The analysis of the AME in NR, TCMR and CR identi ed 41 potential features of which 14 immune and 5 non-immune differed in patients who developed TCMR from NR, 9 immune and 4 non-immune features distinguished TCMR from CR, and 10 immune and 1 non-immune features separated out NR from CR, thereby highlighting complex network of different cell phenotypes speci c for those three AMEs. To determine whether these immune phenotypes could predict patient outcomes, the logistic LASSO regression algorithm was applied to the entire dataset. To improve model robustness, 5-fold crossvalidation was used to determine model parameters, and a 5,000 iterative bootstrapping technique was used to perform feature selection by determining feature importance based on frequency. To ensure that predicted results and model performance are derived from patients not included when training the model, a 1,000 iterative random training/validation split was implemented; performance metrics were averaged and the median prediction of each patient was used for further model evaluation. Examination of the most important features, which present a frequency ≥ 50%, revealed that 8 cell subpopulations contributed the most in generating a model that can accurately differentiate TCMR vs NR (mean ± SD; accuracy 0.89 ± 0.07 and mean area under the curve (AUC) 0.96 ± 0.04) demonstrating a high correlation between median prediction and actual clinical outcome (Spearman correlation coe cient R = 0.77, p = 7.206x10 − 10 ) (Fig. 7a,c and Extended Data Fig. 9a). The highest-ranking immune phenotype was resident memory CD4 + T-cells as a predictor of NR, corresponding to pairwise analysis demonstrating that this immune subset was strongly associated with NR (Fig. 2). In addition, intermediate monocytes, cholangiocytes, and CD16 + M2 macrophages were predictors of NR, whereas PD1 + CD4 + T-cells, HLADR + M2 macrophages, non-classical monocytes, and proliferating hepatocytes were positive predictors of TCMR (Fig. 7b). Application of this modelling approach to differentiate TCMR from CR resulted in 9 highly ranked features that can accurately distinguish these two alloimmune states, with a high sensitivity, speci city, and accuracy, and a mean AUC of 0.96 ± 0.04 (Spearman correlation coe cient between prediction and actual outcome of 0.82, p = 1.782x10 − 9 ) (Fig. 7d,f). Among these features, proliferating and CD16 + M1 macrophages, proliferating and PD1 + CD8 + T-cells, plasma cells, CD3 + and CD16 + CD4 + T-cells predicted TCMR while the CD3 + CD8 + T-cell phenotype was a predictor of CR ( Fig. 7e and Extended Data Fig. 9b). Modelling for differentiating NR from CR was also highly sensitive and speci c, however the sample size for this comparison (38 samples total) may be too small to generate a conclusion (Extended Data Fig. 9). Overall, these results indicate that rare but speci c cell subpopulations identi ed in the present study can potentially harbor high diagnostic value in biopsies obtained across the spectrum of rejection in clinical LT.

DISCUSSION
Our study provides a comprehensive single cell, spatially resolved analysis of the AME in clinical LT, revealing the complexity of alloimmunity in solid organ transplant recipients. Unlike the cancer tumor microenvironment, which has remarkable phenotypic variability between patients and even within the same specimen, our analysis con rms central tenants of transplant immunology, namely that the pathologic features within the AME are similar across individuals, despite differences in patient demographics, underlying etiology of liver disease, features of the donor organ, and timing of rejection episodes. Thus, study of the AME offers an ideal application and proof-of-concept for further development of spatial proteomics immunologic analyses using archival biopsy specimens. Further, exploration of discrete immune subpopulations within the AME of core needle liver biopsies has identi ed immune subsets with exhaustion phenotypes that are enriched in TCMR and largely absent in CR, providing new insights into the mechanistic underpinnings and evolution of liver allograft rejection.
Finally, these data were harnessed to create a predictive model of TCMR and CR using a subset of cell types, which offers potential for clinical use to diagnose rejection states more accurately when compared to the current standard of care that relies on subjective pathologic evaluation using the RAI.
Single cell analysis of the AME has uncovered substantial complexity in allograft rejection, involving at least 32 distinct immune subpopulations. Nearly all prior studies of LT rejection have focused on one or few immune cell types 47,48 . In clinical TCMR, our data demonstrate that diverse cell populations contribute to the underlying pathophysiology. It is increasingly recognized that spatial context is important to completely describe disease phenotypes and that these multiplexed spatial techniques will have critical clinical implications 49 . A recent study in LT recipients examined immune cell type pairs at high resolution to evaluate immune synapse formation and used these data to predict success of immunosuppression withdrawal 11 . These data, taken together with our results, suggest that further characterization of important features of the AME will provide valuable insights into predicting clinical outcomes with greater precision than is currently possible.
Arising from the complex microenvironment was a central theme of immune exhaustion. Our study design captures clinical specimens prior to initiation of rejection treatment, the mainstay of which are calcineurin inhibitors (CNI), which are designed to indiscriminately prevent T-cell proliferation and come with potentially severe side effects. Fortunately, the immune system has physiologic mechanisms to dampen this immune response through PD-1 ligand, a molecule that when knocked out in mice results in auto-immunity, and is also important in chronic in ammatory states 15,50,51 . Our results show PD1 + CD4 + T-cells and PD1 + CD8 + T-cells are expanded in TCMR and spatially interacting with other CD4 + and CD8 + subpopulations, as well as M2 + HLADR + macrophages. Predictive modelling classi ed PD1 + CD4 + T-cells as a feature distinguishing TCMR from NR, while PD1 + CD8 + T-cells were identi ed as a feature distinguishing TCMR from CR. Taken together, these data indicate that immune exhaustion may be a key mediator of TCMR in the AME.
Dysregulated exhaustion states are increasingly recognized as a pathway cancer cells manipulate to mediate immune escape, leading to development of immune checkpoint inhibitor therapies to enhance anti-tumor adaptive immune responses. Our data suggest the opposite therapeutic approach should be explored to counteract pro-in ammatory responses during acute TCMR, via augmentation of physiologic exhaustion. Indeed, prior work using PD-1 agonists in other pro-in ammatory states, including neutrophilic asthma, has demonstrated that therapies designed to promote T-cell exhaustion can mitigate in ammation 16 .
For the past 60 years, pathologic detection of allograft rejection has been conducted using Hematoxylin and Eosin staining. The Banff RAI is then used to characterize rejection by evaluation of portal and/or perivenular in ammatory (immune) in ltrates 52 . Our study suggests that spatial relationships between immune cells and cholangiocytes or endothelial cells may be less important mechanistically. Rather, investigation into the presence or absence of certain immune subpopulations may better inform important considerations of TCMR such as steroid-resistant disease, a disease that often progresses to CR and drives late graft failure. Predictive modelling combined with IMC may bring value in this regard. A recent study used IMC combined with deep learning to predict lung adenocarcinoma progression and patient survival post-surgery with high accuracy 53 . Harnessing multiplexed data together with emerging arti cial intellegnce tools such as deep learning may have profound diagnostic and prognostic value, both in clinical practice and to monitoring responses to treatment in clinical trials. A major bene t of IMC over transcriptomic platforms is that IMC can be employed on archival FFPE samples from a small core biopsy without concern for RNA degradation. Furthermore, proteomics-based platforms like IMC more accurately re ect single cell phenotypes given that RNA is not always linearly correlated with protein translation 54 . In our study, clinically diagnostic areas selected from core biopsies resulted in over 5000 cells per ROI, which is comparable to cell counts obtained from scRNA seq experiments in human liver, without bias caused by pre-analytical variables, including cellular damage and loss, inherent to tissue dissociation and processing techniques necessary to create a single cell suspension 55 .
Our analysis is limited by the inability to conduct complementary transcriptomics or cell-culture based assays as our study was performed on a retrospective set of tissue formalin-xed and para n embedded specimens collected during routine clinical care. Furthermore, there is the possibility of cell classi cation errors within our IMC dataset, particularly within the hepatocyte metacluster that was classi ed based on exclusion of other cell types. To minimize this risk, we conducted an extensive review all annotated tissue specimens alongside raw marker expressions to ensure that tissue labelling was optimized. Due to the unexpected importance of exhaustion phenotypes in our results, our IMC dataset lacked PD-L1 and other exhaustion markers, and these will be incorporated for future studies. Finally, our patient sample size was somewhat limiting for our results, particularly with predictive modelling in CR; however, our analyses support our study being adequately powered for the evaluation of important cell subpopulations and in the modelling of TCMR.
Herein, we provide a novel, detailed, and spatially resolved atlas of clinical liver allograft rejection. Highly        In depth molecular characterization of tissue using bulk RNA sequencing (nCounter) panel con rms a mix of phenotypes enriched in T-cells and macrophages in TCMR. a. T-cell populations predicted from a set of 4 NR tissue specimens used in our dataset and 4 TCMR tissue specimens. Genes correspond to generic Tcells, Th1, cytotoxic and exhausted phenotypes. Activated, cytotoxic as well as exhausted T-cell genes showed a greater expression in TCMR compared to NR. b. Heatmap of scaled expression values of macrophage and NK-related genes including M1, M2a, M2b, M2c phenotypes along with NK-associated genes. Genes belonging to both M1 and M2 polarized macrophages showed a greater expression in TCMR than NR; similarly, NK-associated genes were upregulated in TCMR. c. Histograms showing percentages of upregulated TCR alpha (TRAV) genes, TCR beta (TRBV), and TCR gamma delta (TRGV and TRDV) genes respectively across samples. The increased expression of different TCR alpha, beta and gamma-delta associated with TCMR, is in agreement with expansion in the number of T-cell clones.

Figure 7
Indenti cation of cellular features of liver allograft biopsies are highly predictive of discriminating TCMR