https://orcid.org/0000-0003-3067-0139
Abstract
Genome-wide association studies (GWAS) of Crohn’s disease (CD) in European and leprosy in Chinese population have shown that CD and leprosy share genetic risk loci. As these shared loci were identified through cross-comparisons across different ethnic populations, we hypothesized that meta-analysis of GWAS on CD and leprosy in East Asian populations would increase power to identify additional shared loci. We performed a cross-disease meta-analysis of GWAS data from CD (1621 cases and 4419 controls) and leprosy (2901 cases 3801 controls) followed by replication in additional datasets comprising 738 CD cases and 488 controls and 842 leprosy cases and 925 controls. We identified one novel locus at 7p22.3, rs77992257 in intron 2 of ADAP1, shared between CD and leprosy with genome-wide significance (P = 3.80 × 10−11) and confirmed 10 previously established loci in both diseases: IL23R, IL18RAP, IL12B, RIPK2, TNFSF15, ZNF365-EGR2, CCDC88B, LACC1, IL27, NOD2. Phenotype variance explained by the polygenic risk scores derived from Chinese leprosy data explained up to 5.28% of variance of Korean CD, supporting similar genetic structures between the two diseases. Although CD and leprosy shared a substantial number of genetic susceptibility loci in East Asians, the majority of shared susceptibility loci showed allelic effects in the opposite direction. Investigation of the genetic correlation using cross-trait linkage disequilibrium score regression also showed a negative genetic correlation between CD and leprosy (rg [SE] = −0.40[0.13], P = 2.6 × 10−3). These observations implicate the possibility that CD might be caused by hyper-sensitive reactions toward pathogenic stimuli.
Introduction
Crohn’s disease (CD) is one of the two major forms of chronic inflammatory bowel disease (IBD) characterized by recurring episodes of inflammation of gastrointestinal tract. It is thought to arise by dysregulated mucosal immune responses to the gut flora in genetically susceptible individuals (1). Family and twin studies showed that a positive family history is an important risk factor in both Korea and Western countries (2–4). Although leprosy is a chronic infectious disease caused by Mycobacterium leprae, the role of host genetic factors has been well established through epidemiological and molecular genetic studies (5,6). Dozens of genes and loci have been discovered and replicated though familial clustering, studies of twins, complex segregation analysis as well as genome-wide association study (GWAS) (7–10).
GWASs have suggested that CD and leprosy might share a common underlying genetic susceptibility. Based on GWAS of a Chinese population with leprosy (9–13), replication of European CD or ulcerative colitis (UC) susceptibility loci in leprosy samples of Chinese origin (14), and CD GWAS of European populations (15,16), 10 loci (IL23R, IL18RAP, IL12B, RIPK2, TNFSF15, ZNF365-EGR2, CCDC88B, LACC1, SOCS1-CIITA and NOD2) are established as shared between CD and leprosy with genome-wide significance. Of these 10 loci, four (IL23R, ZNF365-EGR2, SOCS1-CIITA and NOD2) are reported to have independent signals for each disease (10). As these shared loci were emerged from simple cross-comparisons, focusing specifically on genome-wide significant loci across published studies of different ethnic populations, a substantial proportion of the genetic architectures underlying these complex diseases could have been missed. There have been no genome-wide scale comparisons between CD and leprosy susceptibilities among Asian populations to date. We hypothesized that systemic analysis of an integrated dataset of East Asian populations would increase the power to identify shared susceptibility loci. We first performed a meta-analysis considering both diseases as a single phenotype by combining previously published GWAS data for each disease under the assumption of the same or opposite allelic effects in the two diseases. In addition, we used a multi-trait meta-analysis approach, PLEIO (17), which accounts for differing heritabilities of the diseases as well as the genetic correlation between them to maximize power. Using these approaches, we identified one novel shared susceptibility locus and confirmed the 10 established loci in CD and leprosy, which point to response to interleukin-18 and regulation of T-helper 1 type immune response as key shared pathways involved in disease pathogenesis.
Results
Cross-disease meta-analysis and replication
The datasets used in this study and the strategy for identifying shared susceptibility loci are presented in Figure 1. We took two different approaches to select candidate shared loci, while minimizing type I error: simple cross-disease meta-analysis and multi-trait joint analysis using PLEIO (17). The former, the widely used standard meta-analysis, is optimized for finding loci with similar effect size magnitude for two diseases. The latter, a recently developed a variance component method, is optimized for finding loci with effect sizes that conform to the heritability and genetic correlation pattern observed from the whole genome. For this reason, PLEIO may find additional loci unidentified by the standard meta-analysis method. Thus, the two approaches may complement each other and may find a different set of significant loci.
Study design to identify susceptibility loci shared between CD and leprosy.
First, we performed a meta-analysis considering both diseases as a single phenotype. Disease-specific fixed-effects meta-analyses of CD (cohorts I and II, comprising 1621 cases and 4419 controls) and leprosy (cohorts IV and V, comprising 2118 cases and 2822 controls) GWAS datasets identified a total of 27 034 overlapping SNPs with disease-specific P < 0.05 in the CD and leprosy GWAS datasets, which was significantly enriched more than expected by chance (18 968 SNPs). To identify shared susceptibility loci with effects in the same or opposite directions, simple cross-disease fixed-effects meta-analysis or the same analysis with flipped effect directions for all SNPs in the CD GWAS datasets were performed. Applying a threshold of Pmeta < 1 × 10−6 in the cross-disease meta-analyses, with disease-specific P < 0.05, a total of 19 loci, including 10 loci established in both diseases (IL23R, IL18RAP, IL12B, RIPK2, TNFSF15, ZNF365-EGR2, CCDC88B, LACC1, IL27 and NOD2), 2 established CD loci (TBC1D1 and UBE2L3) and 7 novel candidate shared loci (ADAP1, PDGFRL, CTSC, MIR656, CSK, LOC105371070 and NF2), were identified (Supplementary Material, Table S1). Of the 10 established loci, only RIPK2 and LACC1 showed the same effect directions, while the remaining 8 loci showed opposite effect directions in CD and leprosy. Second, we used a recently developed multi-trait joint tool, PLEIO (17), to identify shared susceptibility loci. Input data for PLEIO (17) were the summary statistics from disease-specific fixed-effects meta-analyses of CD and leprosy GWAS datasets. Applying a threshold of Pmeta < 1 × 10−6 in the joint analysis with disease-specific P < 0.05, we identified a total of 10 loci including 8 loci established in both diseases (IL23R, IL12B, RIPK2, TNFSF15, ZNF365-EGR2, LACC1, IL27 and NOD2) and 2 novel candidate shared loci (TBC1D1 and ADAP1) between CD and leprosy (Supplementary Material, Table S2). We focused on the 2 loci (TBC1D1 and ADAP1) that were identified using two different approaches. TBC1D1 is an established CD locus, whereas ADAP1 has not been previously confirmed in either disease. We performed replication analyses of the 2 candidate loci for both CD and leprosy in independent datasets: cohort III for CD and cohort VI for leprosy (Supplementary Material, Table S3). Replication of rs77992257 at 7p22 (ADAP1) in leprosy cohort VI showed association P < 0.05 (Supplementary Material, Table S3). As rs77992257 was not present in CD cohort III data, rs113540766 was selected as a proxy SNP (r2 = 0.71 with rs77992257 in EAS). Furthermore, a TaqMan genotyping assay of rs113540766 in an independent cohort of 613 CD cases and 1083 controls showed a significant association (P = 5.73 × 10−3). Of the 2 candidate loci, only rs77992257 reached genome-wide significance in combined analysis of the CD and leprosy discovery and replication cohorts (Pmeta = 3.80 × 10−11, PLEIO P = 7.59 × 10−10) (Table 1), totaling 11 shared loci associated with both CD and leprosy (Supplementary Material, Table S4). Notably, ADAP1, a novel shared locus, showed opposite allelic effects in CD and leprosy.
The association of rs77992257 at shared susceptibility locus to Crohn’s disease and leprosy
| MAF | Fixed-effects meta-analysis | Cross-disease fixed-effects meta-analysisa | Multi-trait joint analysisb | |||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Locus | Candidate gene | SNP | Position (hg19) | Phenotype | Dataset | Minor allele | Case | Control | OR | P c | OR | P | Phet | P | Phet | P |
| 7p22 | ADAP1 | rs77992257 | 989,986 | CD | Cohort I | A | 0.252 | 0.303 | 0.77 | 8.18 × 10−6 | 0.79 | 9.74 × 10−8 | 6.02 × 10−1 | 3.80 × 10−11 | 8.93 × 10−1 | 7.59 × 10−10 |
| Cohort II | 0.258 | 0.284 | 0.86 | 1.67 × 10−1 | ||||||||||||
| TaqMand | 0.292 | 0.337 | 0.81 | 5.73 × 10−3 | ||||||||||||
| Leprosy | Cohort IV | A | 0.287 | 0.250 | 1.15 | 1.82 × 10−1 | 1.21 | 7.33 × 10−5 | 8.96 × 10−1 | |||||||
| Cohort V | 0.260 | 0.222 | 1.22 | 1.76 × 10−3 | ||||||||||||
| Cohort VI | 0.247 | 0.229 | 1.21 | 3.64 × 10−2 | ||||||||||||
| MAF | Fixed-effects meta-analysis | Cross-disease fixed-effects meta-analysisa | Multi-trait joint analysisb | |||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Locus | Candidate gene | SNP | Position (hg19) | Phenotype | Dataset | Minor allele | Case | Control | OR | P c | OR | P | Phet | P | Phet | P |
| 7p22 | ADAP1 | rs77992257 | 989,986 | CD | Cohort I | A | 0.252 | 0.303 | 0.77 | 8.18 × 10−6 | 0.79 | 9.74 × 10−8 | 6.02 × 10−1 | 3.80 × 10−11 | 8.93 × 10−1 | 7.59 × 10−10 |
| Cohort II | 0.258 | 0.284 | 0.86 | 1.67 × 10−1 | ||||||||||||
| TaqMand | 0.292 | 0.337 | 0.81 | 5.73 × 10−3 | ||||||||||||
| Leprosy | Cohort IV | A | 0.287 | 0.250 | 1.15 | 1.82 × 10−1 | 1.21 | 7.33 × 10−5 | 8.96 × 10−1 | |||||||
| Cohort V | 0.260 | 0.222 | 1.22 | 1.76 × 10−3 | ||||||||||||
| Cohort VI | 0.247 | 0.229 | 1.21 | 3.64 × 10−2 | ||||||||||||
CD, Crohn’s disease; hg19, human genome version 19; MAF, minor allele frequency; OR, odds ratio; P, P-value; Phet, P-value for heterogeneity; Position, chromosome position; SNP, single nucleotide polymorphism.
aA meta-analysis software tool based on fixed-effects model (41).
bA multi-trait joint analysis tool maximizes power by systematically accounting for genetic correlations and heritability of the traits using summary statistics (17).
cP-value was calculated using frequentist association test of SNPTEST (40).
dTaqMan genotyping assay of a proxy SNP (rs113540766, r2 = 0.71 with rs77992257 in East Asians) in an independent cohort of 613 CD and 1083 controls.
The association of rs77992257 at shared susceptibility locus to Crohn’s disease and leprosy
| MAF | Fixed-effects meta-analysis | Cross-disease fixed-effects meta-analysisa | Multi-trait joint analysisb | |||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Locus | Candidate gene | SNP | Position (hg19) | Phenotype | Dataset | Minor allele | Case | Control | OR | P c | OR | P | Phet | P | Phet | P |
| 7p22 | ADAP1 | rs77992257 | 989,986 | CD | Cohort I | A | 0.252 | 0.303 | 0.77 | 8.18 × 10−6 | 0.79 | 9.74 × 10−8 | 6.02 × 10−1 | 3.80 × 10−11 | 8.93 × 10−1 | 7.59 × 10−10 |
| Cohort II | 0.258 | 0.284 | 0.86 | 1.67 × 10−1 | ||||||||||||
| TaqMand | 0.292 | 0.337 | 0.81 | 5.73 × 10−3 | ||||||||||||
| Leprosy | Cohort IV | A | 0.287 | 0.250 | 1.15 | 1.82 × 10−1 | 1.21 | 7.33 × 10−5 | 8.96 × 10−1 | |||||||
| Cohort V | 0.260 | 0.222 | 1.22 | 1.76 × 10−3 | ||||||||||||
| Cohort VI | 0.247 | 0.229 | 1.21 | 3.64 × 10−2 | ||||||||||||
| MAF | Fixed-effects meta-analysis | Cross-disease fixed-effects meta-analysisa | Multi-trait joint analysisb | |||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Locus | Candidate gene | SNP | Position (hg19) | Phenotype | Dataset | Minor allele | Case | Control | OR | P c | OR | P | Phet | P | Phet | P |
| 7p22 | ADAP1 | rs77992257 | 989,986 | CD | Cohort I | A | 0.252 | 0.303 | 0.77 | 8.18 × 10−6 | 0.79 | 9.74 × 10−8 | 6.02 × 10−1 | 3.80 × 10−11 | 8.93 × 10−1 | 7.59 × 10−10 |
| Cohort II | 0.258 | 0.284 | 0.86 | 1.67 × 10−1 | ||||||||||||
| TaqMand | 0.292 | 0.337 | 0.81 | 5.73 × 10−3 | ||||||||||||
| Leprosy | Cohort IV | A | 0.287 | 0.250 | 1.15 | 1.82 × 10−1 | 1.21 | 7.33 × 10−5 | 8.96 × 10−1 | |||||||
| Cohort V | 0.260 | 0.222 | 1.22 | 1.76 × 10−3 | ||||||||||||
| Cohort VI | 0.247 | 0.229 | 1.21 | 3.64 × 10−2 | ||||||||||||
CD, Crohn’s disease; hg19, human genome version 19; MAF, minor allele frequency; OR, odds ratio; P, P-value; Phet, P-value for heterogeneity; Position, chromosome position; SNP, single nucleotide polymorphism.
aA meta-analysis software tool based on fixed-effects model (41).
bA multi-trait joint analysis tool maximizes power by systematically accounting for genetic correlations and heritability of the traits using summary statistics (17).
cP-value was calculated using frequentist association test of SNPTEST (40).
dTaqMan genotyping assay of a proxy SNP (rs113540766, r2 = 0.71 with rs77992257 in East Asians) in an independent cohort of 613 CD and 1083 controls.
A newly identified shared susceptibility locus for CD and leprosy: 7p22
At 7p22, the lead SNP was rs79944617 (P = 1.24 × 10−6) for CD and rs77992257 (P = 7.33 × 10−5) for leprosy (Fig. 2A and B). Both rs79944617 and rs77992257 are in intron 2 of ADAP1 (ArfGAP with dual PH domains 1), 4.1 kb apart and in high linkage disequilibrium (LD; r2 = 0.74, the 1000 Genomes East Asian LD reference [JPT + CHB]). The CD and leprosy meta-analysis identified rs77992257 as a shared lead SNP at 7p22. A regional association plot of leprosy data showed an additional signal at rs60609206 (P = 3.36 × 10−5), 9.6 kb upstream of ZFAND2A (zinc finger AN1-type containing 2A) and 219.5 kb away from rs77992257 with low LD (r2 < 0.2) (Fig. 2C). Following analysis conditioned on rs77992257 using GCTA-COJO analysis (18), the rs60609206 signal was not unaffected (P = 5.26 × 10−4, data not shown); however, in combined of analysis of CD, rs60609206 was not significantly associated (P = 0.26), suggesting that rs60609206 may be leprosy-specific. To confirm that rs60609206 is not an additional signal shared between CD and leprosy at the ADAP1 locus, we performed a colocalization analysis using coloc v3.2–1 (https://github.com/cran/coloc). The colocalization analysis using both summary statistics of CD and leprosy after conditioning on rs77992257 showed a PP4 = 0.017, suggesting that rs77992257 was the only signal shared between CD and leprosy at the ADAP1 locus.
Analysis of Genotype-Tissue Expression (GTEx) (19) data showed that rs77992257 acts as an expression quantitative trait locus (eQTL) for ADAP1 in colon tissue (Supplementary Material, Table S5). The minor allele (G) of rs77992257 (the risk allele for CD and protective allele for leprosy) was associated with higher ADAP1 expression levels than allele A. Search of a blood eQTL database derived from Korean patients with CD (20) also showed that rs73046475 (r2 = 0.69 with rs77992257) was a significant eQTL for ADAP1 (P = 1.91 × 10−7). As a binding protein of the lipid second messenger, phosphatidylinositol-3,4,5-trisphosphate (PtdInsP3) (21,22), ADAP1 contributes to the development and function of regulatory T cells (23) and to the internalization of many microbial pathogens, including Salmonella enterica (24) and Yersinia enterocolitica, into host cells (25).
Colocalization analysis
To check whether GWAS of CD and leprosy showed the same causal SNPs at the 11 shared loci, we used coloc v3.2–1 (https://github.com/cran/coloc), Bayesian procedure providing intuitive posterior probabilities. The colocalization analysis showed that a total of 5 loci (TNFSF15, IL12B, RIPK2, IL27 and ADAP1) were over the threshold of posterior probability 4 (PP4) > 0.8, which means strong evidence of a causal SNP associated with both traits (Supplementary Material, Table S6). NOD2 and IL23R showed 0.8 > PP4 > 0.7; however, the remaining 4 loci including LACC1, IL18RAP, ZNF365 and CCDC88B failed to show supportive evidence of colocalization (PP4 < 0.5) (Supplementary Material, Table S6).
Cis-eQTL analysis
We used the most recent GTEx database v.8 (19), and whole blood cis-eQTL databases for eQTLGen (26), Koreans (20) and Japanese (27) to explore whether the lead SNP at each shared locus significantly influenced gene expression in disease-relevant tissues. To check whether GWAS of CD or leprosy have the same causal SNPs with cis-eQTL data, we performed colocalization analyses using coloc v3.2–1 based on the cis-eQTL datasets of disease-relevant tissues and summary statistics of disease-specific meta-analyses of CD or leprosy. Two loci including IL12B and ZNF365 did not show significant cis-eQTL effects. In the colocalization analysis of the remaining 9 loci applying a threshold of PP4 > 0.8, 11 genes at 3 loci (TNFSF15, IL27 and NOD2) in CD and 19 genes at 5 loci (RIPK2, TNFSF15, CCDC88B, IL27 and NOD2) in leprosy were over the threshold (Supplementary Material, Table S7). A total of 9 genes (TNFSF15, SULT1A1, SULT1A2, SGF29, NPIPB7, SBK1, NOD2, RP11-327F22.6 and RP11-327F22.1) at 3 loci (TNFSF15, IL27 and NOD2) were overlapped between colocalized 11 genes of CD and 19 genes of leprosy (Supplementary Material, Table S7).
Gene annotation and protein–protein interaction network analysis
To prioritize causal genes at shared susceptibility loci between CD and leprosy, we performed gene annotation analysis using Multi-marker Analysis of GenoMic Annotation (MAGMA) v.1.07b. Using the summary statistics for the cross disease meta-analysis between CD and leprosy, we annotated a total of 24 genes at 9 loci with P < 2.91 × 10−6 (Supplementary Material, Table S8). The remaining 2 loci (IL23R, IL12B) failed to show genes with P < 2.91 × 10−6. To identify biological pathways enriched for the annotated genes, we also conducted pathway enrichment analysis using the Gene Ontology resource (http://geneontology.org/). Applying a threshold of false discovery rate (FDR) < 0.05, two biological pathways, response to interleukin-18 and regulation of T-helper 1 type immune response, were associated with the 24 annotated genes (Supplementary Material, Table S9).
Among the 24 genes annotated by MAGMA, we selected 9 genes (SLC9A4, ADAP1, RIPK2, DEC1, EGR2, CCDC88B, CCDC122, SULT1A1 and NOD2) with the most significant P-value at 9 loci as an input data for protein–protein interaction (PPI) network analysis. To construct a functional association network for the 9 genes, we conducted PPI network analysis using the STRING database (https://string-db.org/) (28). One PPI network involving 5 of 9 proteins was found to show significant network connectivity compared with a random set of proteins of the same size (P = 5.46 × 10−7) (Fig. 3), indicating that the 5 genes are biologically connected.
Regional association plots of the 7p22 locus. (A) rs79944617 in CD GWAS data. (B) rs77992257 in leprosy GWAS data. (C) rs60609206 in leprosy GWAS data. SNPs are plotted according to their chromosomal positions (NCBI Build 37) with −log10 P values from the CD or leprosy meta-analysis in the region flanking 750 kb on either side of the marker SNP. The most strongly associated SNP in the discovery stage is shown as a small purple circle. LD (r2 values) between the lead SNP and other SNPs is indicated using colors. The relative locations of the annotated genes and their directions of transcription are shown in the lower portion of the figure. The estimated recombination rates of Asian samples from the 1000 Genomes Project (November 2014) are plotted to reflect the local LD structure. Plots were generated using LocusZoom.
PPI network visualized using STRING. Of 24 genes identified at the 9 shared loci using MAGMA, we selected 9 genes with the most significant P-value at the 9 shared loci and performed PPI network analysis using the STRING database. The color of lines between proteins indicates the types of PPIs. Significant PPI enrichment P-values demonstrate that proteins in the PPI network have more interactions among themselves than what would be expected for a random set of proteins of the same size.
Extent of sharing between CD and leprosy
To explore the proportion of overlapping genetic effects between CD and leprosy, we estimated the genetic correlation and P-values using LD score regression (LDSC) (29) by selecting 19 496 common SNPs in the region of ±250 kb flanking the 11 lead SNPs in the CD and leprosy GWAS datasets. Correlation of genetic effect direction was negative (rg [SE] = −0.52 [0.14]) and significant (P = 3.0 × 10−4). Estimation of the genetic correlation using all common 4 938 574 SNPs detected between the CD and leprosy GWAS datasets in the discovery stage also showed a significant negative correlation (rg [SE] = −0.40 [0.13]; P = 2.6 × 10−3).
We also examined the extent of the overlap of the genetic architecture between CD and leprosy by estimating the variance of CD explained by the polygenic risk scores (PRS) derived from leprosy GWAS (PRSleprosy). First, we calculated variance explained in CD by PRSleprosy based on eight genome-wide significant variants using leprosy effect sizes (Supplementary Material, Table S10). PRSleprosy explained up to 4.27% of phenotype variance of CD (Supplementary Material, Table S11). Furthermore, when the direction of the effects of the LACC1 and RIPK2 loci was flipped to account for the original directional relationship between CD and leprosy, PRSleprosy explained up to 8.21%. Applying PRSice-2 (30), a high-resolution PRSice plot for leprosy predicting CD showed that the best threshold for the PRS model was 1.71 × 10−5 (Supplementary Material, Fig. S1). The skewed distribution of PRS observed at the best P-value threshold appeared to be due to a mixture of normal distributions. Indeed, PRS showed normal distributions when the samples were stratified based on the genotype of the most significantly associated SNP in the major histocompatibility locus (MHC) (rs9271011) (Supplementary Material, Fig. S2). High-resolution best-fit PRSleprosy explained 5.28% of CD variance and was based on 39 SNPs (Supplementary Material, Table S11). Using a logistic mixed model adjusted for age and sex to test for association between PRSleprosy as a continuous variable and CD, each standard deviation decrease in the PRSleprosy was associated with a 48% relative risk increase in CD (OR = 0.52, 95% confidence interval = 0.44–0.62; P = 2.10 × 10−13). Thus, higher PRSleprosy was thus associated with risk reduction for CD.
Discussion
In this study, we have identified one novel susceptibility locus shared between CD and leprosy by conducting a meta-analysis considering both diseases as a single phenotype and using PLEIO (17), a multi-trait meta-analysis approach. We also confirmed 10 previously reported shared susceptibility loci. The novel locus, ADAP1, had not been previously associated with either disease at genome-wide significance. We also identified a significant negative genetic correlation between CD and leprosy. Our study represents the first systematic effort to compare the genetic basis of CD and leprosy in samples of East Asian origin in order to identify risk alleles with pleiotropic effects on two clinically distinct diseases.
Previous studies examining well-established European IBD susceptibility loci in a leprosy dataset from a Chinese population identified 10 common susceptibility loci (10,13,14), all of which, except SOCS1-CIITA, were confirmed in this study, and the reported directions of the effects for six loci (IL18RAP, IL12B, RIPK2, TNFSF15, CCDC88B and LACC1) were consistent with those determined in the current study. As CD and leprosy have independent signals at the SOCS1-CIITA locus (10), the cross-disease meta-analysis failed to confirm the shared susceptibility signal at SOCS1-CIITA. In this study, we determined the directions of the effects for the three additional loci (IL23R, ZNF365-EGR2 and NOD2), which were previously thought to be independent from European CD SNPs. Additionally, we identified a new shared locus at IL27, a well-established CD and leprosy locus (13), and confirmed that its SNPs had opposite effect directions in CD and leprosy. Despite its somewhat limited sample size, our study identified one novel shared locus between CD and leprosy (the ADAP1 locus), which had not previously been discovered in either disease, probably due to the discoveries being made in similar ethnic populations. This association has not been discovered in larger European GWAS probably due to a very low allele frequency of rs77992257 outside of East Asia (0.018 in Europeans, the 1000 Genomes project). ADAP1 is a biologically plausible candidate for both CD and leprosy because of its role in establishing S. enterica (24) and Y. enterocolitica (25) intracellular infection. Based on the eQTL effects of rs77992257 for ADAP1 on colon tissue in the GTEx data and a blood eQTL database derived from Korean patients with CD, allele G, which confers risk for CD and protection from leprosy, was associated with higher ADAP1 expression level. The biological relevance of this locus in CD and leprosy warrants further study.
Our data indicate that 9 of 11 common susceptibility loci for CD and leprosy have opposite genetic effects, with only the effects of LACC1 and RIPK2 variants showing consistent effect directions in both CD and leprosy, in line with published data (10). The LACC1 protein is associated with mitochondrial and nicotinamide adenine dinucleotide phosphate oxidase-dependent production of reactive oxygen species (ROS), bactericidal activity and inflammasome activation in macrophages (31). On NOD2 stimulation of human macrophages, LACC1 associates with the NOD2-signaling complex and is critical for optimal NOD2-induced signaling, mitochondrial ROS production, cytokine secretion and bacterial clearance. The LACC1 p.Ile254Val (rs3764147) loss-of-function variant results in decreased amplification of host pattern recognition receptor-induced mitochondrial ROS signaling, cytokine secretion and bacterial clearance (32). The NOD2-RIPK2 pathway has been implicated in the recognition of pathogens, by activating downstream NF-kB and MAPK signaling pathways (33). Mutations in NOD2, intracellular sensor of bacterial peptidoglycan, were the first discovered and most strongly associated with the risk of development of CD (34). Although NOD2 loss-of-function increases the risk of CD in Europeans, the mechanisms by which altered NOD2 function leads to overt intestinal inflammation remain poorly understood. In the present study, NOD2 shows opposite genetic effects in CD and leprosy, while RIPK2 and LACC1 have consistent effects in both CD and leprosy. The minor T allele of rs1981760 in NOD2 is associated with lower NOD2 expression in the GTEx database (19) and is a risk allele for leprosy and a protective allele in CD. It is plausible that lower expression of NOD2 may impair host defense against M. leprae infection. The fact that lower NOD2 and higher ADAP1 expression levels are associated with protection and risk for CD, respectively, implicates that CD might be caused by over-reaction to pathogenic stimuli.
Finally, we aimed to quantify the genetic overlap between CD and leprosy, as this has not previously been systematically examined using genome-wide data. LDSC analysis, using summary data from disease-specific meta-analyses of CD and leprosy GWAS, showed that these disorders are negatively correlated (rg [SE] = −0.40, [0.13], P = 2.6 × 10−3). Similarly, PRS analysis showed that PRSleprosy, based on 39 SNPs at best P-value threshold, can explain 5.28% CD variance, supporting a significant overlap in genetic architecture between CD and leprosy. Recently, we reported that PRS derived from Korean data (cohort I) (35) explained up to 14.3% of CD phenotype variance (cohort II) (36). Although CD and UC are both IBDs with many similarities, PRSleprosy could explain only 0.13% of UC variance, suggesting key differences in pathophysiology between the two conditions.
Our study has several shortcomings. First, due to the limited power of Korean CD and Chinese leprosy GWAS, our findings on the genetic architectures of CD and leprosy are not comprehensive. Larger GWAS may reveal more shared loci between CD and leprosy in the future. In addition, the estimates of genetic correlation between CD and leprosy using LDSC may be less accurate than those calculated using individual-level genotype data with the genomic restricted maximum likelihood method (37). Second, CD and leprosy share multiple susceptibility loci; however, the most strongly associated variant at any given locus frequently differs. Thus, defining causal genes among potential candidates at each locus is a major challenge. Third, leprosy has two distinct clinical manifestations, designated as tuberculoid and lepromatous, and it would be interesting to examine which type shares common genetic susceptibility loci with CD or specific CD subtypes. In addition, according to the previous report on the genetic overlap between hyper inflammation in leprosy and IBD (38), there is a possibility that leprosy reactions might help to explain the genetic overlap; however, due to our modest sample size and lack of detailed clinical information on patients with leprosy, we were unable to perform such an analysis. Further studies employing larger number of patients with detailed clinical and laboratory data might provide better estimates of the overlap between the two diseases and clues on whether CD might be caused by over-reaction toward pathogenic stimuli.
Materials and Methods
Study subjects
Korean CD dataset
The discovery cohorts consisted of two previously published GWAS including cohort I (896 cases and 4041 controls) genotyped using the OmniExpress and Omni1-Quad (Illumina) (35) and cohort II (725 cases and 378 controls) genotyped using the Infinium Asian Screening Array-24 v1.0 (Illumina) (36). For replication, we used Immunochip II data (the Infinium ImmunoArray-24 v2, Illumina) from cohort III (738 cases and 488 controls) (39). The novel candidate locus not present in the cohort III data was genotyped using a TaqMan genotyping assay in an independent cohort of 613 cases and 1083 controls. Patient clinical characteristics are summarized in Supplementary Material, Table S12. All patients with IBD were recruited from the IBD Clinic of Asan Medical Center.
Chinese leprosy datasets
Summary statistics from three previously published GWAS were used, including cohort IV (706 cases and 1225 controls) (9), cohort V (1412 cases and 1597 controls) (12) and cohort VI (842 cases and 925 controls) (10), genotyped using Human610-Quad BeadChip, Omni Zhonghua chips and Human 660 K-Quad BeadChip, respectively. Cohort IV and V data were used as input for discovery analysis, to maximize the sample size, and cohort VI data were used for replication analysis.
Meta-analyses
Disease-specific meta-analyses of CD and leprosy were first conducted using GWAS dataset summary statistics, calculated using the frequentist association test in SNPTEST (40), based on an additive model. Disease-specific and cross-disease meta-analyses were performed using the inverse-variance method based on a fixed-effects model in meta v1.7 (41). All SNPs with heterogeneity P < 0.05 were excluded because of possible heterogeneity across studies. Disease-specific meta-analyses were performed by combining the cohort I and II CD GWAS datasets, comprising 1621 cases and 4419 controls with 5 885 673 shared SNPs, and cohort IV and V leprosy GWAS datasets, comprising 2118 cases and 2822 controls with 5 608 384 shared SNPs. SNPs with P < 0.05 in each disease meta-analysis were used for cross-disease meta-analysis. For the standard cross-disease meta-analysis, we used two approaches: (1) a directional meta-analysis for associations in the same direction and (2) an opposite allelic effect analysis for associations in the opposite direction, with the direction of associations flipped in the summary statistics from the CD GWAS datasets. Although these two approaches assume different (same or opposite) directions of effects, in terms of magnitude of effects, both approaches are optimized for finding loci that have similar magnitude of effects on both diseases. Since a pleiotropic locus may have differing magnitude of effects on the two diseases, the standard approach alone may not be sufficient. Therefore, we performed an additional joint analysis of CD and leprosy was performed using PLEIO (17), a multi-trait joint analysis tool that maximizes power by systematically accounting for genetic correlations and heritability of the traits using summary statistics. PLEIO does not necessarily assume that the magnitude of effects is similar for two diseases but is rather optimized to find loci whose effects conform to the observed heritability and genetic correlation pattern. We calculated genetic and environmental correlation matrices using summary statistics from CD and leprosy disease-specific meta-analyses, based on LDSC (29) and the 1000 Genomes Project reference panel (JPT + CHB). From the LDSC analysis, we observed the heritability of CD to be 0.34 (SE: 0.10), the heritability of leprosy to be 0.57 (SE: 0.24), the genetic correlation to be −0.40 and the environmental correlation to be 0.01. This observed information, along with the known prevalence of CD (0.00011) and leprosy (0.00010), was used as input data for PLEIO.
The MHC region (chromosome 6: 25–34 Mb) was excluded from the main analysis, unless mentioned, due to its unusual genetic architecture and LD patterns. Established CD loci were defined as SNPs in LD (r2 > 0.2) with reported SNPs, or by using the LD windows calculated by de Lange et al. (15). In leprosy, SNPs in LD (r2 > 0.2) with reported SNPs were defined as established loci (12). To identify novel shared susceptibility loci in CD and leprosy, we selected SNPs with Pmeta < 1 × 10−6 in the CD and leprosy meta-analyses, with the same or opposite direction effects, for replication using independent cohorts III and VI. Following replication, we considered signals that showed P < 0.05 in each disease separately in the replication cohorts and Pcombined < 2.5 × 10−8 in the CD and leprosy cross-disease meta-analysis, including both discovery and replication datasets, as statistically significant. To estimate the linear correlation of allelic effects directions for the shared CD and leprosy susceptibility loci, lead SNPs at each shared susceptibility locus were selected and used to calculate PCC and P-values using the R package ggpubr v 0.1.9 (https://CRAN.R-project.org/package=ggpubr). PCC in the range of −1 to +1 has a negative number when OR values lie on the negative slope in the scatter plot, and 0 indicates no linear correlation.
Colocalization analysis
For colocalization analysis of shared susceptibility loci between CD and leprosy, we used the coloc R package v.3.2–1 to estimate the posterior probability of CD and leprosy sharing the same causal variant. The coloc calculated five posterior probabilities: PP0 (no association with either trait), PP1 (association with CD), PP2 (association with leprosy), PP3 (two independent causal SNPs associated with each trait) and PP4 (a single SNP associated with both traits). A threshold of PP4 > 0.8 was considered as a strong evidence of a colocalized signal. Summary statistics of disease-specific meta-analyses of CD (cohort I and II) and leprosy (cohort IV, V, and VI) were used for the colocalization analysis. Immunochip data of Cohort III were excluded because of low SNP coverage. Shared SNPs between CD and leprosy within 1 Mb window from the lead SNP of cross-disease meta-analysis were included in the colocalization analysis.
eQTL analysis
To gain insight into the potential functional roles of the susceptibility loci shared between CD and leprosy, an extensive cis-eQTL analysis was conducted by searching lead SNPs of cross-disease meta-analysis in publicly available data from the GTEx database (v.8) (19), and whole blood cis-eQTL databases from eQTLGen (26), Koreans (20) and Japanese (27). In GTEx (19), we selected the relevant tissues for CD and leprosy, namely whole blood, small intestine, transverse colon, sigmoid colon, skin and nerve. To check whether GWAS of CD or leprosy have the same causal SNPs with cis-eQTL data, we performed colocalization analyses using an extensive cis-eQTLs of disease-related tissues and summary statistics of disease-specific meta-analyses of CD (cohort I and II) or leprosy (cohort IV, V and VI). We used the coloc R package v.3.2–1 to estimate the posterior probability. Shared SNPs between summary statistics of disease-specific meta-analyses and cis-eQTL data within 1 Mb window from the each lead SNP were included in the colocalization analysis.
Gene annotation and PPI network analysis
We performed gene analysis using Multi-marker Analysis of GenoMic Annotation (MAGMA) v.1.07b (http://ctg.cncr.nl/software/magma) to prioritize causal genes at shared susceptibility loci between CD and leprosy. Using the cross disease meta-analysis of CD (cohort I and II) and leprosy (cohort IV, V, and VI), and LD information of East Asian population, all SNPs located between the transcription start and end sites were aggregate to that gene to calculate the gene P-value based on a multiple regression model. Of 19 257 reference genes, 17 180 genes were included in the gene analysis. All genes with P < 2.91 × 10−6 (0.05/17180) were considered as significantly annotated genes. We also conduct pathway enrichment analysis using the Gene Ontology resource (http://geneontology.org/) to identify biological pathways enriched for the annotated genes.
To construct the PPI networks of annotated genes in susceptibility loci in common between CD and leprosy, we used the Search Tool for the Retrieval of Interacting Genes/Proteins database (STRING v11.5) (28). Among the genes annotated by MAGMA, we selected a gene with the most significant P-value in each locus as an input data. STRING predicts functional interactions between the proteins, based on co-expression, text-mining, biochemical data from experiments, and previously curated pathway and protein-complex data from databases. We applied the default settings of full STRING network type, medium confidence of interaction score (0.4) and FDR < 0.05.
Genetic correlation
To understand the genetic relationship between CD and leprosy, we estimated the genetic correlation using LDSC v1.0.0 (29). Summary statistics from CD (cohorts I and II) and leprosy (cohorts IV and V) meta-analyses were used as input data, including 4 938 574 overlapping SNPs between the CD and leprosy datasets. The East Asian data (JPT + CHB) from the 1000 Genomes Project were used as an LD reference panel for all calculations of the genetic correlation.
Polygenic risk scores
PRSs were used to assess the genetic overlap between CD and leprosy. To estimate the variance of CD explained by the PRS derived from leprosy GWAS data (PRSleprosy), PRSleprosy was computed by summing the risk alleles associated with CD, weighted by the effect size estimates generated by a meta-analysis of leprosy GWAS data. To achieve the largest sample size possible, we used the results of meta-analysis of leprosy cohorts IV, V and VI as base data for estimating effect sizes, and the CD cohorts I, II and III as the target data for evaluating PRS. First, the eight lead SNPs reaching the genome-wide significance threshold (P < 5x10−8) in the meta-analysis of leprosy were included in the calculation of PRS. Those SNPs and their corresponding weights are presented in Supplementary Material, Table S10. Additionally, the directions of effect of two susceptibility loci (LACC1 and RIPK2 loci), which showed directionally consistent associations between leprosy and CD, were manually flipped to maximize the variance explained by PRSleprosy for CD. We then compared the full model (including the PRS) with the null model (with the PRS variable excluded) and estimated the variance explained using Nagelkerke’s pseudo-R2.
Next, PRS values were calculated using PRSice-2 (30) across different P-value thresholds to identify the best P-value threshold for variant selection, which allowed inclusion of variants that failed to reach genome-wide significance, and yet contributed to disease risk. To minimize overfitting due to tight LD, we treated the MHC region (chromosome 6: 25 ~ 34 Mb) with additional caution by selecting only the most significant variant in this region. East Asian (CHB + JPT) 1000 Genomes data were used as a reference panel to calculate LD structure. A clumping algorithm was set to create clumps of SNPs spanning 250 kb in LD with an r2 threshold ≥0.1. PRS values at best P-value threshold were z-transformed before separation of all cases with CD into quintiles based on their individual PRS. The association between PRSleprosy and CD was further evaluated in generalized mixed regression models.
Web resources
The URLs for data presented herein are as follows:
METAL, http://www.sph.umich.edu/csg/abecasis/metal/
Online Mendelian Inheritance in Man (OMIM), http://www.omim.org
The 1000 Genomes Project, http://www.1000genomes.org/
UCSC Genome Browser, http://genome.ucsu.edu/
GCTA, cnsgenomics.com/software/gcta/.
IIBDGC, www.ibdgentics.org
RegulomeDB v2, http://www.broadinstitute.org/mammals/haploreg/haploreg.php
Genotype-Tissue Expression (GTEx) project, http://www.gtexportal.org/home
eQTL Blood Browser, http://www.genenetwork.nl/bloodeqtlbrowser/
Geuvadis/1000 Genomes resources, http://www.ebi.ac.uk/Tools/geuvadis-das/
PLEIO, https://github.com/cuelee/pleio
Data Availability
Summary Korean CD data are available from the European Genome-Phenome Archive under accession no. EGAS00001005026. All the summary statistics of leprosy GWAS are available on the dbGaP (Study Accession: phs000217.v3.p1). Summary data from the cross-disease meta-analysis supporting the findings of this study are available from the GWAS catalog under GCST90103738 (CD ASA and GWAS meta), GCST90103739 (CD_ASA_GWAS + CD_Leprosy_L1_L3 meta, non-flip version), GCST90103740 (CD_ASA_GWAS + CD_Leprosy_L1_L3 meta, flip version).
Acknowledgements
This study was supported by PLSI supercomputing resources of the Korea Institute of Science and Technology Information.
Conflict of Interest statement
The authors disclose the following:
B.D.Y.: grant/research support from CELLTRION, Inc.; lecture fees from Abbvie Korea, CELLTRION, Inc., Janssen Korea, Shire Korea, Takeda Korea and IQVIA; consultant fees from Abbvie Korea, Ferring Korea, Janssen Korea, Kangstem Biotech, Kuhnil Pharm., Shire Korea, Takeda Korea, IQVIA, Cornerstones Health and Robarts Clinical Trials Inc.
S.-K.Y.: grant/research support from Janssen Korea Ltd.
B.H. is the CTO of Genealogy Inc.
The remaining authors have no conflicts of interest.
Funding
Mid-career Researcher Program grant through the National Research Foundation of Korea (2020R1A2C2003275 to K.S.), Ministry of Science and ICT, Republic of Korea; Academic promotion program of Shandong First Medical University to F.R.Z.
Author Contributions
K.S. obtained financial support. K.S. conceived and designed the study. S.W.H., S.H.P., S.K.Y. and B.D.Y. collected clinical samples. K.S., F.Z., J.L. and B.H. supervised the data analysis and interpretation. S.J., D.P., Y.K. J.B. Y.S. and H.L. performed data analyses. S.J. and K.S. drafted the manuscript. K.S., B.H., and H.S.L. revised the manuscript.
References
Author notes
Kyuyoung Song, Jianjun Liu and Furen Zhang jointly directed this work.
