A common rejection module (CRM) for acute rejection across multiple organs identifies novel therapeutics for organ transplantation

We downloaded raw data for eight gene expression studies from organ biopsy specimens from kidney, lung, heart, and liver transplant patients, with and without diagnosis of AR (Table S1 A). To reduce the clinical complexity in defining AR and stable (STA) phenotypes, we used the phenotypes as defined in the corresponding original publications. Phenotype definitions and sample composition for each dataset are described in Materials and methods. Notably, none of the studies had any antibody-mediated rejection samples or did not report this information. We filtered each dataset to include only biopsies from patients with AR and patients who were in STA condition. This filtering reduced the number of samples available for analysis to 236 from 392 samples in eight datasets. After reannotating the probes, each dataset was separately normalized using gcRMA (Irizarry et al., 2003a). Datasets GSE2596 and GSE4470 were not normalized because raw data were not available, and the downloaded data were already normalized.

We applied two meta-analysis approaches to the normalized data as described in Materials and methods. In brief, the first approach combines effect sizes from each dataset into a meta-effect size to estimate the amount of change in expression across all datasets. For each gene in each dataset, an effect size was computed using Hedges’ adjusted g. If multiple probes mapped to a gene, the effect size for each gene was summarized using the fixed effect inverse-variance model. We combined study-specific effect sizes to obtain the pooled effect size and its standard error using the random effects inverse-variance technique. We computed z-statistics as a ratio of the pooled effect size to its standard error for each gene and compared the result with a standard normal distribution to obtain a nominal p-value. P-values were corrected for multiple hypotheses testing using false discovery rate (FDR; Storey, 2002). We identified 180 genes that were measured in all datasets and were overexpressed in AR with P < 0.01 (FDR ≤ 20%; Table S1 B).

We used a second nonparametric meta-analysis that combines p-values from individual experiments to identify genes with a large effect size in all datasets. In brief, we calculated a t-statistic for each gene in each study. After computing one-tail p-values for each gene, they were corrected for multiple hypotheses using FDR. Next, we used Fisher’s sum of logs method (Fisher, 1932), which sums the logarithm of corrected p-values across all datasets for each gene and compares the sum against a χ² distribution with 2k degrees of freedom, where k is the number of datasets used in the analysis. This method identified 1,772 overexpressed genes at FDR < 20% (Table S1 C).

102 genes were identified as significantly overexpressed by both methods (Table S1 D). This group contains genes that (a) had overall large effect size across all datasets and (b) were significant across all datasets. Although our selection criteria may have left out genes with varying expression in AR, the method allowed us to develop robust overlapping transcriptional signals in AR across all transplanted organs.

Using BioGPS, we defined a gene as preferentially expressed in a tissue if its expression in a given tissue was at least three times its median expression across all 84 tissues in BioGPS (Su et al., 2004). We found that the 102 genes are highly expressed in one or more blood cell type that participates in the immune response (Table S1 E), suggesting that our meta-analysis removed tissue-specific bias and identified the relevant pathogenic transcriptional signal of activated infiltrating cells in the graft in AR, rather than being affected by various confounding factors such as organ-specific expression bias, treatment protocols by different groups, or different microarray platforms.

Network analysis of the 102 genes using Ingenuity Pathway Analysis (IPA) revealed that 96 of our genes are part of a network involved in cellular movement and immune cell trafficking (Fig. 1). These genes include major histocompatibility complex class I and II molecules, interferon regulatory factors, granzymes, chemokines, interleukins, transcription factors, and the T cell receptors. All have direct relationships to one another, as defined by the IPA, and each relationship has been experimentally verified in the literature. Canonical pathway analysis of these overexpressed genes using IPA and Pathway-Express (Draghici et al., 2007; Khatri et al., 2007) confirmed that they are in many of the pathways known to be related to regulation of the immune response and our current understanding of graft rejection (not depicted).

In a meta-analysis, one experiment with a large number of samples can overwhelm the results. Also, because of an unequal number of datasets for each transplanted organ in our analysis, it is possible that we may have introduced organ-specific bias. Therefore, to avoid (a) the influence of a single large experiment on the meta-analysis results and (b) organ-specific bias caused by unequal number of datasets (and samples) used in the meta-analysis, we performed leave-one-organ-out meta-analysis. We excluded all datasets from individual organs, one organ at a time, and performed meta-analysis on the remaining datasets from three organs. We hypothesized that the minimal set of genes that are significantly overexpressed, irrespective of the set of organs analyzed, would constitute a CRM in solid organ transplant rejection.

We identified 12 overexpressed genes: BASP1, CD6, CD7, CXCL10, CXCL9, INPP5D, ISG20, LCK, NKG7, PSMB9, RUNX3, and TAP1 at FDR ≤ 20% (Fig. 2, A–L). Network analysis using MetaCore showed that 10 out of the 12 genes are connected to each other, where STAT1 and NF-κB form the central axis of regulation (Fig. 2 M). Regulation of expression by STAT1 and NK-κB for many genes in this network has been verified experimentally in the literature (Chatterjee-Kishore et al., 1998; Der et al., 1998; Sharif et al., 2004; Shi et al., 2005; Robertson et al., 2007; Kuznetsov, 2009; Ellis et al., 2010).

We further validated overexpression of the 12 genes set in three independent cohorts consisting of 794 renal allograft biopsies (Table S2 A). Two of these datasets, GSE21374 (282 samples, AR = 76, STA = 206; Einecke et al., 2010) and GSE36059 (411 samples, AR = 122, STA = 289; Reeve et al., 2013) have previously been published. GSE36059 also identified whether the rejection was T cell mediated (35 samples), antibody mediated (65 samples), or mixed (22 samples). We generated the third dataset of 101 renal transplant biopsies (AR = 43, STA = 58) in our laboratory, referred to as the Stanford cohort in this manuscript (GEO accession no. GSE50058). All datasets contained biopsy-proven AR and STA renal allograft samples. Meta-analysis showed that all genes except CD7 were overexpressed in these cohorts (P < 0.001, FDR < 1%; Table S2 B) in renal transplant biopsies during AR (Fig. 3). Consequently, we excluded CD7 from the list of 12 genes. We defined the set of remaining 11 genes as a CRM that is an important transcriptional axis in AR of transplanted solid organs.

We defined geometric mean of expression of the CRM genes in each sample as a CRM score. In each independent dataset, the CRM score was significantly higher in the AR group (P < 1.5 × 10⁻⁷; Fig. 4, A, C, and E). Each unit increment in the CRM score increased the odds of AR by 4.17, 5.45, and 3.63 in GSE21374, GSE36059, and Stanford cohort, respectively. It was also able to distinguish AR and STA samples with high specificity and sensitivity in GSE21374 area under the curve (AUC = 0.83; Fig. 4 B), GSE36059 (AUC = 0.8; Fig. 4 D), and Stanford cohort (AUC = 0.82; Fig. 4 F).

Because the degree of progressive chronic histological damage is associated with long-term graft survival (Naesens et al., 2011), we further correlated the CRM scores with the extent of graft injury. In GSE1563, which also included biopsies from healthy donors and patients with renal dysfunction without graft rejection, the CRM scores were the lowest for healthy donor kidney biopsies with very low variation (mean = 3.31, SD = 0.12), slightly higher for the STA samples (mean = 4.02, SD = 0.76), and the highest for AR (mean = 5.98, SD = 0.85; Fig. 5 A). Furthermore, the CRM score correlated significantly with the serum creatinine level in these samples (Pearson correlation = 0.43, P = 0.01; Fig. 5 A). However, the CRM score identified two groups of STA patients although serum creatinine does not show any difference in the STA group. The samples in the first group have their CRM scores indistinguishable from that of the healthy donor kidney (HC) samples, which is significantly lower than the CRM scores for the second group of STA patients (P = 0.0014). It is possible that the STA samples with higher CRM scores than healthy controls but lower than AR identify patients that may be at risk for increased graft injury although there is no significant difference in their serum creatinine levels.

Based on this observation, we hypothesized that the CRM score could allow better identification of STA transplant patients who may be experiencing subclinical graft injury and should be monitored closely. Therefore, we investigated the relationship between the CRM score and the chronic histological damage in renal allografts as defined by the chronic allograft damage index (CADI; Yilmaz et al., 2003) using publicly available renal transplant gene expression dataset GSE25902. We have previously described this cohort in an earlier publication (Naesens et al., 2011). In brief, the dataset consists of 120 renal biopsies from 72 patients, which include 24 AR biopsies from 24 patients, and 96 protocol biopsies from 48 patients with varying degrees of chronic allograft injury graded by the CADI. Out of the 96 protocol biopsies with chronic allograft injury, 72 biopsies were obtained from 24 patients at three time points: at the time of implantation and at 6 and 24 mo after transplantation. This cohort of 72 biopsies is referred to as the “longitudinal cohort” in the rest of the manuscript. None of the 24 patients in the longitudinal cohort experienced Banff-grade acute T cell–mediated rejection within the first 2 yr after transplantation. Furthermore, to avoid the confounding influence of donor quality and extended ischemia time on gene expression changes (Naesens et al., 2009), the sample set was carefully selected such that kidney graft quality was excellent at implantation, with minimal chronic damage. Despite the pristine tissue quality at engraftment, the effect of increasing posttransplant time, even in the absence of interval AR, results in a significant increase in chronic histological damage in all histological compartments. At 24 mo after transplantation, 50% of the group (12 patients) had significantly greater histological progression, scored by the incremental CADI score. We defined these samples as a “progressor” group, and the remaining samples as a “nonprogressor” group. At 6 mo, the progressor group had higher mean CADI scores than the nonprogressor (3.2 ± 1.8 vs. 0.75 ± 0.75; P = 0.0003), and similar trends were seen for mean CADI scores for the 24-mo biopsies (7.3 ± 1.6 in progressor in contrast to 1.8 ± 1.4 in nonprogressor; P ≤ 0.0001). No other Banff qualifiers were significantly different at any time point (Naesens et al., 2011).

The CRM scores for the three sample groups (nonprogressor, progressor, and AR) in the 120 samples were significantly different (Fig. 5 B) and positively correlated with the extent of injury (Jonckheere-Terpstra [JT] trend test p-value = 7.99 × 10⁻¹⁵). We further correlated the CRM scores with different Banff qualifiers for AR. The CRM score showed significant positive correlation with Banff t-score (JT p-value = 6.04 × 10⁻¹²; Fig. 5 C) and i-score (JT p-value = 2.72 × 10⁻¹²; Fig. 5 D). The CRM score also showed significant positive correlation with Banff ct-score (JT p-value = 1.995 × 10⁻⁵; Fig. 5 E) and ci-score (JT p-value = 6.195 × 10⁻⁷; Fig. 5 F) in the longitudinal cohort.

Finally, we performed a repeated measures analysis of variance for the longitudinal cohort (Fig. 5 G). Our results showed that the CRM scores were significantly different between the progressor and nonprogressor groups (P = 1.86 × 10⁻⁵), and the CRM scores increased significantly over time for both groups (P = 4.03 × 10⁻⁶). Importantly, as time after transplant increased, the CRM scores for the progressor groups increased at a faster rate than the nonprogressor group (P = 0.0276). Therefore, we used the CRM scores for the 6-mo protocol biopsies to test whether it can predict future histological damage to the graft, as defined by CADI score at 24 mo after transplantation. We found that the CRM scores in 6-mo biopsies can predict future histological damage with high specificity and sensitivity (AUC = 0.882; Fig. 5 H). In summary, our results using two additional independent cohorts of 151 renal transplant biopsies show that the CRM score significantly correlates with graft injury and can be used for identifying patients at risk of developing significant histological damage in future using protocol biopsies.

Because of ubiquitous and correlated overexpression of the CRM during AR and chronic injury in 13 independent datasets consisting of 1,164 biopsies from four types of transplanted organs, we hypothesized that irrespective of the transplanted organ, the CRM is a critical axis of gene regulation during AR and graft injury. Therefore, pharmaceutical disruption of the CRM may reduce immune cell infiltration, reduce graft injury, and increase graft survival.

A literature review found that 6 out of the 11 genes are direct or indirect targets of FDA-approved drugs. Bortezomib is an FDA-approved drug that inhibits PSMB9. It can reverse antibody-mediated rejection and eliminate donor-specific anti–human leukocyte antigen antibodies (Walsh et al., 2010). Mycophenolate mofetil, which is also FDA approved and primarily targets IMP dehydrogenase 2 (IMPDH2), reduces expression of INPP5D by more than twofold (van Leuven et al., 2010). It is a potent immunosuppressive drug that reduces the risk of AR (Knight et al., 2009) and has a possible beneficial effect on chronic graft survival (Ojo et al., 2000). BASP1 and CXCL9 are selectively targeted by doxycycline (Hartl et al., 2009) and sulindac (Sakaeda et al., 2006), respectively. Multiple studies using in vitro and in vivo models have shown that dasatinib (BMS-354825, Sprycel; Bristol-Meyers Squibb) binds to and potently inhibits LCK. Dasatinib is an ATP competitor approved for Imatinib-resistant chronic myeloid leukemia. Lee et al. (2010) have shown that dasatinib binds to the ATP-binding site of LCK and freezes the kinase domain in an inactive conformation. Schade et al. (2008) have shown that inhibition of LCK by dasatinib, which is critical for T cell receptor signaling, is responsible for reduction in T cell activation in vitro and in vivo. In addition, Blake et al. (2008) have shown that dasatinib also reduces NK cell cytotoxicity. Atorvastatin (Lipitor) is an HMG-CoA reductase inhibitor that slows the production of cholesterol and is used for treatment of hyperlipidemia. Out of nine chemokines and four endothelial cytokines investigated in plasma samples from patients with Crohn’s disease, atorvastatin reduced CXCL10 plasma levels but did not affect the other chemokines and cytokines (Grip and Janciauskiene, 2009). Atorvastatin treatment has been shown to also significantly reduce plasma levels of CXCL9 in systemic lupus erythematosus patients (Ferreira et al., 2010). Both CXCL9 and CXCL10 are ligands of CXCR3, which is a chemokine receptor most consistently expressed during rejection that promotes localization of CTLs to inflamed tissues (Wang et al., 2008). Furthermore, statin has been shown to demethylate the FOXP3 promoter (Kim et al., 2010). The transcription factor FOXP3 is critical for development and function of regulatory T cells (T reg cells). FOXP3⁺ T reg cells are a unique subset of CD4⁺ T cells that mediate immunosuppression. Baron et al. (2007) have shown that DNA demethylation in the human FOXP3 locus discriminates T reg cells from activated FOXP3⁺ conventional T cells.

Both dasatinib and atorvastatin are FDA-approved drugs for nontransplant conditions. Because atorvastatin and dasatinib have been FDA approved for treating nontransplant conditions, we hypothesized that they could be repositioned in organ transplantation. To test this hypothesis, we used both drugs in an established mouse FVB-to-C57BL/6 heterotopic cardiac transplant model (Corry and Russell, 1973) to investigate the effect of peritransplant drug administration on mitigating cell infiltration in the transplanted graft. We chose a cardiac transplant model in mouse to also illustrate that the CRM is indeed common during AR in multiple organs. We also used cyclosporine as a positive control in this model. After each drug treatment for 7 d, the grafts were evaluated by comparing against untreated AR using standard graft histology, a count of the infiltrating cell subsets in the graft, and by transcriptional analysis of the graft’s Q-PCR.

Gene expression profiling of nontransplanted hearts (FVB mice) and untreated, transplanted hearts showed that a majority of the 102 cross-organ rejection genes were significantly overexpressed in untreated AR (FDR ≤ 2%; Table S3 A), including all of the CRM genes (FDR ≤ 0.1%). Pathway analysis of down-regulated genes in each treatment group against each untreated AR group using IPA showed that only cyclosporine affected the T cell–related pathways. Atorvastatin affected monocyte- and macrophage-related pathways, and dasatinib affected cell cycle–related pathways (Table S3 B). Using Q-PCR, we found that the majority of the CRM genes were down-regulated in each of the treatment groups (Table S3 C). Not all CRM genes were down-regulated by any of the drugs used, which is expected because each drug acts through different mechanisms.

Furthermore, immunohistochemistry showed that there were fewer infiltrating cells in the atorvastatin and dasatinib treatment groups compared with untreated AR (P < 0.005) and were equivalent to treatment with cyclosporine (P > 0.05, i.e., statistically not significant; Fig. 6, A–F; and Table S3 D). For each treatment group, we measured the number of total infiltrating CD45⁺ cells, CD4⁺ T cells, CD8⁺ T cells, B220⁺ B cells, CD11c⁺ dendritic cells, F4/80⁺ macrophages, Gr1⁺ neutrophils, and NK1.1⁺ NK cells (Fig. 6, F–M). Although all three drugs reduced CD4⁺ and CD8⁺ T cells compared with untreated AR, cyclosporine reduced the number of infiltrating CD8⁺ T cells significantly more than atorvastatin and dasatinib (Fig. 6, G and H; and Table S3 D). However, atorvastatin and dasatinib significantly reduced the number of infiltrating B220⁺ B cells compared with cyclosporine (Fig. 6 J). Atorvastatin and dasatinib also reduced the number of infiltrating macrophages, dendritic cells, and NK cells, whereas cyclosporine’s effect was not statistically significant (Fig. 6, K–M). To test whether reduction in graft-infiltrating cells by atorvastatin and dasatinib results in extruded graft survival, we repeated the experiment, in which the mice were treated for up to 30 d. Using Cox proportional hazard analysis, we found that when treated with atorvastatin or dasatinib compared with untreated AR, the hazard ratio for graft survival was 36.33 (P = 0.002) and 66.26 (P = 0.0007), respectively (Fig. 7 A). Median survival for the untreated AR group was 10 d, but was 17 d for atorvastatin and 24.5 d for dasatinib.

Next, we tested whether combination of cyclosporine and atorvastatin can reduce cyclosporine dose and, hence, help avoid drug toxicity. We treated a group of six mice with atorvastatin and 50% reduced dose of cyclosporine (from 20 mg/kg/day to 10 mg/kg/day). The survival curve for this group, followed for up to 30 d, was similar to that of the mice treated with only atorvastatin (not depicted). Using Cox proportional hazard analysis, when treated with atorvastatin and cyclosporine compared with untreated AR, the hazard ratio for graft survival was 11.38 (P = 0.003). These results suggest that atorvastatin and cyclosporine may work through different pathways, each of which results in extended graft survival compared with the untreated group.

To validate the suggested benefits of statin use in a large clinical transplant population, we used electronic medical records from all 2,515 patients that received renal transplant between January 1989 and March 2012 at the University Hospitals Leuven. Out of the 2,515 patients, 1,566 received statin within the first 180 d after transplantation, with graft surviving at least 180 d. None of the patients started or stopped statin because of renal function evolution or intragraft phenomena, suggesting there is no direct bias in the association between statin use and graft outcome. In Cox proportional hazards analysis, after censoring for when a patient stopped taking statin, graft failed, or recipient death, statin use was associated with improved graft survival (HR = 0.701, P = 0.01; Fig. 7 B). This effect was statistically significant after adjusting for donor and recipient age, repeat transplantation, and calendar year (Table S4).

The goal of this study was to ascertain whether there are common immunological and inflammatory alterations during AR in an allotransplant, regardless of tissue source. Identifying a common response could be useful for monitoring rejection and could suggest novel targets for rational immunosuppressive drug discovery and design.

Despite several transcriptional studies in organ transplantation, their limited sample sizes make it difficult to capture the molecular heterogeneity of AR (Chatterjee-Kishore et al., 1998; Der et al., 1998; Sarwal et al., 2003; Sharif et al., 2004; Shi et al., 2005; Robertson et al., 2007; Kuznetsov, 2009; Ellis et al., 2010). Furthermore, experimental confounders such as platform variability and organ-specific profiles instead of a rejection-specific profile can easily overwhelm expression measurements in individual experiments. We performed a meta-analysis of 236 biopsy samples from publicly available transplant tissue microarray datasets using a novel method. We propose that using publicly available data from multiple laboratories implicitly accounts for the underlying molecular heterogeneity of injury, the variability of the host response, differences in treatment protocols, and other confounding factors. By integrating multiple datasets from different transplanted organs in different hospitals, we were able to increase sample size, avoid organ-specific bias, and identify a CRM. We identified 102 immune-related genes that were overexpressed during AR, which reflected signal from activated cellular infiltrates as expected. More importantly, we identified a CRM comprising 11 genes that are overexpressed irrespective of the transplanted organ. These genes are variably expressed in B cells (Drayton et al., 2006), NK cells (Kondo et al., 2000; Obara et al., 2005), and T cells (Zhao et al., 2002; Yilmaz et al., 2003). Using five independent cohorts, consisting of 928 renal transplant biopsy samples, we validated overexpression of the 11 CRM genes during AR and showed that the CRM score for a sample, defined as geometric mean of expression of the 11 genes, correlated with extent of graft injury. Furthermore, using 6-mo protocol biopsy, the CRM score can predict with high specificity and sensitivity patients who are likely to develop severe histological damage within 2 yr after transplantation.

It is important to note that many of the 102 genes have been previously reported in the literature, including MHC class I and II molecules, chemokines, markers of T cell activation, and co-stimulatory molecules (Fig. 1). For instance, chemokine receptor CXCR3 and its ligands, CXCL9, CXCL10, and CXCL11, are most consistently overexpressed during rejection (Karason et al., 2006; Wang et al., 2008; Chen et al., 2010). However, the CRM genes contain only two of the CXCR3 ligands, CXCL9 and CXCL10. CXCR3 was not identified in our analysis because it was down-regulated in the liver dataset GSE13440 (not depicted), whereas CXCL11 was not identified because it was not measured in the liver dataset GSE13440 (not depicted). However, all samples in GSE13440 were HCV positive (Table S1 A). Hence, down-regulation of CXCR3 in the liver dataset does not necessarily indicate down-regulation of CXCR3 during AR in liver transplant because overexpression of CXCR3 during HCV infection could explain down-regulation or no differential expression of CXCR3 in GSE13440.

In contrast, there are several possible explanations for why our list of genes does not include many genes that have been reported in the literature. For instance, we recently proposed a five-gene set model for diagnosis of AR renal transplantation (Li et al., 2012). Although these five genes were validated in an independent randomized cohort from 12 US transplant programs, there are important factors that distinguish it from our current study. First, our AJT study used peripheral blood samples from transplant patients, whereas the current study used only biopsy samples from transplant patients. The set of genes that are differentially expressed during AR in blood and allograft can be very different. Second, the AJT study was performed only on pediatric samples, whereas all datasets in the current study used adult samples. Third, it has been hypothesized previously that the triggers, which lead to allograft rejection, may be tissue specific (Wang et al., 2008). Experiments that use samples from the same type of organ are more likely to identify these tissue-specific triggers than the effector immune response that converge into a common mechanism. Our meta-analysis approach that uses samples from multiple organs in multiple independent experiments favors the common mechanism over the tissue-specific triggers.

The CRM also provides a mechanism for understanding varying biological outcomes after transplantation and AR therapy regardless of the tissue source of the organ. Different immunosuppressive therapies are likely to have different influences on the allogeneic response against distinct cell lineages. Most currently used immunosuppressive drugs prevent AR by inhibiting T cell activation. This is achieved through various means, such as inhibiting antigen-presenting cell development, cytokine production, or co-stimulatory signals for T cell activation (Halloran, 2004; Lechler et al., 2005; Walsh et al., 2010). This effect was also observed in our mouse model. Cyclosporine reduced infiltrating CD4⁺ and CD8⁺ T cells along with Gr1⁺ neutrophils but did not significantly reduce B220⁺ B cells and antigen-presenting cells including F480⁺ macrophages, CD11c⁺ dendritic cells, and NK1.1⁺ NK cells. In contrast, atorvastatin has been shown to suppress interferon-induced neopterin formation in monocytic cell lines (Neurauter et al., 2003). Neopterin is a metabolite of guanosine triphosphate, which is produced by monocyte-derived macrophages. Atorvastatin reduced F480⁺ macrophages the most in the mouse model (Fig. 6 K).

CXCL9, CXCL10, and CXCL11 are ligands of CXCR3 that polarize CD4⁺ T cells toward the Th1 phenotype and promote localization of CTLs to allograft. Interaction between CTLs and B cells leads to the secretion of CXCR3 ligands by B cells, forming a positive feedback loop, which in turn amplifies the inflammatory signal further (Deola et al., 2008). By reducing CXCL9 and CXCL10 expression in serum, atorvastatin potentially inhibits this positive feedback loop, which can explain lower B220⁺ B cells and other innate immune cells infiltrating the graft in the atorvastatin-treated group.

Current immune suppression therapies in transplant have been associated with increased incidence of cancer after transplantation (Vajdic et al., 2006). Recently, atorvastatin was shown to reduce cancer-related mortality in the cancer patients (Nielsen et al., 2012), which in turn suggests that atorvastatin could also be useful for reducing the risk of cancer in transplant patients in addition to increased graft survival. Furthermore, dasatinib is an oral BCR-ABL inhibitor for the treatment of chronic myeloid leukemia (Shah et al., 2008). Because dasatinib is a cancer treatment drug, it may also reduce the risk of various cancers after transplantation. In addition, it has potent inhibitory effects on TNF, IFNγ, IL2, IL6, IL10, IL17, and IL22 and has been show to inhibit T cell activation and proliferation and reduce NK cell cytotoxicity (Blake et al., 2008), suggesting it would likely show efficacy as an immunosuppressive agent (Schade et al., 2008). Our data show that dasatinib could be an adjunctive immunosuppressive agent, given that it has potent inhibitory effects on B cells and components of innate immune response including macrophages, dendritic cells, and NK cells.

Our finding of atorvastatin’s immunomodulatory effect is not a complete surprise. Variants of statins have been shown to inhibit induction of MHC-II expression by IFN-γ and repress MHC-II–mediated T cell activation (Kwak et al., 2000). Indeed, an immunosuppressive role of statins has been suggested in solid organ transplant clinical trials (Navaneethan et al., 2009). Pravastatin and simvastatin have been shown to confer better survival than the statin-naive control group in cardiac transplant recipients (Kobashigawa et al., 1995; Mehra et al., 2002). It is not clear whether there are any differences in immunomodulation among different statins, as pravastatin has been shown to reduce AR incidence significantly in renal transplant patients (Katznelson et al., 1996), whereas simvastatin and fluvastatin have not demonstrated this effect (Holdaas et al., 2001; Kasiske et al., 2001). Atorvastatin was recently studied in renal transplant patients to show immunoregulatory effects on T cells (Guillén et al., 2010).

However, one of the significant limitations of these studies, with the exception of the ALERT study, is that they were underpowered (Navaneethan et al., 2009). Furthermore, almost all randomized clinical trials using statins enrolled patients immediately after transplantation and followed them for ∼3 mo, and none looked at graft survival rates. Although many of these trials have concluded that statins have no effect on AR incidence rates, it may be because the patients are heavily immunosuppressed immediately after transplantation (typically using three immunosuppressive drugs).

The next logical phase of this study would be to find a combination of these drugs and appropriate dosage of each drug that is maximally beneficial. Although our results show that reducing cyclosporine dose by 50% in combination with atorvastatin does not improve graft survival as much as cyclosporine alone, it does not provide a complete picture. We need to consider bilateral pharmacokinetic interaction between cyclosporine and various statins. Cyclosporine and many statins (atorvastatin, lovastatin, and simvastatin) are metabolized by cytochrome P4503A4 (CYP3A4; Asberg et al., 2001). Cyclosporine has a significant impact on plasma concentrations of all statins, which results in a several-fold higher exposure compared with non-cyclosporine–treated controls (from 2-fold increase for fluvastatin to 23-fold increase for pravastatin; Arnadottir et al., 1993; Regazzi et al., 1993; Goldberg and Roth, 1996; Olbricht et al., 1997; Mück et al., 1999). Several studies have also shown that many statins tend to induce increased cyclosporine levels from no change for fluvastatin to 115% increase by lovastatin (Kuo et al., 1989; Cheung et al., 1993; Campana et al., 1995; Xu et al., 1998; Renders et al., 2001). However, atorvastatin has been shown to significantly reduce systematic exposure of cyclosporine by ∼10% in renal transplant patients (Asberg et al., 2001). Hence, bilateral pharmacokinetic interactions of different statins with cyclosporine and other widely used immunosuppressive drugs such as tacrolimus must be explored in more detail to identify an appropriate and effective combination of a statin and cyclosporine to treat transplant patients. It is possible that retrospective analysis of electronic medical records of multiple cohorts, similar to our analysis in Fig. 7 B, could facilitate identification of drug combinations that should be further tested in preclinical models.

To the best of our knowledge, our retrospective analysis of statin usage in more than 2,500 patients is the largest study with patients followed for up to 22 yr. Our analysis in this large kidney transplant patient cohort showed that statin use significantly improved death-censored renal allograft survival. These data suggest that statins may be important adjunctive agents to consider for de novo treatment in all organ transplant recipients because of their propensity to reduce cardiovascular morbidity and their synergistic immunosuppressive potential. Moreover, this analysis provides a proof of concept that identification of relevant molecules and pathways can be tested rapidly in large retrospective clinical datasets, which should then foster subsequent randomized controlled clinical trials.

Most importantly, our success in identifying the CRM in solid organ transplant rejection and the ability to reposition the existing drugs based on their interaction with the CRM targets underscores the importance of exploiting the data in public repositories for downstream clinical applications. Our results demonstrate that by broadening the sources of samples and pooling evidence from multiple studies for AR, we are able to identify a common, biologically relevant, rejection-specific signature in the tissue that can be useful for diagnostic purposes and therapeutic repositioning. NCBI GEO contains data from more than a million samples in experiments related to diseases that contribute to more than a third of human disease-related mortality in the United States (Butte, 2008). However, in recent years, as the amount of data in public repositories has increased, the number of FDA-approved drugs every year has declined, and largely stagnated to ∼20 drugs per year, from a peak of 53 drugs in 1996, amid an increase in the cost of drug discovery from $138 million in 1975 to $1.3 billion in 2006 and more than 15 yr needed, on average, in developing a single drug (Dudley et al., 2010). Our method can provide a rapid, economical, and targeted approach toward repositioning drugs with known human safety, for many unrelated diseases, quite different from their originally desired application.

We downloaded eight transplant gene expression datasets from four solid organs from GEO (Table S1 A). Each dataset was manually curated to select only tissue biopsy samples from AR and STA patients. For each study, we used the sample phenotypes as defined by the corresponding original published study. For all heart transplant studies, the samples were graded according to ISHLT guidelines. GDS1684, GSE2596, and GSE4470 used samples with ISHLT grade 0 as STA samples, and grade 3A or higher as AR, suggesting these studies used samples with extreme phenotype. GSE9377 also used ISHLT grade 0 as STA samples, which did not have any infections. Furthermore, GSE9377 used additional stringent criteria for sample selection (Table S1 A). GDS724 had a heterogeneous mix of STA samples that consisted of living donor controls and protocol biopsies >1 yr after transplant with good renal function and normal histology. In the case of GSE9493 (renal transplant) and GSE13440 (liver transplant), two independent pathologists classified the biopsies. Interestingly, all samples in GSE13440 were Hepatitis C virus positive. None of the studies had any antibody-mediated rejection samples or did not report this information.

Oligonucleotide arrays were checked for quality to ensure that they were free of experimental artifacts. Microarrays from a cDNA-based platform (GSE2596 and GSE4470) were not checked, as raw data were not available from GEO. Each oligonucleotide dataset was normalized using gcRMA (Irizarry et al., 2003b). Microarray probes in each dataset were mapped to Entrez Gene identifiers (IDs) to facilitate meta-analysis. If a probe matched more than one gene, the expression data for the probe were expanded to add one record for each mapped gene (Ramasamy et al., 2008).

The eight solid organ transplant datasets were analyzed using two different meta-analysis methods: (1) combining effect size and (2) combining p-values. We estimated the effect size for each gene in each dataset as Hedges’ adjusted g, which accounts for small sample bias. If multiple probes mapped to a gene, the effect size for each gene was summarized using the fixed effect inverse-variance model.

The study-specific effect sizes for each gene were then combined into a single meta-effect size using a linear combination of study-specific effect sizes, f_i, where each study-specific effect size was weighted by inverse of the variance in the corresponding study (Eq. 1). After computing meta-effect size, significant genes were identified using z-statistic, and p-values were corrected for multiple hypotheses testing using Benjamini-Hochberg FDR correction.

We used Fisher’s sum of logs method (Fisher, 1934) for meta-analysis by combining p-values. For each gene, we summed the logarithm of the one-sided hypothesis testing p-values across k studies and compared the result to a χ² distribution with 2k degrees of freedom. This process allowed us to identify significant genes (Eq. 2).

We selected 102 genes that satisfied the following criteria: (a) meta-effect size > 0 (i.e., overexpressed genes), (b) FDR ≤ 20% across all datasets when combining effect size, (c) measured in all eight datasets, and (d) when combining p-values using Fisher’s test, p-value < 0.2 for a gene to be up-regulated.

To account for the unequal number of datasets for each organ as well as to find a set of genes that was overexpressed independently of the source organ, we performed meta-analysis by removing organ-specific datasets one organ at a time. Thus, in the first iteration, we removed the datasets from heart transplants (GDS1684, GSE2596, GSE4470, and GSE9377) and performed meta-analysis on the remaining datasets, which were from kidney, lung, and liver. In the second iteration, we removed the datasets from kidney and performed meta-analysis on the datasets from heart, lung, and liver. At each iteration, we performed meta-analysis by combining effect sizes and by combining p-values. Using FDR ≤ 20% as a threshold, at each iteration, we identified 12 genes that were overexpressed in all organs used in the given iteration.

We performed functional pathway analysis using Pathway-Express (Draghici et al., 2007; Khatri et al., 2007; Tarca et al., 2009). Meta-effect size was used as fold change in Pathway-Express to identify significant pathways. We used FDR ≤ 10% as a threshold for identifying significant pathways. We performed network analysis using IPA with an option to include only “direct relationship” to avoid spurious connections caused by “indirect relations.” Direct relationships in IPA result from publications citing experimental evidence for an interaction.

We used the JT trend test to correlate the CRM score with various Banff criteria using GSE25902. It tests for a monotone trend in terms of the class parameter. In our case, the class is defined as various values of Banff criteria (e.g., t-score, i-score, etc.). The test uses the number of times that a sample with higher Banff criterion value has higher CRM score for trend test (Flandre and O’Quigley, 2007).

RT-PCR for the CRM genes in mouse allografts was performed using a high-throughput RT-PCR instrument (BioMark; Fluidigm). Total RNA was extracted from flash-frozen apical graft portions (one third of the allograft) using TRIzol reagent (Invitrogen) according to standard protocols. cDNA was generated using Superscript II (Invitrogen) and was preamplified on an Eppendorf Thermocycler. Preamplified cDNA was mixed with TaqMan Universal PCR Master Mix (Applied Biosystems) and Sample Loading Reagent (Fluidigm) and pipetted into the sample inlets of a Dynamic Array 96.96 chip (Fluidigm). TaqMan gene expression assays (Applied Biosystems) for the 12 genes plus 18S as endogenous control gene were diluted with Assay Loading Reagent (1:2; Fluidigm) and pipetted into the assay inlets of the same Dynamic Array 96.96 chip. After distributing assays and samples into the reaction wells of the chip in the NanoFlex controller (Fluidigm), a total of 1,876 qRT-PCR reactions were performed in the BioMark RT-PCR system in a total of 40 cycles. Data were analyzed using BioMark RT-PCR Analysis Software Version 2.0. Gene expression in each sample was calculated relative to the expression in a universal RNA sample (human universal RNA; Agilent Technologies), using the ΔΔCt method. The IDs of the assays used in PCR are as follows: 18S (Hs99999901_s1), BASP1 (Mm0234432_s1), CD6 (Mm01208285_m1), CD7 (Mm00438111_m1), CXCL9 (Mm00434946_m1), CXCL10 (Mm00445235_m1), INPP5D (Mm00494987_m1), ISG20 (Mm_00469585_m1), LCK (Mm0080297_m1), NKG7 (Mm00452524_g1), PSMB9 (Mm_00479004_m1), RUNX3 (Mm00490666_m1), and TAP1 (Mm00443188_m1).

For each kidney allograft biopsy, a separate core was stored in RNAlater (Ambion) and stored at −20°C until RNA extraction. Total RNA was extracted from each biopsy using TRIzol reagent. RNA integrity was checked with an RNA 6000 Nano LabChip kit (Agilent Technologies) on a 2100 Bioanalyzer (Agilent Technologies). RNA was amplified to cDNA and biotin labeled using the Ovation Biotin System (NuGEN Technologies). cDNA fragments were hybridized onto GeneChip Human Genome U133 Plus 2.0 Arrays (Affymetrix) according to the manufacturer’s protocol. The microarrays were scanned using a GeneChip Scanner 3000 (Affymetrix). Raw CEL files and gcRMA-normalized expression data for these 101 samples are available from the NCBI GEO series GSE50058.

For whole mouse genome expression analysis, we used Agilent Whole Mouse Genome 4 × 44K 60mer oligonucleotide arrays (G2519F; Agilent Technologies). A total of 100 ng of total RNA was used in the Agilent LIRAK PLUS, two-color Low RNA input Linear Amplification method, according to the manufacturer’s instructions. In brief, total RNA was reverse transcribed into cDNA using a T7-promotor primer and MMLV RT. The cDNA was transcribed into cRNA, during which it was fluorescently labeled by incorporation of cyanine (Cy)5-CTP (exposed samples) or Cy3-CTP (negative control samples). cRNA was purified with an RNeasy mini kit, and cRNA yield and Cy incorporation efficiency (specific activity) into the cRNA was measured with a NanoDrop Spectrophotometer (NanoDrop Technologies). cRNAs with yields > 825 ng and specific activity of 8–20 pmol/µg were selected for further processing.

Equal amounts (825 ng) of the exposed and negative control samples were competitively hybridized onto Agilent Whole Mouse GE 4 × 44K microarrays in a hybridization oven at 65°C for 17 h. Slides were washed according to the manufacturer’s instructions and dipped in Stabilization and Drying Solution (Agilent Technologies) to protect them from environmental ozone. They were scanned on an Agilent scanner and further processed using Agilent Feature Extraction Software. Agilent universal mouse reference RNA (#740100; Agilent Technologies) was used as a reference sample.

C57BL/6J (H2^b) and FVB (H2^q) mice were purchased from the Jackson Laboratory. We used male mice aged 6–10 wk with an average body weight of 25 g. Animals were maintained in the animal care facility at Stanford University. All experiments were approved by the Stanford University Institutional Animal Care and Use Committee, and were performed in accordance with the Guide for the Care and Use of Laboratory Animals (1996).

FVB donor hearts were implanted into the abdomen of C57BL/6 WT mice representing a complete MHC mismatch, as described previously (Corry and Russell, 1973; Fischbein et al., 2000). Animals were divided into three treatment groups (atorvastatin, dasatinib, and cyclosporine) and one nontreated control group. There were six animals in each group. Animal activity, body weight, and graft viability (abdominal palpation) were assessed daily.

Commercially available atorvastatin (PZ0001; Sigma-Aldrich) and dasatinib (S1021; Selleck Chemicals LLC) were suspended in sterile PBS (AccuGENE; Lonza) at concentrations of 9 mg/ml and 13.5 mg/ml for low- and high-dose atorvastatin, respectively, and at 4.5 mg/ml for dasatinib. Drug suspensions were aliquoted and stored at 4°C (atorvastatin) or −20°C (dasatinib). We used a cyclosporine solution of 250 mg/ml in ethanol and polyoxyethylated castor oil USP grade for i.v. injection (Bedford Labs). Cyclosporine was ordered through Stanford Hospital Pharmacy and diluted to 1 mg/ml in sterile saline solution.

Atorvastatin (75 mg/kg body weight/day) and dasatinib (25 mg/kg body weight/day) were administered daily by oral gavage. Cyclosporine (20 mg/kg body weight/day) was administered daily intraperitoneally. Before oral gavages of atorvastatin and dasatinib, aliquots were mixed thoroughly. Unused formulations were discarded. Treatment started the day before transplantation and lasted until the day before sacrifice by exsanguination at postoperative day 7. Grafts were harvested and divided into three equal parts for downstream analyses.

For the survival study, we used the same protocol described above with another group of 6 mice for each of the three treatments (total: 18 mice) and 6 mice for untreated AR group. Graft viability (abdominal palpation) was assessed daily for these animals until postoperative day 30.

One third of the explanted allograft heart (postoperative day 7) was immediately fixed in 20% buffered formalin for x min, embedded in paraffin, and subsequently stained with hematoxylin and eosin for histological assessment of tissue. Standard protocols were used. Pictures of the graft tissue were taken at 10 magnification using an E600 light microscope (Nikon) and SPOT version 4.6 imaging software (SPOT Imaging Solutions).

FITC, PE, or allophycocyanin (APC)-conjugated mAbs specific for mouse CD4 (GK1.5), CD8a (53-6.7), F4/80 (BM8), B220 (RA3-6B2), Gr1 (RB6-8C5), CD11c (N418), NK1.1 (PK136), CD45 (30-F11), and their isotype controls were purchased from BD, eBioscience, or BioLegend. Immediately after graft explantation at day 7 after transplantation, one third of each cardiac allograft was homogenized in RPMI 1640 media with 2 mg/ml Collagenase D (Worthington Biochemical Corporation) and 10% FCS for 2 h at room temperature. Cells were incubated with normal hamster serum, normal mouse serum (Jackson ImmunoResearch Laboratories, Inc.), and 5 µg/ml anti-CD16/32 mAb (2.4G2; BD) and then stained with FITC-, PE-, and APC-conjugated mAbs for 30 min at 4°C. To exclude dead cells, we added 7-Amino-Actinomycin D (BD) and incubated cells for 10 min immediately before analysis. Expression of markers was determined with a FACSCalibur flow cytometer (BD) and FlowJo software (Tree Star).

Table S1, included as a separate Excel file, shows details of the gene expression datasets used for discovery of the CRM across all transplanted organs, statistically significant genes identified as differentially expressed across all datasets, and the types of blood cells these genes are expressed in. Table S2, included as a separate Excel file, shows details of the independent gene expression datasets used for validation of the CRM and effect size and p-values for the CRM genes in the validation datasets. Table S3, included as a separate Excel file, shows differential expression, pathway, and graft-infiltrating cell analysis of treated cardiac mouse transplants. Table S4, included as a separate Excel files, shows multivariate Cox proportional hazards analysis of predictors of death-censored graft survival in cohort 1 (2,664 renal allograft recipients) and cohort 2 (2,515 renal allograft recipients); statin use was significantly associated with graft survival, censored for recipient death and for statin stop, independent of recipient and donor age, repeat transplantation, and calendar year.

P. Khatri and A.J. Butte were funded by the US National Cancer Institute (grant R01 CA138256). A.J. Butte was also funded by the Lucile Packard Foundation for Children’s Health and National Institute of General Medical Sciences (grant R01 GM079719).

The authors have no conflicting financial interests.