Visual Overview

Abstract

Background and objectives At least 20 susceptibility loci of IgA nephropathy have been identified by genome-wide association studies to date. Whether these loci were associated with disease progression is unclear.

Design, setting, participants, & measurements We enrolled 613 adult patients with IgA nephropathy for a follow-up of ≥12 months. All 20 IgA nephropathy susceptibility loci were selected and their tag single nucleotide polymorphisms (SNPs) were genotyped. After strict quality control, 16 SNPs and 517 patients with IgA nephropathy were eligible for subsequent analysis. Progression was defined as ESKD or 50% decrease in eGFR. A stepwise Cox regression analysis of all SNPs on Akaike information criterion was performed to select the best model.

Introduction

IgA nephropathy is the most common primary GN worldwide. The formation of galactose-deficient IgA1 and antibodies specific for its hinge-region glycans were thought to be key steps of IgA nephropathy development (1–3). Familial aggregation and various reports on prevalence and prognosis in patients with different origins suggested a genetic contribution (4–6). To date, there have been five genome-wide association studies (GWAS) carried out in patients with IgA nephropathy with different ancestries, and at least 20 susceptibility loci have been identified (7–11).

The prognosis of IgA nephropathy is highly variable and it is challenging for clinicians to accurately predict kidney outcome at the time of biopsy. Studies on risk factors were carried out and risk models were established by combining independent clinical or pathologic risk factors together to enhance the predictable power of models (12–19). Most recently, an extended multicenter study developed and validated two risk models by using clinical and Oxford predictors on the basis of Chinese patients with IgA nephropathy. The clinical model from this study contained age, sex, eGFR, hemoglobin, and urine protein excretion. The combined model (clinical–pathologic model) included age, eGFR, and Oxford M and T scores (19). However, no risk prediction models tested the contributions of genetic variants. It is possible that incorporating genetic risk factors such as susceptibility loci identified in GWAS to current risk models could enhance their performance because the multihit pathogenesis model suggests these loci may potentially influence the pathogenesis and severity of IgA nephropathy (3). A post-GWAS revealed that △CFHR3,1 (a deletion variant encompassing CFHR1 and CFHR3) underlining the top signal on chromosome 1 might reduce tubulointerstitial injury (20), and tag single nucleotide polymorphisms (SNPs) in HLA-DQ/DR region might accelerate IgA nephropathy progression (10,21). Further, a genetic risk score including nine susceptibility loci identified in GWAS could predict ESKD in IgA nephropathy (22). However, these results need to be validated and whether genetic risk factors could add predictive value to clinical–pathologic risk equations needs to be further evaluated. Thus, we genotyped 20 SNPs representing all IgA nephropathy susceptibility loci from previous GWAS and select the best genetic model by using a stepwise Cox regression analysis of SNPs passing quality control process on Akaike information criterion (AIC). The combined effect of candidate SNPs was tested and a genetic risk score was established. We then compared the genetic risk score derived from this study and two previously reported genetic risk scores (10,22). Finally, we tested whether the genetic risk score in this study could add predictable value to the most recently published IgA nephropathy risk prediction models (19).

Materials and Methods

Study Design

The enrolled patients were Chinese Han and were biopsied in Ruijin Hospital between 1989 and 2014. IgA nephropathy was defined by kidney biopsy (23). Inclusion criteria were age at biopsy of ≥18 years old, a follow-up time of ≥12 months, and available baseline data. Exclusion criteria were eGFR<15 ml/min per 1.73 m2 at the time of biopsy or with any evidence of Alport syndrome, thin basement membrane disease, SLE, and Henoch–Schonlein purpura or other systemic diseases. A total of 613 patients with IgA nephropathy were enrolled for genotyping 20 candidate SNPs. A total of 517 patients and 16 SNPs were recruited for the subsequent analysis after quality control. Among all patients, Oxford-MESTC (M, mesangial hypercellularity; E, presence of endocapillary proliferation; S, segmental glomerulosclerosis/adhesion; T, severity of tubular atrophy/interstitial fibrosis; C, presence of crescent) score was available in 401 patients.

Detailed baseline data were collected from all patients at the time of kidney biopsy. eGFR was calculated using the CKD Epidemiology Collaboration equation (24). Histologic changes were scored using the Oxford-MESTC score (23,25,26). The start of follow-up was taken as kidney biopsy. The primary outcome was defined as a combination of ESKD and a 50% reduction of eGFR from baseline. In addition to supportive therapy, patients with hypertension and/or proteinuria >0.5 g/24 h were treated with renin-angiotensin system blocker. Immunosuppressive agents were added if patients had massive proteinuria or no response to renin-angiotensin system blocker, presented with rapidly progressive GN, or pathology showed massive crescents, severe inflammatory cell infiltration, or vascular necrosis.

This study was approved by the Institutional Review Board of Ruijin Hospital, Shanghai Jiao Tong University School of Medicine and was conducted in accordance with the principle of the Helsinki Declaration. Written informed consent was obtained from all involved participants.

SNP Genotyping

Twenty SNPs independently associated with IgA nephropathy in previous GWAS (7–11) were selected. Genomic DNA was isolated from peripheral blood sample using the Fujifilm Blood DNA Kit and stored at −80°C for genotyping. All SNPs were genotyped using Sequenom MassArray platform, which combined multiple PCR, Mass Array iPLEX single base extension chemistry, and matrix-assisted laser desorption/ionization-time of flight mass spectrometry. PCR primers and extension probes were designed with the Sequenom online primer design system. The Spectro-CHIP bioarray with purified products was detected using a matrix-assisted laser desorption/ionization-time of flight mass spectrometer. The final results were analyzed using MassArray Workstation software (Typer 4.0).

Statistical Analyses

Continuous variables were summarized as means±SDs (normally distributed) or median (range) (non-normally distributed). Categorical variables were presented as frequencies (percent). We tested for Hardy–Weinberg equilibrium and pairwise linkage disequilibrium (LD) for all SNPs. Haplotypes were phased using PLINK v1.07 and their frequencies were estimated between progressive group and nonprogressive group on the basis of whether the patient progressed to combined outcome within 5 years. Associations between SNPs and kidney function progression were estimated by Cox proportional hazard model, assuming underlying dominant, recessive, and log-additive models. A backward stepdown selection process using AIC was used to select the best outcome prediction model. The selection was stopped with AIC not decreasing after removing the last SNP and we finally chose the model with the lowest AIC. A genetic risk score was calculated based on genotypes of each SNP, which was simply scored using the risk genotypes with a dominant model of SNPs included in the best model. A genetic risk stratification was established on the basis of the tertiles of genetic risk score and separated all patients to low genetic risk group (genetic risk score 4–5), middle genetic risk group (genetic risk score 6), and high genetic risk group (genetic risk score 7–8). Progression rates among different genetic risk groups were compared by Kaplan–Meier analysis. The contributions of genetic risk score to a clinical risk model and a clinical–pathologic risk model (19) were evaluated. Reclassification improvement was assessed using net reclassification improvement (NRI) according to the 5-year risk of combined outcome. The NRI was defined as ([number of events with increased predicted risk−number of events with decreased predicted risk]/number of events+[number of nonevents with decreased predicted risk−number of nonevents with increased predicted risk]/number of nonevents) and estimated according to the method by Pencina et al. (27,28), which was calculated by R package “PredictABEL” (29) with our specific data in this study. We compared this genetic risk score with two established genetic risk scores from prior studies on 5-year progression risk. The first was the standardized genetic risk equation (standardized genetic risk score) by Kiryluk et al. (10), using β-coefficients of the allelic risk. The second used an unweighted approach of nine SNPs (unweighted nine-SNP genetic risk score) by Zhou et al. (22), which was the simple counts of the total number of risk alleles. Four SNPs failed to pass our quality control and were removed from the above two equations, and rs9357155 was replaced by rs2071543, which is in its LD segment (r2=1.00, calculated by SNAP server). (30) The c-statistic was estimated as an area under the receiver operating characteristic curve to predict the 5-year progression. All statistical analyses were performed using SNPassoc package version 1.9–2, survival package version 2.40–1, ggkm package version 1.0, PredictABEL package version 1.2–2, and pROC package version 1.11.0, with R version 3.2.2.

Results

Study Cohort and SNPs Genotyping

Twenty susceptibility loci from five GWAS were genotyped in 613 patients with IgA nephropathy (7–11). Four SNPs and 96 patients were excluded from this study because they did not pass the quality control filters, which included a per-SNP call rate >95%, per-individual call rate of 100%, and minor allele frequency >0.10. In detail, call rates of rs6677604 (CFH c.1696+2019G>A), rs4077515 (CARD9 p.Ser12Asn), and rs10086568 (DEFA3A/DEFA5A intergene region) were <95%, and the minor allele frequency of rs11574637 (ITGAX c.430+189T>C) was <10%. Individual call rates of 96 patients were <100%. Finally, 16 SNPs and 517 patients with IgA nephropathy were eligible for further analysis (Figure 1). Of these patients, 266 were men and 251 were women. The mean age at biopsy was 38±12 years old, and mean eGFR was 72±32 ml/min per 1.73 m2. Median baseline proteinuria was 1.27 (range, 0.02–12.70) g/24 h. Ninety-six (19%) patients reached the combined outcome during a median follow-up time of 50 (range, 12–238) months (Table 1). Information of the 16 SNPs were shown in Table 2.

Pairwise Interaction Screen among Enrolled SNPs

All 16 SNPs were tested for Hardy–Weinberg equilibrium and inter-SNPs LD. Most of the SNPs agreed with the Hardy–Weinberg equilibrium (Table 2). LD analysis showed that HLA-DQ/DR rs2856717 was in high LD with rs9275224 and rs9275596 (r2=0.38 and r2=0.53, respectively). Two variants around DEFA region rs12716641 and rs2738048 were also in partial LD (r2=0.21) (Supplemental Table 1).

Establish the Best Genetic Predictive Model and Genetic Risk Score

Backward stepwise Cox regression analysis were performed to choose the best genetic predictive model. The best model contained rs11150612, rs7634389, rs2412971, and rs2856717 (Table 3). The genetic risk score using the four SNPs was associated with disease progression before (hazard ratio [HR], 1.65; 95% confidence interval [95% CI], 1.29 to 2.12) and after adjustment by the clinical risk model (HR, 1.29; 95% CI, 1.03 to 1.62) and the clinical–pathologic risk model (HR, 1.35; 95% CI, 1.03 to 1.77) (Table 3). Then, we stratified all patients with IgA nephropathy into three groups by the tertile of genetic risk score. The median outcome-free time in low-, middle-, and high-risk groups were 191, 111, and 87 months, respectively (log-rank test, P=0.001) (Figure 2A). Compared with the low-risk group, patients in the middle-risk group had a 2.12-fold (95% CI, 1.33 to 3.40) increase in the risk of disease progression, whereas patients in the high-risk group had a 3.61-fold (95% CI, 2.00 to 6.52) risk increase (Table 3).

Cox Regression Analysis on Haplotypes

We carried out haplotype analysis of four SNPs (rs7763262, rs9275224, rs2856717, and rs9275596) in the HLA-DQ/DR region. We used the haplotype CGCT (0.65 in progressive group and 0.66 in nonprogressive group; P=0.78), which was most common and equally distributed in two groups as a reference haplotype. The result showed no specific haplotype had increasing progression risks (Supplemental Table 2).

We compared genetic risk score established in our study with a previously published standardized genetic risk score and unweighted nine-SNP genetic risk score on progression prediction. No significant association was seen between combined outcome and the two published equations (Supplemental Table 3).

To test if the performance could be improved, we added genetic risk score to the clinical risk model and the clinical–pathologic risk model to predict 5-year progression. An increase in r2 (from 17.4% to 18.5% in the clinical model and from 15.1% to 16.8% in the clinical–pathologic model) and c-statistic (from 0.83 to 0.86 in the clinical model and from 0.82 to 0.85 in the clinical–pathologic model) showed a better model fit (Figure 2, B and C, Table 4). The significant increase of continuous NRI (0.42; 95% CI, 0.12 to 0.72) and categorical NRI (0.14; 95% CI 0.05 to 0.22) after adding genetic risk score to the clinical risk model indicated improvement in predictive discrimination (Supplemental Table 4, Table 4).

The discrimination performance of different models predicting the risk of combined outcome at 60 months of follow-up

Discussion

The high heterogeneity of prognosis in patients with IgA nephropathy brings difficulties to patient management and treatment. Risk factors were identified, including baseline kidney function, BP, proteinuria, serum albumin, hemoglobin, blood type, body weight, Oxford-MESTC score, and follow-up proteinuria and BP (13–19). To enhance the predicting accuracy, risk equations were established that showed better performance in predicting prognosis (12–15,17,19,31). Xie et al. (19) developed and validated two risk equations on the basis of a multicenter study: the clinical risk model used five clinical variables, including age at biopsy, sex, baseline eGFR, hemoglobin, and proteinuria; the clinical–pathologic risk model used two clinical and two pathologic variables, including age, baseline eGFR, and Oxford M and T score. The two risk equations showed a good performance for predicting disease progression, but the accuracy still needs to be further enhanced. Variants on candidate genes including HLA-DQB, ACE, GPIa, ICAM-1, MYH9, TNFSF13, DEFA, CFHR, C1GALT1, and ST6GALNAC2 were associated with phenotypes and progression of patients with IgA nephropathy (20,32–36). Here we systemically studied 20 SNPs derived from all five GWAS in an extended Chinese cohort and investigated associations of these loci with IgA nephropathy progression. We found the best genetic model, using rs11150612, rs7634389, rs2412971, and rs2856717, was independently associated with IgA nephropathy progression after adjusted by clinical and pathologic variables. Notably, we determined that a genetic risk score using the four-SNP model could increase the performance of the most recently reported clinical and clinical–pathologic risk models of IgA nephropathy from a large-scale risk assessment study (19).

Zhou et al. (22) established the unweighted nine-SNP genetic risk score by nine SNPs to predict ESKD in a Chinese cohort of 297 patients with IgA nephropathy. Kaplan–Meier plots showed patients with the unweighted nine-SNP genetic risk score ≥16 had a poorer prognosis than patients with the unweighted nine-SNP genetic risk score <16. Kiryluk et al. (10) developed a 15-SNP IgA nephropathy genetic risk calculator, a good tool to predict IgA nephropathy onset risk. The above two equations were validated in our cohort, yet showed no significant association with progression. This may be because of racial differences and the four removed SNPs. Besides, risk alleles of disease onset and progression may be different.

We identified a four-SNP risk model, rs11150612 (ITGAM-ITGAX), rs7634389 (ST6GAL1), rs2412971 (HORMAD2), and rs2856717 (HLA-DQ/DR), that is associated with IgA nephropathy progression. The underlying pathogenic mechanism might be autoantigen formation, antigen-presenting, and thereafter formation and modification of autoantibody regulated by ITGAM-ITGAX, HORMAD2, HLA-DQ/DR, and ST6GAL1 taking part in IgA nephropathy development and progression. HLA-DR/DQ, coding for MHC II molecules and involved in antigen presentation, was detected by all five GWAS regardless of racial origins (7–11). HLA-DRB1 classic alleles have been reported to be related to disease progression in a Chinese IgA nephropathy cohort (37). In this study, we determined patients carrying rs2856717-T had an increased risk for disease progression. The result replicated a recent study (21) from southwest China showing associations between rs2856717/TT or TC and increased risk for IgA nephropathy progression. To further dissect the HLA region, we next phased four SNPs (rs7763262, rs9275224, rs2856717, and rs9275596) in the HLA-DQ/DR region and tested their associations with disease progression. No specific haplotype showed an increasing risk for disease progression, which might be because of the limited sample size of our study. ST6GAL1 encodes ST6Gal-1 sialyltransferase, which involves adding terminal α2,6-linked sialic acids to N-glycans within the Fc domain of IgG (38). Crosslinking reaction between IgG and Fc receptors on the surface of immune cells is the prominent factor in autoimmune diseases like SLE, ANCA-associated vasculitis, and rheumatoid arthritis (39,40). Data showed that the mRNA level of ST6GAL1 increases significantly in the glomeruli of patients with IgA nephropathy compared with normal kidney tissue. Notably, a higher mRNA level of ST6GAL1 in patients with IgA nephropathy with normal kidney function (eGFR>90 ml/min per 1.73 m2) than that in patients with IgA nephropathy with mild kidney dysfunction (eGFR 60–90 ml/min per 1.73 m2) was reported (https://www.nephroseq.org/resource/main.html). Rs2412971 is located in the intronic region of HORMAD2, which belongs to a HORMA domain-containing gene family and plays a role in DNA repair (41). The HORMAD2 locus also encodes several cytokines that interact with complement pathways (42). Kiryluk et al. (6) detected a multiplicative interaction between rs6677604 (CFHR3/R1) and rs2412971 (HORMAD2) loci. Rs2412971-A allele is reversed among CFHR3,1-del homozygotes. The biologic basis of this statistical finding remains unknown. ITGAM-ITGAX encodes integrin αM and αX, which mark intestinal dendritic cells that maintain the balance between inflammation and tolerance. ITGAM and ITGAX also combine with the integrin β2 chain to form leukocyte-specific complement receptors 3 and 4. ITGAM was found to take part in the regulation of intestinal IgA-producing plasma cells in mice (43).

As we know, impaired kidney function, increased BP, proteinuria, and MESTC score are risk factors for a poor outcome among patients with IgA nephropathy (12–15,17,25,26,30,44,45). We proposed that the genetic risk score could add predictive value to the current clinical and pathologic equations. After incorporating the genetic risk score to Xie et al.’s clinical and clinical–pathologic equations, the new equations showed an increased discrimination in predicting IgA nephropathy prognosis. Our study determined that genetic factors could help clinicians to predict kidney outcomes of IgA nephropathy more precisely in addition to clinical and pathologic risk models.

There were some limitations in our study. First, it was a retrospective study from a single center. Second, because of the sample size of our cohort, we did not have enough power to analyze rare SNPs or SNPs with low effect size. Third, we only included SNPs identified by GWAS, other genetic variants might contribute to disease progression and need to be further evaluated. At last, the four SNPs (rs6677604, rs4077515, rs10086568, and rs11574637) excluded from this study during quality control were still possibly associated with IgA nephropathy progression. Larger prospective designed study with different ethnic populations will be helpful to overcome these limitations.

In conclusion, we established a genetic risk score using rs11150612 (ITGAM-ITGAX), rs7634389 (ST6GAL1), rs2412971 (HORMAD2), and rs2856717 (HLA-DQ/DR), which was independently associated with IgA nephropathy progression. The genetic risk score could enhance the performance of clinical or clinical–pathologic risk models of IgA nephropathy.

Disclosures

None.

Acknowledgments

This work was supported by grants from the National Key Research and Development Program of China (2016YFC0904100), National Natural Science Foundation of China (81570598, 81370015, and 81000295), Science and Technology Innovation Action Plan of Shanghai Science and Technology Committee (17441902200), Shanghai Municipal Education Commission, Gaofeng Clinical Medicine grant (20152207), the Chinese Medical Association clinical research special fund (13030280413), and Shanghai Jiao Tong University School of Medicine multicenter clinical research project (DLY201510).

Footnotes

Published online ahead of print. Publication date available at www.cjasn.org.