Antibody Mediated Rejection and T-cell Mediated Rejection Molecular Signatures Using Next-Generation Sequencing in Kidney Transplant Biopsies

Recently, interest in transcriptomic assessment of kidney biopsies has been growing. This study investigates the use of NGS to identify gene expression changes and analyse the pathways involved in rejection. An Illumina bulk RNA sequencing on the polyadenylated RNA of 770 kidney biopsies was conducted. Differentially-expressed genes (DEGs) were determined for AMR and TCMR using DESeq2. Genes were segregated according to their previous descriptions in known panels (microarray or the Banff Human Organ Transplant (B-HOT) panel) to obtain NGS-specific genes. Pathway enrichment analysis was performed using the Reactome and Kyoto Encyclopaedia of Genes and Genomes (KEGG) public repositories. The differential gene expression using NGS analysis identified 6,141 and 8,478 transcripts associated with AMR and TCMR. While most of the genes identified were included in the microarray and the B-HOT panels, NGS analysis identified 603 (9.8%) and 1,186 (14%) new specific genes. Pathways analysis showed that the B-HOT panel was associated with the main immunological processes involved during AMR and TCMR. The microarrays specifically integrated metabolic functions and cell cycle progression processes. Novel NGS-specific based transcripts associated with AMR and TCMR were discovered, which might represent a novel source of targets for drug designing and repurposing.

Recently, interest in transcriptomic assessment of kidney biopsies has been growing.This study investigates the use of NGS to identify gene expression changes and analyse the pathways involved in rejection.An Illumina bulk RNA sequencing on the polyadenylated RNA of 770 kidney biopsies was conducted.Differentially-expressed genes (DEGs) were determined for AMR and TCMR using DESeq2.Genes were segregated according to their previous descriptions in known panels (microarray or the Banff Human Organ Transplant

INTRODUCTION
Long-term kidney allograft survival is mainly limited by the occurrence of rejections [1,2].To improve kidney injury detection, the biennial revision of the Banff classification emerged as the gold standard for the diagnosis of rejection during the past 3 decades [3,4].From histology assessment of kidney biopsies, combined with clinical and immunological parameters, the classification is now encompassing molecular and digital biomarkers to improve its sensitivity and provide new diagnostic tools for the clinicians.Recently, transcriptome analysis has shown its capacity to accurately detect injuries and the degree of activity in solid organ transplant biopsies [5].Previous studies focusing on the implementation of microarrays paved the way for the molecular understanding of rejection and allowed the development of gene expression-based classifiers [6].However, this technology suffers from its necessity to design probes, limiting the past studies to the coding transcriptome only.
While lacking protein-coding ability, long noncoding RNAs (lncRNAs) act as functional RNA molecules, regulating protein-coding gene expression through interactions with generegulatory proteins and microRNAs.Growing evidences in the literature showed the pivotal role played by lncRNAs in the establishment and maintenance of the immune response [7][8][9].Therefore, they represent a complete novel source of biomarkers for the diagnosis of various cancers [10][11][12].However, lncRNAs implication in the solid organ transplantation field remains poorly investigated.Combining the non-coding transcriptome on top of the coding might help our understanding of the pathophysiological mechanisms involved during kidney allograft rejection, could improve the molecular classifiers for its detection and prediction and provide new and unknown targets for drug designing and repurposing.
In the present study we investigated the discovery capability of Next-Generation Sequencing (NGS) to unravel both coding and non-coding transcriptome differentially expressed during rejection.For that purpose, we built a real-world, multicentric and extensively phenotyped cohort of 540 patients (770 biopsies) from two clinical studies: EU-TRAIN (NCT03652402) and KTD-Innov (NCT03582436).We performed an Illumina sequencing, analyzed the samples with differential gene expression analysis, identified known genes according to published gene panels (microarray or the Banff Human Organ Transplant) to identify new transcripts and implemented pathway enrichment analysis on the different subgroups.

Study Population and Biopsy Cohort
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.

Kidney Allograft Phenotypes
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.

Detection and Characterization of Circulating Donor-specific anti-HLA Antibodies
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.

Experimental Procedures
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.with the STAR algorithm (v2.7.4a) [17] and the Hg38.p13reference 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.p13GTF annotation file and the BAM files resulting from the alignment.Quality controls results can be found in the Supplementary Table S1.

Differential Gene Expression Analysis
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-cellmediated 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 2 fold change was applied and all significant (adjusted p-value <0.05) differentiallyexpressed genes were considered during the analysis.
The complete description of differentially expressed gene symbols, mean expressions, log 2 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 Enrichment Analysis
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.

Statistical Analysis
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.

Characteristics of the Population
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, 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.73m 2 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).

Molecular Landscape for T-cell Mediated Rejection
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).

DISCUSSION
In this study we aimed at defining and describing the molecular profiles and biological functions associated with antibodymediated and T-cell mediated rejection, combining a deeply phenotyped cohort of kidney allograft biopsies and nextgeneration sequencing analyses.For this purpose, we used the histological labels and the gene expressions as inputs for a differential expression analysis and ranked the significant genes according to the adjusted p-values.We, then, queried publicly available biological databases to understand the pathophysiological mechanisms derived from the upregulated DEGs.In this study, an emphasis was made to discriminate genes from known gene panels (B-HOT and microarray) already validated and used in clinical practice [23][24][25][26] and new genes discovered using the NGS technology.
In the present study, active antibody-mediated rejection was found in 9.4% of the analyzed samples.This incidence aligns well with the most recently reported incidence of AMR ranging from 3% to 12% in a recent systematic review including 28 studies [27].Its molecular signature included features of macrophages activation (CD40, CD58, IDO1), NK cells activation (GNLY, FGFBP2, CD16a), cytotoxic T cells activation (CD8), helper T cells activation (CD4), endothelial cells activation (ICAM1, PECAM1, VCAM1, CDH5), and B cells activation (CD22, CD40, CD86), which showed great consistency with the microarray studies [28,29].From both innate and adaptive immune systems, the enrichment analysis confirmed the ability of the B-HOT gene panel to capture both components occurring during rejection but presented a lack of metabolic functions, such as SUMOylation processes and cell cycle progression and checkpoint that are specifically present in the microarrays.Regarding the NGS-specific gene panel, 603 new genes (comprising 264 upregulation) were found associated with AMR but no annotation was available in the public repositories.They were mainly composed of long non-coding RNAs that are poorly described in the literature.PELATON, for instance, was part of the new NGS-specific top ranked genes and was found to be a regulator specifically located in macrophages and monocytes nucleus, for which the downregulation is associated to decreased phagocytosis functions [30].In this study, we found that PELATON was upregulated during AMR (log 2 FC = 1.56), in line with a probable increased phagocytosis function occurring in the microcirculation inflammation and, consequently, potentially leading to increased differentiation into antigen-presenting cells, T-cell recruitment and activation and, ultimately, B-cell proliferation and transformation into plasma cells.
Compared to the AMR signature, the TCMR signal presented a similar profile compared to the published studies in terms of genes (CD72, LAG3, CD8A, CD28, ANKRD family) and activated cell types and functions [31].However, a key difference existed in the repertoire of inhibited cell adhesion molecules, showing strong inhibition of the endothelial and epithelial cells receptors, emphasizing the cell infiltration observed at the histological level.Regarding the different gene panels, the microarray was specific of mRNA maturation processes and nonsense mediated decay, which could be due to a lack of annotation of the different repositories.In our study, the B-HOT panel was enriched by the main immunological functions but did not include the top adjusted p-value ranked genes, potentially limiting its ability to accurately diagnose TCMR.The addition of new genes, for example, from the microarray or discovered with the NGS technology, could potentially help the molecular classifiers.Finally, for the NGSspecific markers, they were composed of lncRNA which lacked annotation in the current repositories.Few of them are described in the literature such as MIR3142HG which was shown to be a positive regulator of IL-8 and CCL2 [32].
The main advantage of the present study is that the cohort's diverse phenotypes encompass most of the clinical scenarios encountered in routine practice.It also gathered samples and patients representing a real-life setting in terms of population demographics, rejection prevalence and immunosuppression therapies.Lastly, this is, to our knowledge, the first RNA-seq experiment applied in such cohort characteristics (size, heterogeneity, description) to study the molecular signature of rejection in kidney allograft biopsies.A literature review on PubMed comprising the key words "NGS," "transplantation," "kidney" and "rejection" resulted in 46 articles published over the last 5 years: 11 (23.4%) were related to cell-free DNA, 9 (19.6%) were related to infections (comprising also BK virus), 4 (8.7%) were focusing on cell subpopulations, 5 (10.9%) were related to response to treatment and 5 (10.9%) were related to HLA matching.Five references mentioned either the use of NGS or the B-HOT gene panels but showed limitations in the number of patients/biopsies, number of genes under study, in their design (sick vs. well, single centre), or in the representativity of the different diagnoses [33][34][35][36][37].
Regarding the study limitations, one of the main issues is the sampling bias regarding the technique requiring an extra core.The sequenced core might be different from the one analyzed by the pathologist both in terms of quality (i.e., number of glomeruli and arteries) and severity of the disease, which is not the case for the Nanostring technology and the B-HOT gene panel where an extra core is not needed.Second, while NGS might help to discover new genes and physio-pathological pathways, its use in clinical practice is limited in terms of access to the technology and its cost.In our study, most of the genes associated with AMR and TCMR were included in the microarray and the B-HOT-gene panels, validating the relevance and the accuracy of the genes included.Finally, from a clinical aspect, our cohort was mainly treated with corticosteroids, mycophenolic acid and tacrolimus, which might have an impact on the observed molecular expressions.The presented results should be validated on patients treated with different types of immunosuppressive therapies including mTOR inhibitors or Belatacept.

CONCLUSION
We discovered 9.8% and 14.0% novel transcripts associated with antibody-mediated rejection and T-cell mediated rejection, respectively.The main immunological functions were positively captured by both the microarray and B-HOT gene panels.Those new NGS specific transcripts might represent a novel source of targets for drug designing and repurposing.

FIGURE 1 |
FIGURE 1 | Antibody-mediated rejection molecular signature.Volcano plot of the significant differentially expressed genes associated with antibody-mediated rejection.Dots are related to each gene.The significant transcripts are displayed according to a 0.05 threshold on the adjusted p-value (vertical grey line).NGS-specific transcripts are highlighted in red, B-HOT-related in yellow and microarray-related in blue.X-axis represents the -log 10 of the adjusted p-value (the higher, the smaller is the p-value) and the y-axis represents the log 2 fold change.Differences in gene expression between the AMR and non-AMR group are marked with positive (negative) values correspond to up-(down-)regulated transcripts in the AMR group.In total, 6,141 genes were differentially expressed showing the following distribution: 358 included in the B-HOT gene panel, 5,180 included in the microarray gene panel and 603 NGS-specific.Abbreviations: AMR: antibody-mediated rejection; B-HOT: Banff Human Organ Transplant; NGS: next-generation sequencing.

FIGURE 3 |
FIGURE 3 | T-cell mediated rejection molecular signature.Volcano plot of the significant differentially expressed genes associated with T-cell mediated rejection.Dots are related to each gene.The significant transcripts are displayed according to a 0.05 threshold on the adjusted p-value (vertical grey line).NGS-specific transcripts are highlighted in red, B-HOT-related in yellow and microarray-related in blue.X-axis represents the -log 10 of the adjusted p-value (the higher, the smaller is the p-value) and the y-axis represents the log 2 fold change.Differences in gene expression between the TCMR and non-TCMR group are marked with positive (negative) values correspond to up-(down-)regulated transcripts in the TCMR group.In total, 8,478 genes were differentially expressed showing the following distribution: 439 included in the B-HOT, 6,853 included in the microarray and 1,186 NGS-specific.Abbreviations: B-HOT: Banff Human Organ Transplant; NGS: next-generation sequencing; TCMR: T-cell mediated rejection.

TABLE 1 |
Baseline patient characteristics.Delayed graft function was defined as the use of dialysis in the first postoperative week.

TABLE 2 |
Histological, immunological and functional characteristics at the time of biopsy.