This study is part of the FP7-funded BIOMArkers of Renal Graft INjuries research program (BIOMARGIN, https://cordis.europa.eu/project/id/305499, ClinicalTrials.gov, number NCT02832661), led by a European research consortium in search of biomarkers of allograft immune lesions in kidney transplant recipients, with diagnostic and/or prognostic value. BIOMARGIN is a multicenter, prospective, multiphase study, performed at four kidney transplant centers in Europe (Necker Hospital, Paris, France; University Hospitals, Leuven, Belgium; Hannover Medical School, Hannover, Germany; and University Hospital, Limoges, France). Samples were prospectively collected at the time of screening and indication biopsies between April 2013 and June 2015 ( Figure 1A ). All transplantations were performed with negative complement-dependent cytotoxicity crossmatches. In these four clinical centers, protocol biopsies were performed at 3, 12, and sometimes at 24 months after transplantation, according to local center practice, in addition to the indication biopsies. The institutional review board at each site approved the study, and all the patients provided written informed consent.

The study was divided into three phases ( Figure 1A ). In the first case-control discovery phase, N=88 biopsy samples were used for global miRNA expression analysis (N=754 miRNAs, excluding potential normalizing miRNAs). We selected samples based on availability and histological criteria (excluding cases with a diagnosis of glomerulonephritis or polyomavirus-associated nephropathy and cases with an unclear diagnosis). Based on local biopsy readings, a first selection was made, which was further refined through central reading by a group of pathologists, independent from the original center. Statistical analysis (see miRNA selection) was used to select the extended list of miRNAs, most differentially expressed between histological phenotypes, that were subsequently quantified in the 2nd phase.

A similar case-control study design was used for the second (selection) phase, where targeted miRNA expression analysis (N=143) was performed on N=99 biopsy samples. The same statistical pipeline was applied to select an even more restricted list of miRNAs, to be subsequently quantified in the 3rd phase.

Finally, in the third (validation) cohort, all samples prospectively and consecutively collected according to the BIOMARGIN protocol between June 24, 2014 and July 2, 2015 and with adequate kidney allograft biopsy quality, were assessed for the 38 miRNAs selected as a result of the second selection phase. In this real-life setting cohort, no selection was made based on histology, demographics, time or any factor other than sample availability and quality control criteria (adequacy of biopsy and RNA yield).

A robust statistical pipeline was specifically developed for the BIOMARGIN study in R (R Core Team (10), version 1.0.44), as an extension of the biosigner R package (11). Briefly, in the discovery phase, we carried out a feature selection strategy including univariate and multivariate approaches. For univariate selection, we used a nonparametric test (Wilcoxon-Mann-Whitney test) and a parametric test (Student’s t-test) and computed the results into an “univariate score”. To rank the variables, we considered their false discovery rate (FDR)-adjusted P-value with a significance cutoff of 0.05. For multivariate selection, we carried out five approaches comprising sparse partial least squares (SPLS) (12), support vector machine-recursive feature elimination (SVM-RFE) (13), random forest-recursive feature elimination (RF-RFE) (14), elastic net and shrunken centroids (15). Principal component analysis was used to identify statistical anomalies, resulting in the exclusion of two outliers in the Discovery cohort (BIOS-03-0023, BIOS-06-0238). The prediction scores of each miRNA transcript in each multivariate model were integrated to yield a “multivariate score”. By combining the results of these univariate and multivariate analyses, we were able to rank and select a subset of miRNA candidates to be part of the extended list to be quantified during the next phase.

For this post-hoc analysis focusing on MVI pathways, the median normalized expression of each miRNA was compared between cases and controls using the Wilcoxon test in all three cohorts. We controlled for multiple testing using the Benjamini & Hochberg false discovery rate (FDR) method.

The individual histological lesion scores were semi-quantitatively assessed according to the Banff 2019 criteria (16). The term “microvascular inflammation” (MVI) was used to refer to any degree of combination of glomerulitis (g) and/or peritubular capillaritis (ptc). All ABMR cases met the first 2 criteria of the Banff 2015 or Banff 2017 classification (histologic evidence of acute tissue injury and evidence of current/recent antibody interaction with vascular endothelium), but not all met the third criterion [serologic evidence of donor-specific antibodies (DSA) and/or C4d staining]. DSAs after transplantation were determined per local center practice, with DSA positivity defined as detectable donor-specific serum anti-HLA antibodies with a mean fluorescence intensity (MFI) value of 500 at the time of biopsy or any time before.

Two needle cores were taken at each kidney allograft biopsy. One was used for conventional histological grading and at least half of the other was immediately stored in Allprotect Tissue Reagent (Qiagen Benelux BV, Venlo, The Netherlands). The Allprotect tubes were stored at 4°C (minimum 24 hours to a maximum of 72 hours), and then stored at -20°C until RNA extraction. Total RNA was isolated from the kidney allograft biopsy specimens using the Allprep DNA/RNA/miRNA Universal Kit (Qiagen Benelux BV) on a QIAcube instrument (Qiagen Benelux BV). The quantity (absorbance at 260 nm) and purity (ratio of the absorbance at 230, 260, and 280 nm) of the isolated RNA were measured using the NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific/Life Technologies Europe BV, Ghent, Belgium). RNA integrity number (RIN) was evaluated with the Eukaryote nano/pico RNA Kit (Agilent Technologies Belgium NV, Diegem, Belgium) on the Bioanalyzer 2100 instrument (Agilent Technologies BelgiumNV). Quality control indices were not significantly different between the MVI and the No MVI groups [RIN (mean ± SD): 5.6 ± 1.96 vs 5.4 ± 1.97, P=0.66; 260/280 ratio (mean ± SD): 1.87 ± 0.06 vs 1.87 ± 0.11, P=0.40]. For RNA samples with a low RIN (<6.5), we excluded samples with a 260/280 ratio <1.6. The extracted RNA was subsequently split and stored at -80°C. Half of the RNA extract was used for miRNA profiling and the other half for mRNA transcriptomic analysis. The associated data are available from the Gene Expression Omnibus (GEO, GSE179772, https://www.ncbi.nlm.nih.gov/geo).

Specific reverse transcription of 150 ng of total RNA was performed using Megaplex™ RT Primers Human Pool A v2.1 (Step 1) or Custom RT Primers Human Pool (Step 2 and 3) (Thermo Fisher Scientific, Les Ulis, France), on a Veriti Thermal Cycler (Applied Biosystems™, Thermo Fisher Scientific). No complementary (c)DNA preamplification was required. miRNA expression was assessed by qPCR, using TaqMan^® Array Cards (Applied Biosystems™, Thermo Fisher Scientific), where the primers and probe of each miRNA were spotted and dried in duplicate wells of a microfluidic card. We used TaqMan^® Array Human MicroRNA A+B Card Sets v3.0 for the discovery cohort (a two-card set enabling quantitation of 754 unselected human miRNAs) and Custom TaqMan^® Array MicroRNA Cards for the selection and validation cohorts. Data analysis was performed by using Expression Suite software version 1.0.3. Assays with a cycle threshold (CT) values >35 were considered to be not expressed. RNU44 and RNU48 showed consistent expression in all the samples [mean (RNU44+RNU48) coefficient of variation was 3.73% in Array Card A and 3.79% in Array Card B] and were used as housekeeping genes for data normalization across all sub-studies.

Transcriptomic analysis was performed by microarray as already described (17). Briefly, total RNA extracted from the biopsy samples was first amplified and biotinylated to complementary RNA (cRNA) using the GeneChip 39 IVT PLUS Reagent Kit (Affymetrix) and hybridized to the Affymetrix GeneChip Human Genome U133 Plus 2.0 Arrays (Affymetrix), which comprised 54,675 probe sets covering the whole genome. We used a Robust Multichip Average (RMA) normalization. The median normalized expression of each mRNA was compared between cases and controls using the Wilcoxon test. We controlled the inflation of risk alpha due to multiple testing by applying the Benjamini & Hochberg false discovery rate (FDR) method. The microarray data were handled in accordance with the Minimum Information About a Microarray Experiment guidelines.

Biological Interpretation Of Multiomics EXperiments (BIOMEX) software was used for transcriptomic data differential analysis, using Ensembl ID annotation (18). Ingenuity Pathway Analysis (IPA, Build: 478438M Content version: 44691306, Qiagen) and OMICSnet (19) were used to determine the sets of genes regulated by the selected miRNAs. OMICSnet was used to perform Gene Ontology analysis with the Reactome pathway database, using the complete pathway (20). Gene networks building was fully unsupervised and relied on the number of reported molecular interactions of each single gene with all other genes in the network. Genes presenting the highest number of interactions were automatically positioned at the center of the network, while the genes with a lower number of interactions were spread into the periphery. The cellular origin of selected miRNAs was investigated using the Fantom5 project (21), which provides a cap analysis of gene expression in human cells and tissues. For miRNAs expression analysis, we considered all kidney-related structural cells and all circulating hematopoietic cells, with no further selection.

Previously published human single-cell RNA sequencing (scRNA-seq) data from two ABMR biopsies and four healthy references corresponding to transplant surveillance biopsies were used. The associated raw counts or matrices were downloaded respectively from the GEO (GSE145927, https://www.ncbi.nlm.nih.gov/geo) (22) and Kidney Precision Medicine Project (KPMP, https://atlas.kpmp.org/repository). Raw gene expression matrices generated per sample were merged and analyzed using Seurat V4 package (23). The parameters applied to create and filter Seurat objects were as follows: inclusion of cells expressing between 500 and 10000 detected genes and less than 25% mitochondrial transcripts. After filtering, all objects were integrated using the SCTransform integration workflow on Seurat (24). Briefly, 3000 features were used for integration using PrepSCTIntegration function, and the FindIntegrationAnchors function was used to limit batch effect. A full Uniform Manifold Approximation and Projection (UMAP) data set was generated using the DimPlot function and the top 16 principal components. Clusters were built using the FindNeighbors and FindClusters functions with a 0.6 resolution. Cluster identification was performed using the FeaturePlot function by evaluating the expression of specific markers in each cluster. Dot plots were generated using the DotPlot and FeaturePlot plotting functions, with normalized counts in the RNA assay as input data.

Patient and donor characteristics were described by numbers, percentages and frequencies for categorical variables. We reported continuous variables using means and standard deviations (SD) if normally distributed, medians and interquartile ranges (IQR) for variables with a skewed distribution. We compared their characteristics using the Fisher or the Wilcoxon test when appropriate.

We used the Kruskal-Wallis test to compare the median miRNA expression according to MVI intensity, and the Spearman test for correlation analysis. We considered statistical significance for P-values less than 0.05. We used GraphPad Prism (version 8; GraphPad Software, San Diego, CA) and RStudio (version 4.0.3) for statistical analysis and data presentation. The following R packages were used: heatmap3 (heatmap3_1.1.9) for heatmap analyses, fmsb (fmsb_0.7.0) and scales (scales_1.1.1) packages for radar plots, ggplot2 (ggplot2_3.3.5) for bar graphs, and corrplot (corrplot_0.90) for correlation matrix analysis.

Between April 2013 and July 2015, 646 patients enrolled in the BIOMARGIN study provided 716 indication or surveillance biopsies. Inadequate biopsies (N=49), biopsies with no paired sample for research use (N=135) or samples not passing molecular biology QC (N=49) were excluded, leaving 483 samples from 441 individuals available for this study ( Figure 1A ). Baseline and clinical characteristics of the study groups are provided in Supplementary Tables S1-3 , and the histological characteristics of the biopsies are provided in Supplementary Tables S4-6 .

In the discovery cohort, a total of 27 patients met the MVI composite endpoint (MVI+), defined as clinical ABMR diagnosis and/or any level of MVI (g and/or ptc >0). The other 59 patients (MVI-) presented various diagnoses comprising T-cell mediated rejection (TCMR, N=8), interstitial fibrosis and tubular atrophy (IFTA, N=20), and normal biopsies (N=31), but no MVI. From the 754 miRNA candidates, 18 miRNAs were differently expressed between MVI+ and MVI- groups ( Figure 1B ). In the selection cohort, a restricted list of 143 miRNAs were quantified in 99 samples. In this selection cohort, 14 miRNAs were differently expressed between cases with (N=28) and without (N=71) ABMR ( Figure 1C ). Finally, in the largest trans-sectional cohort, 38 miRNAs were quantified in 298 samples, of which 12 miRNAs were differentially expressed between ABMR (N=29) vs no ABMR (N=269) biopsies ( Figure 1D ).

Next, we intersected the 3 lists of differentially expressed miRNAs, and identified 6 miRNAs which were consistently associated with MVI or ABMR across the 3 independent cohorts ( Figure 2A ). miR-142-3p, miR-150-5p, miR-155-5p, miR-222-3p and miR-223-3p were overexpressed in MVI/ABMR, while miR-139-5p expression was decreased. Radar plots ( Figure 2B ) display for each step the differential expression between case and controls of the 6 miRNAs and support the robustness of the miRNA profiles across the three independent cohorts.

Next, we studied how the 6 miRNAs’ expression is associated to individual histological lesions across all 483 allograft biopsies taken together. The expression of the 5 up-regulated miRNAs gradually increased with MVI lesions intensity, while the expression of miR-139-5p was negatively correlated ( Figure 2C ). MiR-139-5p/142-3p/150-5p/155-5p/223-3p were not only highly associated with ABMR-related features (g, ptc, C4d deposition), but also with TCMR-related lesions (i, t) and IFTA lesions (ct, ci), suggesting their involvement in common, non ABMR-specific, inflammatory and scarring processes. Alternately, some intrinsic collinearity might by nature exist between the lesions, resulting in our inability to identify true specificity for one lesion. MiR-222-3p did not correlate with histological lesions of TCMR or IFTA ( Figure 2D ).

Next, we evaluated the differentially expressed genes between MVI+ and MVI- cases, in the previously published discovery set (25). Of the 54,675 probes, a total of 2755 differentially expressed genes (DEGs) were identified between the MVI+ and MVI- groups, comprising 1085 down-regulated and 1670 up-regulated genes ( Figure 3A ). The most differentially expressed mRNAs in MVI+ samples were CXCL9 and CXCL10, in line with previous reports (17, 26) that these genes are the most overexpressed in ABMR. Next, we used IPA and OMICSnet to identify in silico the putative target genes for each miRNA. A total of 4504 mRNAs were predicted to be targeted by at least one of the 6 miRNAs. Then, we integrated these 4504 predicted targets with the DEGs in allograft biopsies. Integrated miRNA/mRNA analysis identified 61 DEGs targeted by miR-139-5p that were significantly increased in MVI+ biopsies. It also identified 28, 58, 52, 43 and 27 DEGs significantly decreased in MVI+ biopsies that were respectively targeted by miR-142-3p, miR-150-5p, miR-155-5p, miR-222-3p, and miR-223-3p ( Figure 3B ).

The cellular origin of the six miRNAs was characterized in various renal cell and hematopoietic cell subtypes using an miRNA atlas available online [Fantom database (27)]. Unsupervised clustering discriminated renal cell-derived miRNAs from immune cell-derived miRNAs ( Figure 4A ). It also revealed one or two dominant cell types for each miRNA. MiR-139-5p appeared to be mainly expressed by glomerular endothelial cells (EC). MiR-222-3p, although its expression was not correlated with tubular inflammation (Banff lesions “t”, Figure 2D ), was primarily related to renal epithelial cells. The four remaining miRNAs were primarily expressed in peripheral blood mononuclear cells (PBMC) or granulocytes, with miR-142-3p enriched in NK cells and monocytes, miR-150-5p in CD4+ and CD8+ T cells, miR-155-5p in CD19+ B cells and macrophages, and miR-223-3p in neutrophils.

To characterize the potential role of each MVI-associated miRNAs, we analyzed publicly available sc-RNAseq data from two renal allograft biopsies with ABMR and high MVI scores (g2, ptc3) (22). These datasets were integrated with 4 healthy reference biopsies with no MVI. After quality control and filtering, 31545 cells were detected and unsupervised clustering revealed 12 clusters ( Figure 4B ). We used well established markers (28) to identify all the main resident and infiltrating cell subtypes and label them ( Supplementary Figure S2 ). Next, we focused on the cell populations from which our six miRNAs were thought to originate from. In this aim, we stratified the cells into 4 major clusters corresponding to endothelium, epithelium, lymphocytes and APCs ( Figure 4C ). In line with previous analysis of these data (22), neither neutrophilic nor NK cell populations were detected in the sc-RNASeq datasets.

MiR-139-5p was the only downregulated miRNA in our study in cases with MVI/ABMR ( Figure 2A ) and was found to be mainly originating from glomerular EC ( Figure 4A ). Sixty-one of its target genes were confirmed to be upregulated in bulk microarray from MVI+ cases ( Figure 3B ) and their expression was further investigated at single-cell resolution in the 3 previously identified EC subclusters ( Figures 4B, C ). As depicted in Figure 5A , the activated EC cluster was only detectable in MVI+ samples. Besides, half of the 61 genes (N=30) were consistently increased in sc-RNAseq MVI+ biopsies, regardless of the three EC subclusters. Among those, GBP5 was the most increased gene. Next, gene ontology analysis using these 30 endothelial-derived transcripts identified UBE2D1 and YWHAG as central players in the gene network ( Figure 5B ). More specifically, immune-related pathways and major histocompatibility complex (MHC)-mediated antigen processing and presentation pathways were enriched ( Figure 5C ). Pseudotime analysis was performed to better characterize gene expression changes in the course of EC activation. Across cells state trajectory ( Figures 5D, E ), UBE2D1, YWHAG and GBP5 continuously increased during MVI associated endothelial activation. Very similar trajectories were found for MHC-related genes HLA-A, -B and -DRA but also for bona-fide endothelial activation and inflammation markers such as VCAM1, CXCL9 and CXCL10, suggesting that all these genes are concomitantly expressed during endothelial activation. Altogether, these results suggest that miR-139-5p is decreased during MVI lesions and no longer repress the MHC antigen presentation pathway in activated EC.

Next, we focused on miR-222-3p, an upregulated miRNA in MVI samples with a renal epithelial cell origin and a more exclusive correlation to ABMR-related histological lesions ( Figures 2A, D , 4A ). In bulk allograft tissue, the expression of 43 of its target genes was downregulated ( Figure 3B ), and further evaluated at a single-cell resolution in the different renal epithelial populations we previously identified ( Figures 4B, C ). A massive decrease of 41 of these 43 (95%) genes across all identified renal epithelial cell clusters suggests a gene-silencing profile in MVI ( Figure 6A ). Interestingly, 4 genes involved in renal physiological pathways (ERBB4, ATP1B1, TIMM50 and KIF3B), but also immune-related genes (TRAP1 and PEBP1) were identified as central genes in the network ( Figure 6B ). Gene ontology of the 43 genes revealed an enrichment in ErbB signaling pathways but also in immune-related pathways and cytokine response ( Figure 6C ) and ERBB4 expression was further confirmed across all renal cells ( Figure 6D ). Altogether, these results suggest that during MVI, a high level of miR-222-3p in renal epithelial cells associates with strong repression of genes involved both in renal electrolyte homeostasis pathways and immune-related pathways in the epithelium.

Similarly, we investigated the regulatory role of (increased) miR-150-5p during MVI in the T-cell cluster, wherein it was most dominantly expressed ( Figure 4A ). We plotted the 58 decreased MVI-associated DEGs identified in bulk allograft biopsy samples ( Figure 3B ) and we confronted it to gene expression in the sc-RNAseq T-cell cluster ( Figure 7A ). The proportion of genes both targeted by miR150-5p and specifically underexpressed in T-cells during MVI was surprisingly low (15/58, 25.9%) suggesting that the global decrease of these genes in vivo was not exclusively due to the T lymphocyte compartment. Nevertheless, from this highly selected list of 15 genes and using OMICSnet platform, we constructed a gene network ( Figure 7B ) and performed gene ontology analysis ( Figure 7C ). We observed that miR-150-5p overexpression was associated with NF-κB regulation in T lymphocytes during MVI.

Finally, the exact same strategy was applied for the 52 miR-155-5p-target genes, underexpressed during MVI ( Figure 3B ). In contrast with miR-150-5p in T lymphocytes, we observed that almost two-third of the miR-155-5p targeted genes were consistently decreased in the APC sc-RNAseq cluster (31/52, 59.6%, Figure 8A ). Gene network analysis revealed that AIFM1 and CIAO1 were central genes ( Figure 8B ). Moreover, gene ontology revealed a strong enrichment in pathways related to mRNA splicing ( Figure 8C ).