EU-TRAIN (NCT03652402) and KTD-Innov (NCT03582436) studies are large, prospective multicenter cohorts that follow adult kidney transplant recipients for 1 year after transplantation. They involve collaboration between transplant centers, analytical platforms, and industrial partners across France and Europe.

The studies focus on adult patients (18 years or older) receiving a living or deceased kidney transplant. Participants must be willing to comply with study procedures and signed an informed consent. Patients with a history of multi-organ transplants, language barriers hindering participation, or vulnerability (minors, pregnant women, etc.) were excluded.

Both EU-TRAIN and KTD-Innov involve baseline visits at the time of transplant, followed by checkups at 3- and 12-months post-transplantation. Additional visits may be scheduled if a patient’s kidney function deteriorated or protein levels raised. KTD-Innov recruited participants between July 2018 and December 2019 and focused on seven French transplant centers (Paris-Necker, Paris-Saint-Louis, Nantes, Bordeaux, Toulouse, Lyon, and Montpellier). The EU-TRAIN study, with a slightly broader enrollment window (November 2018–June 2020), encompasses nine centers across Europe (Paris-Saint-Louis, Paris-Necker, Nantes, Barcelona-Bellvitge, Barcelona-Vall d’Hebron, Berlin-Charité Mitte, Berlin-Charité Virchow, Geneva, Paris-Kremlin-Bicêtre). 770 renal biopsies were collected from 540 patients from the two prospective studies as well as two retrospective cohorts from Necker and St Louis hospitals (Paris, France) between 2006 and 2021. This study was approved by local institutional review boards and written informed consent was obtained from all patients.

Lesions from biopsies were graded by local renal specialist from 0 to 3 according to the 2019 international Banff classification [13], and comprised: glomerulitis (g), peritubular capillary inflammation (ptc), interstitial inflammation (i), tubulitis (t), total inflammation (ti), endarteritis (v), transplant glomerulopathy (cg), interstitial fibrosis (ci), tubular atrophy (ct), vascular fibrous intimal thickening (cv), arteriolar hyalinosis (ah). C4d staining was performed by immunohistochemistry on paraffin sections using the human C4d polyclonal antibody. C4d staining was graded from 0 to 3 by the percentage of peritubular capillaries with linear staining. Earlier biopsies were reclassified to take into account the evolution of the classification.

The presence of circulating donor-specific anti-HLA-A, -B, -Cw, -DR, -DQ and -DP antibodies was analyzed using single-antigen bead assays (One Lambda, Inc., Canoga Park, CA, United States) on a Luminex platform on serum samples collected at the time of transplantation and at the time of biopsy. For each patient, we recorded the number, class, specificities and mean fluorescence intensity (MFI) of all donor-specific HLA antibodies. Positiveness of a DSA was defined by a threshold of 500 for the mean fluorescence intensity. The maximum MFI for the immunodominant DSA (Anti-HLA iDSA MFI) was defined as the highest ranked donor-specific bead. HLA typing of donors and recipients was performed using DNA typing.

After collection, all biopsies were stored in the RNAlater^® solution at −20°C. They were then centralized at the Paris Cardiovascular Research Center (PARCC) in order to be processed by the Paris Transplant Group Precision Pathology Platform for total RNA extraction using the Promega^® Maxwell^® RSC miRNA Tissue Kit [14]. All samples were selected according to a minimal concentration of RNA of 20 ng/μL and an RNA integrity number superior or equal to 7. Purified RNAs were, then, stored in a −80°C fridge while waiting to be sent and sequenced by the GENOM’IC platform at Cochin hospital where the library was prepared according to the Illumina^® Stranded mRNA Prep Ligation protocol [15] with a capture of the polyadenylated RNAs using oligo (dT) magnetic beads. Finally, an Illumina sequencing has been performed in order to obtain 2 × 30 millions paired-end reads on average.

After the sequencing, we used FastQC (v0.11.9) [16] to assess the pre-alignment quality controls. We performed the alignment with the STAR algorithm (v2.7.4a) [17] and the Hg38.p13 reference genome. We finally verified its quality with STAR, Picard tools (v 2.22.9) [18] and RSeQC (v3.0.1) [19] metrics. Raw counts have been generated using the featureCounts program with the GC_000001405.39_GRCH38.p13 GTF annotation file and the BAM files resulting from the alignment. Quality controls results can be found in the Supplementary Table S1.

The identification of differentially expressed genes (DEGs) was performed using the DESeq2 method (v1.30.1) [20]. Gene expression count matrix has been pre-filtered by removing lowly-expressed genes using a threshold of at least 1 Fragment Per Kilobase Million (FPKM) in 20% of the total samples for each gene. The number of filtered genes reduced from 44,613 to 15,563. Fold changes (FC) and Wald statistics were computed for each comparison of interest with a correction for multiple hypothesis testing (Benjamini-Hochberg) and genes were ranked according to increasing adjusted p-values.

Two differential gene expression analysis were conducted including antibody-mediated rejection (AMR) and T-cell-mediated rejection (TCMR). Each diagnosis was tested against all histopathological diagnoses available in the cohort to obtain a molecular signature specific for the diagnosis of interest. This control group included all biopsies diagnosed with either TCMR or AMR (according to the design), isolated interstitial fibrosis and tubular atrophy, acute tubular necrosis, polyomavirus-associated nephropathy, thrombotic microangiopathy, recurrent or de novo glomerulonephritis, calcineurin-inhibitor toxicity or biopsies with no evidence of specific lesions (i.e., normal biopsies). Missing information, borderline (N = 69), mixed (N = 20) and suspicious rejection (N = 12) samples have been excluded from both designs. No threshold on the log₂ fold change was applied and all significant (adjusted p-value <0.05) differentially-expressed genes were considered during the analysis.

The complete description of differentially expressed gene symbols, mean expressions, log₂ fold changes, standard errors and Wald statistics as well as descriptions of the genes previously described in gene panels (B-HOT or microarrays) are shown in the Supplementary Material.

Pathways analysis was performed using both Reactome and the Kyoto Encyclopaedia of Genes and Genomes (KEGG) repositories with ReactomePA (v1.34.0) [21] and clusterProfiler (v3.18.1) [22], respectively, by either choosing as an input the entire list of upregulated genes (Reactome and KEGG) or a subset consisting of the upregulated transcripts included in the B-HOT or the microarray gene panel (Reactome only). Raw p-values were corrected for multiple hypothesis testing with the Benjamini-Hochberg FDR controlling technique and two cut-offs were applied to filter non-significant results: threshold of 0.05 on the adjusted p-value and 0.2 on the q-value. Q-values correspond to the proportion of false positive results in a set of signaling pathways that are at least as significant (adjusted p-value) as the signaling pathway under consideration. While the adjusted p-value gives the expected false positive rate, the q-value gives the expected positive false discovery rate. Pathway names, annotations and statistics are reported in the Supplementary Material.

Continuous variables were described by using means and standard deviations or medians and interquartile ranges. All analyses were performed using R (version 4.0.5, R Foundation for Statistical Computing). Values of p < 0.05 were considered significant, and all tests were 2-tailed.

The study cohort comprised a total of 770 kidney allograft biopsies from 540 patients collected between 2006 and 2021 from 11 international European centers (See Supplementary Material). Baseline characteristics including recipient and donor characteristics are summarized in Table 1. The population was mainly composed of men (n = 336, 63.2%) with a mean age of 51.1 ± 15.9 years at the time of transplantation, a history of glomerulonephritis as end stage renal disease (n = 125, 23.2%) and no history of a prior kidney transplant (n = 435, 80.9%). The majority of the transplantations were performed from deceased donors (n = 439, 81.5%) with 185 (42.1%) exhibiting expanded criteria. In total, 141 (27.9%) patients had pre-existing anti-HLA DSA.

The median time from transplantation to the biopsy was 3.9 months (IQR: 3.0–12.1) (Supplementary Figure S1) with 460 (60%) protocol biopsies and 370 (40%) for cause biopsies. The mean number of biopsy per patient was 1.4 (median = 1), with a maximum of 5 biopsies per patient. The mean eGFR at the time of the biopsy was 42.8 ± 19.1 mL/min/11.73 m² with a mean proteinuria of 0.53 ± 1.58 g/g. One-third of the patients (n = 226, 31.3%) had positive anti-HLA DSA with the immunodominant DSA belonging mainly to the class II (n = 148, 67.6%) (Table 2).

Kidney allograft biopsies were either classified as normal (n = 152, 20.6%) or had histological evidence for one or multiple of the following diagnoses: T-cell mediated rejection (n = 72, 9.9%), antibody-mediated rejection (n = 88, 12.0%), mixed rejection 14 (10.6%), borderline rejection (n = 69, 9.4%), interstitial fibrosis and tubular atrophy (n = 365, 49.6%), recurrence of the initial nephropathy (n = 15, 2.0%), de novo or recurrent glomerulonephritis (n = 18, 2.4%), acute tubular necrosis (n = 68, 9.1%), polyomavirus-associated nephropathy (n = 30, 4.0%), calcineurin inhibitors-related toxicity (n = 65, 8.7%), and thrombotic microangiopathy (n = 34, 4.6%) (Table 2).

60 AMR were compared to 576 non-AMR samples, resulting in 6,141 differentially expressed genes (DEGs). 358 (5.8%) were included in the Banff Human Organ Transplant (B-HOT) gene panel, 5,180 (84.4%) were included in the microarray gene panel, and 603 (9.8%) were new and defined as NGS-specific transcripts (Figure 1 and Supplementary Material). Genes included in the microarrays were highly represented throughout the entire molecular signature (from 0.0% to 84.4% among the increasing top ranked genes and stabilizing at top 2,500 genes), while the B-HOT-related genes were mainly ranked in the top genes (from 100.0% to 5.8%%, reaching 47.0% among the top 100 genes, 25.4% among the top 500 genes and 17.0% among the top 1,000 genes), and the new NGS-specific genes were constantly comprised across the signature between 5.2% and 12.5% of the total (Supplementary Figure S2). Among the top 30 ranked genes, 20 genes (66.7%) were included in the B-HOT gene panel including PLA1A, GBP1/4, GNLY, CCL4, IL15, IDO1, CXCL10/11, 7 (23.3%) genes were included in the microarrays panel (WARS1, GJD3, CLEC1A, CHN1, APOL3, SQLE, LILRA1), and 3 genes (10%) were specific to the NGS gene panel with CCL4L2, PELATON (a long non-coding RNA) and GBP1P1 (Supplementary Table S2).

The list of upregulated and differentially expressed genes was composed of 2,876 genes from which 2,299 (79.9%) were included in the microarrays, 313 (10.9%) were included in the B-HOT panel, and 264 (9.2%) were NGS-specific. Pathway analysis was performed using the Reactome repository. Top ranked (adj.p-value<0.05) known pathological categories were related to immune response: interferon signaling (q-value = 1.78e⁻¹¹), neutrophil degranulation (q-value = 4.57e⁻¹¹), signaling by interleukins (q-value = 4.11e⁻¹⁰), Toll-like receptors cascades (q-value = 6.80e⁻⁰⁶), class I MHC mediated antigen processing and presentation (q-value = 1.37e⁻⁰⁵), Fc Gamma/Fc Epsilon receptors (q-value = 2.18e⁻⁰⁴ and q-value = 2.15e⁻⁰⁵, respectively), signaling by the BCR (q-value = 5.40e⁻⁰⁵), cell surface interactions at the vascular wall (q-value = 2.11e⁻⁰⁴) and PECAM1 interactions (q-value = 8.14e⁻⁰³). In addition, the TCR signaling (q-value = 1.78e⁻¹¹), the PD-1 signaling (q-value = 3.50e⁻⁰⁸), the CD28 co-stimulation (q-value = 7.43e⁻⁰³) and the CTLA4 inhibitory signaling (q-value = 1.78e⁻⁰²) pathways were found significantly enriched in the AMR signature (Figure 2; Supplementary Figure S3 and Supplementary Material). Enrichment analysis derived from the KEGG database demonstrated additional significant pathways including the NK cell mediated cytotoxicity (q-value = 2.44e⁻¹⁰), Th17 cell differentiation (q-value = 2.98e⁻¹⁰) and Th1 and Th2 cell differentiation (q-value = 4.97e⁻⁰⁷), and provided access to the entire set of cytokine and receptors (CCL4, CCL11, CXCL5/6/9/10/11, XCL2, IL15/16/27/34/35, TNF, TGFβ) and cell adhesion and endothelium-related molecules (CD58, MHCI/II, CD40, ITGA, CD2, CD4, PD-L1, CDH5, PECAM1) involved during antibody-mediated rejection (Supplementary Figures S4 and S5 and Supplementary Material).

The analysis of enriched pathways restrained to the different panels highlighted specific functions. The B-HOT panel captured all the above-mentioned significant functions with a total of 191 entries in Reactome (Supplementary Figure S6 and Supplementary Material) while the microarray panel was specifically enriched in SUMOylation processes, RHO/RAC GTPase cycle, cell cycle progression and FCGRIIIA-mediated phagocytosis with only 44 entries (Supplementary Figure S7 and Supplementary Material). Finally, despite its 264 upregulated DEGs, the NGS-specific genes were only enriched in 4 non-specific metabolic functions (Supplementary Figure S8).

48 TCMR were compared to 589 non-TCMR samples and the molecular signature was defined by 8,478 genes. 439 (5.2%) were included in the B-HOT panel, 6,853 (80.8%) were included in the microarrays and 1,186 (14.0%) were NGS-specific (Figure 3). After ranking genes by their adjusted p-value, the proportions of transcripts included with each gene panel were mostly in favor of the microarray panel (from 100.0% to 80.8% with a local minimum of 60.8% among the top 265 genes). The proportion of genes included in the B-HOT first increased to reach a maximum of 30.5% among the top 118 genes, before decreasing to reach a minimum of 5.1%. Except among the top 5 genes, the newly discovered NGS genes were relatively stable (between 7.6% and 18.1%) (Supplementary Figure S9). Among the top 30 ranked genes, 22 (73.4%) were comprised in the microarray gene panel, 4 (13.3%) were comprised in the B-HOT panel (CD72, LAG3, CD8A, CD28) and 4 (13.3%) were NGS-specific (ANKRD23, TSPOAP1-AS1, LOC374443, MIR3142HG) (Supplementary Table S3).

The list of upregulated differentially expressed genes was composed of 4,482 genes from which 3,612 (80.6%) were included in the microarrays, 367 (8.2%) were included in the B-HOT panel and 503 (11.2%) were NGS-specific. Using the entire list of upregulated genes, significantly immunological Reactome enriched pathways comprised pathways related to: interferon signaling (q-value = 1.76e⁻²²), signaling by ROBO receptors (q-value = 2.66e⁻²²), TCR signaling (q-value = 2.40e⁻¹⁷), class I MHC mediated antigen processing and presentation (q-value = 1.88e⁻¹⁶), signaling by interleukins (q-value = 2.42e⁻¹⁶), signaling by the BCR (q-value = 8.04e⁻¹³), Fc Epsilon receptor signaling (q-value = 1.79e⁻¹¹) and Fc Gamma receptor dependent phagocytosis (q-value = 1.09e⁻⁰⁴), TLR cascades (q-value = 9.18e⁻¹⁰), co-stimulation by the CD28 family (q-value = 3.88e⁻⁰⁹), PD-1 signaling (q-value = 4.14e⁻⁰⁹) and neutrophil degranulation (q-value = 6.69e⁻⁰⁹). Out of the 466 Reactome entries, emphasis was also given on nonsense mediated decay and maturation of mRNA functions, SUMOylation processes, metabolism of non-coding RNA and cell cycle progression (Figure 4 and Supplementary Material). The KEGG repository significantly presented enrichment of the Th17 cell differentiation (q-value = 4.72e⁻¹⁴), Th1 and Th2 cell differentiation (q-value = 8.13e⁻¹²), and NK cell mediated cytotoxicity (q-value = 1.10e⁻⁰⁸). A wider range of activation/inhibition of cell adhesion molecules was also presented with lower/higher fold changes (min = −1.8, max = 3.9) compared to the AMR signature. The TCMR signature included the addition of CD22, PDCD1 and SELL and the inhibition of a multitude of molecules at the surface of the endothelial cells (Supplementary Figure S10 and Supplementary Material).

Focusing on the different gene panels, the genes included in the B-HOT panel captured all the immunological functions described previously (Supplementary Figure S11 and Supplementary Material) while the microarray genes specifically captured the nonsense mediated decay, SUMOylation, translation and mRNA maturation processes and the cell cycle progression (Supplementary Figure S12 and Supplementary Material). Finally, with 503 upregulated genes, no enriched pathways were annotated for the new NGS-specific genes.