For baseline patient data (Supplementary Table S1), we included patients (n=55) who underwent liver transplantation at Organ Transplantation Center, the Affiliated Hospital of Qingdao University, among November 2020 and October 2021. Data related to patients consisted of gender, age, ABO blood group, and etiology while donor also included organ type, gender, age, ABO blood group, and laboratory data. All liver grafts were voluntarily donated after cardiac death or by living donors and we also included several healthy people as control, which were approved by the Ethics Committee of the Affiliated Hospital of Qingdao University (QYFYWZLL26550). For scRNA-seq, histological analysis, and flow cytometry, all human liver (n=4) and blood samples (n=4) were collected from patients after liver transplantation.

After liver puncture tissues were obtained, we placed them into tissue preservation solution immediately. Then, the tissues were transported to be dissociated routinely. Briefly, tissues were cut into approximately 1-2 mm³ pieces in the RPMI-1640 medium (Gibco) and digested into single-cell suspensions with the Tumor Dissociation Kit, human (Miltenyi) at 37°C for 60 min. For blood, we collected about 5 ml peripheral blood from one person using anticoagulant tube, which was then centrifuged to prepare PBMC (13). We prepare PBMC using the density gradient centrifugation technique. This is done using an isotonic solution (lymphocyte separator, HISTOPAQUE) with a specific gravity between red blood cells, multinucleated white blood cells, and lymphocytes. The specific gravity of lymphocytes and monocytes is less than or equal to that of the stratified solution, and they float on the surface of the stratified solution after centrifugation. The dissociated cells were subsequently passed through a 70-µm SmartStrainer and centrifuged at 400 g for 5 min. After the supernatant was removed, the pelleted cells were suspended in red blood cell lysis buffer (Miltenyi) to lyse red blood cells. After washing twice with 1× PBS (Gibco), the cell pellets were re-suspended in sorting buffer [PBS supplemented with 2% fetal bovine serum (FBS, Gibco)].

Cells were loaded approximately 17,400 cells/chip position using the 10x Chromium Single Cell 3’ Library, Gel Bead & Multiplex Kit and Chip Kit (10x Genomics, V3.1 barcoding chemistry) according to the manufacturer’s instructions. Purified libraries were analyzed by an Illumina nova-seq 6000 sequencer with 150-bp paired-end reads. We applied Cell Ranger (version 3.0.1, 10x Genomics) to generate single-cell data, which processed chromium single-cell RNA-seq outputs to align reads and generate the feature barcode unique molecular identifier (UMI) matrices.

We used the Seurat R package (version 3.2.0) importing raw data to quality control, normalize, scale, and further process data (14). Low quality cells were removed with the following criterion: cells had either fewer than 501 expressed genes, 2001 UMIs, or more than 25% mitochondrial counts. The remaining high-quality cells were normalized and scaled with the default parameters. Highly variable features were identified using FindVariableFeatures function. Then we performed PCA analysis on the scaled data with the determined variable features, dimension reduction, and cluster using FindNeighbors (dims = 1:10) and FindClusters (resolution = 0.5) function. Then we run non-linear dimensional reduction (tSNE) to explore and visualize data.

We used FindMarkers and FindAllMarkers function to find differentially expressed features of each cluster. Cell types were annotated to known biological types with the canonical marker genes (Supplementary Table S2) according to CellMarker database (http://bio-bigdata.hrbmu.edu.cn/CellMarker/) (15) and published articles (16–19). Furthermore, we also used the SingleR package (version 1.2.4) to help identify the cell types (20).

After the annotation of each cell type, we employed the functional enrichment analysis to differentially expressed genes between different clusters for Gene Ontology (GO) and KEGG (Kyoto Encyclopedia of Genes and Genomes) analysis to illustrate the biological process and potential function of different cells using clusterProfiler package (version 3.17.0) (21) and org.Hs.eg.db package (version 3.11.4). The p value cutoffs of GO and KEGG were both 0.05. The top 10 terms of results were visualized by barplot or dotplot.

Trajectory analysis was performed using monocle package (version 2.17.0) (22). Including T cells, NK cells, B cells, and myeloid cells, we performed analysis of each group with the following parameters: lowerDetectionLimit=0.5, min_expr=0.1, num_cells_expressed>= 10. For visualization, plot_cell_trajectory function was used to plot the potential trajectory according to pseudotime, seurat clusters, and meta data.

Cell-to-cell interactions were done using CellChat package (version 1.1.2) (23). The majority of ligand-receptor interactions were mainly on the basis of KEGG signaling pathway database and recent peer-reviewed experimental studies. The main steps of inference of intercellular communications are as follows: (1) Identification of differentially expressed signaling genes. (2) Calculation of ensemble average expression. (3) Calculation of intercellular communication probability.

The transcriptome data for validation were downloaded and processed from the NCBI Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) with accession number GSE145780 (7), which measured gene expression by microarrays in 235 liver transplant biopsies, including normal, rejection, and fibrosis samples. We performed further statistic analysis using R software (version 4.0.2).

We used CIBERSORT, a method for characterizing complex immune infiltration of tissues according to gene expression profiles, which was developed by Newman et al. (24) from the Alizadeh Lab and was popular in indication of immune infiltration.

The fresh liver tissues were fixed with 4% paraformaldehyde and embedded in paraffin, then sectioned in 4 μm in thickness and stained with hematoxylin and eosin (HE), reticular fiber for morphological evaluation, and Masson’s trichrome staining for fibrosis detection.

Briefly, slides were rehydrated with a series of graded ethanol solutions in deionized water. Antigen retrieval was performed, then slides were serially stained with the following antibodies: CD4, CD8, CD11b, CD14, CD15, and FOXP3. Subsequently, Opal IHC Detection Kit (Akoya Biosciences) was applied as a secondary label and antibody signals. Image acquisitions were performed using the Vectra Polaris multispectral imaging platform (Akoya Biosciences), with the entire slide image being scanned and 3-5 representative regions of interest chosen by the pathologist (25).

Paraffin sections are routinely dewaxed to hydration, and washed with distilled water. Following incubation in 3% H₂O₂ for 10 min, antibodies anti-CD4, CD8, CD11b, CD14, and CD15 were added and incubated at 4°C for overnight. The specimens were incubated with secondary antibodies at 37°Cfor 1 h, followed by diaminobenzidine staining (26).

Cells were stained using fluorochrome-conjugated antibodies anti-CD3, CD4, CD8a, CD14, CD15, CD45, HLA-DR, and CD11b (BioLegend) according to the protocols provided by the manufacturer. Samples were analyzed using FlowJo software.

For the experimental data, GraphPad Prism 8 (GraphPad Software) was used to perform statistical analyses and graphics production. Data represent mean ± SEM. Results were considered significant when P < 0.05. Differences among the three groups were analyzed using a one-way analysis of variance and Newman-Keuls test for post-hoc comparisons.

To reveal the cellular diversity and gene signatures, we performed single cell sequencing from 4 in human liver penetrating tissues and 4 PBMCs from patients after liver transplantation using 10X genomics sequencing, and to discriminate between liver-resident and circulating leucocytes (Figure 1A). The clinical meta data and HE staining of samples are shown in Supplementary Table S1 and Figure 1G. We also merged 3 healthy liver data of our previous work. We applied strict filtering on data, including cells expressed genes, unique molecular identifiers (UMIs), and mitochondrial counts (Supplementary Figures S1A, B), and further performed normalization and scaling.

We conducted dimensional reduction and clustering including principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE), of the combined liver tissues and PBMCs with 68,172 high quality cells (Supplementary Figure S1C). Compared with PBMCs, cells from liver tissues expressed higher levels of CXCR4 and CCL3L1 genes (Supplementary Figures S1D–F), which were liver-resident markers. These clusters were annotated using canonical marker genes (Supplementary Table S2 and Figures 1B, D), across 6 cell lineages, each containing cells from both health and liver transplantation tissues: T cells (27,976 cells, 41.0%, marked with CD3D, CD3E and CD3G), NK cells (14,477 cells, 21.2%, marked with CD7, FGFBP2, and KLRF1), B cells (4534 cells, 6.7%, marked with CD79A, CD79B, and MS4A1), myeloid cells (15,687 cells, 23.0%, marked with CD68, CD14, and CD163), endothelial cells (3944 cells, 5.8%, marked with PECAM1, ICAM2, and ERG), and fibroblasts (1554 cells, 2.3%, marked with COL1A2, COL3A1, and ACTA2). It was found that T, NK, and myeloid cells were dominant in the liver transplantation landscape across samples (Figures 1C, E). The differentially expressed genes of each type are shown in Figure 1F and Supplementary Table S3.

To delineate the intrinsic landscape and potential function of the whole T, NK cell of liver tissues, we first re-clustered the 31,300 liver-resident T, NK cells from 7 livers (n=3 health and n=4 liver transplantation), revealing 15 clusters (Figures 2A–C). It emerged 7 clusters for CD8+ T cells, 3 clusters for CD4+ T cells, 1 cluster for double negative T cells, and 2 clusters for NK cells (Figure 2D and Supplementary Figures S2A–D), each cluster containing cells from both healthy and transplanted livers (Figure 2B). Cells of the first CD8+ cluster, C3-GZMK-CD8T, specifically expressed marker genes including GZMK, DUSP4, and COTL1 (Figure 2E and Supplementary Table S4), which was thought to be cytotoxic CD8+ T cells. The second cluster, C4-CRTAM-CD8T, was dominant in healthy tissue significantly (Figure 2F), expressing the high levels of CRTAM, HLA-DQA2, and MLF1. The CRTAM gene was involved in regulation of CD8+ T cell regulation, so this cluster was identified as regulatory-like CD8+ T cells. Enrichment analysis also confirmed the phenomenon called “regulation of mRNA metabolic process” (Supplementary Figure S2E and Supplementary Table S4). However, it did not express traditional Treg markers such as FOXP3 and IL2RA. Moreover, this cluster showed a myeloid phenotype due to the expression of HLA-DQA2 and MLF1. The third cluster, C5-GZMH-CD8T, also showed a cytotoxic phenotype, characterized by specific expression of killing markers including GZMH and GZMA (Figure 2G). Furthermore, EGR1 expression pointed to the transcriptional regulator function of this cluster. Particularly, we found 3 CD8+ T clusters with higher proportion in liver transplantation tissues. C9-HBB-CD8T holds plenty of unique genes consisting of HBB, HBA, and ALB, pointing out its specific functions, which needs to be studied further. C10-CTLA4-exhausted CD8T, holding 1495 cells, predominantly composed of cells from liver transplantation tissues, expressed high levels of exhaustion markers CTLA4, PDCD1, and LAG3. C11-MKI67-proliferative CD8T was labelled as the proliferative cluster for the high level of related factors including MKI67 and TOP2A (Figure 2G).

CD8+ T clusters showed comparable distribution among patients (Figures 2A, C, F). C4- regulatory-CD8T with CRTAM expression and C5-cytotoxic-CD8T with GZMH expression were prevalent in health tissue, while C10-exhausted-CD8T and C11-proliferative-CD8T cells were dominant in liver transplantation tissues. Interestingly, we performed pseudotime analysis to order CD8+ T cells in pseudotime to project the developmental trajectories and found that clusters showed relative time-process, which began with C11-proliferative-CD8T and ended with C10-exhausted-CD8T (Supplementary Figures S2F, G). Thus, exhausted CD8+ T cells were highly enriched in the period of the last phase, indicating the T cell state went from proliferation and exhaustion.

Similarly, we identified 3 major CD4+ T clusters (Figure 2A), which all expressed high levels of IL7R and LTB (Supplementary Table S4). The first cluster, C1-IL7R-CD4T, was significantly dominant in heath tissue, expressing the highest level of IL7R gene, which was identified as T helper cells (Figure 2F). Moreover, it especially expressed markers such as TNF and CD40LG, which indicated the activation vs. Type1 immunity. The second cluster, C6-CCR7-CD4T, showed high levels of LEF1, CCR7, and CD62L. CCR7 and CD62L have been shown to mediate the migration of memory T cells, thus this cluster was identified as central memory CD4+ T cells. It also showed the enhancer function for the LEF1. While the last CD4+ cluster, called C8-CCR6-CD4T, similarly expressed CCR7, was identified as another memory T cluster for the special expression of CCR7 and CCR6. With poor expression of CD4 and CD8, we found double negative T cell clusters in the graft tissues, C12-CCL3-DNT, whose marker genes were CCL3L1, IFIT2, and DNAJB4. We identified it as exhausted DNT cluster (Figures 2A, E). Two NK clusters (C0-XCL1-NK, C2-FGFBP2-NK) showed different distributions among tissue. The former, C0-XCL1-NK, accounting for the highest proportion of T, NK clusters, showed high levels of XCL1, AREG, and IRF8. While the latter cluster, C2-FGFBP2-NK, expressing FGFBP2, MYOM2, and FCGR3A, was dominant in liver transplantation tissues (Figures 2A, F and Supplementary Table S4).

B cells also play great roles in the process of liver transplant rejection. There were 2091 cells further re-clustered and showed functional profile in healthy and liver transplantation tissues. We matched 11 clusters and acquired the differentially expressing genes of each cluster (Supplementary Figures S2H, I and Supplementary Table S5). C0-RASSF6 expressed high levels of RASSF6, SSPN, and CD82. RASSF6 functions to encode an RAS effector protein which induces apoptosis. Thus, we named this cluster as C0-apoptosis-cluster. C1-ZBTB16, C2-FCRL3, C3-IL6, and C5 might show similar functions for the common high levels of the list of genes, including IGHD, YBX3, IL4R, and FCER2, of which FCER2 is a B cell specific antigen while both IL4R and FCER2 involve the regulation of IgE production. C1-ZBTB16 expressed genes related to cell cycle and cell division, including ZBTB16 and CENPM, simultaneously. C2-FCRL3 holds plenty of unique genes such as FCRL5, FCRL3, and ARL4D, which implicates B cell development. Furthermore, C4-IGHG showed a lot of IGHG family members: IGHG1, IGHG2, IGHG3, and IGHG4. The group of C5-HRK significantly differentially expressed HRK, CCDC191, NEIL1, and PCDH9. Notably, HRK encodes a family of BCL2 proteins, members of which are involved in promoting or inhibiting apoptosis. Thus, C5-HRK was a separate group of apoptosis from the C0-RASSF6. However, the proportion of these two groups of cells in liver transplantation tissue was decreased (Supplementary Figure S2J). Both C6-IL32 and C7-IGHA2 expressed T and NK cell related genes, including IL32, NKG7, CD7, GZMK, and GZMA, which are involved in cytotoxicity.

Compared with liver biopsy, peripheral blood is easier to obtain clinically, which has certain potential value in the occurrence and development of transplant rejection and tolerance (27). Thus, we similarly built the liver transplantation atlas of PBMC at the single cell level. By merging 4 healthy data (28), we obtained 46,808 high quality cells after quality control. Dimensionality reduction, cell cluster, and annotation turned out 22 clusters, including 30,465 (65.09%) T/NK cells, 12,557(26.83%) myeloid cells, and 3786 (8.09%) B cells (Supplementary Figures S3A–C). We performed further analysis by re-clustering each cell type and detailed annotation (Supplementary Figure S3D).

Re-clustering T and NK cells obtained 5 CD4+ T cell clusters, 8 CD8+ T cell clusters, and 2 NK clusters (Figures 3A, B, Supplementary Figure S3E, and Supplementary Table S6). We found CD8+ T cells were much more enriched in liver transplantation tissues while CD4+ T cells were relatively lower than those in normal tissue (Figure 3C), which was further confirmed by flow cytometry (Figure 3D). Similarly, we also identified cytotoxic CD8+ T cells, covering C7-GZMH-CD8T, C8-GZMK-CD8T, C9-KLRB1-CD8T, C10-CCL5-CD8T, and C14-GZMA-CD8T, using canonical markers including GZMA, GZMK, and GZMH (Supplementary Figure S3F). Interestingly, we found C7-GZMH-CD8T was significantly dominant in liver transplantation tissue, and it showed NK cell-related genes including NKG7 and FGFBP2. In addition, LAG3 was highly expressed in this cluster, pointing out the potential exhausted state (Supplementary Figure S3F). However, we did not find a proliferative CD8+ T cluster. We found C1-TSHZ2-CD4T (TSHZ2 and FHIT highly expressing) and C4-ITGB1-CD4T (ITGB1 and AQP3 highly expressing) in liver transplantation tissue were significantly lower than those in healthy tissue (Figure 3E). Importantly, C5-GZMB-NK (GZMB+ SPON2+) is much more prominent in liver transplantation tissue, pointing out that this cluster may play an important role in the process of liver transplant rejection and tolerance.

A total of 3786 B cells were further analyzed (Figures 3F–H). We also found several clusters had different proportions between liver transplantation and normal tissue (Figure 3I). C2-IL4R-B cell expressed high levels of IL4R, TCL1A, and H1FX, while C4-NEIL1-B cell holds increased expressing level of ALI39020.1, PCDH9, and NEIL1. The proportion of both clusters in liver transplantation tissues was significantly higher than that in healthy tissue. On the contrary, C3-AIM2-B cell and C8-IGLL5-B cell were dominant in healthy tissue. The former expressed S100A11, AIM2, and TLR10 genes while the later expressed IGLL5, IGJ, and HSPA6 genes.

The percentage of C10-exhausted-CD8T cells were increased significantly in liver transplantation graft tissues (Figure 2F), especially in rejection tissues (Figure 4D), which indicated that it may play a significant role in the period of immune tolerance and transplant rejection. Thus, we further explored this cluster, which showed extra markers, such as HAVR2, TIGIT, and DUSP4 (Figure 4A). In addition, we validated the high expression of PD-1 using immunohistology of liver transplantation tissues (Figure 4B). With identified differentially expressed genes of this cluster, we conducted pathway enrichment analysis. The results showed that these genes were involved in the regulation of T cell activation, T cell receptor signaling pathway, PD-L1 expression, and PD-1 checkpoint pathway (Figure 4C, Supplementary Figure S4A and Supplementary Table S4). The group of C10-exhausted-CD8T cells represented the highest proportion in the rejection group (Figure 4D). We confirmed this finding in the bulk RNA-seq, as the expression of exhausted CD8T cells marker genes such as CTLA4 and LAG3 in the rejection group were significantly higher than the non-rejection group (Figure 4E, Supplementary Figure S4B), confirming the exhausted CD8T cells increased in the rejection group. Moreover, CIBERSORT deconvolution algorithm was used to analyze immune infiltration according to bulk RNA-seq data, which showed that CD8+ T cells were significantly higher in rejection samples than those in stable samples (Supplementary Figure S4E).

To further delineate the structure of exhausted CD8+ T cells, we re-clustered this group and acquired 3 detailed clusters (Figure 4F). The top 5 marker genes of each cluster are shown in Figure 4G. C10-exhausted-CD8T were divided into three groups and cluster 3 expressed higher CD4 and FOXP3 than the others (Figure 4H). In addition, we confirmed this finding using multiplex immunofluorescence staining (Figure 4I). According to pseudotime analysis, we mentioned above that C10-exhausted-CD8T was highly enriched in the period of the last phase (Supplementary Figure S2F). Curious about cluster 2, called CD4+CD8+FOXP3+ T cells, we further speculated on the internal evolution of exhausted CD8 cells. Interestingly, we found CD4+CD8+FOXP3+ T were at the beginning of differentiation (Supplementary Figures S4C, D).

We identified 14 clusters by clustering the whole myeloid cells in liver tissue, which were annotated as dendritic cells (DCs), Kupffer cells (KCs), macrophage cells (TMs), and tissue monocytes (Monos) (Figures 5A–D and Supplementary Tables S8, 9). Clusters 7, 9, and 13 were identified as KCs because of the high expression of known markers including MARCO, TIMD4, CD5L, and VCAM1 (Figures 5A, E). C7-CD163-Kupffer expressed high levels of classical markers of macrophage such as CD163, which was dominant in health tissue. C9-APOE-Kupffer also holds the larger proportion in healthy tissue with expression of complement C1q family, including C1QA, C1QB, and C1QC. C13-GZMA-Kupffer showed potential cytolytic function, whose markers were GZMA and NKG7, similar to cytotoxic NK cells. Then, we performed pseudotime analysis to 3 KCs, which implied that C13-GZMA-Kupffer was at the end of differentiation (Figure 5F). KCs showed similar results to exhausted CD8+ T cells (Supplementary Figure S2F). We also found 3 clusters of conventional macrophage cells: C2-THBS1-Macro, C5-CDKN1C-Macro, and C6-FOLR3-Macro. It is worth noting that C2-THBS1-macro was expanded in normal liver significantly.

Clustering of myeloid cells acquired 5 clusters of DCs (Figures 5A, G). C0-CD1C-DC showed the largest expression of CD1C, and other markers were FCER1A, CD1E, and SPIB. The second cluster was C4-SDS-DC, with special markers including SDS and A2M. C10-CADM-DC and C11-CD141-DC had similar markers, such as CLEC9A, CADM1, and IDO1, while C10-CADM-DC was enriched in the expression of EGLN3, TACSTD2, and the late state (Figure 5H and Supplementary Figure S5C).

We found C3-LDLR-macro was expanded in liver transplantation tissues (Supplementary Figure S5A), so we further explored this cluster. Interestingly, we named it as a group of MDSC with high expression of TMEM176B, S100A8, and S100A9 (Figure 6A). The flow cytometry result for PBMC and tissue staining confirmed the result (Figures 6B, C). Enrichment analysis indicated that this cluster was involved in neutrophil activation and response to lipopolysaccharide (Figure 6C, Supplementary Figure S6A), which means plenty of marker genes of this cluster participated in lipid metabolism. In addition, this cluster expressed high levels of RETN and LDLR, which are thought to be strong lipid metabolism related genes. The percentage of LDLR+ MDSC was increased significantly in liver transplantation tissues, especially in rejection tissues (Figure 6D). MDSC was recognized as a key factor in inducing immune tolerance and inhibiting transplant rejection. Bulk RNA-seq confirmed MDSC markers, S100A8 and S100A9, were significantly higher in liver transplant rejection tissues than those in stable tissues (Figure 6D). We demonstrated these results using immunohistochemistry and fluorescence (Figure 6E).

After defining the population of T, NK cells, and myeloid cells, we performed cellular communication analysis within liver tissues using CellChat, to infer the intercellular communication network (Figure 7A). We constructed the immune populations’ interactions with outgoing and incoming pathways, between ligands, receptors, and cofactors (Supplementary Table S10). We visualized the dominant senders and receivers using scatter plot, which unfolded that KC_1, DC_4, and KC_2 showed the top 3 outgoing interaction strengths while CD8T_1, CD8T_2, and CD8T_6 showed the top 3 incoming interaction strengths (Figure 7B). Further, to explore how multiple cell types and signaling pathways coordinate, we identified global communication patterns (Supplementary Figure S7A). Two patterns of secreting cells basically distinguished T cells and myeloid cells, which were involved in different pathways individually. Three patterns of target cells emerged, noting that DC_2, KC_1, and KC_2 gathered as the same pattern with specific signals.

Having proved that C10-CTLA4-exhausted CD8T was dominant in liver transplantation tissue, with special characteristics, we further explored the interactions between C10-CTLA4-exhausted CD8T and the whole myeloid cells (Figure 7C). Importantly, we found the TIGIT-NECTIN2 pathway particularly existed in C10-CTLA4-exhausted CD8T and C3-LDLR-MDSC (Figure 7C). C10-CTLA4-exhausted CD8T showed the most powerful outgoing strength in the TIGIT signaling pathway network while LDLR+MDSC placed second to receivers, and TIGIT-NECTIN2 contributed most of TIGIT signal’s ligand-receptor pair (Figure 7D). Meanwhile, two clusters showed similar results in the NECTIN signaling pathway network (Supplementary Figure S7B). Moreover, we visualized the signaling gene expression distributed among clusters, showing TIGIT and NECTIN2 expressed high levels in C10-CTLA4-exhausted CD8T and C3-LDLR+MDSC separately (Figure 7E).

Next, we found the TIGIT gene was highly expressed in T cells from liver transplantation tissue while NECTIN2 genes were dominant in myeloid cells from liver transplantation tissue (Figure 7F). In addition, we validated this finding using the immunohistochemical technique (Figure 7G). Bulk RNA-seq confirmed TIGIT and NECTIN2 are significantly expressed higher in liver transplant rejection tissues than that in stable tissues (Figure 7H). In summary, our exploration of the key ligand-receptor interactions between C10-CTLA4-exhausted CD8T and C3-LDLR+MDSC mentioned the TIGIT-NECTIN2 pathway as an important regulator in the process of transplant rejection and induced immune tolerance.