Visual Abstract
Abstract
Background and objectives IgA nephropathy is the most common form of primary GN worldwide. The evidence of geographic and ethnic differences, as well as familial aggregation of the disease, supports a strong genetic contribution to IgA nephropathy. Evidence for genetic factors in IgA nephropathy comes also from genome-wide association patient-control studies. However, few studies have systematically evaluated the contribution of coding variation in IgA nephropathy.
Design, setting, participants, & measurements We performed a two-stage exome chip–based association study in 13,242 samples, including 3363 patients with IgA nephropathy and 9879 healthy controls of Han Chinese ancestry. Common variant functional annotation, gene-based low-frequency variants analysis, differential mRNA expression, and gene network integration were also explored.
Results We identified three non-HLA gene regions (FBXL21, CCR6, and STAT3) and one HLA gene region (GABBR1) with suggestive significance (Pmeta<5×10−5) in single-variant associations. These novel non-HLA variants were annotated as expression-associated single-nucleotide polymorphisms and were located in enhancer regions enriched in histone marks H3K4me1 in primary B cells. Gene-based low-frequency variants analysis suggests CFB as another potential susceptibility gene. Further combined expression and network integration suggested that the five novel susceptibility genes, TGFBI, CCR6, STAT3, GABBR1, and CFB, were involved in IgA nephropathy.
Conclusions Five novel gene regions with suggestive significance for IgA nephropathy were identified and shed new light for further mechanism investigation.
Introduction
IgA nephropathy is the most common form of primary GN worldwide, and as a progressive disease, it is an important cause of kidney failure (1⇓–3). The evidence of geographic and ethnic differences and familial aggregation of the disease support a strong genetic contribution to IgA nephropathy. Of the many pathogenic mechanisms reported, defects in IgA1 glycosylation that lead to formation of immune complexes have been consistently indicated (1⇓–3). It has been demonstrated by several studies that the estimated heritability of aberrant IgA1 glycosylation is 54%–80% (4⇓–6). Evidence for genetic factors in IgA nephropathy also comes from genome-wide association patient-control studies (GWAS) (7⇓⇓⇓–11). Recent work has identified a number of genetic factors, mostly associated with mechanisms of defense against infection (12,13). While emphasizing the important new insights that have emerged from GWAS of IgA nephropathy, it should not be forgotten that even the latest and largest GWAS meta-analysis attributes only 7.6% of disease risk to genetic factors in Chinese populations and 6.2% in Europeans.
GWAS-implicated single-nucleotide polymorphisms (SNPs) usually confer relatively small effects and lie within nonprotein-coding sequence, suggesting additional genetic variants yet to be discovered (14⇓–16). Few studies have systematically evaluated the contribution of coding variation in IgA nephropathy. However, a study for coding variants primarily relies on gene sequencing that is typically expensive, especially to systematically scale the whole genome in a number of samples. Recently, Illumina HumanExome BeadChips, also referred to as exome chips, were designed as a comprehensive and cost-effective alternative to explore the role of coding variants in human complex diseases, which made it applicable in large cohorts. Since then, the exome arrays have been applied to studies for complex diseases and have detected several significant coding variants; therefore, they provide great opportunities to further disentangle disease biology (14⇓–16).
To provide insights into the biologic mechanisms of IgA nephropathy, we performed a two-stage exome chip association study in 13,242 samples, including 3363 patients with IgA nephropathy and 9879 healthy controls of Han Chinese ancestry.
Materials and Methods
Participants
All patients with IgA nephropathy were recruited from Peking University First Hospital, and all controls were healthy blood donors. The study was approved by the Institutional Review Board of Peking University First Hospital, China, and it was conducted according to the Declaration of Helsinki principles (institutional review board nos. 2013–548 and IRB00001052–17089, Chinese Human Genetic Resources 2008–014).
The cohort in total was composed of 3363 patients with IgA nephropathy and 9879 healthy controls (Table 1, Supplemental Figure 1, Supplemental Table 1). The diagnosis of IgA nephropathy was defined by kidney biopsy according to histopathologic detection of immunofluorescence showing at least 2+ (scale =0–3+) mesangial deposition of IgA, with IgA comprising the dominant Ig deposited in the glomeruli and excluding individuals with cirrhosis, IgA vasculitis, hepatitis B–associated GN, HIV infection, Alport syndrome, thin basement membrane disease, SLE, or other systemic diseases.
Sample information for different arrays involved in this study
Genotyping and Quality Controls in the Discovery Stage
There were 640 patients and 4295 controls genotyped by Illumina Human Exome V1.2 BeadChip in the discovery stage. This array contains about 240,000 variants (http://genome.sph.umich.edu/wiki/Exome_Chip_Design), 93.5% of which had minor allele frequency (MAF) ≤5%. The genotype calling and clustering were conducted using GenomeStudio V.2011.1. Illumina human exome genotyping array clustering and quality control were undertaken as reported (17,18).
Samples with call rates <95% were excluded for further analysis. Potential genetic relatedness in terms of pairwise identity by state was examined using Plink V.1.90 (http://pngu.mgh.harvard.edu/purcell/plink/). Close relatives were identified as having an estimated identity by descent (IBD) proportion of genome-wide shared alleles >0.1875, which was halfway between the expected IBD for third- and second-degree relatives. Sample outliers were checked on the basis of the principal components analysis method using Eigensoft 7.2.1 tools. There were 258 subjects excluded in the further analysis due to high heterozygosity rate (44 subjects; over ±3 SD from the mean), failure of IBD test (162 subjects), and principal component outliers (52 subjects; continental outliers were checked, and individual outliers were defined by over ±3 SD for the first two principal components) (Supplemental Figure 2).
SNPs were excluded if they showed call rate <95%, MAF<0.001, or significant departure from Hardy–Weinberg equilibrium in controls (P<10−4). At the end of quality control, genotype data for 33,655 autosomal biallelic variants in 601 patients and 4076 controls qualified for association tests in the discovery stage. All of the analyses were accomplished in Plink v1.90.
Single-Nucleotide Polymorphisms Selection and Genotyping in the Replication Stage
The flow chart of the study design is shown in Supplemental Figure 1. A total of 144 non-HLA SNPs of MAF≥1% with association P<0.005 in the discovery study were genotyped for replication in an independent cohort of 8565 samples (2762 patients and 5803 controls; Xu-jie Zhou et al. unpublished data with the same criteria of quality control) using Illumina Human 610-Quad BeadChip (610), Illumina Infinium OmniZhongHua-8 v1.3 BeadChip (Zhonghua8), or Illumina Infinium Global Screening Array-24 v1.0 BeadChip (GSA).
Histocompatibility Leukocyte Antigen Imputation
To identify association at the MHC region, we extracted the SNPs within the extended MHC region (chr6: 25–34 Mb) and imputed classic alleles, untyped SNPs, and amino acid residues using SNP2HLA with the Asian HLA imputation reference panel (19,20). Imputation of genotypes for patients and controls was performed together. In order to reduce the uncertainty of imputation, we restricted our analysis to variants with high imputation quality (r2>0.8). Imputed dosage results were used in all subsequent analyses.
Functional Annotation of Novel Variants
The regulatory effects for new variants were annotated in Haploreg V4.1. B cells have been widely implicated in the risk of IgA nephropathy, and thus, the effects of newly identified SNPs were evaluated in primary B cells from the National Institutes of Health Roadmap Epigenomics Consortium (2015). The expression quantitative trait locus (eQTL) effects for each novel SNP were obtained from public databases. The newly identified genes and other genes that are the closest or the ones reported previously in IgA nephropathy GWAS were used in STRING v11.0 to explore gene networks.
Quantification of Genes Expression
A total of 10 ml whole blood was collected in an EDTA vacutainer and centrifuged at 600×g for 30 minutes to remove plasma. Then, blood cells were diluted with PBS and layered on top of Ficoll-Paque (GE Healthcare). After centrifugation at 400×g for 30 minutes, the PBMC fraction was collected and washed in PBS to remove excess Ficoll-Paque. Platelets were removed following a second centrifugation at 300×g for 10 minutes. The extraction of RNA from PBMCs (obtained from 59 patients with IgA nephropathy and 28 healthy donors who were newly enrolled independent individuals, age: 40±13 versus 34±10 years and men: 56% versus 50%, respectively) was collected for total RNA extraction using TRIzol Reagent (Invitrogen) following the manufacturer’s protocol. The quantification of genes expression of PBMC was tested by RNA microarrays, which were performed by Compass Biotechnology Co., Ltd. (Beijing, China) using PrimeView Human Gene Expression Array.
Statistical Analyses
We adopted joint analysis for this two-stage GWAS (Supplemental Figure 1). For single-variant association, logistic regression was applied to test in both the discovery and replication stages, assuming an additive model. The first two principal components were included as covariates. Meta-analysis was conducted using METAL (https://genome.sph.umich.edu/wiki/METAL_Documentation) by default “SAMPLESIZE” analysis scheme with P value and direction of effect weighted according to sample size. The widely accepted thresholds P<5×10−8 and P<5×10−5 were used for genome-wide significance and suggestive significance, respectively. We adopted a nominal P<0.005 in the discovery stage and P<0.05 in the replication stage with the same direction of effect, as reported in the Chinese (14). The variance explained by the identified SNPs was estimated by GCTA software (https://cnsgenomics.com/software/gcta). For gene-based association, low-frequency variants (defined as MAF≤0.05) were analyzed for their potential differential clustering between patients and controls within certain genes (i.e., gene burden analysis) using a sequence kernel association test (SKAT-O) according to the EPACTS-3.3.0 tool (http://www.sph.umich.edu/csg/kang/epacts/) with the same significance threshold of 5×10−3. As normal distributions for the relative gene expressions of the candidates, we adopted unpaired t tests in checking the group difference.
Results
Summary of Exome Chip Association Study
In the discovery stage, after stringent quality control as introduced in Materials and Methods, we evaluated the association of 33,655 autosomal biallelic variants in 601 patients and 4076 controls. There was no substantial population structure in the qualified subjects for the exome chip association study (Supplemental Figure 2). The genomic inflation factor (λ) for single-variant associations in the discovery stage was 1.09. After removing SNPs within the MHC region, λ decreased to 1.05. The Manhattan plot and the quantile-quantile plot of the observed association P values are shown in Supplemental Figure 3.
Only one HLA variant surpassed the genome-wide significance of 5×10−8. For the extended HLA region (chr6: 25–34 Mb) in the exome-wide data, among 2018 variants included in the array 67 HLA SNPs associated with IgA nephropathy with P<0.005. The intronic variant rs29243 in GABBR1 showed the most significance (P=1.82 × 10−8). We observed 144 non-HLA SNPs with a possible association (P<0.005) in the array (Supplemental Table 2).
We checked the 28 reported associated variants associated with IgA nephropathy (Supplemental Table 3). In spite of limited data of exome array after quality control, ten variants could be retrieved. We successfully replicated associations of rs9275596 (HLA-DQA1-DQB1-DRB1, P=0.02) and rs2412971 (HORMAD2, P=0.001) with IgA nephropathy. By further checking the reported loci separately, rs17366568 (ST6GAL1, P=0.03), rs75858201 (ODF1-KLF10-UBR5, P=0.03), rs33952257 (ACCS, P=0.03), rs9901675 (TNFSF13, P=0.03), and rs2412973 (HORMAD2, P<0.001) showed P<0.05. Because these associated variants were the same or in LD with reported GWAS lead SNPs (Supplemental Figure 4, Supplemental Table 3), it suggested that this exome array could replicate associations for six reported loci.
Associations of Histocompatibility Leukocyte Antigen Variants
HLA (GABBR1) was the only region associated with IgA nephropathy in the discovery exome chip (Figure 1). With genotyped data, after stepwise conditional analysis, only missense variant rs12722039 in HLA-DQA1 showed additional independent association (P=1.85×10−5) with P threshold <5×10−5. Because the top association was somewhat different compared with previous report, we conducted HLA imputation and checked the contribution of specific HLA alleles and amino acid residues to the overall risk of IgA nephropathy. After quality controls, 5722 variants had high imputation quality (r2>0.8). Single-variant rs16895058 (3′ UTR of UBD), which was in high LD (r2=1.0) with rs29243 in GABBR1, showed the most significant association (P=6.22×10−9). However, no classic alleles or amino acid residue surpassed genome-wide significance in this exome-wide variant panel. HLA-DQB1*05:03 showed the most significance (P=3.27×10−4) in classic class 2 HLA alleles. DRB1–86 showed the most significance (P=8.88 × 10−5) in amino acid positions.
HLA gene variants associated with IgA nephropathy. The regional plots for HLA in single-variant associations are shown. The data for HLA region (chr6: 25–34 Mb) before imputation are shown and are plotted by LocusZoom. rs29243 in GABBR1 showed the most significance (P<5×10−8), and a second signal was from HLA-DQA1. SNP, single-nucleotide polymorphism.
Nonhistocompatibility Leukocyte Antigen Variants with Additional Replication Data and Meta-Analysis
Data for 133 non-HLA variants of the 144 could be retrieved in an independent cohort of 8565 samples. We found 35 SNPs exhibiting a significant association in the replication stage (P<0.05) (Supplemental Table 4).
Among these associated variants, two SNPs located in HORMAD2, which has been consistently associated with IgA nephropathy in previous GWAS (9⇓–11), showed the genome-wide significance (rs2412971, Pmeta=1.19×10−11; rs2412973, Pmeta=1.56×10−11). They were in high LD (r2=0.96) and were the lead SNPs in reported GWAS (9,11).
We found five loci with suggestive associations defined with P<0.001, namely, C5orf55 (chr5, rs139173233, Pmeta=4.72×10−5), FBXL21 (chr5, rs40986, Pmeta=1.13×10−5), CCR6 (chr6, rs3093024, Pmeta=1.39×10−6; rs3093023, Pmeta=1.05×10−6), NAALAD2 (chr11, rs11018879, Pmeta=3.02×10−5), and STAT3 (chr17, rs1026916, Pmeta= 1.75×10−5; rs744166, Pmeta=2.41×10−5)-HSD17B1 (chr17, rs605059, Pmeta =5.74×10−6)-COASY (chr17, rs615942, Pmeta=1.11×10−5)-PSMC3IP (chr17, rs2292751, Pmeta=1.33×10−5) (Figure 2).
Five novel non-HLA loci associated with IgA nephropathy with suggestive significance. The regional plots for novel non-HLA loci in single-variant associations are shown. The visualization of association results was plotted by LocusZoom (http://locuszoom.sph.umich.edu/). It shows the P values obtained in the exome-wide discovery samples. The –log10 P values (y axes) of SNPs are presented against their chromosomal positions (x axes). SNPs highlighted are the validated SNPs with P value from combined analysis of discovery and validation samples, and their linkage disequilibrium (LD) (r2) values with nearby SNPs are indicated by different colors. The recombination rate (centimorgans per megabase; estimated using the 1000 Genomes November 2014 American Society of Nephrology samples) is shown with blue lines, and the genes within the region are shown as black arrows. All coordinates are on the basis of National Center for Biotechnology Information build 37.
For the CCR6 and STAT3-HSD17B1-COASY-PSMC3IP gene regions, as multiple variants have been identified to be associated with IgA nephropathy, we further checked the LD among them. In both regions, the associated variants were in high LD (r2=0.99–1.00). In GWAS reported as associated with other autoimmune diseases (the NHGRI-EBI catalog of published GWAS, https://www.ebi.ac.uk/gwas/), rs3093023 and rs744166 have been shown as the lead SNPs, with CCR6 and STAT3 as the candidates. We thus further focused on these two variants.
In the FBXL21, CCR6, and STAT3 gene regions (up- and downstream 400 kb of the gene), apart from the above associated SNPs, there were additional variants in partial LD showing significant associations (P<0.005), further supporting their genetic associations.
Common Variant Functional Annotation
SNPs including rs40986, rs3093023, and rs744166 were further checked for functional annotation. rs40986 is a missense variant (p.Pro209Leu) of FBXL21, and rs3093023 and rs744166 are intronic variants of CCR6 and STAT3, respectively.
In the eQTLGen Consortium data (an eQTL meta-analysis performed in blood samples from 31,684 individuals) (21,22), the associated SNPs had significant cis-eQTL effects: for example, rs40986 on TGFBI (P=1.09×10−19), rs3093023 on CCR6 (P=3.24×10−145), and rs744166 on STAT3 (P=6.91×10−59). By HaploReg v4.1 (Supplemental Table 5), rs3093023 and rs744166 were observed to be located in enhancer regions enriched in histone marks H3K4me1 in primary B cells. As enhancers are distal cis-regulatory elements that modulate gene expression, it may partially account for the underlying mechanism of the eQTL association.
Gene-Based Group of Low-Frequency Variants
It has been suggested that a gene-based test provides better power by aggregating low-frequency variants compared with single-variant analyses (23,24). To test for accumulation of low-frequency variants within a gene, we performed a gene burden analysis with SKAT-O (25). Differences between patients and controls did not survive correction for multiple testing for 8255 genes analyzed (P< 6.06×10−6). Detailed information of the SKAT-O gene burden analysis is provided in Supplemental Table 6. However, of interest, as a key factor that triggers alternative pathways of complement, CFB showed the lowest P value (3.98× 10−4).
Differential Gene Expression of Novel Genes in Immunoglobulin A Nephropathy Kidney Tissue and Blood
For the genes identified above, we interrogated the open-access Nephroseq Classic (v4) database (www.nephroseq.org), which integrates a total of 36 datasets from 2015 samples, including two studies from human IgA nephropathy kidneys microdissected into glomerular and tubulointerstitial compartments. As shown in Supplemental Figure 5 and in both datasets, TGFBI and CFB showed consistent significantly higher expression in the glomeruli and tubulointerstitium of patients with IgA nephropathy. Compared with healthy living donors, TFGBI mRNA expression was at least 2.7- and 1.4-fold higher in glomeruli and tubulointerstitium, respectively (P<0.05). CFB transcripts were at least 1.5- and 2.3-fold higher in glomeruli and tubulointerstitium, respectively (P<0.05). Glomerular STAT3 and tubulointerstitial CCR6 mRNAs showed higher expressions in patients with IgA nephropathy in one of the two studies. One study suggested that GABBR1 was marginally expressed higher in IgA nephropathy glomeruli.
As these associated genes exert important functions in immunity regulation, to further check the gene expressions in blood, we enrolled a larger sample size (59 patients and 28 controls) compared with previous reports. As can be seen in Figure 3, TGFBI, STAT3, and GABBR1 were significantly upregulated in PBMCs from patients with IgA nephropathy (P<0.05). Differential mRNA levels were confirmed by detecting two or more probes.
TGFBI, STAT3, and GABBR1 showed differential expressions from PBMC of patients with IgA nephropathy (IgAN) compared with controls. Among the five novel genes observed, only showed the positive findings of differential gene expressions confirmed by two or more probes in microarray. The scatterplots are presented with mean ± SD and P values. (A) TGFBI was detected by three probes, 11734549_s_at (IgAN versus controls, 8.85±0.42 versus 8.63±0.51, P=0.03), 11755281_a_at (9.01±0.46 versus 8.73±0.88, P=0.05), and 11734550_x_at (9.37±0.40 versus 9.15±0.59, P=0.05). (B) STAT3 was detected by three probes, 11720755_at (IgAN versus controls, 8.09±0.22 versus 7.89±0.17, P<0.001), 11720756_at (8.39±0.27 versus 8.09 ±0.47, P<0.001), and 11750674_a_at (8.19±0.27 versus 7.93±0.42, P<0.001). (C) GABBR1 was detected by two probes, 11730725_a_at (IgAN versus controls, 6.19±0.32 versus 6.02±0.45, P=0.04) and 11721553_a_at (7.31±0.28 versus 7.19±0.31, P=0.07).
Gene Network Analysis
The above blood eQTL analysis and differential gene expression data suggested that TGFBI rather than FBXL21 might be the effector gene. We thus further investigated TGFBI as a candidate. Using the Protein-Protein Interaction Network tool STRING v11.0 (https://string-db.org/), the novel suggestive genes identified involved in IgA nephropathy, including CCR6, STAT3, TGFBI, GABBR1, and CFB, were found connected with known IgA nephropathy genes (CFH, VAV3, ST6GAL1, HLA-DQA1, DEFA1, KLF10, CARD9, ACCS, ITGAM, TNFSF13, and HORMAD2) with a protein-protein interaction enrichment P=1.13×10−5 (Figure 3); they jointly suggested several pathways relevant to immune response (Supplemental Table 7). The most significant biologic process was regulation of immune system process (GO:0002682; PFDR=3.76×10−15), and the most significant KEGG pathway was Staphylococcus aureus infection (hsa05150, PFDR=1.95×10−12).
Lastly, we evaluated the contribution of the associated SNPs to overall IgA nephropathy risk. Cumulatively, by GCTA software, these novel association signals explained about 1.67%±1.21% (P =1.67×10−16) of disease variance (Table 2).
Meta-analysis results of the three novel non-HLA IgA nephropathy susceptibility loci
Discussion
In this study, we identified five suggestive susceptible genes for IgA nephropathy, namely, CCR6, STAT3, TGFBI, GABBR1, and CFB, by a multitiered exome-wide study in 13,242 Han Chinese samples. We suggest that the novel implicated genes might modulate the immune response and share a similar pathway with infections in IgA nephropathy. To the best of our knowledge, this is the first large-scale exome array study for IgA nephropathy, which provides insights into the biologic mechanism of IgA nephropathy.
Although our understanding of the genetic susceptibility of IgA nephropathy has been advanced dramatically recently, the total explained heritability is still limited. Importantly, elucidating the genetic mechanism of IgA nephropathy underlying each susceptibility signal is still a great challenge, which can partially be attributed to the property of GWAS SNPs within nonprotein-coding or intergenic regions (26,27). Exome-wide association studies have been validated to be a powerful approach to identify disease-causing genes (28). In this study, we first confirmed the strong genetic associations between variants in HLA and HORMAD2 in IgA nephropathy, both of which passed genome-wide significance, highlighting their etiologic significance (10,11). More importantly, we identified five genes with suggestive evidence for an association with IgA nephropathy. These associated variants were annotated as functional of expression-related SNPs, and SNP-containing genomic regions may be possible exonic or intronic transcriptional regulatory elements affecting gene expression in B cells. This possibility aids prioritization of their effector genes. The above identified loci have not been discovered in previous GWAS in IgA nephropathy, and the overall evidence did not reach genome-wide significance for the newly identified gene loci. However, in our multistage genetic study, combinations of clues from integrative statistical, bioinformatics, expression, and pathway analyses support the causality of the genes and their biologic relevance in the development of IgA nephropathy. Furthermore, these important findings are concordant with some more recent statistical and functional evidence supporting the involvement of these genes in IgA nephropathy (29⇓–31) (Table 3).
Novel IgA nephropathy loci and gene functional prioritization
CFB, GABBR1, TGFBI, and CCR6 are all involved in complement regulation pathways, which have been considered a dominant driver of kidney injury in IgA nephropathy (32), whereas CCR6 signaling may induce the JAK-STAT3 pathway in T cells, which has been suggested to play an important role in Th17 signaling and IgA1 production (29) (Figure 4). Using network analysis, CCR6, STAT3, TGFBI, and CFB are involved in two pathways, namely, response to stimulus (GO:0050896, P<0.001) and regulation of biologic process (GO:0050789, P=0.002) (Supplemental Table 7). Integrated with known IgA nephropathy loci, the top pathways were all involved in immune regulation and defense response, supporting the idea that initiating pathogenic insult in IgA nephropathy must arise from aberrant immunity (9). Five of the top ten KEGG pathways are consistent with infections (S. aureus infection, tuberculosis, leishmaniasis, herpes simplex infection, and toxoplasmosis). Of particular interest, all these infections have been reported to be associated with IgA nephropathy. The top suggested Staphylococcus is also by far the most frequent etiologic agent as IgA-dominant infection-associated GN (33). The remaining top KEGG pathways include known IgA nephropathy pathways: intestinal immune network for IgA production and inflammatory bowel disease. Furthermore, our novel data in combination with previous reports support the likely biologic significance for the Jak-STAT signaling pathway and the Th17 cell differentiation in IgA nephropathy. These pathways are interconnected by these associated genes, and the novel ones warrant more attention.
Interactive network analysis of the proteins encoded by novel and known IgA nephropathy genome-wide association patient-control study (GWAS)–implicated genes. The functional interaction network in this figure was constructed for the genes/loci identified in previous GWAS and this study using STRING v11.0. Inputted the newly identified genes in this study and the 11 closest or reported genes implicated in previous IgA nephropathy GWAS and found that all of the new genes interacted with known IgA nephropathy genes. The five novel genes were highlighted with circles.
Limitations of the study should also be noted. Because this is an exome-wide study, the number of qualified variants is relatively small (33,655 autosomal biallelic variants) and not optimal in imputation, thus limiting the discovery ability. Further studies in large independent samples will be needed to replicate these findings. The discovered variants may act as susceptibility biomarkers rather than causal variants. Future studies involving fine mapping of targeted regions and deep functional studies are required to delineate the molecular mechanisms underlying the risk variants identified.
In summary, we performed the first large-scale exome-wide study for IgA nephropathy in a Han Chinese population and identified five novel suggestive genes for IgA nephropathy. The findings provide novel clues into the biologic mechanisms of IgA nephropathy.
Disclosures
All authors have nothing to disclose.
Funding
Support was provided by Beijing Natural Science Foundation grant Z190023; Beijing Nova Program Interdisciplinary Cooperation Project grant Z191100001119004; Beijing Youth Top-notch Talent Support Program grant 2017000021223ZK31; Chinese Academy of Medical Sciences Innovation Fund for Medical Sciences grant 2019-I2M-5-046; Clinical Medicine Plus X-Young Scholars Project of Peking University grant PKU2020LCXQ003; Fok Ying Tung Education Foundation grant 171030; National Natural Science Foundation of China grants 82022010, 81970613, and 82070733; Natural Science Foundation for Innovation Research Group of China grant 81621092; the Training Program of the Major Research Plan of the National Natural Science Foundation of China grant 91642120; and the University of Michigan Health System–Peking University Health Science Center Joint Institute for Translational and Clinical Research grant BMU2017JI007.
Acknowledgements
The authors thank their collaborators, Dr. Ali G. Gharavi, Dr. Krzysztof Kiryluk, and Dr. Elena Sanchez-Rodrıguez (Division of Nephrology, Department of Medicine, Columbia University, New York), for providing previous GWAS data, assisting statistical analysis, and giving advice for the manuscript. The authors thank Dr. Kirsten Herold for elaborative language editing on our manuscript and the Human Genetic Resources Core Facility and Center for Life Sciences of Peking University for enrolling controls. They also thank the patients, their families, and healthy donors for their cooperation and for giving consent to participate in this study. The authors apologize to the scientists whose important works are not cited here due to space limitation.
The funders had no role in study design, data collection and analysis, or preparation of the manuscript.
H. Zhang and X.-j. Zhou designed the study; Y. Hu, L. C. Tsoi, and X.-j. Zhou carried out experiments and data analysis; T. Gan, P. Hou, Yang Li, Yanming Li, L.-j. Liu, J.-c. Lv, Y.-y. Qi, S.-f. Shi, Y.-n. Wang, H.-j. Xu, and Y.-m. Zhang contributed to sample collection and data analysis; C. C. Berthier, K. He, Yanming Li, M. T. Patrick, H. Zhang, and X.-j. Zhou drafted and revised the paper; and all authors approved the final version of the manuscript.
Supplemental Material
This article contains the following supplemental material online at http://cjasn.asnjournals.org/lookup/suppl/doi:10.2215/CJN.06910520/-/DCSupplemental.
Supplemental Figure 1. Flow chart of the study design of this study.
Supplemental Figure 2. Principal components analysis plot in the exome array data.
Supplemental Figure 3. Manhattan plot and quantile-quantile plot in the discovery stage.
Supplemental Figure 4. The regional plots for GWAS-reported loci for IgA nephropathy in single-variant associations.
Supplemental Figure 5. Differential gene expression for novel loci using Nephroseq.
Supplemental Table 1. Some clinical parameters for patients with IgA nephropathy in the exome array.
Supplemental Table 2. Summary statistics of the associated genetic variants with P<0.005 in the discovery samples.
Supplemental Table 3. Replication data for reported variants/loci in the exome array.
Supplemental Table 4. Summary statistics of the associated genetic variants with P<0.05 in the validation samples.
Supplemental Table 5. Functional annotation for four SNPs within suggestive genes.
Supplemental Table 6. Detailed information of the SKAT-O gene burden analysis (data shown for significant findings).
Supplemental Table 7. Gene ontology and KEGG pathways implicated by newly implicated genes and known IgA nephropathy risk genes using STRING (http://string-db.org/).
Footnotes
Published online ahead of print. Publication date available at www.cjasn.org.
See related editorial, “Assessing Genetic Risk for IgA Nephropathy: State of the Art,” on pages 182–184.
- Received May 7, 2020.
- Accepted December 11, 2020.
- Copyright © 2021 by the American Society of Nephrology
References
In this issue
