Peripheral Blood RNA Sequencing Unravels a Differential Signature of Coding and Noncoding Genes by Types of Kidney Allograft Rejection.

INTRODUCTION: Peripheral blood (PB) molecular patterns characterizing the different effector immune pathways driving distinct kidney rejection types remain to be fully elucidated. We hypothesized that transcriptome analysis using RNA sequencing (RNAseq) in samples of kidney transplant patients would enable the identification of unique protein-coding and noncoding genes that may be able to segregate different rejection phenotypes.
METHODS: We evaluated 37 biopsy-paired PB samples from the discovery cohort, with stable (STA), antibody-mediated rejection (AMR), and T cell-mediated rejection (TCMR) by RNAseq. Advanced machine learning tools were used to perform 3-way differential gene expression analysis to identify gene signatures associated with rejection. We then performed functional in silico analysis and validation by Fluidigm (San Francisco, CA) in 62 samples from 2 independent kidney transplant cohorts.
RESULTS: We found 102 genes (63 coding genes and 39 noncoding genes) associated with AMR (54 upregulated), TCMR (23 upregulated), and STA (25 upregulated) perfectly clustered with each rejection phenotype and highly correlated with main histologic lesions (ρ = 0.91). For the genes associated with AMR, we found enrichment in regulation of endoplasmic reticulum stress, adaptive immunity, and Ig class-switching. In the validation, we found that the SIGLEC17P pseudogene and 9 SIGLEC17P-related coding genes were highly expressed among AMR but not in TCMR and STA samples.
CONCLUSIONS: This analysis identifies a critical gene signature in PB in kidney transplant patients undergoing AMR, sufficient to differentiate them from patients with TCMR and immunologically quiescent kidney allografts. Our findings provide the basis for new studies dissecting the role of noncoding genes in the pathophysiology of kidney allograft rejection and their potential value as noninvasive biomarkers of the rejection process.

Kidney transplantation remains the preferred treatment for end-stage renal disease., However, chronic immune-mediated allograft rejection, fundamentally driven by the humoral effector pathway of adaptive immunity, remains the main cause of accelerated graft loss., Recent clinical and experimental evidence suggest that chronic AMR progressively appears as a continuum process, preceded by T-cell immune activation, which can lead to clinical or subclinical TCMR.5, 6, 7, 8 While the downstream effector mechanisms of rejection seem to be compartmentalized and driven by specific effector immune pathways, such as alloreactive T cells and donor-specific antibodies (DSAs), a broader activation of the alloimmune response is more likely to occur during the rejection process.9, 10, 11

Currently, while different rejection phenotypes have been classified according to defined histologic changes in kidney transplant biopsy specimens, such pathologic lesions are not highly specific as no histology-independent biopsy diagnostic system has been fully validated yet.4, 5, 6 Therefore, the study of the underlying molecular mechanisms of different histologic injury patterns using genome-wide gene expression on kidney allograft biopsy specimens has gained considerable interest to more precisely characterize the different pathologic processes of transplant rejection.,12, 13, 14, 15, 16 Notably, the evaluation of the different molecular mechanisms of allograft rejection has been also investigated in PB, in an attempt to better characterize the distinct effector immune pathways leading to allograft rejection while avoiding an invasive tissue allograft assessment.

The majority of the blood-based differential gene expression studies have found common PB signatures of rejection, but transcriptional profiling with cDNA and oligo-based microarrays,, as well as other unbiased approaches, such as proteomics,,, protein arrays, and transcriptional profiling,,,,, have not been able to discern specific molecular variations distinguishing between different rejection phenotypes, because many rejections show combined pathologic features of both TCMR and AMR, the so-called mixed rejections. For instance, AMR transcripts have been identified from biopsy specimens to be strongly regulated by interferon-γ, but a similar axis of interferon-γ regulation has also been identified in TCMR and mixed rejections in PB. Therefore, while previous published transcriptional studies have allowed for a better understanding of the overall biology of kidney transplant rejection, distinct differences in molecular pathways regulating TCMR and AMR have been difficult to dissect in any biologically meaningful manner.

We hypothesized that some of the difficulty in identifying specific PB features of AMR, distinct from TCMR, may relate to the specific differences, possibly resulting in changes in microRNAs or long noncoding RNAs, which are not usually measured by standard microarray technologies. With this in mind, we conducted RNAseq on unique biopsy-matched PB samples with the aim to conduct transcriptome analysis at much greater resolution because RNAseq allows for the differentiation between protein-coding genes and noncoding genes, including pseudogenes, microRNA, and long noncoding RNA.23, 24, 25 The aim of this approach was to gain a more granular view of the possible mechanisms of the main gene pathways of the adaptive and innate immune responses in PB during an active rejection episode in the graft, specific to either TCMR or AMR, compared with immunologically quiescent stable kidney transplant patients.

Methods

Patients

One hundred kidney transplant patients were evaluated in this study. A first cohort of 37 consecutive adult kidney transplant recipients undergoing either for-cause or surveillance biopsy procedures between January 2014 and December 2015 at Bellvitge University Hospital (Barcelona, Spain) were enrolled in the study. Three different groups of patients were categorized according to the histologic diagnosis made by a blinded expert transplant pathologist following the updated Banff 2017 classification: patients with pure acute TCMR (n = 13), patients with AMR (n = 12), and patients with a normal/preserved kidney allograft parenchyma (STA, n = 12). All blood samples were obtained at the time of the kidney allograft biopsy procedure and before any immunosuppression rescue therapy was given. Subsequently, 63 additional kidney transplant patients from 2 distinct transplant centers undergoing consecutive for-cause or surveillance biopsy procedures between January 2017 and October 2018 were evaluated for targeted mRNA expression of most representative noncoding and coding-related genes observed in the first discovery cohort of patients. From this second cohort, 33 patients were from Bellvitge University Hospital, Barcelona, Spain (validation cohort I) and 29 patients were from University California San Francisco, San Francisco, USA (validation cohort II). In the first validation cohort (I), the main histologic diagnoses in for-cause biopsy procedures were TCMR (n = 10) and AMR (n = 9), whereas 8 patients were evaluated in surveillance biopsy procedures and were classified as STA (n = 15). In the second validation cohort (II), the main histologic diagnoses were STA (n = 12), AMR (n = 7), and mixed rejections (n = 10). Each biopsy specimen was scored by a site’s specialist pathologist blinded for any result of the transcriptional study. The main clinical, immunologic, and histologic characteristics are shown in Supplementary Table S1. All patients in the study were receiving the same triple immunosuppressive-based therapy with tacrolimus, mycophenolate mofetil, and prednisone. All patients evaluated in the study gave written informed consent to participate and the institutional review boards at Bellvitge University Hospital (PR228/13) and the University of California–San Francisco approved the study (14-13573).

Renal Allograft Histology

Kidney allograft biopsy procedures were performed for cause because of kidney allograft dysfunction in all patients with TCMR and AMR rejection, whereas STA patients showed stable allograft function and a preserved graft parenchyma in 6-month surveillance biopsy procedures. All renal biopsy specimens were analyzed following the Banff 2017 score classification and the histologic analysis was blindly evaluated by an expert renal pathologist at each center before submitting samples for molecular evaluation. Details are shown in Supplementary Methods.

Statistical Analysis

We performed differential expression (DE) analysis to find those genes associated with the outcome. We applied DESeq2 with a binomial distribution in all pairwise combinations (AMR vs. STA, TCMR vs. STA, and AMR vs. TCMR) and Elastic Net (ENET) with a multinomial distribution for the 3-way comparison (AMR vs. TCMR vs. STA). With the results from ENET we performed a network analysis to show the association between each gene and each of the clinical variables using a linear regression model and we represented the final network using the Fruchterman-Reingold layout, and finally we performed a linear regression model to find those protein-coding genes associated with each noncoding gene and gene ontology (GO) annotations to predict the potential function of the noncoding RNA genes.

For comparison with previously published microarray analysis, we also applied ENET to find genes differential expressed between rejectors (AMR plus TCMR) and STA. For this validation, we compared our results with 2 lists of previously published differentially expressed genes between STA and rejections (Supplementary Table S2). One is the list of 2382 differentially expressed genes (false discovery rate [FDR] < 0.05) identified in 4 different datasets using 3 different microarray platforms measured in either whole blood or PB mononuclear cells and the other is the list of 2977 identified probesets belonging to 1907 genes in PB in adult kidney transplant patients associated with kidney rejection.

For the validation cohorts, we used 1-way analysis of variance to determine significant differences between groups and we used the Tukey multiple comparison test to compare the difference between each pair of means. Further details of the statistical analyses are found in the Supplementary Material and Methods.

Results

Clinical and Histologic Characteristics of the Discovery Cohort

To understand the fundamental molecular mechanisms underlying each histologic phenotype, we used RNAseq to profile unique blood samples from 37 consecutive kidney transplant patients, classified in 3 specific clinical phenotypes: STA (n = 12), AMR (n = 12), and TCMR (n = 13). The baseline clinical characteristics of this patient cohort is shown in Table 1. Patients with TCMR and AMR showed worse allograft function than STA patients who showed a good and stable kidney allograft function. Only patients with AMR displayed DSA. TCMR Banff scores ranged between IA and IIA (6 IA, 6 IB, and 3 IIA).

Table 1

Main clinical and demographic characteristics

Main clinical variables	Clinical Phenotypes
	TCMR, n = 13	AMR, n = 12	STA, n = 12	P value
Donor age, yr ± SD	61.5 ± 14.09	43.2 ± 21.92	50.6 ± 17.24	0.048
Recipient age, yr ± SD	62.2 ± 12.25	45.7 ± 15.26	55.3 ± 11.98	0.013
Recipient sex, female, n (%)	5 (38.5)	4 (33.3)	5 (41.7)	0.91
Donor sex, female, n (%)	8 (61.5)	5 (41.7)	5 (41.7)	0.51
Cause of ESRD, n (%)				0.55
Unknown	6 (46)	3 (25)	5 (36)
Glomerular	4 (31)	5 (42)	2 (16.7)
Interstitial	0 (0)	2 (17)	1 (8.3)
Vascular	1 (7.7)	1 (8.3)	2 (16.7)
Diabetes	2 (15.4)	0 (0)	1 (8.3)
APKD	0 (0)	0 (0)	1 (8.3)
Others	0 (0)	1 (8.3)	0 (0)
Type of transplant, died, n (%)	10 (77)	11 (91.7)	11 (91.7)	0.46
No. of transplants ± SD	1.15 ± 0.38	1.6 ± 0.79	1.17 ± 0.39	0.1
1 vs. >1, n (%)	11 (85)	7 (58)	10 (83)	0.23
No. HLA antigen mismatch ± SD	3.0 ± 0.95	3.58 ± 0.9	3.25 ± 1.3	0.41
Induction type, n (%)				0.064
None	0 (0)	2 (16.7)	1 (8.3)
Anti-CD25 mAb	10 (76.9)	4 (33.3)	10 (83.3)
rATG	3 (23.1)	6 (50.0)	1 (8.3)
DSA at biopsy (yes), n (%)	0	12 (100)	0	<0.001
Class I		3
ClassII		7
Class I and II		2
eGFR at biopsy, ml/min ± SD	30.1 ± 20.42	28.8 ± 20.9	49.7 ± 14.42	0.016
Proteinuria at biopsy, g/24 h ± SD	0.76 ± 0.9	1.79 ± 1.37	0.19 ± 0.23	0.001
Graft loss after biopsy assessment, n (%)	4 (30.8)	8 (66.7)	1 (8.3)	0.01
Time to biopsy, months ± SD	4.8 ± 3.8	91.6 ± 83.1	6.8 ± 2.4	<0.001

AMR, antibody-mediated rejection; APKD, autosomic polycystic disease; DSA, donor-specific antibody; eGFR, estimated glomerular filtration rate; ESRD, end-stage renal disease; HLA, human leukocyte antigen; mAb, monoclonal antibody; rATG, rabbit antithymocyte globulin; STA, stable; TCMR, T cell–mediated rejection.

In Table 2, all semiquantitative histology scores defining the distinct clinicopathologic groups of the study are described. Each of these histologic measures and time posttransplant were assessed as correlates in the context of RNAseq data.

Table 2

Main histologic lesions of the patients of the study

Mean Banff scores in all kidney graft compartments	Histologic phenotypes
	TCMR, n = 13	AMR, n = 12	STA, n = 12	P value
Acute lesions
Ag	0.8 ± 0.8	1.5 ± 0.9	0.2 ± 0.4	<0.001
Ai	2.1 ± 0.8	1.1 ± 0.7	0.2 ± 0.6	<0.001
at	2.3 ± 0.6	0.8 ± 0.5	0.5 ± 0.9	<0.001
ti	2.1 ± 0.9	1.4 ± 0.8	0.2 ± 0.4	<0.001
ptc	0.5 ± 1.0	1.0 ± 0.8	0.09 ± 0.3	0.030
av	0.2 ± 0.8	0.2 ± 0.4	0.08 ± 0.3	0.81
C4d	0.1 ± 0.4	1.7 ± 1.0	0.1 ± 0.6	<0.001
Chronic lesions
cg	0.08 ± 0.3	1.7 ± 1.3	0.08 ± 0.3	<0.001
ci	0.8 ± 0.4	1.5 ± 0.8	0.8 ± 0.8	0.03
ct	0.7 ± 0.5	1.7 ± 0.7	0.8 ± 0.8	0.001
cv	0.3 ± 0.5	0.7 ± 0.9	0.2 ± 0.6	0.15
ah	0.1 ± 0.5	1.1 ± 1.2	0.2 ± 0.6	0.01
cm	0 ± 0	0.9 ± 1.1	0.08 ± 0.3	0.002

ag, acute glomeruli; ah, arterial hyalinosis; ai, acute interstitium; AMR, antibody-mediated rejection; at, acute tubuli; av, acute vascular; cg, chronic glomeruli; ci, chronic interstitium; Cm, chronic mesangial; ct, chronic tubuli; cv, chronic vascular; ptc, peritubular capillaritis; STA, stable; TCMR, T cell–mediated rejection ti, total interstitial inflammation.

TCMR: at, ai, av, and cv.

AMR: ag, ptc, c4d, cg, and cv plus donor-specific antibody/antiHLAb.

Interstitial fibrosis tubular atrophy: ci and ct.

Differentially Expressed Genes Classify Patients into the 3 Clinical Phenotypes

We performed DE analysis between the 2 types of rejection phenotypes and STA patients using RNAseq first on only coding and after adding the noncoding genes. We initially considered both types of rejection as a single outcome (REJ) and then considered AMR and TCMR separately, addressing all possible comparisons. Subsequently, a 3-way comparison analysis was performed. We used both DESeq2 and ENET for these analyses. As shown in Table 3, based on the Jaccard index (JI), ENET, which is a penalized technique that considers the correlation structure between the variables, classified all the samples more accurately than DESeq2, suggesting the importance of considering the correlation structure between the genes. When noncoding genes were added, the results significantly improved not only increasing the JI but also lowering the number of genes selected to classify the outcome, highlighting the importance of noncoding genes in the rejection process. When all the patients who rejected were grouped together (REJ) within 1 category, neither DESeq2 nor ENET were able to properly classify the data, whereas when separated according to the rejection phenotype (AMR or TCMR) the differences were significantly stronger, especially when using the 3-way comparison, where ENET found 102 genes (63 coding and 39 noncoding genes) associated with AMR (54 upregulated), TCMR (23 upregulated), and STA (25 upregulated) showing 3 main clusters with a JI of 1. Figure 1 shows a heat map and a cluster plot of the 3-way comparison analysis. The heat maps showing the results of all comparisons made are shown in Supplementary Figures S1–S10 and Supplementary Table S3. To test for possible overfitting with ENET using a multinomial distribution, we performed a permutation test shuffling the clinical outcome 10 times and applying exactly the same strategy as with the original data. None of the permutations had accurate results with a mean JI of 0.27 (Supplementary Figure S11). Finally, to test whether time to biopsy procedure was a confounding factor for the gene expression signature reported here, we applied the same method (ENET with a multinomial distribution) to the same model, adding this variable. The results remain the same with the same 102 genes selected (Supplementary Table S4).

Table 3

Differentially expressed genes for all the possible comparisons (pairwise and 3-way) using DESeq2 and ENET and considering only coding genes and the addition of noncoding genes

Method	Gene list	Differentially expressed genes (Jaccard indexa)
		REJ vs. STA	AMR vs. STA	TCMR vs. STA	AMR vs. TCMR	AMR vs. TCMR vs. STA
ENET (optimal alpha, lambda by CV)	Coding plus noncoding	36 (0.47)	59 (1)	1 (–)	23 (1)	102 (1)
	Coding	328 (0.43)	50 (1)	1 (–)	1073 (0.64)	131 (0.68)
DESeq2 (FDR < 0.05, |FC| > 1.5)	Coding plus noncoding	1176 (0.40)	4774 (0.55)	0 (–)	2099 (0.50)	—
	Coding	875 (0.40)	3541 (0.55)	0 (–)	1739 (0.50)	—
DESeq2 (FDR < 0.05)	Coding plus noncoding	1391 (0.40)	5482 (0.55)	0 (–)	2244 (0.50)	—
	Coding	1094 (0.40)	4221 (0.55)	0 (–)	2092 (0.50)	—

AMR, antibody-mediated rejection; FC, fold change; FDR, false discovery rate; REJ, rejection (AMR plus TCMR); STA, stable; TCMR, T cell–mediated rejection.

Jaccard index measures the similarity of the samples within the same category. Closer to 1 means more similar, closer to 0 means more different.

Figure 1

(a) Heatmap showing the 102 genes selected by Elastic Net (ENET) using a multinomial distribution. The selection is based on the optimal alpha and lambda parameter using cross-validation. (b) Cluster plot showing the similarity between the samples using the 102 genes selected by ENET. Cluster 1 classifies the T cell–mediated rejection patients (blue), cluster 2 the antibody-mediated rejection patients (red), and cluster 3 the stable patients (green). The color scale shown in the heatmap matrix represents the differential expression normalized by column to show the differences by samples per gene.

Comparison of RNAseq Data and Microarray Data in Kidney Allograft Rejection

We then compared the 1094 coding genes detected by DESeq2 package (Table 3) with those identified in previous published microarray studies using the same FDR (<5%) between rejection and stable patients. Of the 2382 differentially expressed genes identified in 3 different microarray platforms from 4 different data sets evaluated in either peripheral blood mononuclear cells or whole blood, 225 genes were also present in our analysis. Of the 1907 identified genes in PB in adult kidney transplant patients developing kidney allograft rejection, 183 genes were also observed in our work. The list of 1094 genes with the validated ones are shown in Supplementary Table S5. However, even though we showed that by grouping all the patients who rejected into 1 single group does not result in clean differentiation between AMR and TCMR by microarray-based transcriptional profiling, the overlapping rejection-specific genes detected in PB, across multiple datasets, using different platforms and patient samples, supports the quality and biological validity of the current dataset and analytics.

Network Analysis Finds Association Between the Differentially Expressed Genes and the Fundamental Histologic Lesions of Allograft Rejection

The 102 genes found to be associated with AMR, TCMR, and STA were also associated with the fundamental histologic lesions defining each clinical phenotype, time to biopsy procedure, and DSA. Figure 2 shows the correlation matrix among all these parameters. As expected, histologic chronic lesions (chronic vascular, chronic tubuli, chronic interstitium, chronic glomeruli, arterial hyalinosis, and chronic mesangial) correlated between them and with DSA and time to biopsy procedure. Regarding acute histologic lesions, 2 main correlation networks were observed, 1 clustering the acute tubuli, acute interstitium, and total interstitial inflammation lesions and another at the acute vascular, peritubular capillaritis, acute glomeruli, and C4d lesions. The latter net did also correlate with DSA and chronic glomeruli, which altogether defined the clinical phenotype AMR. Also, chronic interstitium and chronic tubuli lesions, defining IFTA, showed the highest correlation (ρ = 0.91).

Figure 2

Correlation matrix for all fundamental histologic lesions, donor-specific antibodies, and time to biopsy procedure. The coefficients belong to a Pearson correlation and are colored if the P value is significant (< 0.05).

For the association analysis, we performed a linear model between each of the 102 genes and all fundamental histologic parameters, time to biopsy assessment, and DSA. Those that passed a FDR threshold of 0.05 were plotted in a network using the Fruchterman-Reingold layout (Figure 3). We observed 3 networks in the correlation matrix: 1 main gene cluster associated with main chronic histologic lesions (chronic glomeruli, chronic interstitium, chronic tubuli, chronic vascular, arterial hyalinosis, and chronic mesangial), time to biopsy assessment, and DSA. Interestingly, we observed that the genes that were upregulated in AMR (AP4S1, ARMCX4, BTD, CACHD1, CDK11B, DMAP1, ELF3, LA16c-360A4.1, LINC00298, LINC01278, LRIG1, MFSD6L, MPP7, NECAB3, P2RX5-TAX1BP3, PARD6G-AS1, RAB30, RN7SL246P, RP11-1008M1.1, RP11-133M8.3, RP5-1107A17.4, RPL4P6, SEPT2, SETDB2, SUPT5H, SYVN1, TADA3, TP53BP1, UBAC2, UC_338, USP21, and ZMYF841) were positively associated with these clinical parameters, whereas those upregulated in TCMR and STA were negatively associated with them. Furthermore, the acute histologic lesions acute glomeruli, peritubular capillaritis, and C4d were also close to this main cluster as also shown in the correlation matrix, but they also shared genes with the other cluster made by acute tubuli, acute interstitium, and total interstitial inflammation lesions. Notably, the genes associated with this cluster were positively associated with those upregulated in TCMR, such as AK9, RAD23B, RPP21, RP11-43N22.1, UCN, S100A5, RP11-392P7.6, MERTK, and ACTR3B, and negatively associated with those upregulated in STA, and few of them were associated with AMR. The parameter acute vascular was not associated with any gene. The set of genes associated with the different clinical and histologic parameters is shown in Supplementary Table S6.

Figure 3

Network plot showing the association between the 102 genes and the fundamental histologic lesions. The vertex represents the genes belonging to each cluster shown in Figure 1 and the histologic lesions and the edges link those that were associated with a false discover rate < 0.05 in a linear regression model. The width of each line represents the statistical significances (log10 P value) and the color whether they are positively associated (gray) or negatively associated (orange).

Functional Annotation of Differentially Expressed Genes Using GO Biological Terms

We annotated the 102 genes DE by the 3 different phenotypes (AMR, TCMR, and STA) with the GO biological terms using the platform FUMA., We included separately the 54 genes upregulated in AMR and the 23 upregulated in TCMR. We found several enriched GO terms (FDR < 0.05) involving 13 overlapping genes that were upregulated in AMR (SYVN1, UBAC2, DMAP1, SETDB2, TP53BP1, TADA3, USP21, SUPT5H, CDK11B, SEPT2, RAB30, STC17, and DACT1; Table 4). All of the GO terms might be grouped by similarity into 3 main groups: regulation of protein stability, Golgi organization, and protein acetylation. We also found significant enriched GO terms (FDR < 0.05) using the upregulated genes in TCMR with 13 overlapping genes (CMTM3, HCST, MERTK, LGALS1, PSBM1, ANGPT1, UCN, AK9, DCTPP1, ARPP19, SNRPF, RAD23B, and ZNF354A; Table 5), which are grouped by similarity in immune system development, regulation of transference activity, protein phosphorylated, cellular macromolecular complex assembly, and response to organic cyclic compound. ANGPT1, which belongs to the angiopoietin family, is shared among the majority of all GO terms in TCMR.

Table 4

GO biological terms for the functionality of the upregulated genes in antibody-mediated rejection

GO biological term	FDR	Overlapping genes
		SYVN1	UBAC2	DMAP1	SETDB2	TP53BP1	TADA3	USP21	SUPT5H	CDK11B	SEPT2	RAB30	STX17	DACT1
Negative regulation of response to endoplasmic reticulum stress (GO 1903573)	0.04	X	x
Peptidyl lysine modification (GO 0018205)	0.04			X	x	x	x
Chromatin modification (GO 0006325)	0.04			X	x		x	x	x
Covalent chromatin modification (GO 0016569)	0.04			X	x		x	x
Mitotic nuclear division (GO 0140014)	0.04				x		x			x	x
Chromatin organization (GO 0006325)	0.04			X	x		x	x	x
Golgi organization (GO 0007030)	0.04											x	x
Regulation of protein stability (GO 0031647)	0.04	X					x							x
Positive regulation of protein import (GO 1904591)	0.04			x										x
Organelle fission (GO 0048285)	0.04				x		x			x	x
Regulation of mRNA metabolic process (GO 1903311)	0.04								x	x
Positive regulation of nucleocytoplasmic transport (GO 0046824)	0.05			x										x
Protein acetylation (GO 0006473)	0.05			x			x

FDR, false discovery rate; GO, gene ontology.

Table 5

GO biological terms for the functionality of the upregulated genes in T cell–medicated rejection

GO biological term	FDR	Overlapping genes
		CMTM3	HCST	MERTK	LGALS1	PSMB1	ANGPT1	UCN	AK9	DCTPP1	ARPP19	SNRPF	RAD23B	ZNF354A
Regulation of immune system process (GO 0002684)	0.01	x	x	x	x	x	x	X
Positive regulation of cell communication (GO 0010647)	0.01	x	x		x	x	x	X
Phosphorylation (GO 0016310)	0.01		x	x		x	x		x
Positive regulation of response to stimulus (GO 0048584)	0.01	x	x		x	x	x	X
Phosphate containing compound metabolic process (GO 0006796)	0.01		x	x		x	x	X	x	x
Regulation of immune response (GO 0050776)	0.01	x	x			x	x
Positive regulation of transport (GO 0051050)	0.01			x			x	x			x
Protein phosphorylation (GO 0006468)	0.01		x	x		x	x
Immune system development (GO 0002520)	0.01			x	x		x
Positive regulation of transferase activity (GO 0051347)	0.01					x	x	x
Negative regulation of protein modification process (GO 0031400)	0.01					x	x	x
Cellular macromolecular complex assembly (GO 0034622)	0.02						x					x	x
Macromolecular complex assembly (GO 0065003)	0.02						x			x		x	x
Regulation of cell death (GO 0010941)	0.03			x	x	x	x	x
Positive regulation of immune system process (GO 0002684)	0.03	x			x	x
Multi organism reproductive process (GO 0044703)	0.03			x				x					x
Negative regulation of cell death (GO 0060548)	0.03			x			x	x
Positive regulation of intracellular signal transduction (GO:1902533)	0.03		x		x		x
Response to organic cyclic compound (GO 0014070)	0.03				x			x						x
Regulation of transferase activity (GO 0051338)	0.03					x	x	x
Negative regulation of multicellular organismal process (GO 0051241)	0.04				x		x	x
Regulation of protein modification process (GO 0031399)	0.04					x	x	x			x
Catabolic process (GO 0009056)	0.04					x				x		x	x
Regulation of transport (GO 0051049)	0.04			x			x	x			x
Negative regulation of protein metabolic process (GO 0051248)	0.04					x	x	x
Response to external stimulus (GO 0009605)	0.05	x					x	x						x
Locomotion (GO 0040011)	0.05	x		x			x
Positive regulation of protein modification process (GO 0031401)	0.05					x	x	x

FDR, false discovery rate; GO, gene ontology.

Functional Annotation Inferred for Noncoding DE Genes

We next investigated whether noncoding RNA also influenced the different clinical phenotype and its association with coding RNAs based on the idea of co-expression matrix as previously done by others.33, 34, 35 As shown in Figure 4 in a Circos plot, we found 529 coding–noncoding associations (linear regression model, FDR < 0.05) among the 102 DE genes associated with the different clinical phenotypes (Supplementary Table S7). The order of the genes shown in the plot represents the cumulative effect size of their relationships; RP5-1107A17.4, TRDC, and LINC01278 are the noncoding genes with the larger effect size and PSMB1, CACHD1, and SEPT2 are the coding genes with larger effect size. Notably, PSMB1 was also shared among many GO terms enriched in the upregulated genes in TCMR. Since coding and noncoding genes were strongly associated, we subsequently evaluated the functionality of the noncoding genes, by inferring the function of the respective associated coding gene through enriched GO biological terms (FDR < 0.05). As shown in Supplementary Table S8, we observed that most noncoding genes with an inferred GO biological term were upregulated in AMR (RP5-1107A17.4, LINC01278, RPL4P6, LINC00299, PARD6G-AS1, LA16c-360A4.1, AP000936.1, SIGLEC17P, LINC00298, AC159540.14, and RP11-1008M1.1), but 1 noncoding gene (RP11-644F5.11) that was upregulated in TCMR and none in STA. All these genes shared the common term regulation of transcription, DNA templated (GO 0006355) and were associated with the same coding genes (CDK11B, DACT1, DMAP1, ELF3, SETDB2, SUPT5H, TADA3, TP53BP1, USP21, ZMYM6, ZNF138, ZNF354A, ZNF436, and ZNF841) from which 5 of 14 belong to zinc finger transcriptional activators.

Figure 4

Circos plot showing the 529 significant associations (false discovery rate < 0.05) selected from the 102 coding and noncoding genes obtained in the differential expression analysis using ENET. The links represent each association between the coding and the noncoding genes. The coding genes (from left to right) and the noncoding genes (from right to left) are ordered by the sum of effect size.

Characterization of Noncoding Genes and Related Coding Genes in 2 Distinct Cohorts of Kidney Transplant Patients

From our RNAseq in biopsy-matched PB samples of kidney transplant patients, we identified unique protein-coding and noncoding genes. We next aimed at characterizing these targeted mRNA of noncoding genes and related coding genes in 2 additional kidney transplant cohorts displaying a variety of kidney graft histologic phenotypes according to the last Banff classification., We investigated the pseudogene sialic acid–recognizing Ig-like lectins 17P (SIGLEC17P) and 9 SIGLEC17P-associated coding genes, the natural killer (NK)–related transcript (NCAM1/CD56) because of its association with SIGLEC17P and a highly upregulated coding gene associated with B-cell Ig class switch recombination (TP53BP1). We selected the noncoding SIGLEC17P gene because though it is a noncoding gene it is still highly expressed in high levels in NK cells, which is widely accepted to mediate AMR. SIGLECs 13 and 17 were rendered nonfunctional during Hominin evolution. While SIGLEC 17P underwent a single base pair deletion in humans, they are still preserved in chimpanzees and other new world monkeys. The SIGLEC 17P’s mRNA are still widely expressed in NK cells. This gene was significantly enriched in the common GO term and is also found to be upregulated in AMR (downregulated in TCMR) in our RNAseq data. Therefore, such an association is less likely to be coincidental, so we selected SIGLEC17P over other noncoding genes. In cohort II (Figure 5a), expression of SIGLEC17P was also significantly upregulated in AMR patients compared with TCMR and STA. Notably, the 9 coding genes associated with SIGLEC17P (AP4S1, ZMYM6, USP21, DMAP1, SUPT5H, TP53BP1, NECAB3, BTD, and DACT1) were also found to be highly upregulated in AMR compared with TCMR and STA. Likewise to SIGLEC17P, the matured NK cell specific marker NCAM1/CD56 as well as the B cell–related protein coding gene TP53BP1 were also found to be significantly overexpressed in AMR compared with TCMR and STA. Likewise, in cohort III (Figure 5b) the expression of SIGLEC17P was found to be upregulated in AMR than STA but not than mixed rejections. Out of the 9 coding genes highly associated with SIGLEC17P, only BTD and ZMYM6 could be assessed in cohort III. As shown, they were also significantly upregulated in AMR compared with STA but not to mixed rejections.

Figure 5

Bar graphs showing normalized mRNA expression of SIGLEC17P and associated coding genes. We used 1-way analysis of variance to determine significant differences between groups and the Tukey multiple comparison test to compare the difference between each pair of means. (a) Differential gene expression in cohort I. SIGLEC17P antibody-mediated rejection (AMR) versus T cell–mediated rejection (TCMR), P = 0.022; AMRversus stable (STA), P = 0.0007; STA 0.07125 ± 0.03705, AMR 0.2323 ± 0.1133, TCMR: 0.1160 ± 0.09869. AP4S1 AMR versus TCMR, P = 0.031; AMR versus STA: P = 0.037; STA: 0.7652 ± 0.4108, AMR: 1.687 ± 1.399, TCMR: 0.6691 ± 0.5628. ZMYM6 AMR versus TCMR, P = 0.0014; AMR versus STA, P = 0.0038; STA: 0.3861 ± 0.1316, AMR: 0.9380 ± 0.6814, TCMR: 0.2818 ± 0.09229. USP21 AMR versus TCMR, P = 0.0043; AMR versus STA, P = 0.0120; STA: 0.8361 ± 0.3075, AMR: 1.832 ± 1.350, TCMR: 0.6185 ± 0.3997. DMAP1 AMR versus TCMR, P = 0.0049; AMR versus STA, P = 0.0195; STA: 1.251 ± 0.7931, AMR: 2.806 ± 2.091, TCMR: 0.7854 ± 0.2629. SUPT5H AMR versus TCMR, P = 0.0056; AMR versus STA, P = 0.0017; STA: 0.6141 ± 0.3482, AMR: 2.453 ± 1.996, TCMR: 0.7128 ± 0.6880. TP53BP1 AMR versus TCMR, P = 0.0037; AMR versus STA, P = 0.0061; STA: 1.254 ± 0.6196, AMR: 2.743 ± 1.650, TCMR: 1.056 ± 0.7488. NECAB3 AMR versus TCMR, P = 0.0305; AMR versus STA, P = 0.0188; STA: 3.483 ± 1.314, AMR: 8.225 ± 6.061, TCMR: 3.541 ± 2.566. BTD AMR versus TCMR, P = 0.0131; AMR versus STA, P = 0.0108; STA: 4.098 ± 2.208; AMR: 15.35 ± 15.89, TCMR: 3.554 ± 2.415. DACT1 AMR versus TCMR, P = 0.0071; AMR versus STA, P = 0.0046; STA: 10.04 ± 6.806, AMR: 69.81 ± 75.63, TCMR: 9.828 ± 6.197. NCAM1 AMR versus TCMR, P = 0.0154; AMR versus STA, P = 0.0013; STA: 0.7785 ± 0.3715, AMR: 4.438 ± 3.340, TCMR: 1.540 ± 1.667. RAB30 AMR versus TCMR, P = 0.0422; AMR versus STA, P = 0.0342; STA: 1.707 ± 0.916, AMR: 3.68 ± 2.823, TCMR: 1.644 ± 1.351. (b) Differential gene expression in cohort II. SIGLEC17P AMR versus STA, P = 0.0296; STA: 7.618 ± 2.855, AMR: 15.91 ± 8.117. BTD AMR versus STA, P = 0.0334; STA: 0.2316 ± 0.1042, AMR: 0.3905 ± 0.1548. ZMYM6 AMR versus STA, P = 0.0416; STA: 1.537 ± 0.4455, AMR: 2.391 ± 0.8787.

Discussion

While an important body of evidence suggests the preponderant role of the alloimmune response driving kidney graft injury and ultimately graft loss, the main molecular mechanisms participating in the rejection process are still poorly understood. Herein, using RNAseq analysis of PB samples timed with allograft kidney transplant biopsy spciemens we show that different biological processes govern the T- and B-cell effector mechanisms triggering allograft rejection, which are otherwise believed to be both clinically and pathologically diverse processes, thus requiring complete different rescue immunosuppressive therapies. Most importantly, using advanced machine learning techniques, we show the relevance of including the noncoding genes in the analysis to find biological pathways discriminating both types of rejections. The robustness of our findings is further sustained by the characterization of these transcripts in 2 distinct cohorts of kidney transplant patients displaying the same histopathologic lesions compatible with AMR.

Previous gene expression studies,, show a slightly overlap with our results. The fact that previous analyses are only focused on coding genes and acute rejection as a whole make it difficult to find robust differences in pathways and biomarkers that clearly distinguish TCMR and AMR by PB analysis alone, as shown here. Moreover, the technologies used are completely different because our analysis is based on RNA sequencing and the previous studies are based on microarray data, but the advent of next-generation sequencing addresses some of the limitations of microarray analysis, such as the requirement of known genome sequence targets, limited dynamic range, inconsistent hybridization probe properties, and the high level of background noise caused by cross-hybridization. Notably, RNAseq allows for higher resolution and dynamic range and serves not only as an additional tool for clinical diagnosis but also for an in-depth biological assessment of the main molecular mechanisms and pathways involved in allograft rejection after kidney transplantation—and even though we are working with a small dataset and genes that have a small effect might be lost, we found a great number of genes previously associated with rejection and therefore we believe in the quality of our data and analysis. However, we used the Tempus Spin RNA Isolation Kit (Applied Biosystems, Foster City, California, USA) for total RNA extraction from the blood sample, and unfortunately small RNAs are lost using this purification method because of the low molecular weight cutoff from the column (approximately 200 bp). This is why only long noncoding RNAs but no microRNAs could be identified. There are other RNA isolation kits that are specific for isolation of microRNA, but we would have lost all the coding and long noncoding RNAs. We chose the kit that isolates most RNA types rather than just one kind.

Using the 3-way comparison with ENET, we found 102 genes (63 coding and 39 noncoding genes) associated with AMR (54 upregulated), TCMR (23 upregulated), and STA (25 upregulated) that perfectly clustered our samples. Using the expression of these 102 genes, we proposed 3 main clusters defined by the genes that were upregulated within each one, and they highly correlated with the main fundamental histologic lesions used to classify the type of rejection. Indeed, those genes upregulated in AMR positively correlated with transplant glomerulopathy, microvascular inflammation (ag and ptc), C4d deposition, and the presence of DSA, whereas they negatively correlated with genes upregulated in TCMR and STA. Likewise, the histologic lesions of tubulitis, interstitial inflammation defining TCMR, positively correlated with genes upregulated in TCMR.

The biological insights under this set of genes revealed interesting and novel data; a number of upregulated genes in AMR are enriched in the negative regulation of response to endoplasmic reticulum stress GO term involving the SYVN1 and UBAC2 genes, which have been depicted to play a role modulating the adaptive immune response in organ transplantation by orchestrating the synthesis, maturation, and degradation of secreted and transmembrane proteins intervening during ischemia–reperfusion or rejection. Importantly, while RNAseq allows for the discovery of noncoding transcripts, the functional annotation is still poor and difficult to assess. Notably, we found 16 RP11 genes from which 11 were upregulated in rejection samples (6 in AMR and 5 in TCMR). Interestingly, RP11 genes measured in urine have been previously associated with TCMR and suggested as novel biomarkers in kidney rejection. By inferring the functionality of the noncoding genes through the associations with the coding genes, we found that some noncoding genes were significantly enriched in a common GO term (regulation of transcription, DNA-templated) and they all shared 14 overlapping coding genes. Out of the different noncoding genes, we highlight ELF3, TADA3, TP53BP1, UP21, and SIGLEC17P that are upregulated in AMR and connected with the AKT1 and mitogen-activated protein kinase/extracellular signal-regulated kinase pathway.39, 40, 41, 42, 43, 44 Moreover, while TP53BP1 has been described to bind to p53 and play a key role in the repair of double-strand DNA breaks, it has also been shown to play a key role in Ig class-switch recombination during antibody genesis,, a key feature for DSA formation. We have inferred the noncoding genes’ functionality based on the idea of co-expression as has been previously done by other authors,33, 34, 35 although the fact that coding and noncoding gene are associated does not imply a common biological function.

Another relevant finding in our study is the high enrichment in the common GO term of the noncoding SIGLEC17P gene in patients with AMR. SIGLEC17P encoded by SIGLEC genes belongs to the family of transmembrane receptors that bind sialic acid–containing ligands. Despite being a noncoding gene, SIGLEC17P has been shown to be highly expressed in high levels in NK cells. Indeed, it is widely suggested that the pathogenic role of NK cells participating in AMR is caused by the enrichment of NK-related transcripts found in kidney allograft tissue samples of patients undergoing AMR., Therefore, we decided to validate the expression pattern of this pseudogene along with its direct downstream coding genes in PB using quantitative reverse transcription polymerase chain reaction studies n 2 additional cohorts of kidney transplant patients with distinct types of allograft rejections (TCMR, AMR, and mixed rejection as well as STA). Interestingly, a significantly increased expression of SIGLEC17P and 9 different direct downstream coding genes together with the NCAM/CD56 transcript was observed in AMR compared with STA and TCMR patients. Moreover, the selectivity of this gene was further observed when patients with mixed rejection were also included in the validation analysis. As shown, this gene expression did also capture the humoral pattern of these mixed phenotypes whereas it did not in pure TCMR patients, suggesting the potential value of evaluating these genes for identifying the type of allograft rejection and guide decision-making regarding the type of rescue therapy.,

Previous studies using RNAseq in PB samples in the field of transplantation are scarce. Dorr et al. first depicted the molecular changes occurring in PB over time in 32 stable kidney transplant patients. Kurian et al. compared the molecular signatures both in PB and tissue biopsy samples of kidney transplant patients undergoing either clinical or subclinical TCMR using either microarray or RNAseq platforms and showed equivalent predictive performance between both technologies. Unlike these previous studies, we aimed to evaluate the molecular differences between the 2 main types of rejection processes. More recently, an 8 coding-gene signature has been described in PB using whole genome transcriptomic arrays, discriminating kidney transplant patients with AMR from TCMR and STA. We do not replicate this signature in our data as per the design nature of our study, since we used a multinomial model, used RNAseq, and included noncoding genes in the analysis. Nevertheless, in the list of genes extracted from the model (AMR vs. STA) using only coding genes and applying the DESeq2 package, we did find that the CXCL10 and IL15 genes were recently reported.

To recognize the distinct molecular changes between STA, TCMR, and AMR in PB, we carefully selected the patients for this analysis, with paired biopsy procedures at the time of the blood draw subjected to blinded histologic diagnosis performed before the molecular evaluation and before treatment intensification. Thus, though these strategies increase the value of the data, it has necessarily limited the sample numbers available for the study. Nevertheless, to overcome this limitation we used 2 strategies; first, we used an advanced machine learning technique, ENET, that allowed us to study the genes in a multivariate way and to consider the 3 clinical phenotypes (STA, AMR, and TCMR) in the same model using a multinomial distribution, thus overcoming the univariate methods for pairwise comparisons frequently used in DE analysis., Second, we performed a permutation test to validate our results and to make sure we were not overfitting the model. Importantly, the characterization of the same noncoding and coding genes in 2 additional kidney transplant cohort confirms the validity of the data obtained. Also, we cannot completely rule out that the DE genes found among AMR patients compared with TCMR and STA may also have captured the chronicity stage of these histologic lesions because AMR patients had a longer time of follow-up after transplantation than TCMR and STA patients. However, it is well-described that the main immune effector mechanism triggering these lesions does not differ between acute and chronic AMR, regardless the time being diagnosed. Moreover, to overcome this issue, time after transplantation was taken into account in the multinominal ENET model and it did not impact the data obtained.

In summary, our study highlights the capacity of RNAseq to identify robust molecular signals that can accurately distinguish between the 2 main rejection phenotypes in PB and highlights the important role of specific noncoding genes in the pathophysiology of the distinct types of kidney allograft rejection. Additional larger studies are highly warranted to confirm our data and eventually use them as novel biological markers of specific kidney rejection types to ultimately provide guided customized immunosuppressive rescue therapies.

Disclosure

All the authors declared no competing interests.