This study aimed to develop a clinical decision-point risk stratification model for early graft dysfunction evaluation rather than an asymptomatic screening tool. This hypothesis-generating exploratory analysis was a single-center retrospective cohort study conducted at Beijing Chaoyang Hospital (August 2021-June 2024). For rejection cases, baseline measurements were intentionally defined at the first clinical suspicion of graft dysfunction to evaluate whether the model could provide superior acute rejection risk quantification using routinely available tests at this earliest point of clinical concern. The study aimed to identify potential immune signatures rather than test predefined hypotheses, thus no formal sample size calculation was performed a priori. The study protocol (IRB No. 2024-s-725) was approved by the Beijing Chaoyang Hospital Ethics Committee. Written informed consent was obtained for all peripheral-blood sampling; a waiver of consent was granted for retrospective collection of anonymized clinical data. All transplant procedures utilized organs from deceased Chinese citizens following post-mortem organ donation. The organ procurement and allocation strictly adhered to the ethical guidelines and laws of the People’s Republic of China.

The study population comprised liver transplant recipients managed at Beijing Chaoyang Hospital’s outpatient clinic. All patients were categorized into: (1) stable controls with stable liver function, and (2) rejection cases hospitalized for suspected acute rejection, confirmed by subsequent liver biopsy and reviewed by two independent hepatopathologists (Banff 2019 criteria) (16). The primary endpoint of this study was biopsy-proven ACR.

To emulate a prospective prediction scenario within our retrospective design, we employed a landmark analysis approach. For all patients, the ‘baseline’ measurements used in the model were defined as those obtained during routine clinical monitoring prior to any clinical suspicion of rejection. This ensures that the model predictions are based on information that would have been available to clinicians at the time of routine follow-up, thus validating its pre-emptive potential.

Exclusion criteria were systematically applied to eliminate confounding factors: (1) incomplete clinical or immunological data; (2) pre-existing immune-related conditions including autoimmune hepatitis, HIV infection, or hematologic malignancies; (3) ABO incompatible transplants; (4) re-transplantation or combined organ transplantation; (5) age < 18 years at the time of transplantation; (6) follow-up period < 6 months.

All patients received intraoperative induction therapy with basiliximab (20 mg on day 0 and day 4 post-transplantation), as per our institutional standard protocol. The immunosuppressive regimen consisted of: (1) tacrolimus with initial target trough levels of 5-8 ng/mL, adjusted based on clinical response; (2) mycophenolate mofetil at 500 mg twice daily, with dose reduction to 250 mg twice daily for white blood cell counts 2.0-3.0×10^9/L and complete discontinuation for counts < 2.0×10^9/L; and (3) prednisone discontinued at 3 months post-transplant per institutional practice in patients with benign liver diseases, while switched to sirolimus at 1 month in those with kidney injury or with malignant indications (17–19).

For acute rejection episodes, tacrolimus dose adjustments followed a standardized protocol (20, 21): (1) primary escalation to target 8-11 ng/mL for rejection; (2) secondary augmentation with a round of steroids: intravenous methylprednisolone 500 mg on day 1, 240 mg on day 2, then daily reduction of 40 mg until 20 mg/day is reached, followed by oral prednisolone 20 mg/day for approximately one month if AST/ALT failed to decline > 50% within 72 hours or if neurotoxicity (tremor ≥ grade 2) occurred. Steroid-resistant cases (defined as < 50% improvement in liver biochemistry within 5 days) received rabbit anti-thymocyte globulin (Thymoglobulin, Genzyme) at 1.5 mg/kg/day for 5 days, preceded by rigorous screening for active infections.

Peripheral blood mononuclear cells (PBMCs) were analyzed following our established protocol (22, 23). The sequential gating strategy for lymphocyte subset identification is illustrated in Supplementary Figure 1. Briefly, PBMCs processed within 2 hours of collection to ensure > 95% viability (Trypan blue exclusion) were isolated from EDTA blood by Ficoll density centrifugation and stained with the following fluorochrome-conjugated antibodies from BD Biosciences: CD3-FITC (clone UCHT1, Cat. 555332), CD4-PerCP-Cy5.5 (clone SK3, Cat. 566923), CD8-PerCP-Cy5.5 (clone RPA-T8, Cat. 560662), CD19-PE (clone HIB19, Cat. 555413), CD56-PE (clone B159, Cat. 555516), CD16-APC (clone B73.1, Cat. 561304), CD11c-APC (clone B-ly6, Cat. 559877), CD123-PerCP-Cy5.5 (clone 7G3, Cat. 560904). Samples were acquired on a NovoCyte D2060R flow cytometer with a minimum of 10,000 viable lymphocytes per tube. Technical duplicates were run for 20% randomly selected samples to ensure reproducibility. Daily calibration was performed with CS&T beads, and compensation was verified using single-stained controls. Data were processed with NovoExpress software using standardized gating strategies.

Baseline tacrolimus levels and lymphocyte subset percentages were defined differently for the two study groups. In the non-rejection group, these values represent levels measured during routine outpatient monitoring, immediately prior to drug administration. In the rejection group, baseline referred to measurements obtained at the earliest indication of graft dysfunction, preceding the diagnostic liver biopsy. This included either the visit when liver function tests first demonstrated unexplained elevation, or the most recent routine monitoring visit before hospitalization for suspected rejection. Patients with biopsy-proven acute rejection underwent additional measurement approximately 2 weeks after the initiation of anti-rejection therapy.

All statistical analyses were conducted using R software (version 4.3.1; R Foundation for Statistical Computing) with the following packages and methodologies: Normality of continuous variables was assessed using Shapiro-Wilk tests with visual confirmation by Q-Q plots. Between-group comparisons employed Student’s t-tests for normally distributed data or Mann-Whitney U tests for non-parametric distributions, while longitudinal analyses used paired t-tests or Wilcoxon signed-rank tests as appropriate. Multiple testing correction was performed using the Benjamini-Hochberg false discovery rate (FDR) method with α=0.05. Primary immunological endpoints were FDR-adjusted (q<0.05); all other comparisons are reported as nominal P values.

High-dimensional flow cytometry data were visualized using Uniform Manifold Approximation and Projection (UMAP; uwot package, version 0.1.16) with default parameters (n_neighbors=15, min_dist=0.1, metric=“euclidean”). Dynamic time warping (DTW) analysis was implemented through the dtw package (version 1.23-1) to quantify temporal patterns in lymphocyte subset changes, using the symmetric2 step pattern with open-end alignment. Euclidean metric was selected for DTW to preserve absolute magnitude differences in lymphocyte counts, as opposed to cosine similarity which would emphasize relative patterns. DTW was chosen over linear mixed models to capture non-linear temporal patterns in immunosuppression response. Distributional differences between groups were evaluated using Wasserstein distance metrics computed via the transport package (version 0.13-0). The Wasserstein distance quantifies the minimum work required to transform one probability distribution into another.

Given the small sample size and rare outcome (18 events among 98 patients), Firth-penalized logistic regression was employed to reduce small-sample bias and obtain more reliable estimates. Internal validation via 1000 bootstrap iterations provided optimism-corrected performance metrics. Multicollinearity among candidate predictors was assessed using variance inflation factors (VIFs) computed from linear models. VIF values < 5.0 were considered acceptable, indicating no substantial multicollinearity. CD3+ parameters were excluded from multivariate models due to collinearity with CD4+/CD8+ subsets (VIFs >5), as confirmed by the car package (version 3.1-2).

Model performance was assessed via the area under a receiver operating characteristic (AUC-ROC) analysis (pROC package, version 1.18.0), calibration curves, and decision curve analysis (rmda package, version 1.6), with 1,000 bootstrap iterations for internal validation. Sensitivity analyses with 20-fold multiple imputation confirmed robustness (ΔAUC <0.01). Lead-time analysis was performed to quantify the early warning capability of the model. For each patient who developed rejection and was identified as high-risk by the model, lead time was calculated as the interval between the date of model positivity (based on the Youden-optimized probability threshold) and the date of first biochemical evidence of graft injury (defined as ALT/AST > 2× upper limit of normal). The study is powered to detect Pearson correlations ≥ 0.30 (two-sided α = 0.05, power = 0.80) as calculated with the ‘pwr’ package; weaker associations may have been missed.

The baseline demographic and clinical characteristics of the study cohort are summarized in Table 1. Among the 98 liver transplant recipients, 18 developed biopsy-proven ACR (AR group) and 80 maintained stable graft function (Non-AR group). The median time to acute rejection was 5 months (range, 1-11 months). No statistically significant differences were observed in age, gender distribution, liver disease, MELD score, or key operative parameters including operating time, warm ischemia time, cold storage time, anhepatic phase duration, operative bleeding, and transfusion requirements between the two groups (all P > 0.05). Similarly, the distribution of immunosuppressive regimens or complications did not differ significantly between groups. This comparability strengthens the subsequent analysis by minimizing the confounding effects of these baseline variables.

However, patients in AR group demonstrated significantly lower pre-treatment tacrolimus trough levels (5.15 ng/mL [3.85-6.30]) compared to non-rejection controls (6.60 ng/mL [4.90-7.80], p=0.0086). Following rejection therapy, tacrolimus levels increased to 8.25 ng/mL [7.15-9.35] (pre- vs post-treatment p<0.001). The rejection cohort showed higher peak ALT (498 U/L [145-820]) and bilirubin (3.4 mg/dL [2.1-5.0]) at diagnosis compared to stable controls (ALT < 60 U/L, bilirubin < 1.5 mg/dL). Post-treatment infection rates included cytomegalovirus viremia (11.1%), bacterial (22.2%), and fungal infections (16.7%).

Comparative analysis revealed significant temporal changes in T cell subsets during post-transplant follow-up (FDR-adjusted q < 0.05, Supplementary Table 1; Figure 1): CD3+T% increased from 44.7 ± 24.5% at ≤ 1 month to 66.7 ± 14.4% at > 6 months (p < 0.001, q = 0.029), while CD4+T% showed early elevation from 18.3 ± 13.4% to 34.8 ± 11.9% at 2-3 months (p < 0.001, q = 0.029). Absolute counts demonstrated parallel increases, with CD3+T cells rising from 388 ± 278 to 951 ± 510 cells/μL (p < 0.001, q = 0.029), CD4+T cells from 170 ± 171 to 404 ± 213 cells/μL (p < 0.001, q = 0.034), and CD8+T cells exhibiting delayed recovery from 220 ± 158 to 549 ± 344 cells/μL (p < 0.001, q = 0.034). No significant differences were observed in any lymphocyte subsets when stratified by sex (male vs female, all q > 0.73), age groups (≥ 60 vs < 60 years, all q > 0.91), or primary disease (malignant vs benign, all q > 0.62) (Supplementary Figures 2-4).

UMAP analysis identified strong correlation between UMAP1 coordinates and CD4+T cell parameters (CD4+T%: r = 0.49, p < 0.001; CD4+T count: r = 0.58, p < 0.001). UMAP2 primarily reflected CD8+T cell distribution (CD8+T%: r = 0.69, p < 0.001). Early post-transplant patients (≤ 1 month) clustered distinctly in UMAP space (Wasserstein distance=3.21, Supplementary Figure 5). The immune landscape remains relatively homogeneous across sex, age, and primary liver disease in stable transplant recipients (Supplementary Figure 6).

These findings suggest distinct temporal reconstitution patterns for major T cell populations, while other lymphocyte subsets remained stable across demographic and etiological subgroups.

To evaluate the immunological alterations following anti-rejection treatment in liver transplant recipients, we compared peripheral blood lymphocyte subsets (both percentages and absolute counts) in 18 patients before and after therapy. CD3+T cells decreased from 75.67 ± 9.53% to 47.20 ± 16.05% (p < 0.001, q < 0.001); CD4+T cells declined from 38.82 ± 11.97% to 20.96 ± 10.18% (p < 0.001, q < 0.001). CD8+T cells showed moderate reduction (36.65 ± 8.00% vs. 26.36 ± 9.27%, p < 0.001, q < 0.001). No significant changes occurred in CD19+B cells, NK cells, NKT cells, or dendritic cell subsets (all q > 0.05) (Figures 2A, B; Supplementary Table 2).

More dramatic depletion was observed in absolute counts. CD3+T cells plummeted from 861 ± 372 to 320 ± 227 cells/μL (p < 0.001, q < 0.001). CD4+T cells dropped from 441 ± 228 to 143 ± 116 cells/μL (p < 0.001, q < 0.001). CD8+T cells decreased from 420 ± 186 to 178 ± 122 cells/μL (p < 0.001, q < 0.001). NKT cells showed marginal reduction (20.7 ± 21.4 vs. 10.3 ± 9.9 cells/μL, p = 0.021, q = 0.045). CD19+B cells and NK cells demonstrated non-significant downward trends, while dendritic cells remained stable (all q > 0.05) (Figures 2C, D; Supplementary Table 2). This synchronous T-cell depletion likely reflects glucocorticoid receptor-mediated apoptosis (24), while preserved innate immunity suggests calcineurin inhibitor-specific effects on nuclear factor of activated T-cells signaling (25).

The data reveal selective depletion of T lymphocyte subsets, particularly CD4+T and CD8+T populations, with relative preservation of innate immune cells after anti-rejection therapy. The magnitude of reduction was more pronounced in absolute counts than percentages, suggesting both proportional redistribution and genuine cell loss.

Anti-rejection therapy triggered coordinated declines in both the percentage and absolute counts of major T cell compartments (Figure 3; Supplementary Table 3). Specifically, CD3+T cell percentage dropped by a median of 25.7% (q < 0.001), while CD3+T, CD4+T and CD8+T absolute counts declined by 543, 284 and 284 cells/μL, respectively (all q < 0.001). A parallel fall in CD4+T cell percentage (median −17.2%, q < 0.001) and CD8+T cell percentage (median −9.4%, q = 0.028) was observed. Among non-T subsets, only myeloid dendritic cells (mDC) and NKT cells exhibited modest absolute reductions (median −0.9 cells/μL and −10.4 cells/μL; p = 0.023 and 0.021, respectively), although these did not withstand FDR correction. CD19+B cell and NK cell parameters remained stable (q > 0.05). Spearman correlation analysis of the significant change scores (Δ) demonstrated an inverse relationship between ΔCD4+T cell percentage and ΔCD8+T cell percentage (r = −0.43, p = 0.046) and a positive correlation between ΔCD4+T cell percentage and ΔCD4+T cell absolute change (r = 0.59, p = 0.012), indicating coupled CD4 dynamics during immune recovery (Supplementary Figures 7, 8).

Then, to quantify the temporal coherence of these changes, DTW analysis revealed pre- and post-treatment time-series of all 16 lymphocyte parameters. DTW distance ranking highlighted CD3+T cell absolute count, CD4+T cell absolute count, and CD8+T cell absolute count as the three most consistently reshaped subsets (mean DTW distances 28.5, 17.9 and 10.3, respectively), whereas CD19+B-cell, NK, pDC and NKT compartments showed low distances, underscoring minimal or inconsistent responses (Supplementary Figure 9; Supplementary Table 4). Heat-map visualization confirmed a uniform red-to-blue shift (decrease) across the T-cell zone, while non-T subsets remained color-neutral (Supplementary Figure 10). Trajectory plots further revealed tight convergence of individual patient lines toward lower T cell values (Supplementary Figure 11), contrasting with the widely scattered paths of B and NK cells (Supplementary Figure 12).

Taken together, DTW analysis corroborates the selective and synchronous suppression of T-cell subsets as the hallmark of effective anti-rejection therapy.

Comparative analysis revealed significant differences in lymphocyte subsets between non-rejection (n=80) and pre-treatment acute rejection (n=18) groups (Figure 4; Supplementary Table 5). CD3+T cell and CD4+ T cell percentages were markedly higher in rejection patients (CD3+T: 78.5% [69.3-82.3] vs 65.5% [54.5-77.0], q = 0.030; CD4+T:38.8 ± 12.0% vs 28.5 ± 13.4%, q = 0.030). The volcano plot highlighted CD3+% and CD4+% as the most differentially distributed subsets (effect sizes >1.2, Supplementary Figure 13A), confirmed by boxplot visualization (Supplementary Figure 13B). No significant differences were observed in CD8+T cells, NK/NKT cells, or DC subsets (all q > 0.05). Clinical parameters (age, sex, follow-up duration) and primary liver disease distribution were comparable between groups (Supplementary Figure 14).

This integrated analysis demonstrates that acute rejection is associated with specific alterations in T cell homeostasis, while innate immunity remains preserved. The combination of traditional statistics with distribution-based metrics (Wasserstein distance) provides robust characterization of rejection-associated immune signatures.

Guided by univariate analysis identifying tacrolimus levels (OR 0.739, 95% CI 0.555-0.935, p=0.009) and CD4+ T cell percentage (OR 1.060, 95% CI 1.020-1.110, p=0.004) as significant correlates of acute rejection, we constructed a parsimonious prediction model. Comprehensive evaluation of other lymphocyte subsets revealed no additional independent predictive value beyond the core CD4+ T cell-tacrolimus axis. In 98 recipients (18 rejection, 80 non-rejection), the final Firth-penalized logistic model retained only baseline tacrolimus (OR per ng/mL 0.732, 95% CI 0.537-0.947, p=0.015) and baseline CD4+ T cell percentage (OR per % 1.072, 95% CI 1.020-1.137, p=0.005) as independent predictors (Supplementary Figure 15). Collinearity diagnostics confirmed statistical independence between the predictors, with variance inflation factors of 1.04 for both tacrolimus level and CD4+ T-cell percentage, and a weak negative correlation (r = -0.205, P = 0.415). The model demonstrated good discrimination with an AUC of 0.774 (95% CI 0.674-0.874; Figure 5), sensitivity 66.7%, specificity 83.8%, and a Brier score of 0.123. Bootstrap internal validation (R = 1,000) confirmed model stability with a median optimism-corrected AUC of 0.785 (95% CI 0.685-0.885). The model exhibited excellent calibration across the risk spectrum (Hosmer-Lemeshow test: χ²=3.166, p=0.924; Supplementary Figure 16), and decision curve analysis demonstrated net clinical benefit across threshold probabilities of 10-50% (Supplementary Figure 17).

For clinical implementation, risk stratification using the optimal probability threshold of 0.178 provided 83.3% sensitivity and 65.0% specificity. Patients with both high CD4+ T% (≥30.1%) and low tacrolimus levels (<6.2 ng/mL) had substantially elevated rejection risk (Fisher’s exact p<0.001), identifying 72.2% (13/18) of rejection cases while maintaining 85.0% specificity.

Critically, lead-time analysis was performed on patients identified as high-risk by the model using the Youden’s index-optimized threshold (probability > 0.178), which corresponded to 15 model-positive patients (83.3% of the rejection cohort), providing early warning signals a median of 8 days (IQR: 3.5 days; range: 6-14 days) before biochemical evidence of graft injury emerged (Figure 6). Among these model-positive patients, all (15/15, 100%) showed positive lead times, confirming consistent early detection when the model identified rejection risk.

The combined model provided significant incremental value over individual components (Supplementary Figure 18), with absolute AUC improvements of 0.080 over tacrolimus monitoring alone (0.774 vs. 0.694) and 0.041 over CD4+ T cell monitoring alone (0.774 vs. 0.733). Likelihood ratio tests confirmed these improvements were statistically significant (p = 0.007 and p = 0.014, respectively).

Collectively, this streamlined two-variable immunopharmacologic model enables pre-emptive risk stratification for acute liver-transplant rejection, offering both early detection capability and clinically meaningful performance improvements over conventional monitoring approaches.