ABSTRACT
Background Organ availability limits kidney transplantation, the best treatment for end-stage kidney disease. Deceased donor acceptance criteria have been relaxed to include older donors with higher risk of inferior posttransplant outcomes. Donor age, although significantly correlates with transplant outcomes, lacks granularity in predicting graft dysfunction. Better characterization of the biological mechanisms associated with deceased donor organ damage and specifically predictive of transplant outcome in recipients is key to developing new assessment criteria for donor kidneys and developing function-preserving interventions.
Methods 185 deceased donor pretransplant biopsies with clinical and demographic donor and recipient metadata were obtained from the Quality in Organ Donation biobank (QUOD), selected on the basis of 12-month paired posttransplant function and deep proteomic profiles acquired by mass spectrometry. Using a 2/3rd:1/3rd training:test data split, sampling equally across posttransplant function, we applied machine learning feature selection followed by protein-wise relaxed LASSO regression modeling, assessing the performance of the final set of protein models on the test data. Western blotting validated protein changes, and the biological relevance of the final set of protein models was externally validated by contextualization against a published dataset of human healthy and disease kidney transcriptomes.
Results Our analysis revealed 144 proteins carrying outcome-predictive information, all of which showed donor-age modulated associations with posttransplant function, as opposed to age and protein/gene effects being independent terms. Observed associations with inflammatory, metabolic, protein processing and cell cycle pathways suggest biological targets for possible interventions pretransplant. Contextualization of our results against external spatial transcriptomic data suggest a sub-nephrotic spatial localization of the predictive signal.
Conclusions Integrating kidney proteome information with clinical metadata enhances the resolution of donor kidney quality stratification, and the highlighted biological mechanisms open new research directions in developing predictive models and novel interventions during donor management or preservation to improve kidney transplant outcome.
SIGNIFICANCE STATEMENT Currently, organ quality assessment pretransplant relies on key factors such as donor age or clinical information, these lack granularity in depicting graft susceptibility and capacity to function posttransplant. A high-resolution proteomic profiling of 185 pretransplant biopsies of kidneys with known posttransplant function and complete metadata was performed. Integration of donor kidney proteomes with 56 clinical metadata variables using regularized regression modelling resulted in enhancing the resolution of donor kidney quality stratification. Immunometabolic and catabolic processes contributed to donor kidney susceptibility and worse transplant outcomes in an age modulated pattern, validated by western blotting. Comparison of kidney proteomes with a recent transcriptomics dataset of healthy and diseased kidneys provide an additional special single cell resolution to the findings of this study.
INTRODUCTION
Kidney transplantation is the optimal treatment for end-stage kidney disease. Compared to dialysis, transplantation increases life-expectancy, improves quality of life and is cost-effective. Limited availability of suitable donor kidneys impedes treatment of chronic kidney disease, and often prolongs dialysis, increasing morbidity and mortality. Deceased donor organ shortages, living donation decline in some countries and emerging ageing populations drive increased utilization of older deceased donor kidneys, now comprising more than half of offered organs1,2.
Ageing associates with time-dependent decline of organ function, evidenced in kidneys by histologic lesions, such as tubular atrophy, interstitial fibrosis, glomerulosclerosis, and arteriosclerosis. Older kidneys demonstrate fewer functioning glomeruli, less renal mass, podocyte dysfunction, and impaired cellular repair3. Glomerular diseases are more common and associated with worse outcomes in older patients4. Age accelerates the transition from Acute Kidney Injury (AKI) to chronic injury5 and is an independent risk factor of graft dysfunction and loss for deceased donor kidneys6; furthermore, older donors are more likely to suffer from additional risk factors such as diabetes, hypertension or cardiovascular disease.
Donor age is incorporated in clinical scoring algorithms to inform kidney allocation decisions7,8, but is insufficient to consistently predict transplant outcomes. Current front-line models incorporating further clinical factors such as terminal serum creatinine, history of hypertension and diabetes8,9 show consistent performance across demographics but lack granular predictive accuracy10.
Molecular analyses of biopsies plausibly offer higher resolution assessment of organ state; but require ‘big picture’ understanding of mechanisms associated specifically with poor outcome, rather than immediate (but potentially recoverable) acute injury. Deceased donors are frequently assessed as having sustained damage (i.e. AKI) based on serum creatinine levels11, however this metric poorly associates with longer term outcomes11–14.
Mass spectrometry (MS) proteomic studies can provide such a ‘big picture’, but have heretofore lacked cohort capacity to represent demographic diversity15. Advances in high-throughput techniques16 now allow sensitivity and depth without sacrificing throughput capacity. Developments in machine learning and nonlinear regression analyses furthermore offer tools to extract maximal knowledge from limited size experimental cohorts, with applications in disease staging, disease recurrence prediction, treatment response monitoring, and biomarker identification17,18.
Integration of deep proteomic profiles with heterogenous clinical and demographic factors using modern statistical tools can empower the next steps toward precision medicine19. Here, we benefit from the granularity provided by our MS-based proteomic profiling to report age– and immunometabolism-related proteomic signatures in pre-implantation kidney biopsies associated with transplant outcomes.
METHODS
Study Design
Deceased donor pre-transplantation kidney biopsies (n=186; 1 sample excluded during data processing) were obtained from the Quality in Organ Donor (QUOD) biobank, a national multi-center UK wide bioresource of deceased donor clinical samples acquired during donor management and organ procurement. Biopsies were obtained from Donation after Brain Death (DBD) donors and Donation after Circulatory Death (DCD) donors at the back table immediately after kidney procurement.
Selection of biopsies was based on paired 12-month post-transplant outcomes. To minimize the impact of recipient factors, we only included kidneys for which the contralateral kidney was transplanted with similar outcome. Kidneys were selected to cover the outcome continuum i.e. the range of estimated Glomerular Filtration Rate (eGFR) in the recipient at 12 months posttransplant (henceforth, ‘eGFR12’), from primary non-function to eGFR>80 ml/min/1.73 m2, excluding pediatric donors. Samples were linked to corresponding donor and recipient demographic and clinical metadata, provided by NHS Blood and Transplant National Registry.
Study Approval and Ethics statement
Informed consent from donor families was obtained prior to sample procurement. Collection of QUOD samples and research ethics approval was provided by QUOD (NW/18/0187).
Protein Extraction from Renal Tissue Specimens
Deceased donor biopsies were procured, handled and stored according to consistent, predefined collection protocols designed to minimize pre-analytical variability. Donor kidney biopsies were collected ex situ immediately after flush-out and procurement of the kidney in the donor hospital on the back-table. Biopsies were obtained from the upper pole of the donor kidney cortex during back table preparation, using a 23mm needle biopsy gun. The obtained biopsies were divided in two, with one half stored in RNAlater (Thermo Scientific, Illinois, USA), then liquid nitrogen and the other half stored in formalin.
RNAlater biopsy samples were homogenized in RIPA buffer (89900, Thermo Scientific, Illinois, USA) with protease and phosphatase inhibitor cocktail (1861280, Thermo Scientific, Illinois, USA) using a bead beater (Biorad, Hertfordshire, UK) at 6500rpm for three cycles of 40 seconds with intermediate cooling on wet ice between cycles. Biopsy protein concentration was determined using a Pierce Bicinchoninic Acid protein assay kit (23227 Thermo Scientific, Illinois, USA).
In-Solution Trypsin Digestion
50 µg of protein homogenates were prepared. Disulfide bonds were reduced by adding 200 mM of DTT (Sigma) to a final concentration of 5 mM for 60 mins at room temperature. Free cysteine residues were alkylated by adding 200 mM of iodoacetamide (Sigma) to a final concentration of 20 mM and incubated for 60 mins at RT in the dark.
The samples were topped up to 200 µl with 6 M urea, 100 mM TrisHCl pH 8.5. Methanol/chloroform protein precipitation was used to remove detergents before tryptic digestion. In brief, 600 µl of Methanol and 150 µl of Chloroform were added and mixed. Then 450 µl of MilliQ-H2O was added and then centrifuged for 1min at 12,000 g. The upper aqueous layer was carefully removed without disturbing formed protein pellets between layers and 450 µl of Methanol was added and then centrifuged at 12,000 g for 5 min. The supernatant was removed and the protein pellets resuspended in 50 µl of 6 M Urea, 100 mM TrisHCl, pH 8.5. Urea concentration was reduced to 1 M by adding 250 µl of MillQ-H2O. Samples were digested at 37 °C overnight with Trypsin added at a 1:50 ratio (trypsin:protein). Tryptic peptides were acidified and purified using Sep-Pak C18 cartridges (WAT020515, Waters, Wilmslow, UK) and dried by Speed Vac centrifugation. Pellets were resuspended in 80 µl of resuspension buffer A (98 % MilliQ-H2O, 2 % acetonitrile 0.1 % formic acid) for LC-MS/MS analysis.
Mass Spectrometry Analysis
Generation of Fractionated Pool for Spectral Library Reference
A highly fractionated spectral library was generated as a standard reference for the subsequent analysis of individual samples. To create this spectral library, a pool sample was prepared by combining 2 µl of each tryptic digest sample prepared above. Fractionation of the pooled sample was performed using offline high-pH reverse-phase HPLC on an XBridge BEH C18 XP column (3 × 150 mm, 2.5 μm pore size, Waters no. 186006710) over a 100-minute gradient (Buffer A: water, pH10 with ammonium hydroxide. Buffer B: 90 % acetonitrile, 10 % water, pH10 with ammonium hydroxide) with fractions collected every 2 minutes. The fractions of the pool sample for the spectral library were analyzed by nanoLC-MS/MS using a Dionex Ultimate 3000 using a 75 µm x 500 mm (2 µm particle size) C18 EASY-Spray column at 250 nL/min (Thermo Scientific), coupled to an Orbitrap Fusion Lumos mass spectrometer. Peptides were separated using a 60-minute linear gradient from 2-35 % buffer B (Buffer A: 5 % DMSO, 0.1 % formic acid in water. Buffer B: 5 % DMSO, 0.1 % formic acid in acetonitrile). The samples were analyzed on the mass spectrometer in Data-Dependent Acquisition mode. MS1 scans were acquired in the Orbitrap with an m/z range of 400 – 1500 m/z at a resolution of 120,000 and an AGC target of 4 x 105. Precursors between charge states 2+ and 7+ were selected for HCD fragmentation using the Advanced Precursor Determination option with an intensity threshold of 2.5 x 104. Selected precursors were isolated using the quadrupole with a 1.6 m/z isolation window and fragmented using HCD with a normalized collision energy of 30%. MS2 spectra were acquired in the Orbitrap using a resolution of 30,000, a maximum injection time of 54 ms and an AGC target of 5 x 104.
LC-MS/MS Analysis of Individual Biological Samples
The individual samples were analyzed on the same nanoLC-MS/MS system as above. A 45-minute linear gradient was used from 2-35 % buffer B with the same buffer composition as above. In contrast to the pool samples, the individual samples were analyzed using the SWATH DIA method 16. The sample analysis order was randomized, and analyses of aliquots of the pooled sample used for library generation were scheduled every 20 runs throughout the sequence as a quality control. MS1 scans were acquired in the Orbitrap with an m/z range of 350 – 1650 m/z at a resolution of 120,000, an AGC target of 4 x 105 and a 3 second cycle time. MS2 scans were then acquired in stepped isolation windows with a 1 Th overlap between each window; first from 350-380 m/z up to 930-960 m/z in increments of 30 Th (i.e. 21 scans with midpoints of 365, 394 … 916, 945), followed by a 100 Th window scan of 959-1059 m/z and a 592 Th scan of 1058-1650 m/z.
Fragmentation of these windows was performed using a normalized collision energy of 25 % with a stepped collision energy of 10 %. MS2 spectra were acquired in the Orbitrap using a resolution of 30,000, a maximum injection time of 54 ms, an AGC target of 5 x 104 and a scan range of 360-1650 m/z.
Western Blot Validation
Samples were compared as 5 Upper Tertile (UT) versus 5 Lower Tertile (LT)within a single age category (younger; total n=20, or older; total n=20) on each gel. Kidney homogenates at 1μg/μL of protein were denatured at 90LIC in Laemmli buffer and separated using 4-12% pre-cast Bis-Tris gels (Thermo Fisher Scientific, Massachusetts, USA) in the case of VTN and APOE western blots and 16% pre-cast Tricine gels (Thermo Fisher Scientific, Massachusetts, USA) in the case of PREB and CST3 western blots, in both cases using MOPS running buffer and followed by transfer onto a PVDF membrane (Thermo Fisher Scientific, Massachusetts, USA). Membranes were then blocked in Intercept (TBS) Blocking Buffer (LiCOR P/N 927-60001) and incubated with primary antibodies overnight at 4LIC. Antibodies used were rabbit polyclonal anti-PREB (Thermo Fisher Scientific PA5-53125), mouse polyclonal anti-Vitronectin (R&D, MAB2349), rabbit monoclonal anti-APOE (16H22L18) (Thermo Fisher Scientific, 701241), rabbit monoclonal anti-CST3 (D6U3E) (Cell Signaling Technology, 24840) and mouse monoclonal anti-β-actin (Sigma, A5316) for protein normalization. The membranes were then incubated in species appropriate secondary antibodies (RDye® 800CW Goat anti-Mouse IgG Secondary Antibody and IRDye® 680RD Goat Anti-Rabbit IgG Secondary Antibody, LiCOR, USA). Finally, blots were imaged using the LI-COR Odyssey system (LiCOR, Nebraska, USA). Analysis and quantification were performed using LiCOR Image Studio (Version 2.2). Protein expression was normalized against β-actin and the fold change was calculated relative to the per-gel mean value of the Upper Tertile group. For each protein, fold changes relative to UT were combined across all gels within each age category, log2 transformed, and the difference between UT and LT calculated by t-testing without assumption of equal variance.
Statistical Analyses
Proteomic Data Processing
Both the DDA fractionated pools and all SWATH samples and pool data were analyzed simultaneously using DIA-NN v1.7.12 20. MS/MS spectra from the DDA pools were searched against UniProt Reference Homo sapiens database (retrieved 15/10/2020) with default settings, including 1 missed cleavage, oxidation of methionine as a variable modification and carbamidomethylation of cysteine as a fixed modification, to generate a spectral library. This was then cross referenced against the SWATH data to quantify peptides, and thus proteins, within each SWATH sample and pool run. Precursor FDR was set to 1%.
Subsequent analysis of the DIA-NN output and integration with the clinical data was performed in R v4.0.2. One sample had extremely low average intensity and more than 15% missing values, so was eliminated from the analysis. Across the remaining 185 samples and pools, proteins with more than 50% missing values were eliminated. Intensity values were transformed by VSN (R package ‘vsn’ 21), and the remaining missing values were imputed by a k-nearest neighbor approach (R package ‘impute’ 22).
Clinical Variable Preprocessing
Clinical variables were subdivided into Donation after Brain Death (DBD)-specific variables, Donation after Cardiac Death (DCD)-specific variables and jointly applicable variables (see annotation, Supplementary Table 1). DBD– and DCD-specific variables were considered for the purposes of assessing correlation with outcome but were not considered for statistical analyses applied to the full dataset.
For significance testing of outcome differences between donor and recipient characteristics, and for all modelling, eGFR values for both donor measurements and recipient measurements at 12-month post-transplantation were stratified into tertiles – Lower Tertile (UT; eGFR < 40), Intermediate Tertile (IT; 40 ≤ eGFR < 60) and Upper Tertile (UT; eGFR ≥ 60) based on their distribution. Missing values in the continuous and tertile-stratified forms of the variables were imputed for each sample by finding the ‘set of neighbors’ by 3-month eGFR, i.e. the 10 samples with the closest recipient eGFR at 3 months posttransplant, then taking the mean eGFR12 (ignoring any missing values) within that set of neighbors.
We curated the non-DBD/DCD specific clinical variables to remove those inappropriate for modelling (single-instance categories, variables which duplicated another variable in different units); the final list is given in Supplementary Table 2. HLA mismatches were modelled as ‘Match Grade’ with levels None; no mismatches, Favorable; no DR and one or fewer B mismatches, and Non-Favorable; at least one DR or two B mismatches (Supplementary Figure S2 and Supplementary Table ST1).
For curated clinical variables, missing values were imputed for each sample in a similar manner. A set of neighbors was built by finding the 10 samples with the smallest Gower’s distance across all variables (multiple data types). We then took a relevant summary statistic of the non-missing values within the set of neighbors for each imputed value. For numeric values, the mean was taken. For categorical values, the most frequent category across the whole dataset was taken, breaking ties by iteratively shrinking the set of neighbors by removing the most distant from the sample to be imputed, until there was a single most frequent category.
Outcome Subgroup Comparison
Differences between outcome subgroups within DBD and DCD donor types (separately) were tested using Fisher’s exact test for categorical variables and t-tests for numeric variables as appropriate.
Clinical variable association was calculated depending on the type of each variable in pairwise comparisons. Numeric-numeric comparisons were by Pearson’s r. Numeric-categorical comparisons were by the square root of the proportion of variance explained (eta) from an ANOVA model where the numeric variable was taken as the response. Categorical-categorical comparisons were by Cramér’s V. Hierarchical clustering was performed using the single-linkage method to avoid composition of comparison types.
The proteomic and clinical data were matched in R by sample identifier. Based on an initial scree plot of the variance distribution, unsupervised clustering of the proteomic data was performed by k-means clustering with k=4. The number of initial random centroids and the maximum number of iterations were both set to 10.
Training and Test Data Split
We split our data into a training and test sets, sampling equally across stratified eGFR. Of the 185 samples, 2/3 were used for training (‘training set’, 118 kidneys), and 1/3 were used for testing (‘test set’, 61 kidneys). Six paired kidneys (3 pairs) were reserved from both sets and used to assess MS and model performance (‘QC set’, 6 kidneys).
Machine Learning Feature Selection
Using sample assigned to the training dataset, we modelled the combined dataset of all protein abundance data plus the curated donor clinical variables (Supplementary Table S2) against ranked 12 month posttransplant eGFR, using Prediction Rule Ensemble (PRE) modelling (R package ‘pre’23). We used the ‘gpe’ function allowing linear fits, decision tree fits and multivariate adaptive regression spline fits (via R package ‘earth’ 24) to be considered for each ensemble. Decision tree learning was set to be ‘random forest’-like, with 500 initial trees generated and boosting disabled to allow parallel processing. Other settings (including tree pruning and the parameters for linear and spline fits) were left as default; tree depth was therefore limited to the default 3 decision levels. Ensemble modelling was performed in an iterative manner. At each stage, modelling was performed on the set of predictors that excluded proteins listed in the final ensemble of all previous models. Clinical variables were never excluded. Modelling was repeated for 2000 iterations and all proteins appearing in any rule ensemble, or with high correlation (Pearson’s r>0.65) to any of these candidates, were selected for further analysis.
Regression Modelling
Based on the outcome of feature selection, we used regularized regression (relaxed LASSO, package ‘glmnet’, using the cv.glmnet function with 20-fold cross validation, selecting as the final model the model that maximized the shrinkage parameter lambda such that it was within 1 standard error of the lambda value that had the smallest cross-validation error, i.e. the ‘lambda.1se’ model) to regress each protein, donor age and protein:age terms against linearized outcome, defined as the proportion of eGFR12 values in the NHS Blood and Transplant National Registry between 2016 and 2021 inclusive (6 years) that were less than the observed eGFR12 (i.e. a normalization to population distribution quantile). An independent age-only model was also fit against the training dataset using the same parameters. For each protein model and the age-only model we calculated the prediction root-mean-square error (RMSE). Proteins with a higher RMSE than the independent donor age model, or whose model did not retain either the protein or protein:age term were discarded.
Functional Network Analysis
The final list of proteins resulting from regression modelling was queried using the Enrichr platform (R package ‘enrichR’25) against the Reactome 2022 database. Annotation term assignments were filtered at 5% FDR. From the set of candidates with at least one assigned term across all three databases, a connection graph was generated (R package ‘igraph’26) with nodes representing proteins and edges representing shared annotation terms. Nodes (proteins) were clustered using the walktrap community detection algorithm27 as implemented in igraph. Within-community enrichment for annotation terms was calculated by Fisher’s exact test.
Age-Protein-Outcome Clustering
For each protein, eGFR12 was predicted using the corresponding relaxed LASSO model for donor ages 20,21,22…80 and protein quantiles 0,0.01,0.02…1, resulting in a 61×101 matrix. Matrix-matrix distances were computed using the Frobenius distance (R package ‘StatPerMeCo’28), and proteins were clustered according to these pairwise distances by hierarchical clustering.
Spatial Correlation Analysis
Processed, normalized AKI and CKD spatial scRNA-seq expression data generated by Lake et al.29 were obtained from CELLxGENE (https://cellxgene.cziscience.com/collections/bcb61471-2a44-4d00-a0af-ff085512674c).
Proteins identified by feature selection and LASSO regression filtering were matched to the expression dataset by the ‘Genes’ column in the DIA-NN output. AKI:CKD correlation in matched versus unmatched genes in each characterized spatial location was compared using the approach described by Fisher30 (R package ‘cocor’31).
RESULTS
Selected Samples were Demographically Balanced Across Outcome Strata
Samples were selected to provide a balanced representation of the UK donor population (Table 1) and reproduced generally observed trends in terms of the correlation between eGFR12 and clinical variables; we found the strongest associations were donor age (Pearson’s r=-0.52), and recipient age (r=-0.28) (Supplementary Figure S1).
Donor kidney associated metadata. Samples are subdivided by donor type and by final assigned outcome tertile. Numerical variables are given ± standard deviation. Categorial variables are given alongside percentage of total cohort
For an interpretable analysis of key clinical factors, we considered our samples by eGFR12 stratum Lower Tertile (LT; eGFR12 < 40), Intermediate Tertile (IT; 40≤eGFR12<60), and Upper Tertile (UT; eGFR12≥60) all units ml/min/1.73 m2 (Figure 1). We investigated associations between clinical variables and stratified eGFR12 subgroups within each donor type, and between donor types within stratified eGFR12 subgroups (Supplementary Table ST1); the only variables with significant difference between outcome groups across donor types were the UK Kidney Donor Risk Index; UKKDRI8 (ANOVA F-test; DBD: p=1.298e-6; DCD: p=3.946e-7), donor age (ANOVA F-test; DBD: p=1.253e-9; DCD: p=1.196e-7), histories of hypertension (ANOVA F-test; DBD: p=0.0020; DCD: p=0.1069). Histories of diabetes (ANOVA F-test; p=0.6188; DCD: p=0.2348) and terminal serum creatinine levels (ANOVA F-test; DBD: p=0.6972; DCD: p=0.6448) were similar across outcome subgroups although the latter was higher in DBD than in DCD in the UT group (t-test; p=0.0443).
We selected biopsies from QUOD biobank taken from one kidney from each donor pair. Donor kidney samples were selected randomly from pairs where both recipients had similar outcomes. The biopsy samples were subjected to proteomic analysis to yield a snapshot of the organ proteome at kidney retrieval. We analyzed donor characteristics and clinical variables and protein abundances in a combined model against recipient eGFR at 12 months posttransplant (eGFR12; units given in ml/min/1.73 m2).
Integration of Kidney Proteomes with Clinical Metadata by Rule Ensemble and Regression Modelling Identifies Outcome-Associated Proteins
Proteomic analysis quantified 2984 protein groups with 50% or less missing values (out of 7790 identified protein groups in total) over 185 samples and 20 interspersed sample pools (Supplementary Figure S2A). Analysis of sample pools showed minimal sample acquisition-related variance (squared mean pairwise Z-corrected Pearson’s r=0.94). The 3 pairs of kidneys showed high correlation of protein intensity values between donor pairs (Pearson’s r=0.71, 0.92 and 0.91; Supplementary Figure S2B).
We initially explored the proteomic data using Principal Component Analysis (PCA) to find underlying linear trends. Sample variance concentrated in the first two principal components (PC1: 20.01%; PC2: 13.38%; Figure S3). K-means clustering identified 4 distinct clusters (Figure S3; Figure S4A) which associated with donor type, with a preponderance of DBD samples towards Cluster 2 and a preponderance of DCD samples towards Cluster 4 (Figure S4B, upper left panel; p=0.0235), but did not associate strongly with recipient eGFR (p=0.4134), donor eGFR (p=0.1684), or donor age (p=0.7907) (Figure S4B, upper middle, upper right and lower left panels). There was a weakly significant association between cluster membership and donor BMI (p=0.0350) and with donor serum creatinine (p=0.0326) (Figure S4B, lower middle and lower right panels).
To assess individual protein relationships with outcome, we adopted a descriptive modelling approach, using a subset of our data (‘training set’; 2/3 of data) to find protein models which predicted eGFR12 better than donor clinical data alone, then assessing model performance against unseen (i.e. held out) test data (‘test set’, 1/3 of data). Six of the kidneys analyzed (3 pairs) were pairs from the same donor; these were reserved from both training and test data for quality control (‘QC set’, 6 kidneys). To create our train/test split, we randomly selected equal numbers of samples for the training set within each of the eGFR12 strata described above, to sample equally across outcome range.
To identify possible outcome-associated proteins, relevant clinical variables and potential protein-clinical variable interactions, we used iterative Prediction Rule Ensembles23 (PRE) learning on our training set to select features among the set of quantified proteins and donor type-independent donor clinical variables. Briefly, PRE identifies a minimal ensemble of explanatory rule (functions of variables) for a given response (here, ranked eGFR12) allowing for nonlinear associations and variable interactions. To build a catalog (rather than minimal set) of explanatory proteins, we repeated our analysis 2000 times, removing proteins identified in the rule ensemble at each iteration from the dataset for future iterations, retaining clinical variables.
This process generated 195 candidate proteins; we further supplemented this list with proteins that had high correlation (Pearson’s r>0.65) with any of those candidates; bringing the final list up to 255 candidates.
The only donor clinical variable term to feature consistently in rule ensembles was donor age, appearing in all 2000 ensembles generated. We individually tested each protein candidate for association with outcome, including protein, donor age, and age:protein interaction terms using LASSO regression32 to aim for the simplest explanatory models (i.e. that include only the minimum necessary terms). To ensure models would generalize beyond our data, for the final modelling we linearized the eGFR12 measurement against UK NHS Blood and Transplant National Registry data. We discarded candidates whose model discarded both the protein and age:protein terms, or had a higher root-mean-square error (RMSE) of prediction than a model built with donor age alone.
After filtering we had identified 144 proteins which associated with outcome (Supplementary Table 2). The mean RMSE of our protein models was 25.76 ml/min/1.73m2 in our training data, and fell slightly to 22.25 ml/min/1.73m2 in our test data, indicating the models were not over-fit (Supplementary Figure S5).
Functional Analysis of Outcome-Associated Proteins Reveals Immuno-Metabolic Pathway Clusters
We performed a network analysis of shared Reactome pathways (Figure 2A). Walktrap clustering revealed 4 major clusters of shared-pathway proteins (Table 2); Immune Regulation & Complement Activation, Metabolism, Protein Metabolism & Modification, and Cell Cycle. To validate our findings, we selected three representative proteins with available antibodies and known biological relevance (Vitronectin (VTN); fibrosis, Apolipoprotein E (APOE); CKD risk,
A: Shared Reactome pathway membership network analysis of filtered features. Nodes are colored by assigned cluster, and the clusters are annotated according to the top three most enriched pathways within each cluster. B: Western blots comparing younger (age ≤ 49; mean age 34) and older (age ≥ 58; mean age 64) donors between Upper Tertile (UT; eGFR12 ≥ 60) and Lower Tertile (LT; eGFR12 < 40) outcomes. As above, Left-Right: VTN, APOE, CST3. Top row: representative western blots (n=5 per group) from comparison of younger donors. Middle row: representative western blots (n=5 per group) from comparison of older donors. Bottom row: result values for all quantified samples relative to the UT mean. UT values in green, LT values in blue. Error bars indicate ±1 standard deviation; the central wider bar indicates mean. Significance stars indicate t-test comparison p-values (***: p < 0.001, *: p < 0.05).
All eGFR12 units given in ml/min/1.73 m2
Proteins in Figure 2A were clustered by pathway membership, forming 4 major clusters (those with more than a single member). We assigned summary labels to each cluster based on the top 3 pathways with shared membership in each cluster.
Cystatin C (CST3); nephron function) and confirmed that all three showed differences between low and high eGFR outcome sample subgroups by western blot (Figure 2B). Furthermore, all three proteins consistently predicted similar outcomes for each pair of kidneys within the triple-pair QC set, and correctly separated the two below-median outcome (eGFR12<50) pairs from the above-median outcome (eGFR12>50) pair (Figure S6).
Association of Proteins with Posttransplant Outcome is Modulated by Donor Age
For all 144 proteins, the minimal model retained an age:protein interaction term where the predictive effect of protein abundance was modulated by age, independent of the effect of age alone or protein abundance alone (Supplementary Table 2). To visualize these effects, we used each model to predict eGFR12 across increasing donor age and protein quantile. The interplay between protein quantile and age could be broadly grouped (by hierarchical clustering) into 3 clusters with protein to eGFR12 associations that showed different patterns across donor age (Figure 3), with about a third (31%) of candidates, including CST3, falling into a cluster with a profound effect of age on protein to eGFR12 association at high donor age.
Our modelling found that the effect of candidate protein levels on outcome changed with age. A: We quantile-normalized protein levels to a consistent scale and predicted outcome (12-month posttransplant eGFR) across a representative donor age range. Using this prediction matrix to compare the interaction with age between protein models, we found 3 major clusters (Left-Right; blue, pink and green dendrogram branches and corresponding plot background color) representing increasing levels of age modulation. Moving Left-Right across each plot, the eGFR12 (color scale, see key on right) contours across protein quantile (y axis) change as donor age (x axis) increases. B: Illustrative protein quantile-GFR12 plots are shown for each three proteins in each cluster. Traces are shown for donor ages between 20 and 80 as labelled on the right end of each trace.
Comparison to Spatial scRNA-Seq Data Reveals Localization of Outcome-Associated Signal
As an independent validation of our findings, we sought to contextualize them in the wider context of kidney damage. To do this, we compared the expression levels of our candidate protein set to a recent spatial scRNA-seq dataset comparing AKI and CKD29. This initial comparison revealed that all 144 of our candidates were matched to transcripts reported in this external dataset, and, consistent with our proteomic data. Two of our highlighted candidates associated with outcomes, APOE and CST3, were also highlighted as associated with fibrosis and inflammation damaged tubules in this dataset.
We then investigated how the expression levels of transcripts corresponding to our set of 144 proteins changed between AKI and CKD. We reasoned that protein and gene expression data that carries outcome-predictive information must reflect differences in organ injury.
Therefore, between different injurious scenarios, we would expect gene expression that is predictive (i.e. discriminatory) to have lower correlation than genes whose expression change was linked to general inflammation and stress response. We compared the AKI:CKD expression correlation of our set of candidates versus the expression correlation of all other transcripts in the dataset. Since the dataset is spatially resolved, we were able to make this comparison across regions of the nephron. We observed that our candidate protein set had stronger discriminatory power (i.e. significantly lower correlation; Fisher test, p=4.04e-11) than background in distal convoluted tubule epithelial cells and in surrounding structures, but not in leukocytes or in globally aggregated data (Figure 4).
We reasoned that proteins/genes which predict outcome should associate with specific facets of organ injury rather than general inflammatory and stress response, and therefore would be expected to have a lower expression correlation between acute and chronic damage scenarios than non-predictive other proteins/genes. We explored this in a spatially resolved public transcriptomic dataset and found that the degree of correlation was significantly lower for the set of 144 candidate predictors identified by our modelling approach, but that the signal appears to be spatially localized, primarily to distal convoluted tubule epithelial cells and the epithelia of neighboring tubules (indicated by red highlighting in right hand diagram).
This result gave us confidence that the 144-candidate set we report represents signatures of specific organ injury, and suggests that proteomic signals in distal convoluted tubule may be particularly representative of the extent to which injury impacts future functionality.
DISCUSSION
This study presents the first large unbiased human kidney tissue proteomics dataset for deceased donor kidneys. We examined the relationship between posttransplant kidney function with pretransplant kidney proteomes by analyzing 185 deceased donor kidneys with complete donor and recipient associated metadata. Although our study was not aiming to report a single biomarker, our large cohort, our unique machine learning approach, and finally alignment of proteomic data with a large spatial single cell transcriptomics dataset on healthy and diseased kidney substantiates our key findings. We describe an integration of kidney proteome with clinical metadata and show that modelling age modulation of proteomic signals provides enhanced resolution of donor kidney quality stratification and highlights biological mechanisms.
Our analysis reveals the importance of integrating across both subclinical and clinical data. Exploring our data using iterative PRE feature selection, a substantial number of proteins were revealed to be relevant, but only one clinical variable, donor age. Donor age is a key contributor in clinical decisions and is a strongly weighted term in extant kidney allocation scoring systems 8,33, but in the case of all reported candidate proteins, prediction models retained a protein:donor age interaction term. This second-order effect has not (to our knowledge) been explored in transplantation, and may be important to fully understanding molecular predictors.
A clinical variable notable by its omission from our association findings was donor type34. It is possible that outcome-predictive information ‘carried’ by donor type is also derivable from the proteome (proteomic data were able to partially distinguish samples based on type even when considering only variance across linear combinations of proteins; Figure S4). Without disputing donor type-specific mechanisms of kidney damage35, our data are consistent with the idea that location and type of damage may be a greater contributor towards recovery potential36 than mechanism of injury.
Kidney metabolism is altered as a result of biological stress occurring during donor management, in both DBDs and DCDs37. Within our final list of 144 proteins associated with outcome there is a common theme of implication in immune response to kidney injury (including both chronic injury, and acute injury) as a result of metabolic disruption, particularly responses associated with ischemia. We highlighted three proteins from our set of 144 outcome-associated candidates as a sanity check; all three are examples of known chronic damage indicators that we find here to have a donor organ pre-transplant predictive association.
Vitronectin (VTN) is a primary component of the extracellular matrix involved in cell adhesion, enhancing the activity of plasminogen activator inhibitor-1 and inhibition of the terminal complement pathway38, and is a potential fibrotic biomarker39.
Apolipoprotein E (APOE) stands out as having previously reported genetic allele age-related associations with disease and organ dysfunction including risk of Alzheimer’s Disease40 (with the strongest effect manifesting around age 6541), macular dysfunction, atherosclerosis and pulmonary scarring42,43, and evidence for shared allele risk across diseases44. In kidneys, APOE plays an important role in lipid metabolism to regulate the growth and survival of mesangial cells and preserve organ function45; it is a marker for outcome in transplant recipients46–48, and there is evidence for APOE genotype association with kidney dysfunction risk49–51, possibly manifested by lipidomic differences between allelic profiles52. We have previously reported small (not statistically significant) increases in APOE due to ischemic reperfusion injuries53 possibly explained by a role in senescence mediation54. There is existing evidence for similar allele dependent transplant outcome effects in another apolipoprotein (APOL1)55, suggesting that both APOE genotype and the broader apolipoprotein allelic profile may play an important role in posttransplant graft function.
Cystatin C (CST3) is particularly noteworthy as, measured in serum, it is a known and effective general biomarker for kidney function and has predictive power for outcomes in transplant recipients56–58. Our evidence indicates a further association between CST3 levels in donor kidney tissue and outcome; moreover, that this effect is age dependent, starting around age 40. Serum CST3 is relatively independent of age in children and young adults59, but may also show age association in later years60.
A potential hurdle in contextualizing our findings in existing orthogonal transcriptomic data was that proteomes and transcriptomes are generally poorly correlated in terms of abundance changes61. In order to avoid this problem, we considered how the pattern of injury scenario-related differences changed between predictive and background signals within the same transcriptomic dataset29. The spatial component of this dataset further allowed us to show that injury-specific signals within our 144-protein set localized most strongly to the distal convoluted tubule (DCT) epithelium. Damage to DCT epithelial cells can be predictive of outcome62; the interplay of proximal and distal injury is complex63, but previous work has suggested that DCT epithelial cells can play an important protective and damage repair role across the nephron, expressing survival and reparative factors in response to injury64. Indications that aging is associated with accumulation of senescent epithelial cells with a maladaptive stress response29 suggest a mechanism whereby aging related senescent cell accumulation could reduce the ability of DCT epithelia cells to protect and promote recovery of injury sustained by the proximal tubules. This would be consistent with our findings (immuno-metabolic and ischemia response functional association, and spatially localized response) and the negative donor age correlation with transplant outcome, and might explain recent successes with senescence inhibitors in treating CKD65, suggesting that routes to stimulate or enhance the DCT regenerative effect could be avenues for organ preservation.
Our list of outcome-associated candidates cannot be exhaustive. Practicalities of sample acquisition limited sampling of a wide range of outcomes outside the 30-60 donor age range, especially limited Upper Tertile outcome events in older donors. Organ allocation algorithms impose a close link between donor and recipient age in the sample cohort, so while we interpret these age-moderated effects in terms of organ resilience in older donors, it could also represent a greater ability to repair a given level of damage in younger recipients. Further, we consider only chronological donor age, rather than a more nuanced representation of the epigenomic biological clock 66, which may to account for some variation observed with respect to both donors and recipients.
It is immediately clear from our results that the strength of the donor age factor is enormous relative to any other protein or clinical effect; this age effect is liable to dominate any prediction weighting and reduce the accuracy of estimated protein contribution. A much larger cohort could mitigate this issue, especially if paired with variant sequencing to understand genetic diversity. Advances in high-throughput proteomics techniques continue to increase feasible cohort sizes67 but fundamental limitations on organ acquisition remain. Archiving at scale of clinical samples in bioresources such as the QUOD biobank to parallel advancements in big data analysis and interpretation platforms is therefore necessary for future development of granular evidence-based decision making.
However, given the limitations of our dataset, our protein models generalized well to samples unseen during model training (Figure S5), and showed promise for outcome-classification prediction (Figure S6) suggesting that, with a larger training cohort to address the caveats acknowledged above, a substantial component of early posttransplant outcome may be predictable solely using subclinical measurements moderated by donor age. The biological themes of the 144 proteins identified herein reinforce known immuno-metabolic mechanisms of kidney injury but raise interesting possibilities for further work, especially with regard to donor genetic background, and also suggest that the possibility of donor age-moderated weighting should be considered as a matter of course in future work.
SUPPLEMENTARY MATERIALS
Supplementary Table ST1: Clinical variable p-values for association with donor type and outcome.
Supplementary Table ST2: Summary of results for all candidate proteins.
Note that coefficient values provided are for linearized eGFR12, i.e. they relate to eGFR12 quantile ranging from 0 to 1.
Supplementary Figure S1: Donor and recipient clinical and demographic data association with recipient 12-month eGFR rank.
Single-linkage hierarchical clustering of curated, imputed clinical variables by relative association strength (taking distance as 1-association). The outcome variable (ranked recipient eGFR at 12 months post-transplantation) is highlighted in red.
Supplementary Figure S2: Protein quantification quality.
A: Missingness comparison: Proteins are shown ranked by the number of missing values across all samples and the twenty standard pools, excluding one run which was removed due to low signal. 2984 proteins had missing values in 50% or less runs.
B: Paired Kidney Comparison: Protein abundance values from paired kidneys (left/right) from 3 individual donors were compared, as these are effectively biological replicates. x axes: value in left kidney. y axes: value in right kidney. Inset: Pearson’s r correlation coefficient.
Supplementary Figure S3: Scree plot of variance represented by the first ten principle components in the proteomic data.
A reasonable number for the cluster parameter (k) supplied for k-means clustering (see Figure S4) lay between 3 and 5 based on the ‘elbow’ method’; we selected k=4.
Supplementary Figure S4: Unbiased analysis of pretransplant kidney proteomes and cluster associations.
A: Unbiased analysis of proteomic data by k-means clustering. Sample separation by Principal Component Analysis. Top Left: Samples were assigned to four clusters by k-means. Bottom & Right: There was a difference in the distribution of DBD and DCD donors across clusters, with the DBD donors being more heavily concentrated in Cluster 2 (‘+’ symbol; orange shading), and DCD in Cluster 4 (‘x’ symbol; pink shading)
B: There were no associations between proteome clusters and most donor and recipient factors, except for mildly significant differences in donor BMI and creatinine (selected comparisons shown; left-right, top-bottom: donor type, recipient 12-month posttransplant eGFR (outcome), donor eGFR, donor age, donor BMI, donor creatinine at retrieval).
Supplementary Figure S5: Performance of individual protein with age models on unseen data. For each protein, the protein term plus age term and protein:age interaction term were modelled against linearized 12-month posttransplant eGFR (eGFR12) using LASSO. For proteins passing filtering, with eGFR root-mean-square error (RMSE) better than an age-only model, we compared the RMSE achieved on the model training data against the RMSE achieved on held-out test data. We observed that the models actually achieved slightly better performance on the test data, reassuring us that no model was over-fit to the training data (i.e. all models generalized well to unseen data).
Supplementary Figure S6: Performance of VTN, APOE and CST3 models on paired kidney data. Six of the analyzed kidneys were paired (i.e. two came from each of three donors). These data were not used for training or testing but reserved for QC of the proteomics (Figure S2). Although our analysis was focused on identification of outcome-associated protein signals rather than prediction of eGFR12 (due to cohort size limitations, we do not claim to have a maximally biologically representative training set), we nevertheless used these paired kidneys to assess the consistency of our modelling in predicting above/below median outcomes in the three proteins we highlighted as particularly interesting. In all three pair cases the outcomes in (different) recipients were similar within each pair, and for all three proteins the model consistently predicted the correct above/below median classification for all six kidneys.
DATA AND MATERIALS AVAILABILITY
The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE 68 partner repository with the dataset identifier PXD033428.
DISCLOSURE AND FUNDING
This study was supported by NHS Blood and Transplant funding awarded to MK & RJP. SF was supported by Kidney Research UK, grant reference KS_RP_002_20210111 awarded to MK. PDC was supported by a Chinese Academy of Medical Sciences 2018-I2M-2-002 awarded to BMK. Authors declare that they have no competing interests.
Author contributions
Conceptualization: MK
Methodology: PDC, SF, RV, SD, RF, BMK, AS, ES, RJP, MK
Investigation: PDC, SF, RV, PJ, SD, IV, KT, AS
Visualization: PDC
Funding acquisition: BMK, RJP, MK
Project oversight: MK
Supervision: RF, AS, MK
Writing – original draft: PDC, ES, MK
Writing – review & editing: All authors
Data Availability
The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE 67 partner repository with the dataset identifier PXD033428. All other data produced are contained in the manuscript or supplemental materials.
ACKNOWLEDGMENTS
We thank the UK QUOD Consortium and NHS Blood and Transplant UK Registry for providing the samples and the associated clinical and demographic metadata.; in particular we thank Sheba Ziyenge, Lewis Simmonds and Dr Sarah Cross, Dr Sergei Maslau and Mr Tomas Surik for their support on the QUOD sample selection.
We thank members of the Discovery Proteomics Facility within the TDI Mass Spectrometry Laboratory for expert help with mass spectrometry analysis, and members of the Lindgren group at the BDI for informative discussions regarding statistical modelling.
Footnotes
Manuscript revised to include a revised analysis workflow reporting an expanded list of identified proteins and orthogonal contextualisation in a recently published independent dataset. Methods, results and discussion expanded to include additional analysis and findings. Figures and supplemental figures and tables revised.
REFERENCES
