We downloaded mRNA expression profiles of human AR from the Gene Expression Omnibus (GEO) database. In our study, GSE138043 was used to construct co-expression networks and identify hub genes related to AR. The microarray dataset provided gene expression profile in the percutaneous allograft biopsy from 15 AR patients and 37 NAR (25). According to the data processing information of GSE138043, each dataset was normalized independently using Robust Multiarray Average (RMA) followed by log2 transformation and quantile normalization. Data from GSE50058 and GSE343 were used for hub genes validation. In the GSE50058 dataset, 42 AR patients and 58 STA individuals were recruited, and the RNA was extracted from their renal allograft biopsy. In the GSE343 dataset, the total RNA was extracted from the kidney tissue of 25AR patients and 15 NAR. Supplementary Table S5 shows the summary of discovery and validation microarray data sets of clinical biopsy samples from kidney transplants.

“limma” R package was utilized to the differentially expressed screen genes (DEGs) between AR and NAR in the expressing data. The genes with adjusted p-value <0.05 were selected as having significant change. “ggplot2” and “pheatmap” were used respectively to paint the volcano plot and heatmap of all DEGs.

The co-expression network of the genes was constructed based on GSE138043 microarray dataset by the R package “WGCNA.” The soft-thresholding power that we chose was 17 when 0.9 was used as the correlation coefficient threshold. We defined 0.25 as the threshold for cut height to merge possible similar modules.

To obtain further insights into the function of the target module most related to AR, we referred to the Database for Annotation, Visualization and Integrated Discovery (DAVID) (https://david.ncifcrf.gov/) to perform the KEGG enrichment analysis. The results were shown graphically by the R package “ggplot2.”

The green module, which was most significantly related to AR, was imported into Cytoscape with their weighted correlations. We identified the hub gene with the following criteria: (1) DEGs in green module; (2) gene significance (GS) > 0.4 and module membership (MM) > 0.9; (3) and top 10 Maximal Clique Centrality (MCC)value calculated by the Cytohubba package in Cytoscape v3.8.2.

A 10X Genomics Chromium machine was used for single-cell capture and cDNA preparation following manufacture’s instruction. Chromium™ Single Cell 3′ Solution was used to perform reverse transcription on gel bead in emulsion, followed by cDNA cleanup and amplification. The cDNA is digested and broken into fragments of about 200–300 bp, followed by the traditional second-generation sequencing library construction process, and PCR amplification is performed to obtain a DNA library. Illumina sequencing platform of paired-end sequencing mode was used to perform high-throughput sequencing on the established library. Sequence data were processed with Cell Ranger V2.1.0 (10X Genomics).

Then, quality control was performed to filter low-quality cells. For 10X-derived datasets, we only retained cells that had (1) genes more than 200 and <6,000, (2) UMIs more than 500 and <40,000, and (3) <15% of reads mapped to mitochondrial genes.

The Nephroseq v5 online database (http://v5.nephroseq.org/), an integrated data-mining platform for gene expression data sets of kidney diseases, was adopted to validate the correlation between the hub genes and clinical manifestations of AR by Spearman rank correlation coefficient analysis. A p-value of <0.05 was considered statistically significant. r>0.6 was considered strong correlation, and 0.4<r<0.6 was considered medium intensity correlation.

The xCell (https://xcell.ucsf.edu/), which is a webtool that performs cell-type enrichment analysis from gene expression data for 64 immune and stroma cell types, was adopted to reveal different infiltrating cell types between AR and NAR, and adjusted p-value < 0.01 was chosen as the cutoff criterion. Cell-type enrichment score is shown in Supplementary Table S3. To further explore the six hub genes expression in 76 single cell types, we obtained RNA expression values per cell types from Human Protein Atlas Dataset (proteinatlas.org). Single cell-type clusters were normalized separately from other transcriptomics datasets using trimmed mean of M values (TMM). To generate expression values per cell type, clusters were aggregated per cell type by first calculating the mean nTPM in all cells with the same cluster annotation within a dataset.

To further explore the potential function of the selected hub genes in AR, we used gene set enrichment analysis (GSEA_4.1.0) for single hub gene. In the dataset GSE138043, samples were divided into two groups according to the median expression level of hub genes. The h.all.v6.2.sytmbols.gmt in Molecular Signatures Database (MSigDB) was selected as the reference gene set, and adjusted p-value <0.05 was considered significantly different. The results were showed graphically by the R package “ggplot2.”

The Raw264.7 cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM) High Glucose (11965084, Gibco) containing 10% fetal bovine serum (FBS) (12103C, Sigma) and 1% penicillin–streptomycin (15070063,Gibco), incubated at 37°C in a humidified atmosphere containing 5% CO₂, and routinely passaged every 1 or 2 days. To induce RAW264.7 cell line to M1/M2 phenotypes, 10⁵ Raw264.7 cell were seeded in six-well plates 24 h before exposed to IFN-γ (20 ng/ml) + lipopolysaccharide (LPS) (10 ng/ml) to M1 phenotype, and IL-4 (10 ng/ml) to M2 phenotype for 24 h.

Isolation and culture of bone-marrow-derived macrophages (BMDMs) were described by Pineda-Torra et al. (26). In brief, bone marrow was flushed out from the femurs and tibias, cultured and differentiated for 7 days in 1640 supplemented with 10% heat-inactivated FBS and 1% penicillin–streptomycin with macrophage-stimulating factor (M-CSF) (50 ng/ml) at 37°C in a humidified incubator with 5% CO₂.

Flow cytometry was used to measure the phenotypical changes in Raw264.7 macrophages. Single-cell suspension is prepared before staining with fluorochrome-labeled anti-CDF4/80 (clone BM8, BioLegend, San Diego, CA, USA), anti-CD80 (clone 16-10A1, eBioscience, San Diego, CA, USA). Data were analyzed using FlowJo v.10 (Treestar, Ashland, OR).

Research involving human participants was reviewed and approved by the Research Ethics Committee of the First Affiliated Hospital of Zhejiang University School of Medicine. Patients/participants (or their close relatives) provided written informed consent to participate in this study. Peripheral blood mononuclear cell (PBMC) was isolated within 3 h after collection. Percutaneous allograft biopsy was collected and fixed with 4% formalin for paraffin embedding. Biopsies were scored by the revised Banff 2019 classification of renal allograft pathology; rejection cases here were TCMR including borderline cases. Supplementary Table S6 showed the characteristics of the samples used in this study.

Total RNA was extracted from cultured cells and peripheral blood mononuclear cells (PBMCs) by Trizol reagent (Invitrogen, CA, USA). cDNA was prepared using the PrimeScriptTM RT Reagent Kit with gDNA Eraser (No. RR047A, Takara, Shiga, Japan) following the manufacturer’s protocol. Real-time PCR was run using SYBR Green and CFX96™ Real-Time PCR Detection Systems (Bio-Rad, CA, USA). The mRNA levels of selected genes were calculated after normalization to β-actin by using the 2^−ΔΔCt method according to the manufacturer’s protocol. All primer sequences used are shown in Supplementary Table S4.

Immunohistochemistry (IHC) was performed following standard protocol. Briefly, after being dewaxed and rehydrated, the 2-μm paraffin-embedded sections were incubated with Anti-BLAME Polyclonal Antibody (bs-2473R, Bioss, Boston, MA, USA). After washing, the sections were incubated with the horseradish-peroxidase-labeled anti-mouse/rabbit IgG polymer (GK500710, GeneTech, South San Francisco, CA, USA) and diaminobenzidine. The sections were then counterstained with hematoxylin, dehydrated, and cleared. Six random fields of each section were photographed, and the staining was semi-quantified using the National Institutes of Health Image J by an investigator blinded to the experimental protocol.

Coverslips containing RAW264.7 or BMDMs were fixed with 4% paraformaldehyde for 15 min and blocked in phosphate-buffered saline (PBS) containing 0.1% Triton X-100 and 3% bovine serum albumin(BSA) for 30 min at room temperature prior to incubation with Anti-BLAME Polyclonal Antibody (bs-2473R, Bioss), TLR4 antibody (sc-293072, Santa Cruz) overnight in a humidified chamber at 4°C. Slips were incubated with secondary antibody Alexa Fluor 594-conjugated anti-rabbit IgG (1:500) and Alexa Fluor 488-conjugated anti-mouse IgG (1:500) 1 h at 37°C. Sections were then examined by immunofluorescence microscopy (Leica DMLB, Wetzlar, Germany).

Total RAW264.7 cells were lysed in radioimmunoprecipitation assay (RIPA) cell lysis buffer following separated by sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE) and transferred to 0.22-μm polyvinylidene fluoride membranes (Millipore, Burlington, MA, USA). Blocked with 5% milk and incubated with the inducible nitric oxide synthase (iNOS) antibody (ab178945, Abcam, Cambridge, UK) at 4°C overnight. Protein was visualized using secondary anti-rabbit IgG or anti-mouse IgG (Sigma) with conjugated horseradish peroxidase and chemiluminescent substrate (Millipore, Billerica, MA, USA).

The statistical significance differences between the two groups were analyzed using non-parametric test or t-test based on data distribution characteristics. All analyses were conducted using software R4.1.0. A p-value < 0.05 was considered statistically significant.

Kidney transcriptome data from GSE138043 containing 15 AR patients (rejection at 12 months post-renal transplant) and 37 NAR (non-rejection at 12 months post-renal transplant) was used for further analysis. We identified 1,561 DEGs (differential genes, adjusted p-values < 0.05) between the AR and NAR groups. Among 1,561 DEGs, 541 genes were upregulated, while 1,020 genes were downregulated. All genes are displayed in volcano plot in Figure 1A, and DEGs are listed in Supplementary Table S1. Additionally, the red plots represent adjusted p < 0.01, orange plots represent log2FC > 1, and green plots that were annotated represent both. Unsupervised clustering hierarchy was used in heatmap (Figure 1B).

WGCNA package was applied to compile the network. Keeping to the scale-free topology criterion, β = 17 was considered in this study for which the fit index curve flattens out upon reaching a high value (>0.9) and the mean connectivity ≤100 (Figure 2A). Hierarchical clustering was used to generate a hierarchical clustering dendrogram of genes; meanwhile, Dynamic Tree Cut R library was used for detecting clusters. As shown in Figures 2B–D, 12 distinct gene modules (M1–12) were defined, as MEDissThres was set to 0.25 to merge similar modules. Genes failing to fit within a distinct group were assigned to the gray module. The interaction relationships of 1,000 randomly selected genes were presented in the network heatmap (Supplementary Figure S1A).The genes in the same module were highly correlated, while they were weakly correlated to those in other modules. Thus, the reliability of the modules was verified.

Module-trait correlations analyses showed that multiple modules were related to AR (Figure 3A) The summary of significance of all genes in each module related to AR is shown in Figure 3B. Clearly, the green module was most significantly related to AR, followed by the turquoise module. Figure 3C shows the significance of these genes in the green module for AR (cor = 0.35, p = 2.4e−22). To investigate the potential biological function of genes in the green modules, we performed KEGG enrichment analysis. The top 20 pathways terms of green module are shown in Figure 3D. We found that genes in green modules played roles in cytokine–cytokine receptor interaction, chemokine signaling pathway, cell adhesion molecules (CAMs), tuberculosis, and other important regulatory processes. To identify hub genes in the target module, we calculated the MCC values via Cytohubba and constructed a network based on the top 10 genes (Supplementary Table S2). The module was visualized using Cytoscape 3.8 software (Figure 3E). The node colors coded from yellow to red (low to high) correspond to the top 10 MCC values from low to high. In intersection with GS > 0.4, MM > 0.9, and DEGs in the green module, six genes were regarded as the hub genes(Figure 3F), including dedicator of cytokinesis 2 (DOCK2), NCK-associated protein 1 like (NCKAP1L), interleukin 2 receptor subunit gamma (IL2RG), SLAM family member 8 (SLAMF8), CD180 molecule (CD180), and protein tyrosine phosphatase receptor type E (PTPRE). Relative mRNA expression of six hub genes in GSE138043 is shown in Figure 1B.

As expected, the expression levels of hub genes including DOCK2, NCKAP1L, IL2RG, SLAMF8, CD180, and PTPRE were significantly upregulated in AR samples from the GSE50058 dataset (Figure 4A). For verifying hub genes, we obtained another dataset, GSE343, and analyzed the expression levels of the above five genes except for SLAMF8, which was not found in this dataset between AR patients and STA patients (Figures 4B). The expression of these genes was also upregulated in AR samples. We used the Nephroseq v5 online database to explore the correlation between the expression of IL2RG and clinical traits of AR. As shown, there was a positive correlation between the expression of IL2RG in AR with the Banff pathological grading of transplanted kidney (r=0.6540, p<0.0001) (Figure 4C).Correlation between the expression of DOCK2 in AR and the Banff pathological grading is shown in Supplementary Figure S1C (r=0.4863, p<0.0035).

To study the transcriptional changes related to AR further, we performed single-cell RNA sequencing in P1 (patient PBMCs of acute rejection post renal transplant) and C1 (control patient of stable kidney function post renal transplant). The patient had a higher percentage of native T, CD1C+_B DC, NKT, NK, and monocytes in P1 compared with C1 (Figures 5A–C). Expression of DOCK2, NCKAP1L, IL2RG, SLAMF8, CD180, and PTPRE (colored single cells) on UMAP plot projecting PBMCs from P1 and C1 are shown in Supplementary Figure S2. In addition, to further validate these six hub genes, we collected PBMCs from eight non-AR and 10 AR patients to perform qRT-PCR. The results showed that compared with NAR group, the mRNA levels of DOCK2, NCKAP1L, IL2RG, SLAMF8, CD180, and PTPRE were all significantly elevated in AR patients (Figures 5D–I).

Cell-type enrichment analysis was performed by xCell from gene expression data for 64 immune and stroma cell types. Our data revealed that the AR group had higher immune scores and microenvironment scores than the NAR group. Among the 64 cell types, scores of 20 types of cells were significantly differently expressed (p<0.01), including immune cells, such as B cells, basophils, CD4+ Tem, CD4+ memory T cells, CD8+ T cells, CD8+ Tcm, class-switched memory B cells, DC, macrophages, M1, mast cells, memory B cells, monocytes, plasma cells, aDC, cDC, naive B cells, pDC, and stroma cells such as endothelial cells, smooth muscle, and MV endothelial cells (Figure 6). We used the Protein Atlas online database to explore the six hub genes RNA expression in 76 single-cell types (Supplementary Figure S3). Except for SLAMF8, which is more singly expressed in Langerhans cell and macrophages, the other five genes are more widely expressed in immune cells such as B cells, DC, T cells, and monocyte consistent with the different cell types in AR groups, meaning that six hub genes exercise certain biological properties in these cell types during AR.

We performed GSEA, and h.all.v7.4.sytmbols.gmt in MSigDB was used as the reference gene set. The full list of gene sets enriched in samples with DOCK2 (Figure 7A), NCKAP1L (Figure 7B), IL2RG (Figure 7C), SLAMF8 (Figure 7D), CD180 (Figure 7E), or PTPRE (Figure 7F) highly expressed is shown. Four gene sets were enriched in samples with highly expressed DOCK2, NCKAP1L, IL2RG, SLAMF8, CD180, and PTPRE, including “allograft rejection”, “interferon γ response”, “interferon α response”, and “inflammatory response”. Moreover, gene sets “IL6–JAK–STAT3 signaling”, “KRAS signaling”, and “TNFα signaling via NFkB” were enriched in the samples with either DOCK2, NCKAP1L, IL2RG, SLAMF8, CD180, and PTPRE highly expressed.

SLAMF8, one of the six hub genes, was also the top rejection-associated transcripts in TCMR versus everything else including ABMR in previous study (27). In addition, SLAMF8 is more specifically expressed in macrophages, which is exactly the cell type that differentially infiltrated in the AR group in our above research results (Figure 6; Supplementary Figure S3D). Immunohistochemical staining of kidney tissue revealed a higher level of SLAMF8 expression in the renal interstitium from TCMR (T-cell-mediated acute rejection after renal transplantation) kidney tissue than from HC (healthy donor control) (Figures 8A, B). A previous study confirmed that SLAMF8 maintained TLR4 expression on macrophages and promoted LPS-induced mitogen-activated protein kinase (MAPK) activation (28). We subsequently determined whether the SLAMF8 participates in macrophage activation via TLR4. Fluorescence detection of SLAMF8 and TLR4 revealed that approximately complete SLAMF8+ cells express TLR4 (Figure 8C). Thus, we speculated that SLAM8 participates in AR progression, likely through upregulating TLR4 expression.

In order to explore in which phenotype of macrophages SLAMF8 is specifically expressed, murine RAW 264.7 macrophages and primary BMDM cells were treated with 10 ng/ml LPS plus 20 ng/ml IFNγ or 10 ng/ml IL-4 for 24 h. CD80 and CD206 are known to be a specific surface marker of the M1 phenotype and the M2 phenotype, respectively. Flow cytometry and RT-qPCR were used to measure the polarization of RAW 264.7 macrophages (Supplementary Figures S4A–D). Western blotting using the M1 phenotype marker iNOS was confirmed to be increased in the LPS plus IFNγ treated group in RAW264.7 (Supplementary Figure S4E). The mRNA expression of SLAMF8 and TLR4 was increased in LPS and IFNγ treated RAW 264.7 macrophages and BMDMs (Figures 9A, B; Supplementary Figures S5A, B) consistent with the phenomenon in immunofluorescence staining (Figures 9C–E; Supplementary Figures S5C–E). Two-color immunofluorescence analysis revealed that the M1 phenotype contained large numbers of TLR4-expressing SLAMF8+ cells in contrast to M0- and M2-type macrophages (Figures 9C–F; Supplementary Figures 5C–E).