Abstract Elucidating the temporal process of immune remodeling under immunosuppressive treatment after liver transplantation (LT) is critical for precise clinical management strategies. Here, we performed a single-cell multi-omics analysis of peripheral blood mononuclear cells (PBMCs) collected from LT patients (with and without acute cellular rejection [ACR]) at 13 time points. Validation was performed in two independent cohorts with additional LT patients and healthy controls. Our study revealed a four-phase recovery process after LT and delineated changes in immune cell composition, expression programs, and interactions along this process. The intensity of the immune response differs between the ACR and non-ACR patients. Notably, the newly identified inflamed NK cells, CD14^+RNASE2^+ monocytes, and FOS-expressing monocytes emerged as predictive indicators of ACR. This study illuminates the longitudinal evolution of the immune cell landscape under tacrolimus-based immunosuppressive treatment during LT recovery, providing a four-phase framework that aids the clinical management of LT patients. Graphical abstract graphic file with name fx1.jpg [59]Open in a new tab Public summary * • A dynamic immune landscape after liver transplantation unveils a four-phase recovery process, facilitating the precise management of liver transplantation patients. * • T cell activity deficiency does not completely suppress the allogeneic cell rejection. * • The elevated proportion of inflamed NK cell subset is a potential driver for graft rejection in liver transplantation patients. * • Monocytes serve as the major signal source to initialize immune responses after liver transplantation. * • CD14^+RNASE2^+ monocytes and FOS-expressing monocytes are potential early predictive indicators for graft rejection. Introduction Liver transplantation (LT) is the ultimate therapy for patients with end-stage liver disease. Tacrolimus-based immunosuppressive regimens have been widely used in the clinic and have achieved great success in LT,[60]^1 yet graft rejection and severe drug side effects still occur in some LT patients. Consequently, therapeutic drug monitoring (TDM) is essential for individualizing and optimizing tacrolimus therapy. However, achieving precision in immunosuppression management purely based on TDM is a formidable challenge because the change in plasma concentration of tacrolimus only reflects its pharmacokinetics in vivo and cannot directly reflect the immunosuppressive status of LT patients. A comprehensive understanding of immune status dynamics during the LT recovery trajectory and between patients with different clinical outcomes would facilitate the development of precise clinical management strategies (e.g., preventing graft rejection and/or reducing drug toxicity for the right patients at the optimal time). Recently, single-cell mRNA sequencing (scRNA-seq) and cytometry by time of flight (CyTOF) have brought many new insights to the immune status of LT.[61]^2^,[62]^3^,[63]^4^,[64]^5 For instance, Huang et al. mapped a single-cell atlas of the transplanted liver, profiling the physiological process of graft reconstruction and the interaction between the graft and the recipient.[65]^3 In a follow-up study, they used CyTOF to explore the temporal evolution of hepatic macrophage populations.[66]^5 However, these studies were conducted in animal LT models and may not reflect the dynamics in LT patients. Several single-cell studies of intrahepatic immune cells in LT patients have been published,[67]^5^,[68]^6^,[69]^7 but data were only gathered for limited time points. Other studies used CyTOF or T cell receptor repertoire analyses to address the alterations in immune cells among LT patients but with limited cell markers.[70]^2^,[71]^8 Longitudinal, unbiased fate mapping of the global immune cell populations across the LT recovery trajectory in humans remains unexplored. The analysis of peripheral blood mononuclear cells (PBMCs) offers a minimally invasive, rapid, and convenient approach to dynamically reflect systemic immune states, which could remedy the weakness of TDM and allow for a more accurate evaluation of the immune status of patients. In this study, we profiled the immune cell dynamics by scRNA-seq, bulk mRNA-seq, CyTOF, and flow cytometry in PBMCs from LT patients with and without acute cellular rejection (ACR) at multiple time points up to 365 days after LT. We charted a four-phase response during LT recovery, characterized cell heterogeneity between patients with different clinical outcomes, and further validated the presence and clinical relevance of newly discovered immune cell subsets in additional cohorts. Results Study design To assess the dynamics of immune responses after LT, we collected PBMCs from six carefully selected LT patients. Serial blood samples were obtained on days 0, 1, 3, 5, 7, 11, 15, 19, 23, 60, 90, 180, and 365. Within this cohort, two patients experienced ACR, and additional blood samples were collected at the time of rejection (days 31 and 60 after LT, respectively) and postrejection (days 67 and 98, respectively). PBMCs from five healthy donors were collected as controls (HCs). As a result, 69 samples were subjected to scRNA-seq, with a subset also examined via CyTOF and cytokine assays ([72]Figures 1A and 1B). Following this discovery cohort (cohort 1), longitudinal PBMCs from a subset of 37 non-ACR LT patients and 4 HCs (cohort 2) were collected and subjected to bulk RNA-seq, with the relevant clinical data recorded. Longitudinal PBMCs from an additional 24 LT patients (14 without ACR and 10 with ACR) and 7 HCs (cohort 3) were collected and analyzed using flow cytometry and ELISA ([73]Figure 1A). The demographics and clinical features of the patients are detailed in [74]Table S1. Figure 1. [75]Figure 1 [76]Open in a new tab Study design (A) Overview of the sample collection and analysis workflow. (Created with [77]BioRender.com.). (B) Schematic diagram of time points for blood sampling and analysis (top). Line graphs showing longitudinal changes in clinical parameters (mean) in the ACR and non-ACR groups (bottom). MARS-seq, massively parallel RNA single-cell sequencing. Comprehensive analysis of single-cell data, bulk transcriptomics, and clinical parameters reveals a four-phase recovery process after LT During the LT recovery process, the patients’ immune responses exhibited large interindividual differences, yet most ultimately returned to a stable state. The non-ACR patients with long-term follow-up successfully recovered from LT, enabling us to observe the entire recovery course. Initially, we scrutinized the time-based gene expression patterns in cohort 1. Specifically, we compared the combined transcriptomes of immune cells at individual time points to those at other time points to pinpoint time feature genes (TFGs) ([78]Figure S1A; [79]Table S2). The number of TFGs sharply increased on day 1 after LT and then decreased and stabilized at a relatively low level between days 19 and 90. From days 90 to 365, the number of TFGs increased and was maintained at a relatively high level comparable to that in HCs. Gene Ontology analysis[80]^9 of upregulated genes at each time point showed that biological processes and pathways enriched at each time point demonstrated a phased transition ([81]Figure 2A). Notably, pathways enriched on days 90, 180, and 365 demonstrated significant similarity with those in HCs (hypergeometric distribution test, p < 0.0001). These observations suggested that the immune components underwent stepwise remodeling and eventually reached a relatively steady state at approximately day 90 after LT. Figure 2. [82]Figure 2 [83]Open in a new tab Definition of a four-phase recovery process in LT (A) Dot plot illustrating the top 10 selected enriched terms of TFGs at each time point in the non-ACR group of cohort 1. (B) Unsupervised hierarchical clustering analyses (ward.D2 method) of immune-related gene expression profiles (left) and GSVA results (right) in the non-ACR group of cohort 1. (C) Principal-component analysis (PCA) of non-ACR samples from cohort 1 based on the immune-related gene expression. (D) PLSDA plot of samples from cohort 2 based on the immune-related gene expression, with the pairwise Adonis test results on the right. (E) Boxplots depicting the temporal changes in clinical parameters across different phases of LT in cohort 2. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, and ∗∗∗∗p < 0.0001 in the Wilcoxon test. ns, not significant. (F) Results of time-series clustering analysis by Mfuzz (left) based on the non-ACR group of cohort 1. Heatmap of longitudinal DEGs across LT phases with selected enriched pathways labeled at right. Furthermore, we applied unsupervised hierarchical clustering analysis to the immune-related gene expression, the gene set variation analysis (GSVA) results, and the immune cell populations identified in the scRNA-seq data ([84]Figures 2B and [85]S1B). Combining the results, we successfully identified 4 sequential but separable phases as follows: the first 7 days as phase 1, from days 7 to 19 as phase 2, from days 19 to 90 as phase 3, and from days 90 to 365 after LT as phase 4. The immune-related gene expression changes across these phases presented a progressive trend from phases 1 to 4 and gradually approached a healthy state ([86]Figure 2C). To validate the generalizability of this four-phase model, we used bulk RNA-seq on 71 samples from 14 LT patients at multiple time points and 4 HCs in cohort 2. Partial least squares discriminant analysis (PLSDA) of 64 post-LT samples highlighted a clear distinction from phases 1 to 3 (pairwise Adonis, p < 0.01). In contrast, phase 4 and HCs showed negligible differences (pairwise Adonis, p = 0.111) ([87]Figure 2D). Correspondingly, the number of differentially expressed genes (DEGs) between HCs and LT patients decreased from phases 1 to 4 in both the scRNA-seq (from n = 1,073 to n = 459) and bulk RNA-seq (from n = 5,297 to n = 1,216) data ([88]Figures S1C and S1D). Finally, to evaluate the clinical implications of these four-phase recovery processes, we retrospectively analyzed the shift in liver function parameters during these distinct phases in cohort 2 (n = 37), including alanine transaminase (ALT), aspartate transaminase (AST), total bilirubin (TB), direct bilirubin (DB), prealbumin (PA), and prothrombin time (PT) ([89]Figure 2E). These parameters showed constant and phased changes matching the transcriptome-defined phases. Using Mfuzz[90]^10 for the scRNA-seq data, we identified nine distinct transcriptional patterns associated with major immune events across the four recovery phases ([91]Figure 2F). Phase 1 (cluster 6) showed a surge in transcripts for cellular oxidant detoxification (e.g., GPX1, NCF2) and Toll-like receptor binding (e.g., TLR8, S100A12), indicating an early response to the traumatic injury-induced stress. Phase 2 (cluster 3) was characterized by an increased expression of type I interferon (IFN) stimulated genes (ISGs) (e.g., ISG15 and MX1) and genes involved in ubiquitination (e.g., VCP, WSB1). Since the ISGs are important for antigen presentation[92]^11 and facilitation of the alloresponse,[93]^12^,[94]^13 and ubiquitination is required for eukaryotic cells to recover from oxidative stress by decomposing stress granules,[95]^14^,[96]^15^,[97]^16 there could be an ongoing alloresponse and posttraumatic stress reparative response in immune cells during phase 2. In phases 2 and 3 (cluster 2), there was an increase in RNA splicing and chromatin remodeling-related transcripts (e.g., SF3A3, UPF1, CHD2), along with a peak in cell-cycle-related transcripts (e.g., CDK4, TLK1) in phase 3 (cluster 8). These results indicated that immune cells in phase 3 continuously recovered from stress-induced transcriptional inhibition[98]^17^,[99]^18 of phase 1 and enhanced cell proliferation for homeostasis remodeling. By phase 4 (clusters 7 and 9), immune response genes (e.g., IL-2RB, GZMB, IL-6ST) aligned closely with HCs, despite some persistent disparities (e.g., SELL, NOTCH2) (cluster 1). In addition, differential expression and pathway enrichment analysis across 71 bulk RNA-seq samples reaffirmed these phase-specific immune signatures ([100]Figure S1E). Also, we applied weighted gene coexpression network analysis to the bulk RNA-seq data and identified 7 gene modules (M1–M7) with distinct expression patterns throughout the LT recovery phases ([101]Figure S1F). Consistent with the findings from the scRNA-seq data, we found that M2 and M3, which reflected T cell and natural killer (NK) cell responses, exhibited the lowest levels in phase 1 and then progressively increased in later phases. M6, related to the response to IFN, continued to increase during phases 2 and 3. M7, associated with the innate immune response, was enriched in phase 1 after LT. Taken together, these results indicated that each LT recovery phase may have a unique immune signature. Global changes in immune cells across the four-phase recovery process Although transcriptomic analyses could objectively reflect the temporal changes in major immune events, many cell-type-specific changes were still masked. We next examined temporal alterations in cell composition over the four-phase recovery process based on the scRNA-seq data to address this ([102]Figure S2A). After the removal of neutrophils and erythrocytes, 7 major cell clusters were identified ([103]Figures 3A and [104]S2B), including monocytes, CD4^+ T cells, CD8^+ T cells, B cells, NK cells, megakaryocytes, and dendritic cells. As expected, these cell clusters showed different enrichment across the LT recovery trajectory ([105]Figure 3B), with the compositions of immune cells gradually altering from phases 1 to 4 and approaching a healthy state ([106]Figure 3C). Figure 3. [107]Figure 3 [108]Open in a new tab Global changes in immune cells across the four-phase recovery process (A and B) Uniform manifold approximation and projection (UMAP) plots of all merged patient samples colored by cell type (A) and time points (B) in the non-ACR group. (C) PCA result of the distribution of immune cells in the non-ACR group. (D) Heatmap showing the module score of immune cells from the non-ACR group in the distinct phase. (E) Line graphs showing longitudinal changes in cell proportions identified in the scRNA-seq data (mean ± SEM). Comparisons were made between the ACR and non-ACR groups and between phase 4 and HCs. Wilcoxon test. The sample sizes (n) for each group were as follows: non-ACR group: baseline (n = 4), phase 1 (n = 11), phase 2 (n = 12), phase 3 (n = 9), phase 4 (n = 7); ACR group: baseline (n = 2), phase 1 (n = 5), phase 2 (n = 6), phase 3 (n = 4), during rejection (n = 2), postrejection (n = 2); and HCs (n = 5). (F and G) Correlation heatmaps between cell proportions and clinical parameters or cytokines (F) and between cell proportions and gene coexpression modules (G). Spearman’s correlation. (H) t-distributed stochastic neighbor embedding plots of aggregated samples and distinct groups colored by annotated cell types derived from CyTOF data. ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001; ∗∗∗∗p < 0.0001. To determine the cell types that dominate the recovery phases of the non-ACR group, we examined the phase-feature cell types (PFCs), defined as immune cells with higher module scores of phase-feature genes from the scRNA-seq data ([109]Figure 3D; [110]Table S3). In phase 1, the innate immune cells (monocytes and megakaryocytes) were PFCs. For phase 2, PFCs included megakaryocytes, monocytes, CD8^+ T cells, and NK cells. The PFCs in phase 3 were NK cells, CD8^+ T cells, and monocytes. The PFCs of phase 4 were CD8^+ T cells, CD4^+ T cells, and NK cells, aligned with those observed in HCs. Overall, the dominant immune cells during the normal recovery trajectory tended to transition from innate immune cells to adaptive immune cells, which was also confirmed in both bulk RNA-seq data ([111]Figure S2C) and cellular composition ([112]Figure 3E). Correlation analyses of the clinical parameters, cytokines, and the gene modules with the cell proportions highlighted monocytes, T cells, and NK cells as central to LT recovery, which were significantly correlated with one or more of the clinical parameters (e.g., ALT, AST, TB, DB) and gene modules ([113]Figures 3F and 3G). Among these, monocyte proportions showed correlations with the levels of inflammatory response parameters (e.g., IL-6, IL-8, IL-13) and tissue repair parameters (e.g., SCGF and HGF). T cells from LT patients display a weak immune response during LT recovery To elucidate the potential contributions of these immune cells to the LT recovery trajectory and clinical outcomes, we conducted a longitudinal comparison of immune cell subpopulations and their functional changes between the non-ACR and ACR groups and validated our findings with a 42-parameter CyTOF assay on cohort 1 at crucial time points ([114]Figures 3H and [115]S2D). As primary mediators of allogeneic immune responses,[116]^19 T cells present complex dynamics post-LT and are affected by immunosuppressants.[117]^1 We identified 7 T cell subsets, namely naive CD4^+ T cells, central memory CD4^+ T cells, T regulatory cells (Tregs), naive CD8^+ T cells, anergic CD8^+ T cells, effector memory CD8^+ T cells, and terminal effector CD8^+ T cells (CTLs) ([118]Figures 4A–4C). Figure 4. [119]Figure 4 [120]Open in a new tab Dynamic changes in subpopulations of T cells after LT (A) UMAP plot showing T cell subpopulations of all of the samples, with the percentages of each cell type shown. (B) Dot plot showing scaled expression levels of cell-type-specific genes in the non-ACR, ACR, and HCs. (C) UMAP plots showing all merged T cells colored by groups and phases. (D) Bar plot showing the T cell subsets’ frequency change across the phases of LT in different groups. post, postrejection. (E) Line graphs showing longitudinal changes in the relative frequency of the T cell subpopulations in the non-ACR, ACR, and HCs. (F) Dot plot showing scaled expression levels of function-specific genes of all T cells in the non-ACR, ACR, and HCs. (G) Line graphs demonstrating temporal changes of the relative function-specific gene set scores in T cells. (H) Line graph demonstrating temporal changes of the cytotoxicity score of CTLs in the non-ACR, the ACR, and HCs. (I) Dot plot showing scaled expression levels of cytotoxicity-related genes of CTLs in the non-ACR, the ACR, and HCs. In (E), (G), and (H), data are shown as mean ± SEM, and comparisons were made between ACR and non-ACR groups and between phase 4 and HCs. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ∗∗∗∗p < 0.0001. ns, not significant in the Wilcoxon test. We investigated the T cell dynamics in the non-ACR group ([121]Figures 4D and 4E). Consistent with other reports,[122]^2^,[123]^8 we found that the proportions of CD4^+ and CD8^+ T cell subsets exhibited a rapid decline in phase 1 compared to the baseline (except for anergic CD8^+ T cells and Tregs), indicating that T cells were suppressed in phase 1, probably due to the strong immunosuppressant treatment (tacrolimus + mycophenolate mofetil + basiliximab) and surgical stress[124]^20 early after LT. Interestingly, all of the diminished T cell subsets recovered and continuously increased from phases 2 to 4, particularly in the memory subsets, which could be partially ascribed to homeostatic proliferation induced by long-term immunosuppressive therapy[125]^21 and persistent antigen stimulation.[126]^22 A continuous increase in the CD8^+ T cell proportion after phase 1 was also observed in the CyTOF analysis ([127]Figure S3A). The proportions of Tregs and anergic CD8^+ T cells peaked in phase 2 and gradually decreased to a level similar to that of the baseline after that. This pattern may result from immunosuppressive therapy because tacrolimus has been reported to be associated with a decrease in the frequency of the circulating Tregs. In contrast, the frequency of effector T cells is not reported to be affected.[128]^23^,[129]^24 Although we expected significantly different T cell dynamics in the ACR group during the corresponding phases (phases 1–3), this was not observed for most T cell subsets ([130]Figures 4D and 4E). Surprisingly, more anergic CD8^+ T cells and naive CD4^+ T cells were observed in the ACR group. In line with this, T cells in the ACR group displayed higher levels of co-inhibitory molecules and lower levels of co-stimulatory molecules compared to those in the non-ACR group ([131]Figures 4F and 4G), a condition of T cell dysfunction initiation.[132]^25^,[133]^26 Taken together, these results indicated that T cells did not generate an intensive immune response to liver alloantigen stimulation during the LT recovery process in the ACR group. To understand why increased CTLs do not lead to graft rejection in non-ACR patients, we analyzed the cytotoxicity of CTLs ([134]Figures 4H and 4I), one of the major functions causing graft rejection.[135]^27^,[136]^28 As expected, during the LT recovery, the overall cytotoxicity function of CTLs in the non-ACR group remained lower than HCs and gradually decreased. In addition, the cytotoxicity of CTLs (e.g., PRF1, GNLY, GZMB) was significantly lower in the ACR group than in the non-ACR group. This confirmed the impaired T cell functionality in the ACR group. Taken together, these findings indicated that immune cell types other than T cells mediate rejection under tacrolimus-based immunosuppressive treatment after LT. NK cells from ACR patients exhibit a distinct activation of “inflamed” transcriptional programs Unlike T cells, NK cells can recognize and lyse target cells without the need for prior antigen-specific sensitization but via integrating signals from both killer cell inhibitory receptors (KIRs) and killer cell activatory receptors (KARs).[137]^29 Surgical stress[138]^30^,[139]^31^,[140]^32 and tacrolimus[141]^33^,[142]^34 have been reported to cause phenotypic and functional defects in human NK cells. To date, the function of NK cells in LT has been less studied, or even with controversial results.[143]^35 Here, we identified 2 NK cell clusters, including CD56^brightCD16^− NK cells, which may support adaptive immunity through the secretion of cytokines, and CD56^dimCD16^+ NK cells, serving as a cytotoxic population with elevated expression of the effector molecules CX3CR1, GNLY, GZMH, and PRF1 ([144]Figures 5A and 5B).[145]^36^,[146]^37 Figure 5. [147]Figure 5 [148]Open in a new tab Heterogeneity of NK cells across the LT recovery process (A) UMAP plot showing NK cell subpopulations of all of the samples, with the percentages of each cell type shown. (B) Dot plot showing scaled expression levels of NK cell-type-specific genes in non-ACR, ACR, and HCs. (C) Bar plot showing the frequency of NK cell subsets in different groups across the phases of LT. (D) Dot plot showing scaled expression levels of function-specific genes of CD56^dimCD16^+ NK cells in non-ACR, ACR, and HCs. Selected genes are shown. (E) Line graphs demonstrating temporal changes of the relative function-specific gene set scores in CD56^dimCD16^+ NK cells (mean ± SEM). Comparisons were made between ACR and non-ACR groups and between phase 4 and HCs. Wilcoxon test. (F) UMAP plots showing the group information of NK cells. (G) Dot plot showing average expression levels of the top 10 DEGs across non-ACR, ACR, and HCs. (H) UMAP plots depicting 6 subsets of CD56^dimCD16^+ NK cells (left), with cell trajectory (Monocle3) and group distribution (right). (I) Dot plot showing average expression levels of the top 10 cluster-specific genes in CD56^dimCD16^+ NK cells. (J) Boxplots showing the cell-type-specific gene set scores for each CD56^dimCD16^+ NK cell subset (mean ± SEM). (K) Bar plot showing the relative abundance of CD56^dimCD16^+ NK cell subclusters in different groups across phases of LT. (L) Representative flow cytometry analysis of inflamed NK cells from non-ACR and ACR patients. (M) Bar graphs showing the proportion of inflamed NK cells from non-ACR (n = 7) and ACR (n = 7) patients in baseline (left) and from non-ACR (n = 11) and ACR (n = 9) patients in phase 3 (right) (mean ± SD). Student’s t test. (N) Boxplots showing the percentage of inflamed NK cells treated with different concentrations of MePDN in the presence or absence of tacrolimus (n = 6, per group). In vitro experiment. Paired t test. ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001; ∗∗∗∗p < 0.0001. In non-ACR patients, the proportion of NK subsets decreased post-LT but rebounded following surgical stress recovery and tacrolimus adaptation from phase 2 ([149]Figure 5C). Similar results were observed in the CyTOF data ([150]Figure S4A). Given that CD56^dimCD16^+ NK cells were the major effector NK subset following LT, we next focused on their functional dynamics throughout the LT recovery process ([151]Figures 5D and 5E). Phase 1 was characterized by an increase in KIRs and a decrease in KARs, indicative of suppressed NK cell functions. From phases 2 to 4, KARs exhibited a gradual increase, whereas KIRs decreased in phase 2, followed by a subsequent increase. Correspondingly, the cytotoxicity (e.g., GNLY, GZMA) and proliferation of CD56^dimCD16^+ NK cells increased in phase 2 but gradually stabilized thereafter. Moreover, CyTOF analysis revealed that the proportion, proliferation (Ki67), cytotoxicity (NKp46), and activation (CD38) of CD16^+ NK cells were significantly suppressed on day 1 and restored after day 7, but the killing effect downregulator (TIGIT) of NK cells[152]^38 gradually increased after day 7 ([153]Figure S4B). These observations indicated the CD56^dimCD16^+ NK cell activation in the non-ACR group started predominantly in phase 2 and was well controlled. The changes in the NK cell proportion were parallel between the ACR and non-ACR groups ([154]Figure 5C), but the functional dynamics varied ([155]Figures 5D and 5E). Specifically, a significant upregulation of KARs in phases 2 and 3 were observed in the ACR group, indicating that CD56^dimCD16^+ NK cells were continuously overactivated starting from phase 2. Despite comparable levels of cytotoxicity and proliferation between the ACR and non-ACR groups, the IFN-γ response signaling pathway of CD56^dimCD16^+ NK cells was highly upregulated and longer lasting in the ACR group, indicating a robust and prolonged inflammatory response of NK cells in ACR patients. In addition, we found different distributions of CD56^dimCD16^+ NK cells among the non-ACR, ACR, and HCs groups ([156]Figure 5F). CD56^dimCD16^+ NK cells from ACR patients upregulated the transcription of ISGs (e.g., STAT1, XAF1) and ETS1 ([157]Figure 5G), an important regulator for human NK cell development and terminal differentiation.[158]^39 Subsequent reclustering of CD56^dimCD16^+ NK cells identified six distinct subclusters ([159]Figures 5H and 5I). Pseudotime trajectory analysis of these cells revealed that the main path along the pseudotime trajectory started from HCs sets, followed by the ACR sets, and ultimately ended with the non-ACR sets ([160]Figures 5H and [161]S4C), indicating different states of CD56^dimCD16^+ NK cells in the ACR group and the non-ACR group. By comparing our NK subsets to the published ones ([162]Figure 5J; [163]Table S4),[164]^40^,[165]^41 we identified cluster 0 as adaptive NK cells characterized by the expression of CD52, CD3E, and IL-32. Cluster 1 was similar to inflamed NK cells, expressing IFN-responding markers such as STAT1 and XAF1. Clusters 2 and 5 shared similar expression patterns, but cluster 2 downregulated metabolic pathways and had a decreased inflammatory response compared to cluster 5 ([166]Figure S4D), indicating that cluster 2 was similar to terminal NK cells,[167]^40^,[168]^42 whereas cluster 5 was close to mature NK cells. Pseudotime trajectory analysis revealed that cluster 3 was predominantly located at the trajectory’s root and had a higher G1 score ([169]Figure S4E), identifying it as resting NK cells. Cluster 4, highly expressed erythrocyte-related genes (e.g., HBB, HBA1), was defined as erythrocyte-like NK cells, which may be generated by the adhesion of erythrocyte fragments. Significant differences in these NK cell subsets were observed between the ACR and the non-ACR groups ([170]Figure 5K). Adaptive NK cells were more enriched in LT patients, as confirmed by flow cytometry in cohort 3 ([171]Figures S4F and S4G). Most important, the ACR group had more inflamed NK cells than did the non-ACR groups, especially in phase 3 ([172]Figure S4H). Flow cytometry confirmed that inflamed NK cells produced more IFN-γ ([173]Figures S4I and S4J) and were more enriched in the ACR group of cohort 3, both at baseline and in phase 3 ([174]Figures 5L, 5M, and [175]S4K). However, among ACR patients, the abundance of inflamed NK cells tended to decrease after rejection following a high dose of methylprednisolone (MePDN) treatment ([176]Figure 5K). Evidence of the suppressive impact of MePDN on NK cell activity has been reported.[177]^43 Accordingly, we analyzed the influence of MePDN therapy on the inflamed NK cells of LT patients in vitro. We found a significant reduction in the proportion of inflamed NK cells after MePDN treatment, both in the presence and absence of tacrolimus ([178]Figure 5N), indicating that inflamed NK cells may be a potential indicator for early steroid intervention. The aforementioned observations suggested a relatively well-controlled activation of NK cell subsets in the non-ACR group but stronger and longer lasting inflammatory responses of NK cells in the ACR group, which may be attributed to the different distributions of inflamed NK cells. Monocytic subsets play an important role in immune remodeling after LT Circulating monocytes play a key role in acute rejection through diverse pathways, such as innate pattern recognition receptors (PRRs)[179]^44 and antigen presentation.[180]^45^,[181]^46 Although the fate of human monocyte subsets in a steady state and systemic inflammation has been studied,[182]^47 their dynamics post-LT are poorly understood. Here, we identified seven monocyte clusters ([183]Figures 6A and 6B), including the classical monocyte subsets and novel subsets with distinct functions. CD14^+RNASE2^+ monocytes were characterized by the elevated expression of RNASE2 and antioxidant activity-related genes (e.g., ALOX5AP, MGST1) and were involved in the toll-like receptor binding and the RAGE receptor binding pathways ([184]Figure S5A). Both CD14^+FOS^+ monocytes and CD14^+AP-1^+ monocytes expressed FOS, which is believed to mediate acute inflammatory responses.[185]^48^,[186]^49 Figure 6. [187]Figure 6 [188]Open in a new tab Dynamics of monocytic subpopulations across LT (A) UMAP plot showing monocyte subpopulations of all of the samples, with the percentages of each cell type shown. (B) Dot plot showing scaled expression levels of monocyte-type-specific genes in non-ACR, ACR, and HCs. (C) UMAP plots showing all merged monocytes colored by groups and phases. (D) Bar plot showing the frequency of monocyte subsets in the different groups across the phases of LT. (E) Line graphs showing longitudinal changes in the relative frequency of the monocyte subpopulations in the non-ACR, ACR, and HCs. (F) Dot plot showing scaled expression levels of function-specific genes in non-ACR, ACR, and HCs. (G) Line graphs demonstrating temporal changes of the corresponding function-specific scores in all of the monocytes. (H) Flow cytometry analysis of delta median fluorescence intensity (ΔMFI) of RNASE2 in CD14^+CD16^− monocytes from patient samples treated with or without Golgi Stop for 4 h. Representative images (left) and quantification data from paired groups (right, n = 40 per group) are shown. ΔMFI was determined by the subtraction of isotype control from antibody staining. Paired t test. (I and J) Boxplots showing the plasma concentration of RNASE2 at different time points after LT. Paired t tests along the time (I) and unpaired Wilcoxon tests between ACR (n = 9) and non-ACR (n = 10) groups (J) were performed. (K and L) ROC curves for RNASE2 plasma concentration and ACR on day 1 (K) and day 3 (L) after LT. In (E) and (G), data are shown as mean ± SEM, and comparisons were completed between the ACR and non-ACR groups and between phase 4 and HCs. Wilcoxon test. ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001; ∗∗∗∗p < 0.0001. In addition to the diverse cell composition, a high degree of interphase and intergroup heterogeneity was observed in the monocytic subpopulations ([189]Figures 6C and 6D). First, we quantified the abundance of monocyte subsets across each phase after LT ([190]Figures 6D and 6E). Among non-ACR patients, CD14^+ monocytes and CD14^+RNASE2^+ monocytes were the top two most abundant populations in phase 1 and subsequently decreased, along with which the CD14^+CD16^+ monocytes increased sharply in phase 2 and then decreased. Meanwhile, CD16^+ monocytes and CD16^+C1QA1^+ monocytes began to increase and peaked in phase 3. In phase 4, we observed a general trend toward normalizing abundance in the monocyte subset population. Second, we analyzed the functional transition of monocytic subsets during the LT recovery process ([191]Figures 6F and 6G). We found that PRRs rapidly peaked in phase 1, paralleling the proportion change in CD14^+RNASE2^+ and CD14^+ monocytes, both expressing high levels of PRRs (e.g., TLR8, S100A8, CLEC7A). Meanwhile, the proliferation score peaked, whereas the differentiation score dropped to the bottom. These observations indicated that in phase 1, the preexisting monocyte subsets amplified in response to damage-associated molecular patterns (DAMPs) and dominated innate immune activation. From phases 2 to 3, PRRs reduced, accompanied by an increase in antigen presentation and differentiation score, indicating an ongoing functional transition among monocytes. The enhanced antigen presentation was confirmed by the elevated HLA-DR of monocytes from day 7 in the CyTOF data ([192]Figure S5D). Surprisingly, although the scores of PRR, proliferation, and differentiation in phase 4 were almost comparable to those in HCs, the antigen presentation function of the non-ACR group was still higher than that of HCs ([193]Figure 6G). This pattern may indicate that even with allogeneic antigens continuously recognized by recipients’ innate immune systems, novel homeostasis had been reestablished in non-ACR patients, and thus no rejection was triggered. For the proportion of monocyte subsets in the ACR group ([194]Figures 6D and 6E), the change of monocyte subsets in phase 1 was consistent with that in the non-ACR group, but with a relatively lower elevation of CD14^+RNASE2^+ monocytes. Intriguingly, FOS-expressing monocytes, including CD14^+FOS^+ monocytes and CD14^+AP-1^+ monocytes, increased and became the major monocytic subtypes apart from CD14^+ monocytes in the ACR group in phase 2. In phase 3, CD14^+FOS^+ monocytes continued this trend, emerging as the most abundant subtype in the ACR group. This distinct enrichment of FOS-expressing monocytes in the ACR group was supported by flow cytometry data in cohort 3 ([195]Figures S5B and S5C), indicating their potential role in modulating acute rejection. A higher differentiation score of monocytes was observed in the ACR group ([196]Figures 6F and 6G), aligning with the increased diversity of monocyte subsets in this group. In addition, we found that the ACR group had a lower PRR score in phase 1 than did the non-ACR group, indicating a weaker innate immune response in this group. This corresponded with the lower proportion of CD14^+RNASE2^+ monocytes of the ACR group in phase 1. These findings strongly suggest that differences in CD14^+RNASE2^+ monocytes may underlie the attenuated innate immune response in phase 1 of the ACR group. Our data revealed that RNASE2 was highly expressed by monocytes ([197]Figure 6B). In addition, the increased RNASE2 expression during phase 1 was validated in cohort 3 ([198]Figure S5E). Given that RNASE2 is usually secreted by eosinophils,[199]^50 we posited that it could also be secreted by monocytes, a hypothesis that we validated through experiments in cohort 3 ([200]Figures 6H and [201]S5F). ELISA analysis of patients in cohort 3 revealed an increase in plasma RNASE2 levels in phase 1 post-LT, followed by a gradual decrease ([202]Figure 6I). Moreover, the non-ACR group exhibited higher plasma RNASE2 levels than the ACR group on days 1 and 3 ([203]Figure 6J). Receiver operating characteristic (ROC) analysis revealed that plasma RNASE2 levels on days 1 and 3 may be effective early biomarkers for ACR incidence, with an area under the ROC curve values of 0.889 and 0.800 ([204]Figures 6K and 6L), respectively. Differences in the cell-to-cell communication network between the non-ACR and ACR groups Following the identification of 22 distinct cell clusters ([205]Figure S6A), we assessed the potential intercellular communications after LT. Although the total number of communications in the non-ACR group was comparable to that in the ACR group across all of the phases, the non-ACR group demonstrated notably stronger interaction strength, especially during phase 1 ([206]Figures S6B and S6C). Next, we analyzed changes in communication patterns after LT in these two groups ([207]Figure 7A). The communications in the non-ACR group were generally regular and stable, with CD14^+ monocytes serving as the primary signal source across most phases. CD14^+RNASE2^+ monocytes in phase 1 and CD14^+CD16^+ monocytes in phase 2 served as an additional signal source, suggesting a transition in innate immune responses from activation by DAMP signals to antigen presentation. Consistently, starting from phase 2, central memory CD4^+ T cells and B cells were maintained as additional signal sources, suggesting the activation of adaptive immune responses after innate immune activation. Notably, the strength of the above-mentioned signal sources was extremely high in phase 1 but declined from phases 2 to 4 and gradually approached that of HCs, suggesting a well-regulated immune activation in the non-ACR group. Conversely, the cellular communications in the ACR group were maintained at a low level, but with irregularity. CD14^+ monocytes remained the major signal source across most phases, but with lower strength. Notably, the communication strength of CD14^+RNASE2^+ monocytes in phase 1 was much lower than that in the non-ACR group, in line with their proportions and RNASE2 levels. In addition, FOS-expressing monocyte subsets persisted as additional signal sources from phase 2 to the time of rejection, underscoring their potential role in mediating inflammatory responses and graft rejection. Figure 7. [208]Figure 7 [209]Open in a new tab Multicellular interaction analysis after LT (A) Scatterplots comparing the outgoing and incoming interaction strength in the 2-dimensional space across phases in the non-ACR and ACR groups. (B) Heatmaps showing the relative importance of each cell type as the sender, receiver, mediator, and influencer, based on the computed 4 network centrality measures of galectin signaling of ACR and non-ACR based on the baseline and phase 3 datasets. A comparative analysis of information flow further specified the shared and unique signaling pathways in the intercellular communication network between the non-ACR and the ACR groups ([210]Figure S6D). Pathways with shared flow include LCK, MIF, MHC class II, SELPLG, CD22, and GALECTIN, which would be involved in general immune responses. Specifically, the GALECTIN pathway, predominantly driven by monocytes ([211]Figures S6E and S6F), was more responsive in CD56^dimCD16^+ NK cells in the ACR group in both baseline and phase 3 ([212]Figure 7B). As mentioned, inflamed NK cells characterized by enhanced IFN signaling were enriched both at baseline and in phase 3 of the ACR group. Given that galectin-9, a member of the galectin family, is known to diminish NK cell cytotoxicity while augmenting IFN-γ production,[213]^51^,[214]^52^,[215]^53 it is plausible to speculate that monocyte-derived GALECTIN signaling triggers the “inflamed” transcriptional programs in NK cells. Further studies are needed to verify these hypotheses. Discussion Although multiple studies have elucidated pathophysiological elements involved in the modulation of immune tolerance or graft rejection after LT, the longitudinal dynamics of immune cells under immunosuppression conditions remains unknown. This knowledge is key to improving the clinical management of LT patients, such as precision immunosuppression strategies, early identification, and intervention for rejection. To date, immunosuppressant dose adjustment after LT relies merely on the TDM strategy, lacking prognostic biomarkers and assessments, which is far from precision management, often leading to multiple side effects and rejection. Furthermore, achieving target blood concentrations does not always guarantee efficacy or diminish adverse events. Therefore, individualized decisions are required considering the postoperative recovery process and immune status. Our analyses for the first time delineate a four-phase reestablishment process of immune homeostasis under immunosuppressive treatment after LT and provide a longitudinal atlas of peripheral immune cells and their transcriptomic profiles during this process. In addition, our analyses suggest the critical role of the inflamed NK cell subsets, CD14^+RNASE2^+ monocytes, and FOS-expressing monocytes in the development of ACR. These findings have significant implications for understanding the mechanisms underlying the immune remodeling after LT and provide potential targets for the management of post-LT ACR. As the transplanted liver progressively adapts to the host, with immunogenicity and rejection risk diminishing,[216]^54 the necessary dose of immunosuppressive drugs gradually decreases.[217]^55^,[218]^56^,[219]^57 The phases of liver function restoration post-LT have been demarcated into the convalescence period, followed by a stabilizing period.[220]^58 Similarly, Huang et al. identified an acute and stable phase during perioperative graft recovery using a syngeneic mouse LT model.[221]^3 These findings together suggest that a staged immune state exists in LT patients but remains undefined. Immune cells are pivotal in the emergence of postoperative complications such as infections[222]^59 and graft rejections.[223]^60^,[224]^61 Here, meticulously selected homogeneous study subjects without other postoperative complications minimize various confounding factors and allow us to compare the main immunological processes that mediate LT recovery. Our findings suggest and validate a four-phase LT recovery process characterized by the sequential activation of innate and adaptive immune compartments. The chronological changes in cellular constituents and associated immune characteristics elucidate peripheral immune responses in each phase, providing valuable insights into proposing a reasonable immunosuppressive utilization or withdrawal scheme ([225]Figure 8). Figure 8. [226]Figure 8 [227]Open in a new tab A schematic illustration of a four-phase framework that aids in the clinical management of LT patients In phase 1 (the first 7 days), in the scenario of the most potent immunosuppressant prescription, together with surgical trauma, stress responses, allorecognition, and immune responses after LT, we observed significant enrichment and activation of innate immune cells (e.g., CD14^+ monocytes, CD14^+RNASE2^+ monocytes) and inhibition of the activity and proportion of T cells and NK cells. In addition, we identified and validated a higher activation of CD14^+RNASE2^+ monocytes in the non-ACR group in this phase, manifested by the secretion of RNASE2. The strength of innate immune activation together with the restricted adaptive immune responses in phase 1 may guide precise immunosuppressant administration. Indeed, surgical stress and initial treatment of immunosuppressants post-LT can induce excessive immunosuppression, leading to severe adverse effects, such as infection and tumor recurrence.[228]^62^,[229]^63^,[230]^64 Minimizing exposure to calcineurin inhibitors (<10 ng/mL for tacrolimus) early after LT prevents tumor recurrence without increasing ACR incidence.[231]^62 Furthermore, a separate study examining 493 LT recipients under tacrolimus treatment demonstrated that a higher dosage of tacrolimus (10–15 ng/mL) correlated with increased mortality rates compared to a more moderate dosage (7–10 ng/mL), and yet the rates of rejection remained comparably consistent.[232]^65 These findings and our data together indicate that the recommended high-dose immunosuppressant treatment (>10 ng/mL for tacrolimus) may not be necessary during the first week after LT, and monitoring the level of CD14^+RNASE2^+ monocytes/RNASE2 during this period could facilitate the early identification and intervention for patients at risk of ACR. With the relief of injuries and stresses, the graft function is restored. The main immune features switched to alloresponse and posttraumatic stress reparative responses, including an enhancement of monocyte antigen-presenting capacity and the activation of both NK cells and T cells, which peaked at phase 2 (from days 7 to 19) and lasted until phase 3 (from days 19 to 90). Moreover, a longer lasting inflammatory response of NK cells was observed in the ACR group in phase 3, as indicated by the increased number of inflamed NK cells, which can be decreased by high-dose glucocorticoid treatment. Interestingly, the monocytes from ACR patients also showed an enhanced differentiation capacity in phase 2, as indicated by the increased diversity of monocytes, such as the increased proportions of CD14^+FOS^+ monocytes and CD14^+AP-1^+ monocytes. The relatively strong immune responses within 90 days after LT may explain at least in part why the tacrolimus medication guideline recommends potent target tacrolimus trough concentration[233]^56^,[234]^57 and why acute rejections usually occur within the first 3 months after LT.[235]^66^,[236]^67 All of these findings in phases 2 and 3 indicate that LT patients are more susceptible to acute rejection in these phases, and the FOS-expressing monocytes in phase 2 and inflamed NK cells spanning from phases 2 to 3 may be predictive biomarkers and potential therapeutic targets for clinical management and graft rejection. Liver function is normalized and the rejection risk is reduced after 90 days of LT.[237]^67 Importantly, we observed that the immune status of LT patients from phase 4 (from days 90 to 365) shared many similarities to that of HCs, but disparities still persisted, indicating that LT and immunosuppressive therapy result in durable immune status alterations. Nevertheless, our data suggest that 90 days after LT would represent a time point of immune steady state for non-ACR patients and a critical juncture for minimizing tacrolimus exposure to decrease the complications of long-term immunosuppression.[238]^55^,[239]^68 Following the identification of the four-phase immune remodeling process after LT, we further revealed critical immune cells for immune tolerance or rejection dynamically modulated across different phases. Notably, we observed that T cells in the ACR group displayed an unexpectedly weak response during LT recovery, characterized by an increased presence of anergic CD8^+ T cells and reduced CTLs activity compared to the non-ACR group, suggesting more potent T cell suppression in the ACR group. Although the T cell composition in peripheral blood may not mirror that in the liver,[240]^69 this finding somehow challenged the conventional view of inhibiting T cell function for graft tolerance in LT. Blocking T cell activity does not universally prevent graft rejection,[241]^70^,[242]^71^,[243]^72 indicating that existing immunosuppressive strategies do not entirely specifically suppress the alloimmune responses that lead to rejection. In the context of tacrolimus-based immunosuppressive treatment, we identified rejection-related cellular features in NK cells and monocytes through extensive analysis of our longitudinal data. First, we newly identified and validated that inflamed NK cells with enhanced IFN signaling were more enriched in the ACR group in phase 3 after LT. As reported, NK cells in the circulation would swiftly infiltrate the liver graft after revascularization to replace the graft NK cells[244]^73^,[245]^74 and serve as a major source of IFN-γ in grafts.[246]^12 In addition, IFN-γ is a well-characterized factor for ACR after transplantation and is considered a predictive biomarker of acute rejection in both liver and kidney transplants.[247]^75^,[248]^76 We thus reasoned that heightened recruitment inflamed NK cells to the liver in the ACR group triggers or mediates rejection through an IFN-γ-mediated response. The current immunosuppressants for LT do not specifically inhibit monocyte-macrophage lineage cells.[249]^45 Second, we discovered an upregulation in the proportion and communication of FOS-expressing monocytes in the ACR group. Consistently, Cui et al. reported that inflammatory monocytes overexpressing FOS in the liver exhibit the most significant similarity to PBMCs, which are increased during liver inflammation and decreased after inflammation subsided.[250]^77 FOS and JUN co-regulate interleukins to modulate adaptive immunity. We thus propose that FOS-expressing monocytes mediate graft rejection by activating the adaptive immune response, which requires further exploration. Together, these findings suggest the significant roles of both NK cells and monocytes in immune remodeling after LT, and targeting these cells could refine tacrolimus-based immunosuppression strategies for LT patients. The present study has certain limitations, notably the small sample size of cohort 1, particularly within the ACR group. This prompted us to validate the critical results in two independent validation cohorts. In addition, future studies with larger and more diverse cohorts are needed to comprehensively delineate the complexity of immune responses post-LT and to broaden our findings. Also, our study only illustrated the immune state of PBMCs after LT, providing dynamic insights into the systemic immune landscape but still insufficient to clarify the local immune environment in the liver graft itself. Integrating our data with that from the liver graft would undoubtedly be more systematic and informative, leading to more comprehensive conclusions. Such work is ongoing in our group. Collectively, our study has revealed a temporal landscape of immune cells throughout the LT recovery process and provided key insights into the dynamic immune responses associated with ACR. These results represent the molecular and cellular signatures of the immune cell populations during the LT recovery process and generate a four-phase framework that aids in the clinical management of patients. Materials and methods See the [251]supplemental information. Acknowledgments