This article is mentioned in:

Introduction

Kidney transplantation is the preferred renal replacement therapy for patients with chronic renal failure (1). The development of immunosuppressive medications targeting T cell-mediated immunity has led to fewer incidences of and a more effective treatment for acute cellular rejection of kidney allografts (2). However, the currently available immunosuppressive medications do not target humoral adaptive immunity. Therefore, antibody-mediated rejection (AMR) is now one of the most challenging complications, often resulting in the deterioration of renal allograft function and graft loss in a large majority of cases (3–5). Further improvements in long-term graft survival will require a clear understanding of the mechanisms of tissue injury and the identification of biomarkers that can be used to predict AMR (6).

Allograft rejection following transplantation is currently diagnosed by histological features in biopsies. However, there remains considerable inter-observer disagreement (7,8). Therefore, a more accurate assessment of rejection would aid the attempt to reduce the failure of transplants (3,9). The limitations of donor-specific antibody (DSA) assessments and histopathological evaluations have led to increased focus on the investigation of various biomarkers to diagnose transplantation rejections. A novel approach, which combines a biopsy histopathology approach with the gene expression profiling of kidney allografts, provides a more accurate prediction for graft loss. Using a series of investigations, several studies documented the fact that the key transcripts upregulated in AMR reflect endothelial changes in the renal microcirculation characteristics and in natural killer (NK) cells (10–12).

Several genes implicated in AMR have been disclosed. A total of 23 endothelial DSA-selective transcripts (DSASTs) and AMR molecular scores have reflected changes in the micro-circulatory endothelium, which were not previously detected during routine DSA histopathology-based assessments. These provide independent values in terms of risk stratification and prognosis (11,13). Additionally, studies recently demonstrated that microvasculature injury (MVI) scores of two or more were significantly associated with a histological diagnosis of AMR, with increased DSASTs providing plausibility to the Banff MVI threshold for a complement component 4d negative AMR diagnosis (14,15).

Due to the difficulty associated with obtaining samples, particularly from human tissues, and the associated costs involved, small samples sizes are often used for microarray experiments. Integration of multiple microarray dataset has been advocated to improve the gene signature selection process (16). Meta-analyses are the most typically applied in order to detect differentially expressed genes (DEGs) (17), which may serve as candidate gene signatures for the further refinement of clinically useful biomarkers or gene signatures (18).

One previous study showed that integrating gene expression data from a number of sources or meta-analyses could lead to an increase in the statistical power of DEG detection while allowing for heterogeneity assessments. These analysis types may result in accurate, robust and reproducible predictions (19).

In the present study, a meta-analysis of array-based gene expression datasets from kidney transplantation studies was conducted to determine gene expression changes associated with AMR.

Materials and methods

Datasets

Microarray datasets from kidney transplant patients with AMR were identified by searching the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo). The search was conducted in September 2017. The key word used in the search was ‘kidney transplantation’ and the Homo sapiens search filter was applied. Each dataset was manually curated to select biopsy or peripheral blood samples from Homo sapiens. Only original experimental articles that compared the expression levels of mRNAs between patients with AMR and those without AMR (controls) were retained. All eligible publications met the following inclusion criteria: i) The studies were associated with the diagnostic value of mRNA for a diagnosis of AMR; and ii) the studies provided sufficient data with which to assess the diagnostic value of mRNA in AMR. The exclusion criteria included: i) Duplicate publications; ii) studies without sufficient data; and iii) letters, reviews, editorials, meeting abstracts and case reports. Two researchers independently screened the list of publications and evaluated the possibility of inclusion. Inter-rater agreement was assessed with Cohen's κ statistic.

Data extraction

The following information was extracted from each identified study: The GEO accession number, platform, sample type, number of cases and controls, references and expression data. Two independent reviewers extracted data from the original studies and a consensus was reached, or a third reviewer would resolve any discrepancies between the two reviewers. The studies used the gene expression platforms HG-U133 Plus 2.0, HuGene-1_0-st (both Affymetrix; Thermo Fisher Scientific, Inc., Waltham, MA, USA) and Agilent-014850 (Agilent Technologies, Inc., Santa Clara, CA, USA). In total, 9 studies were initially identified. The remaining 6 studies (GSE36059, GSE44131, GSE50084, GSE51675, GSE64261 and GSE93658) were included in the final analysis. The detailed information about the downloaded datasets is summarized in Table I. GSE36059, GSE44131, GSE50084 and GSE51675 were used for the discovery set, while GSE64261 and GSE93658 were used for the replication set. After analyzing each set, the combined dataset was used for the analysis.

Table I Information of the gene expression datasets from GEO.

Table I

Information of the gene expression datasets from GEO.

GSE	Platform	Sources	Control, n	AMR, n
Discovery set
GSE36059	HG-U133_Plus_2	Biopsy	281	65
GSE44131	HuGene-1_0-st	Biopsy	12	11
GSE50084	HuGene-1_0-st	Biopsy, whole blood	25	46
GSE51675	Agilent-014850	PBMCs	10	10
Replication set
GSE64261	Agilent-014850	PBMCs	5	5
GSE93658	HuGene-1_0-st	Biopsy	16	33

[i] GEO, Gene Expression Omnibus; GSE, GEO Series Experiments; HG-U133_Plus_2, Affymetrix Human Genome U133 Plus 2.0 Array; HuGene-1_0-st, Affymetrix Human Gene 1.0 ST Array; Agilent-014850, Whole Human Genome Microarray 4×44K G4112F; PBMCs, peripheral blood mononuclear cells; AMR, antibody-mediated rejection.

Gene expression analysis

Each individual dataset was preprocessed using the log2 transformation and normalization approach. To combine the results of the individual studies and to obtain a list of more robust DEGs between the control and AMR cases, the guidelines provided by Ramasamy et al (2008) (19) for a meta-analysis of gene expression microarray datasets were followed. The R packages MetaQC (20) and MetaDE (17,21–23) were used for quality control (QC) and for the identification of DEGs, respectively. MetaQC implements the 6 quantitative QC measures of internal QC (IQC), external QC (EQC), accuracy QC of the featured genes (AQCg), accuracy QC of the pathway (AQCp), consistency QC in the ranking of featured genes (CQCg) and consistency QC in the ranking of the pathway (CQCp). In addition, the mean rank of all QC measures in each dataset was computed as a quantitative summary score by calculating the ranks of each QC measure among all included datasets. All probe sets on the three different platforms were re-annotated to the most recent National Center for Biotechnology Information Entrez Gene Identifiers (Gene IDs), and the Gene IDs were used to cross-map genes among the three different platforms. When multiple probes matched the same gene symbol, probes presenting the greatest inter-quartile range were selected. Only genes present in all of the selected platforms were considered. The moderated t-statistic was used to calculate the P-values in each dataset, and a meta-analysis was conducted with the MetaDE package to identify DEGs using the Fisher's (24), maximum P-value (maxP) (25), rth ordered P-value (roP) (26), Stouffer (27) and naive rank summation (SR) (28) methods.

Gene set analysis (GSA) of DEGs

In order to select genes and pathways associated with AMR, genes were annotated using GSA, as conducted by GSA-SNP software version 1.0 (29). The GSA-SNP analysis uses Gene Ontology (GO) (30), the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (31,32) and the Molecular Signatures Database (MSigDB) (33). Genes that showed nominal significance levels of P<0.05 were selected.

Construction of a disease-related gene co-expression network

To identify the phenotype-related modules and protein-protein interactions (PPIs), the PPI databases BioGRID (34) and Agile Protein Interactomes Data Server (35) were downloaded, and the list of the identified genes from the present study was imported into Cytoscape version 3.5.1 (36). The molecular complex detection (MCODE) clustering algorithm was used to identify sub-network modules (37). A network score was calculated based on the complexity and density of each sub-graph. A module with an MCODE score >2 was considered significant. Post-filtering was performed to remove low-quality modules. During the filtering process, the parts of each module that showed consistent expression and high connectivity levels were selected to constitute the final module through a manual review.

Results

Short overview of included studies

The study selection process is presented in Fig. 1. Following the search and selection process, the 6 studies of GEO36059, GEO44131, GSE50084, GSE51675, GSE64261 and GSE93658 from whole blood or peripheral blood mononuclear cells and biopsy samples of renal transplant patients met the inclusion criteria. The value of Cohen's κ coefficient was 0.81, indicating a sufficient level of agreement (coefficient value >0.80). All probe sets were re-annotated with the most recent Gene IDs and then mapped, yielding 27,047 common genes across three different platforms. The resulting dataset contained 328 samples from the control and 132 samples from the AMR cases for the discovery set, and 21 samples from the control and 38 samples from the AMR cases for the replication set. Principal component analysis plots for the discovery set were plotted with MetaQC to visualize the quality of the studies via a systematic analysis (Fig. 2). The first two PCs captured 99% of the variance.

Figure 1 Study workflow. Identification of eligible gene expression datasets for a meta-analysis of AMR in KT. GEO, Gene Expression Omnibus; AMR, antibody-mediated rejection; KT, kidney transplantation.

Figure 2 Principal component analysis biplot of the quality control set, including IQC, EQC, AQCg, AQCp, CQCg and CQCp measures in the four datasets of the discovery set. QC, quality control; IQC, internal QC; EQC, external QC; AQCg, accuracy QC of the featured genes; AQCp, accuracy QC of the pathway; CQCg, consistency QC in the ranking of the featured genes; CQCp, consistency QC in the ranking of the pathway.

DEGs in antibody-mediated rejection

Five systematic analysis methods were performed by combining the P-values in the MetaDE package. These were the Fisher's, maxP, roP, Stouffer and SR methods. A total of 608 DEGs were selected by the maxP and roP methods (P<0.05) in the discovery set (Fig. 3); 291 genes were upregulated in AMR cases compared with those in the control cases, while 317 genes were downregu-lated in AMR cases compared with those in the control cases. C-X-C motif chemokine ligand 10 (CXCL10), CXCL9 and guanylate binding protein 1 (GBP1) were most significantly upregulated among the AMR cases in the discovery set. A total of 317 DEGs were selected by the maxP and roP methods (P<0.05) in the replication set; 219 genes were upregulated in AMR cases compared with those in the control cases, while 98 genes were downregulated in AMR cases compared with those in the control cases. CXCL10, CXCL9 and GBP1 were most significantly associated with AMR in the replication set. The 20 top-ranked differentially downregulated genes associated with AMR in the replication set are listed in Table II. In total, 608 DEGs were identified by the maxP and roP methods (P<0.05) in the combined discovery and replication sets; 248 genes were upregulated in AMR cases compared with those in the control cases, while 360 genes were downregulated in AMR cases compared with those in the control cases. Among the genes that showed a significant association with AMR in the present study, CXCL10, CXCL9, GBP1, complement C1q A chain (C1QA), complement C2 and interferon-induced protein with tetratricopeptide repeats 2 (IFIT2) were most significantly upregulated, while renin binding protein (RENBP), kynurenine 3-monooxygenase (KMO), thyroid hormone receptor interactor 6 (TRIP6) and glutathione S-transferase θ1 (GSTT1) were most significantly downregulated (Table III).

Figure 3 Heat map identifying differentially expressed genes in antibody-mediated rejection cases and controls of the discovery set subjected to a maxP systematic analysis when the false discovery rate was <0.05. A, antibody-mediated rejection; C, control.

Table II Top-ranked differentially upregulated genes associated with antibody-mediated rejection in kidney transplant cases in the combined discovery set with the replication set.

Table II

Top-ranked differentially upregulated genes associated with antibody-mediated rejection in kidney transplant cases in the combined discovery set with the replication set.

Gene symbol	Fisher's	maxP, roP	Stouffer	SR
CXCL10	1.00×10−20	1.00×10−20	1.00×10−20	1.00×10−20
CXCL9	1.00×10−20	1.00×10−20	1.00×10−20	1.00×10−20
GBP1	1.00×10−20	1.00×10−20	1.00×10−20	1.00×10−20
C1QA	1.00×10−20	1.00×10−20	1.00×10−20	1.00×10−20
C2	1.00×10−20	1.00×10−20	1.00×10−20	1.21×10−5
IFIT2	1.00×10−20	1.00×10−20	1.00×10−20	2.17×10−5
STAT1	1.00×10−20	3.22×10−6	1.00×10−20	1.61×10−6
GNG11	1.00×10−20	3.22×10−6	1.00×10−20	4.11×10−5
CD180	1.00×10−20	3.22×10−6	1.00×10−20	7.81×10−5
TBX21	1.00×10−20	3.22×10−6	4.03×10−6	3.27×10−4
WARS	1.00×10−20	4.03×10−6	1.00×10−20	8.06×10−7
ANKRD22	1.00×10−20	4.03×10−6	1.00×10−20	3.22×10−6
OAS3	1.00×10−20	4.83×10−6	1.00×10−20	3.96×10−4
DOCK4	1.00×10−20	7.25×10−6	8.06×10−20	5.94×10−4
PARP14	1.00×10−20	8.86×10−6	1.00×10−20	5.64×10−6
PIK3AP1	1.61×10−6	1.21×10−5	1.00×10−20	2.82×10−5
IFI44L	1.00×10−20	1.29×10−5	1.00×10−20	1.35×10−4
MX2	1.00×10−20	1.69×10−5	1.00×10−20	1.37×10−5
BST2	1.00×10−20	1.85×10−5	1.85×10−5	4.03×10−6
DDX60	7.41×10−5	2.17×10−5	1.61×10−5	2.13×10−3

[i] maxP, maximum P-value method; roP, rth ordered P-value method; SR, naive rank summation method.

Table III Top-ranked differentially downregulated genes associated with antibody-mediated rejection in kidney transplant cases in the combined discovery set with the replication set.

Table III

Top-ranked differentially downregulated genes associated with antibody-mediated rejection in kidney transplant cases in the combined discovery set with the replication set.

Gene symbol	Fisher's	maxP, roP	Stouffer	SR
RENBP	8.94×10−4	8.86×10−6	1.53×10−4	4.72×10−3
KMO	3.54×10−5	1.61×10−5	6.44×10−6	2.18×10−4
TRIP6	2.21×10−5	1.85×10−5	6.59×10−4	1.02×10−2
GSTT1	2.50×10−5	3.46×10−5	7.25×10−6	8.63×10−4
KCNJ10	9.51×10−5	5.40×10−5	3.38×10−5	2.34×10−3
CA2	1.86×10−4	7.97×10−5	3.30×10−5	2.43×10−3
PINK1	4.91×10−3	1.10×10−4	1.46×10−3	1.10×10−2
MCAT	5.78×10−4	1.18×10−4	1.28×10−4	4.36×10−3
SLC10A2	5.70×10−3	1.46×10−4	1.69×10−3	1.21×10−2
MYO7A	1.85×10−3	1.89×10−4	4.54×10−4	6.55×10−3
SDC2	3.83×10−4	1.76×10−4	1.14×10−4	6.19×10−3
MYO7A	1.85×10−3	1.89×10−4	4.54×10−4	6.55×10−3
AGAP1	1.42×10−3	1.93×10−4	3.85×10−4	7.91×10−3
ERBB2	9.26×10−5	3.25×10−4	4.11×10−5	2.92×10−3
JUP	1.36×10−4	3.38×10−4	3.95×10−5	1.19×10−2
ZNF358	6.66×10−3	3.83×10−4	2.85×10−3	2.09×10−2
AHCY	6.33×10−4	4.14×10−4	2.75×10−4	5.04×10−3
TMEM37	1.19×10−3	4.22×10−4	3.79×10−4	5.27×10−3
UCHL1	2.32×10−3	5.09×10−4	7.55×10−4	1.46×10−2
XYLB	1.02×10−3	5.36×10−4	3.75×10−4	6.67×10−3

[i] maxP, maximum P-value method; roP, rth ordered P-value method; SR, naive rank summation method.

GSA of DEGs

In the combined discovery and replication sets, the GO algorithm identified >231 gene sets associated with AMR (all P<0.05). The highly ranked gene sets are shown in Table IV. The KEGG pathway identified 31 gene sets associated with AMR (all P<0.05). The MSigDB identified >1,720 gene sets associated with AMR (all P<0.05). The regulation of the immune response gene set, in this case β2-microglobulin (B2M), baculoviral IAP repeat-containing 3 (BIRC3), C1QA, CD3e molecule (CD3E), CD79A, major histocompatibility complex class II DM α (HLA-DMA), HLA-DOA, HLA-DOB, HLA-DPB1, HLA-DQA1 and HLA-DRA, was significantly associated with AMR (P<1×10−5). Allograft rejection genes, in this case HLA-DMA, HLA-DOA, HLA-DOB, HLA-DPB1, HLA-DQA1, HLA-DRA, HLA-DMB, tumor necrosis factor (TNF) and granzyme B, were also significantly associated with AMR (P<1×10−5). Hallmark interferon-γ response genes, specifically transporter 1, ATP binding cassette subfamily B member (TAP1), MX dynamin-like GTPase 2, CXCL10, GBP4, CXCL9, interleukin 10 receptor subunit α, SAM and HD domain containing deoxynucleoside triphosphate triphosphohydrolase 1 and intercellular adhesion molecule 1, were also significantly associated with AMR (P<1×10−5).

Table IV Top-ranked gene sets associated with antibody-mediated rejection traits in kidney transplant cases.

Table IV

Top-ranked gene sets associated with antibody-mediated rejection traits in kidney transplant cases.

Set name	z-score
Gene Ontology
MHC protein complex	14.003
Regulation of immune response	13.938
Immune response-activating cell surface receptor signaling pathway	13.358
Antigen receptor-mediated signaling pathway	12.653
Positive regulation of immune response	12.517
Interferon-γ-mediated signaling pathway	12.211
KEGG pathway
hsa04612_Antigen processing and presentation	12.096
hsa05330_Allograft rejection	12.087
hsa04145_Phagosome	11.901
hsa04514_Cell adhesion molecules (CAMs)	10.889
hsa04650_Natural killer cell-mediated cytotoxicity	7.916
hsa04620 _Toll-like receptor signaling pathway	7.598
MsigDB
Hallmark_allograft_rejection	18.105
Ichiba_graft_versus_host_disease_35d_up	15.368
Rodwell_aging_kidney_up	14.291
Flechner_biopsy_kidney_transplant_rejected_vs_ok_up	14.205
Hallmark_interferon_γ_response	13.654
GSE10325_CD4_Tcell_vs_Bcell_DN	12.701

[i] KEGG, Kyoto Encyclopedia of Genes and Genomes; MSigDB, Molecular Signatures Database.

Construction of a disease-related gene co-expression network

There were two modules with a MCODE score of >2 in AMR cases. Fig. 4 shows the sub-network modules. One module consisted of immunoglobulin heavy constant γ1 (G1m marker) (IGHG1), Fc fragment of IgG receptor IIIa (FCGR3A), FCGR3B, amyloid P component serum (APCS) and glycerol kinase (GK), and the other module consisted of CD74, HLA-DMA, HLA-DMB, HLA-DRA, HLA-DRB1, GBAA0374, proteasome subunit β5 (PSMB5), PSMB8, PSMB9, TAP1 and TAP2.

Figure 4 AMR-related network modules. The solid lines indicate the interactions between a protein and a protein. The oval-shaped nodes are proteins that interact. (A) Modules showing proteins of CD74, HLA-DMA, HLA-DMB, HLA-DRA, HLA-DRB1, GBAA0374, PSMB5, PSMB8, PSMB9, TAP1 and TAP2, which are associated with AMR in kidney transplant cases. (B) Modules showing proteins encoded by IGHG1, FCGR3A, FCGR3B, APCS and GK, which are associated with AMR. AMR, antibody-mediated rejection.

Discussion

Biopsy histopathology combined with gene expression profiling in kidney allografts can provide more accurate predictions of graft loss compared with histopathology alone. A meta-analysis of a gene expression profile is a powerful tool for elucidating gene signatures, and this method has been widely applied to augment statistical power and provide validated conclusions (17). In the present study, DEGs were identified that were associated with AMR in kidney transplant cases using a meta-analysis from GEO data.

The results showed that the upregulated expressed genes of CXCL10, CXCL9, GBP1 and C1QA were associated with the occurrence of AMR in kidney transplant cases. The CXCR3 ligand chemokines, CXCL9 and CXCL10, which are interferon-γ (IFNG)-induced small cytokines, can act as biomarkers and are being increasingly investigated as screening tools for the early diagnosis of renal transplantation dysfunction (38,39). Once the antibody recruits NK cells by virtue of IgG-Fc and complement receptors, NK cells can release IFNG. In the present study, the meta-analysis revealed that the upregulated expressions of the IFNG response genes of signal transducer and activator of transcription 1, IFI30 lysosomal thiol reductase, PSMB 8-9-10 and GBP1 were associated with acute rejection, which is consistent with the findings of a previous study (40).

Genes that were downregulated, including RENBP, KMO, TRIP6, and GSTT1 were associated with acute rejection in the present study. The dominant negative mutant or RNA-mediated interference of TRIP6 reportedly inhibits nuclear factor-κB activation by TNF, IL-1, Toll-like receptor 2 or Nucleotide-binding oligomerization domain-containing protein 1 (41). Detection of anti-GSTT1 antibodies in a recipient with a GSTT1-null genotype prior to transplantation was reported to be predictive of graft rejection in the event of a GSTT1-positive donor (42).

The GSA method has been used for the identification of meaningful associations between markers and diseases or traits of interest in a large set of genes or proteins (29). When a GSA analysis of disease conditions was conducted in the present study, genes or pathways commonly involved with the immune response were found. A network analysis also identified the HLA family. DSA against the HLA antigens is critical with regard to the diagnosis of AMR (43). Additionally, IGHG1, FCGR, APCS and GK were noted. Mice that were deficient in molecules essential for the recognition, internalization or lysosomal DNA degradation of apoptotic cells, including serum amyloid P, were previously reported to develop systemic autoimmune disorders (44).

In conclusion, the present study identified DEGs associated with AMR in kidney transplant cases using a meta-analysis from publicly available datasets. Although additional challenges are encountered when attempting to define the role of genes that affect the pathophysiological mechanism of AMR, it is hoped that these results will lead to a more in-depth understanding of the molecular mechanism of AMR. Promising genes or pathways can be utilized as drug targets. Future studies are required to validate the identified DEGs and pathways.

Acknowledgments

Not applicable.

Funding

This study was supported by the Seoul National University Research Grant in 2017 (no. 370C-20170060) and the National Research Foundation of Korea grant funded by the Korea government (MSIT) (no. NRF-2018R1A2B6001859).

Availability of data and materials

The datasets analyzed during the current study are available from the corresponding author on reasonable request.

Authors' contributions

IWK, NH, SK, YSK and JMO designed the study. IWK, SK, and YSK planned the statistical analysis. IWK and JHK analyzed the study results and wrote and revised manuscript. NH, SK, YSK and JMO reviewed the manuscript. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References