Gene expression changes are associated with loss of graft function and Interstitial Fibrosis and tubular atrophy: Diagnosis versus Prediction

Patient demographics and clinical variables

Average recipient age was 44 ± 16.4 years for the NA group and 45 ± 9.9 years for the IF/TA group. Donor age was 37 ± 14.3 years and 33 ± 12.4 years for the NA and IF/TA groups respectively. No statistically significant changes were identified between the two groups in regard to cold ischemia time, warm ischemia time, incidence of delayed graft function (DGF), and acute rejection (AR). A summary of the patient demographic and clinical data is found in Table 1A.

Gene expression profiling of IF/TA and NA samples

Differential gene expression analyses were carried out using samples (needle biopsy) collected from kidney grafts in order to identify gene expression changes occurring in IF/TA allografts compared to NA that could be used as early detection or diagnostic biomarkers of IF/TA. We have identified 2,223 probesets (1,729 unique genes). A volcano plot and a histogram of the p-values from the two-sample t-test are displayed in Figure 1. When applying partitioning around mediods to the 80 probe sets that remained after filtering, k=2 yielded the largest average silhouette width of 0.507 indicating a borderline reasonable structure. Hierarchical clustering was then performed; the resulting heatmap and dendograms are shown in Figure 2. For the two large clusters identified in the heatmap from the hierarchial clustering procedure, Cluster 1 included 15 NA samples and one IF/TA sample and Cluster 2 included 23 IF/TA and 2 NA samples. Six samples within Cluster 2 clustered apart from the other pathological cases of IF/TA (Cluster 2a). Of these, four were IF/TA samples, which in addition to signs of IF/TA also showed histological signs of subclinical acute rejection (N= 3) or BK virus (N= 1).

Interaction networks and functional analysis

Gene ontology analyses were performed to identify biological processes and pathways over-represented by the differentially expressed genes. Functional analyses indicate that up-regulated genes are primarily involved in pathways related to the immune response, T-cell activation, natural killer cell mediated cytotoxicity, and programmed cell death (p-value <0.05). Down-regulated genes were associated with various general metabolic processes such as carboxylic acid metabolism, fatty acid and amine metabolism, and oxidative phosphorylation (p-value <0.05). The top 15 up-regulated and 15 down-regulated biological processes are shown in Table 2. A complete list can be found in Supplemental Table 1.

Pathway analyses were also performed. The top 30 canonical pathways (p-value <0.05) identified are shown in Figure 3. A complete list can be found in the Supplemental Table 2. Interestingly, two of the pathways identified were allograft rejection and graft-versus-host signaling. Natural killer cell mediated cytotoxicity and antigen presentation were also identified. Genes in the allograft rejection pathway included: FAS, GZMB, GZMA, CD40, CD28, CD86, PRF1, IFNG and several HLA genes (i.e. HLA-A, HLA-B, HLA-C, HLA-DRA, HLA-DRB1, HLA-DQB1, HLA-E, HLA-F, and HLA-G). Additionally, our analysis also identified a number of genes associated with natural killer cells including FAS, MICB, BID, CASP3, ITGB2, ITGAL, LCK, and LAT. HLA genes were also identified as part of the antigen processing and presentation pathway along with CD8B, CD8A, CD4, CD74, KLRD1, KLRC3, and KLRC1. As with the biological processes, down-regulated pathways suggest a general decline in metabolism with oxidative phosphorylation, valine, leucine and isoleucine degradation, and TCA.

LASSO model and N-fold cross-validation

A least absolute shrinkage and selection operator (LASSO) model was derived from the differentially expressed probesets, resulting in a final model of 10 probesets with 100% training set accuracy. Because training set accuracy is always an optimistic estimate, to obtain an unbiased estimate of classification error, N-fold cross-validation was used. Using N-fold cross-validation, the 41 different LASSO models included an average of 8 probesets (range 2 to 12). A total of 36 different probesets were included in the models derived using N-fold cross-validation (Supplemental Table 3). The N-fold cross-validated error was 2.4% with only 1 IF/TA sample being misclassified (sensitivity 95.8%, specificity 100%). However, when the final LASSO model was applied to gene expression data from an independently collected set of protocol biopsies (34 biopsies taken at 3-months post-transplantation) to test the predictive ability of the detected gene expression changes, all samples were classified as NA.

Sample profiling at 3 months post-kidney transplantation

A summary of the patient demographic and clinical data is found in Table 1B. Microarray data for the thirty-four protocol biopsies taken at 3-months were segregated into 2 groups based on eGFR at 9-months post-transplantation: Group 1 had 24 samples with eGFR >45 mL/min and normal histology, and Group 2 included 10 samples with eGFR ≤45 mL/min and IF/TA (IF (ci≥1) and TA (ct≥1) involves more than 25% of the cortical area). No probe sets were significant using a Bonferroni correction, though eleven probesets were differentially expressed between the two groups at the P<0.001 level. Moreover, 143 probe sets were significant when a P<0.005 was used. HGF and FGF signaling and THBS1 expression were up regulated in the Group 2 samples. IL10 signaling was up regulated in patients from Group 1. A LASSO model was also derived and N-fold cross-validation conducted. The final LASSO model included 17 probesets and had 100% training set accuracy (Table 2B). A total of 93 different probesets were included in the 33 different models derived using the N-fold cross-validation procedure. The 93 probe sets are listed in the Supplemental Table 4.

Validation of final LASSO model generated from 3 month transplantation biopsies to predict progression to IF/TA in an independent set of samples collected after 6 month-post-transplantation

When the LASSO model derived using the 3 month-post-transplantation biopsies was applied to the independent set of biopsies collected after 6-months-post transplantation (mean=9.25, range 6–12). These samples were classified as NA (N=17) or allograft with IF/TA (N=7) using histological findings and graft function at the biopsy time. The classifier’s accuracy was 58.3%. The positive predictive value was 71%.

QPCR validation

Microarray results were validated for FOS, CCL5, GZMA, CD86, IL10RA, ITGB2, CXCL9, ABAT, and ACAT1. All genes tested showed a similar trend to the microarray data with a Pearson’s correlation of 0.881 (Supplemental Figure 3).

Kidney samples

One-hundred-one allograft kidney samples from adult deceased donor kidney recipients (KTRs) were included. Written informed consent was obtained from all patients. The Institutional Review Board (IRB) at VCU approved the study protocol. No living donors, HIV positive or re-transplantation patients were included. Renal allograft tissue was obtained through an 18-gauge biopsy needle. Samples were placed in RNAlater (Ambion) immediately after collection. The study design is resumed in Figure 1. Samples collected at different time points post-transplantation were studied. First, the analysis included normal allografts (NA, n=17) and IF/TA (n=26) samples (classified based on histopathology using the Banff criteria (34,35)). NA samples were defined as allograft biopsies from unique recipients collected at 9-months or later post-transplantation with normal histology and the estimated glomerular filtration rate (eGFR) values at the time of collection were >45 mL/min. Estimated GFR was calculated using the Modification of Diet in Renal Disease formula (35). Second, samples collected at 3-months post-implantation (n=34) were analyzed. This last set of samples was also divided into 2 groups based on histological findings and eGFR at 9 months post-KT. Group 1 (n=24 with normal histology and eGFR >45 mL/min) and Group 2 (n=10, TA (ct≥1) and IF (ci≥1) involving more than 25% of the cortical area, as previously defined (36,37)) and eGFR ≤45 mL/min. Finally, an independent set of samples collected after 6-month post-transplantation (n=24) was used for validation.

Microarray data statistical analysis

Gene expression microarray data for forty-one samples were available (24 IF/TA and 17 NA). Probe level data were read into the R programming environment using the affy Bioconductor package and the robust multiarray average method was used to obtain probeset expression summaries. Quality assessment of each GeneChip was performed as it was previously published (38, 39). Prior to performing analyses, all control probe sets and probe sets declared absent for all samples were removed, leaving 17,641 probe sets for statistical analyses.

A two-sample t-test was performed and probe sets were considered significant using α=0.05 with a Bonferroni correction. Additionally, P-values from the t-test were subsequently used in obtaining the false discovery rate (FDR) using the q-value method. To identify whether there were any structures in the gene expression data, the data were filtered by retaining probe sets called present or marginally present for 80% or more of the samples and having a gene expression standard deviation among the 99^th percentile or higher. After filtering, partitioning around medoids was applied allowing the number of clusters to range from two to six and the average silhouette width was estimated as a measure of the strength of clustering. Hierarchical clustering using Ward’s method was also applied using 1- |Pearson’s correlation| as the dissimilarity measure.

Probe sets declared absent by the detection call algorithm for all 34 arrays (3-month biopsies) were removed leaving 16,300 probe sets for statistical analysis. Samples were segregated into 2 groups based on eGFR histological findings at 9-months post-transplantation and: Group 1 (G1) had 24 samples with eGFR >45 mL/min and normal histology, and Group 2 (G2) included 10 samples with eGFR ≤45 mL/min and IF/TA (IF (ci≥1) and TA (ct≥1) involves more than 25% of the cortical area). To identify probe sets exhibiting differential expression between Group 1 and Group 2, a two-sample t-test was performed and probe sets were considered significant using α=0.05 with a Bonferroni correction.

The 24 samples collected after 6 months post-transplantation were also classified as with IF/TA vs. no-IF/TA based on histological findings at the time of collection and graft function.

Prediction modeling

The least absolute shrinkage and selection operator (LASSO) model, which maximizes the likelihood subject to the constraint that the sum of the absolute values of the regression coefficients is less than some tuning parameter t, was fit where the dependent variable predicted was IF/TA vs. NA or Group1 vs. Group2 (40). The final model was selected as the smallest multi-degree of freedom model that minimized the training error. The predicted class was IF/TA or Group 2 if the fitted probability was ≥0.50 and normal allograft or Group 1 otherwise. To obtain an unbiased estimate of classification error, N-fold cross-validation was used (41,42). To identify the predictive ability of the identified GE signature, the final LASSO model from the IF/TA vs. NA analysis was also applied to samples collected at 3-months post-transplantation. The LASSO models were fit using the glmpath package (42). All analyses were conducted in the R programming environment using appropriate Bioconductor packages (43).

Interaction networks and functional analysis

Gene ontology and gene interaction analyses were executed using ToppGene (http://toppgene.cchmc.org/) (44). Gene lists containing Entrez GeneID numbers were used as input.

(Supplemental information about the used methods in: Supplemental methods)