NFKB1 and MANBA Confer Disease Susceptibility to Primary Biliary Cholangitis via Independent Putative Primary Functional Variants.

BACKGROUND & AIMS: Primary biliary cholangitis (PBC) is a chronic and cholestatic liver disease that eventually leads to cirrhosis and hepatic failure. We recently identified several susceptibility genes included NFKB1 and MANBA for PBC in the Japanese population by genome-wide association study. However, the primary functional variants in the NFKB1/MANBA region and the molecular mechanism for conferring disease susceptibility to PBC have not yet been clarified.
METHODS: We performed high-density association mapping based on a single-nucleotide polymorphism (SNP) imputation analysis, using data from a whole-genome sequence reference panel of 1070 Japanese individuals and the previous genome-wide association study (1389 PBC patients, 1508 healthy controls). Among SNPs (P < 5.0 × 10-7) in the NFKB1/MANBA region, putative primary functional variants and the molecular mechanism for conferring disease susceptibility to PBC were identified by in silico/in vitro functional analysis.
RESULTS: Among the SNPs in the NFKB1/MANBA region, rs17032850 and rs227361, which changed the binding of transcription factors lymphoid enhancer-binding factor 1 (LEF-1) and retinoid X receptor α (RXRα), respectively, were identified as putative primary functional variants that regulate gene expression. In addition, expression-quantitative trait locus data and gene editing using a clustered regularly interspaced short palindromic repeat (CRISPR)/Cas9 system supported the potential role of rs17032850 and rs227361 in regulating NFKB1 and MANBA expression, respectively.
CONCLUSIONS: We identified independent putative primary functional variants in NFKB1/MANBA and showed the distinct molecular mechanism by which each putative primary functional variant conferred susceptibility to PBC. Our approach was useful to dissect the pathogenesis not only of PBC, but also other digestive diseases in which NFKB1/MANBA has been reported as a susceptibility locus.

We have identified novel putative primary functional variants for primary biliary cholangitis (PBC) in NFKB1/MANBA, and shown the molecular mechanism by which the variants confers susceptibility to PBC. This study indicates that NFKB1/lymphoid enhancer-binding factor 1 and MANBA/retinoid X receptor α have important roles in the development of PBC.

Primary biliary cholangitis (PBC) is a chronic and progressive cholestatic liver disease that is characterized histologically by chronic non-suppurative destructive cholangitis, ductopenia, interface hepatitis, fibrosis, and, eventually, biliary cirrhosis. Because of the existence of well-defined autoantigens (such as the E2 component of the pyruvate dehydrogenase complex [PDC-E2]), homogeneous clinical expression, and an overwhelming predominance in females, PBC is considered an organ-specific autoimmune disease.1, 2 The destruction of small bile ducts is mediated not only by adaptive immune responses, including CD4+ T cells, CD8+ T cells, and B cells, but also innate immune responses, including natural killer cells against biliary epithelial cells.2, 3, 4

The higher concordance rate for PBC in monozygotic twins compared with dizygotic twins, and the higher estimated relative sibling risk (λs), indicates the involvement of strong genetic factors in the development of PBC.5, 6 Previous genome-wide association studies (GWASs) and subsequent meta-analyses identified HLA and 31 non-HLA susceptibility regions (IL12A, IL12RB2, STAT4, IRF5, MMEL1, SPIB, DENND1B, CD80, IL7R, CXCR5, TNFRSF1A, CLEC16A, NFKB1, RAD51L1, MAP3K7IP1, PLCL2, RPS6KA4, TNFAIP2, ELMO1, IRF8, TNFSF11, SH2B3, CRHR1, TYK2, IL1RL2, CCL20, DGKQ, C5orf30, LOC285626, TNFAIP3, and chromosome 17q12-21) as susceptibility genes for PBC in patients of European descent.7, 8, 9, 10, 11, 12, 13, 14 In the Japanese population, we previously identified HLA, TNFSF15, POU2AF1, PRKCB, NFKB1, MANBA, IL7R, and the chromosome 17q12-21 locus as susceptibility regions for PBC by GWAS.15, 16

The NFKB1/MANBA gene region, which is located on human chromosome 4q24, has been identified as a susceptibility region in several digestive disorders, such as ulcerative colitis (rs3774937, P = 5 × 10-14; rs3774959, P = 4 × 10-12)17, 18 and Crohn's disease (rs13126515, P = 4 × 10-8), as well as PBC. Nuclear factor-κB (NF-κB) is a well-known transcription factor complex that includes p50 (NF-κB1 and its precursor p105), p52 (NF-κB2 and its precursor p100), Rel (c-Rel), RelA (p65), and RelB, and has important roles in apoptosis, lymphoid organogenesis, innate and adaptive immunity, and central and peripheral tolerance. In Nfkb1-/- mice, proliferation of T cells and production of T-helper type 2 (Th2) cytokines reportedly are reduced. In addition, these Nfkb1-/- mice showed susceptibility to experimental autoimmune encephalomyelitis, typhlocolitis, and infection with Leishmania major, but resistance to asthma.21, 22, 23, 24 The protein product of MANBA on the other hand, β-mannosidase, is the final exoglycosidase in the N-linked glycoprotein oligosaccharide catabolism pathway, and is localized to the lysosomal compartment. β-mannosidosis is a rare lysosomal storage disease caused by the deficiency of β-mannosidase. The most severe cases of human β-mannosidosis show not only mental retardation, developmental delay, and dysmorphology, but also frequent infection.26, 27, 28, 29, 30, 31 Taken together, the evidence indicates that both NFKB1 and MANBA may play important, but different, roles in the regulation of the immune system. However, largely because of the existence of a shared linkage disequilibrium (LD) block in this region (Figure 1A), primary functional variants in the NFKB1/MANBA loci and the molecular mechanisms by which these variants confer susceptibility to PBC have yet to be clarified.

Figure 1

LD structure around (A) SNPs in this locus for the analysis of LD were filtered by Minor allele frequencies (MAF) > 0.05 and has publications in the NCBI variation viewer, and 95 SNPs were selected. LDs of these SNPs were analyzed by LDlink using the 1000genome Japanese in Tokyo (JPT) data. (B and C) r2 of each SNP with (B) rs230534, which was the NFKB1 top-hit SNP with susceptibility to PBC and (C) rs2866413, which was the MANBA top-hit SNP with susceptibility to PBC, using the present genotyping data. (B) rs2866413 and (C) rs230534 are shown as red dots. Stronger colors indicate higher r2 scores. Chr, chromosome.

In the present study, to identify putative primary functional variants in the NFKB1/MANBA region, we performed high-density association mapping based on single-nucleotide polymorphism (SNP) imputation analysis using the data from a whole-genome sequence reference panel of 1070 Japanese individuals and our previous GWAS.15, 16 We next performed in silico and in vitro functional analyses to identify putative primary functional variants that alter gene expression. Finally, we attempted to identify the molecular mechanisms underlying disease susceptibility to PBC associated with these functional variants by expression-quantitative trait locus analyses and gene editing using a clustered regularly interspaced short palindromic repeat (CRISPR)/Cas9 system.

Results

SNP Imputation Analysis and High-Density Association Mapping

Only 16 SNPs located in NFKB1/MANBA had been genotyped in our previous GWAS using the Affymetrix (Santa Clara, CA) Axiom Genome-Wide ASI1 Array.15, 16 To analyze all SNPs located in this region, SNP imputation analysis and high-density association mapping for disease susceptibility to PBC were performed using a reference panel of 1070 Japanese individuals (1KJPN) and genotype data from our previous GWAS.16, 32

Besides the SNPs which were significant to a level of P < 5.0 × 10-7 in our previous GWAS (rs230534, rs2866413, rs223489, and rs2903283), an additional 131 SNPs were found to be significant to this level by SNP imputation analysis, probably because of the existence of strong LD in this region (Figure 1 and Supplementary Table 1). Of these SNPs, 16 were located in the 5’ intergenic region of NFKB1, 67 were located in NFKB1 introns, 1 was located in the 5’ intergenic region of MANBA, 2 were located in the coding region of MANBA, and 49 were located in MANBA introns (Supplementary Table 1 and Figure 2A and B).

Figure 2

High-density association mapping for the selection of candidate putative primary functional variants in (A) Results of high-density association mapping in NFKB1/MANBA based on SNP imputation analysis. A whole-genome sequence reference panel of 1KJPN and genotype data from our previous GWAS (1389 PBC patients, 1508 healthy controls) were used. Dots show the P value of each SNP, and red dots show experimentally genotyped SNPs from GWAS. (B) Schematic of the selection process for candidate putative primary functional variants. There were 724 and 135 SNPs whose P values were less than 1.0 × 10-3 and 5.0 × 10-7, respectively, that were located in NFKB1/MANBA. Of the 135 SNPs, rs230534, which showed the most significant association with susceptibility to PBC in the NFKB1/MANBA region, and the candidate putative primary functional variants located in a transcription regulatory element (NFKB1: rs2272676, rs230504, rs17032850, rs3774959, and rs7678055; MANBA: rs227361) were chosen. (C) The positions of the 7 candidate functional variants in NFKB1/MANBA are shown. Chr, chromosome.

Although there was a strong LD block in the NFKB1/MANBA locus, LDs were not as strong between NFKB1 SNPs and MANBA SNPs, which showed relatively stronger association with susceptibility to PBC (Figure 1B and C). Therefore, primary functional variants seemed to be located in both the NFKB1 locus and the MANBA locus.

Identification of the Putative Primary Functional Variant in NFKB1

From the total of 83 NFKB1 SNPs, we used the RegulomeDB database (http://www.regulomedb.org/index) to identify those that were located in DNase hypersensitivity clusters and histone acetylation regions, and that were predicted to affect the binding of transcription factors. Twelve SNPs with RegulomeDB scores higher than 2c, which were supported by transcription factor binding data and DNase peaks, were selected as potential candidates (Table 1) (rs747559 in the 5’ intergenic region; rs2272676 in intron 1; rs230532 in intron 2; rs230493, rs230504, rs230519, and rs230520 in intron 5; rs227361 in intron 10; rs1598859, rs3774959, and rs17032850 in intron 11; and rs4648055 in intron 12). Among these 12 SNPs, 5 SNPs (Table 1) (rs2272676, rs230504, rs17032850, rs3774959, and rs4648055), which were located in DNase I hypersensitivity clusters and H3K27Ac markers in at least 1 cell type in the UCSC genome browser (http://genome.ucsc.edu/index.html), were selected as the final candidates. In addition, the most significant SNP in the NFKB1/MANBA region identified from high-density association mapping (Table 1) (rs230534) also was selected as the final candidate in the NFKB1 locus for the conferral of susceptibility to PBC (Figure 2B and C). The PICS algorithm (https://pubs.broadinstitute.org/pubs/finemapping/pics.php), which calculates the most likely causal SNPs given the observed association signal at a locus, indicated that all of the final candidates showed the PICS values of greater than 0.001 were indexed by rs230534 (ie, the top-hit SNP in the NFKB1 locus) (Supplementary Table 2).

Table 1

SNPs Associated With Susceptibility to PBC in the Japanese Population in NFKB1/MANBA by High-Density Association Mapping

SNP_IDa	GWAS or imputationb	Position (chromosome 4)c	Pd	ORe	RegulomeDBf	UCSCg	Location
rs230534	GWAS	103449041	1.34E-08	1.37	4	△	NFKB1 intron 2
rs35185023	Imputation	103448080	2.61E-08	1.37	6	×	NFKB1 intron 2
rs2978641	Imputation	103575280	1.17E-07	1.33	6	×	MANBA intron 12
rs227277	Imputation	103594938	1.23E-07	1.34	6	×	MANBA intron 8
rs230531	Imputation	103450377	1.24E-07	1.34	6	×	NFKB1 intron 2
rs11097786	Imputation	103406511	1.27E-07	1.34	6	×	5' of NFKB1
rs141936164	Imputation	103401723	1.37E-07	1.34	6	×	5' of NFKB1
rs228613	Imputation	103577465	1.42E-07	1.33	No data	×	MANBA intron 12
rs228612	Imputation	103575385	1.56E-07	1.33	No data	×	MANBA intron 12
rs62328536	Imputation	103399049	1.59E-07	1.34	4	×	5' of NFKB1
rs2272676	Imputation	103423326	1.87E-07	1.34	2b	○	NFKB1 intron 1
rs2866413	GWAS	103557077	2.40E-07	1.32	No data	Nonsynonymous	MANBA exon 14 (T701M)
rs227368	Imputation	103611845	2.80E-07	1.33	5	Nonsynonymous	MANBA exon 6 (V253I)
rs230504	Imputation	103481561	3.58E-07	1.33	1d	△	NFKB1 intron 5
rs3774959	Imputation	103511114	3.74E-07	1.32	1b	△	NFKB1 intron 11
rs17032850	Imputation	103507703	3.74E-07	1.32	2b	△	NFKB1 intron 11
rs4648055	Imputation	103515313	4.27E-07	1.32	1f	△	NFKB1 intron 12
rs227361	Imputation	103586977	4.32E-07	1.32	1b	○	MANBA intron 10

SNPs with underlines were the final candidate primary functional variants.

Genotyped by our previous GWAS or the imputed genotypes by the high-density association mapping in the present study.

Position of the SNPs in hg19.

P values calculated by the Pearson chi-square test for the allelic model.

Odds ratio (OR) of minor allele from the 2-by-2 allele frequency table.

Functinoal prediction scores of each SNP by the RegulomeDB database.

Probability of the functional damages checked by the UCSC genome browser.

To evaluate the effect of final-candidate SNPs on the binding affinity of transcription factors, electrophoretic mobility shift assays (EMSAs) were performed using biotin-labeled probes corresponding to the different alleles of each SNP and nuclear extracts of either the HepG2 (human liver carcinoma) or Jurkat (human T lymphocyte) cell lines. Among the final candidates, a difference in mobility shift between the major allele (G allele) and the minor allele (C allele, PBC susceptible) of rs17032850 was detected (Figure 3A and B). In HepG2 cells, the thickness of the shifted band was different between major alleles and minor alleles (Figure 3A), whereas in Jurkat cells, the position of the shifted bands was different (Figure 3B). These band shifts were abrogated by incubation with a 200× concentration of a nonlabeled competitor probe (Figure 3C). However, no difference in mobility shift was detected for the other candidate SNPs (rs230534, rs2272676, rs230504, rs3774959, and rs7678055) (Figure 3A and B). These results indicated that the putative primary functional variant in NFKB1 was rs17032850.

Figure 3

In vitro functional analysis of each candidate variant of (A and B) EMSA of each of the 6 candidate primary variants using biotin-labeled probes corresponding to the major and the minor alleles, and nuclear extracts of (A) HepG2 and (B) Jurkat. rs17032850 was the only variant to show a difference in mobility shift between the 2 alleles. (C) Competitor assay, using HepG2 and Jurkat nuclear extracts and a 200× concentration of unlabeled probe corresponding to either the C (ie, PBC susceptibility) or G alleles of rs17032850. (D) Outline of reporter plasmid constructs. PCR fragments of intron 11 of NFKB1, containing rs17032850, were subcloned into the pGL4.23 vector. (E and F) Transcription was measured by cellular luciferase activity, 24 hours after transfection of (E) HepG2 cells and (F) Jurkat cells. The luciferase activity of cells transfected with the PBC susceptibility allele (C allele) of rs17032850 was decreased compared with those transfected with the G allele. Three independent experiments with triplicate measurements were performed for each assay, and data represent means ± SD. *P < .05 (Student t test).

Molecular Features Associated With the PBC Susceptibility Allele of rs17032850

To determine differences in transcription efficiency between the major and PBC susceptibility alleles of rs17032850, luciferase reporter assays were performed in HepG2 and Jurkat cells. In HepG2, the luciferase activity of cells 24 hours after transfection with a reporter construct containing the C allele (ie, the PBC susceptibility allele) of rs17032850 was reduced significantly compared with those transfected with a plasmid containing the G allele (Figure 3D and E). However, in Jurkat cells, a difference in luciferase activity between the rs17032850 C allele and G allele was not observed (Figure 3F). This likely is owing to differences in transcription factor binding to these alleles between the different cell types.

By using the TRANSFAC database (Qiagen, Valencia, CA; http://www.gene-regulation.com/pub/databases.html), we predicted that the G allele of rs17032850, but not the C allele, constituted a lymphoid enhancer-binding factor 1 (LEF-1) binding motif (Figure 4A–C). We first confirmed that LEF-1 is expressed in both HepG2 and Jurkat cells (Figure 4D). Subsequently, the binding of LEF-1 was investigated in a supershift assay using nuclear extracts of HepG2. Consistent with the prediction of transcription factor binding, the mobility shift associated with the G allele of rs17032850 was supershifted by pre-incubation with an anti–LEF-1 antibody before electrophoresis (Figure 4E).

Figure 4

Identification of transcription factors targeting the G allele of rs17032850. (A and B) Prediction of transcription factor binding using the TRANSFAC database. (A and B) The G allele of rs17032850 produces a (A) LEF-1 binding site, however, this sequence is disrupted in the (B) C allele. (C) The binding motif of LEF-1, with the position of rs17032850 shown as a rectangle. (D) LEF-1 expression levels in cancer cell lines. Both Jurkat and HepG2 cells showed abundant expression of LEF-1. Data were obtained using GeneCards from the Weizmann Institute of Science (http://www.genecards.org/). (E) Supershift was observed after incubation of HepG2 cell nuclear extracts with an anti–LEF-1 antibody. Three independent experiments were performed in each assay. Ab, antibody.

To assess the influence of rs17032850 on endogenous gene expression, the expression levels of 11 genes located within 1 Mb of NFKB1/MANBA were compared (using the GTEx portal database; http://gtexportal.org/home/) for the different genotypes of rs17032850 in every organ in which the genes’ expression levels were above the threshold for detection. Among these 11 genes, NFKB1 (muscle_skeletal: P = .000062) and UBE2D3 (esophagus_mucosa: P = .000013) showed effects on expression that were associated significantly with the rs17032850 genotypes (Figure 5). These results indicated that NFKB1 and UBE2D3 were candidate effector genes of rs17032850, a PBC susceptibility allele that alters a putative LEF-1 binding site.

Figure 5

rs17032850 genotypes are associated with differences in endogenous Comparison between rs17032850 genotype and the endogenous expressions of (A) NFKB1 and (B) UBE2D3 in all tissues registered in the GTEX database. Statistical significance levels after Bonferroni multiple comparison correction: P = .000108 (0.05/[11 genes × 45 organs]).

By expression-quantitative trait locus analysis, not only primary functional variants, but also other variations that show strong LD with the primary functional variants, are associated with gene expression levels. Therefore, the effect of rs17032850 on gene expression levels was assessed using rs17032850 genotype knock-in versions of the Jurkat cell lines constructed using a CRISPR/Cas9 system. Among the candidate effector genes, only NFKB1 showed an expression level that differed significantly between rs17032850-C/C clones (n = 5) and rs17032850-G/G clones (n = 5) (P = .016; Mann–Whitney U test) (Figure 6A and B). In addition, the difference in NFKB1 expression level between rs17032850-C/C clones and rs17032850-G/G clones was eliminated by knock-down of LEF-1 (Figure 6C). These results confirmed that rs17032850 is a primary functional variant of the NFKB1 locus, such that the rs17032850 SNP alters LEF-1 binding to, and regulation of, the NFKB1 gene.

Figure 6

The effects of rs17032850 on gene expression in the knock-in clones of the Jurkat cell line generated using a CRISPR/Cas9 system. (A and B) Expression levels of (A) NFKB1 and (B) UBE2D3 when comparing between rs17032850-C/C knock-in clones and rs17032850-G/G knock-in clones. (C) NFKB1 expression levels in each knock-in clone after the knock-down of LEF-1. *P < .05 (Mann–Whitney U test). miRNA, microRNA.

Identification of the Putative Primary Functional Variant in MANBA

From the 52 total MANBA SNPs that showed P values less than 5.0 × 10-7 by SNP imputation analysis, candidate functional variants that were predicted to be located in transcription regulatory elements were selected. One SNP with a RegulomeDB score higher than 2c was selected as a potential candidate (Table 1) (rs227361 in intron 10). We confirmed that this SNP is located in a DNase I hypersensitivity cluster and a H3K27Ac marker region for at least 1 cell type, using the UCSC genome browser, and this SNP consequently was selected as the final MANBA locus candidate for the conferral of PBC susceptibility (Figure 2B and C). The PICS algorithm indicated that rs227361 showed the highest PICS values indexed by rs2978641 (ie, the top-hit SNP in the MANBA locus) (Supplementary Table 3).

By EMSA, a difference in mobility shift between the major allele (A allele) and the minor allele (G allele, PBC susceptible) of rs227361 was detected using nuclear extracts from HepG2 and Jurkat (Figure 7A and B). The bands shifted by the major and minor alleles were specifically diminished by incubation with a 200× concentration of a nonlabeled competitor probe (Figure 7A and B).

Figure 7

In vitro functional analysis of (A and B) EMSA and competitor assays of rs227361 using biotin-labeled probes corresponding to the A and G alleles, and nuclear extracts of (A) HepG2 and (B) Jurkat. There was a difference in mobility shift between the 2 alleles. (C) Outline of reporter plasmid constructs. PCR fragments of intron 10 of MANBA, containing rs227361, were subcloned into pGL4.23. (D and E) Transcription was measured by cellular luciferase activity 24 hours after transfection of (D) HepG2 cells and (E) Jurkat cells. The luciferase activity of cells transfected with the PBC susceptibility allele (G allele) of rs227361 was decreased compared with those transfected with the A allele. Three independent experiments with triplicate measurements were performed for each assay, and data represent means ± SD. *P < .05 (Student t test).

Among the MANBA SNPs which showed P values less than 5.0 × 10-7 by SNP imputation analysis, rs227368 (V253I) and rs2866413 (T701M), located in the 6th and 14th exons, respectively, also had the potential to be selected as candidate functional variants (Table 1). However, after SWISS model (https://swissmodel.expasy.org/) and Polyphen2 analysis (http://genetics.bwh.harvard.edu/pph2/), neither nonsynonymous substitution was predicted to have a deleterious effect on MANBA protein folding and function (Figure 8).38, 39 These results indicated that the only putative primary functional variant in MANBA was rs227361.

Figure 8

Prediction of the effect of nonsynonymous substitutions on MANBA protein functionality. (A–C) By both (A) SWISS model and (B and C) Polyphen2, the minor alleles of the indicated nonsynonymous substitutions were predicted to have no damaging effect on MANBA protein folding or function.

Molecular Features Associated With the PBC Susceptibility Allele of rs227361

To assess differences in transcription efficiency between the major and the PBC susceptibility alleles of rs227361, luciferase reporter assays were performed in HepG2 and Jurkat cells. The luciferase activity of cells 24 hours after transfection of a reporter construct containing the G allele (ie, PBC susceptibility allele) of rs227361 was significantly lower compared with cells transfected with a plasmid containing the A allele in both cell lines (Figure 7C–E).

The G allele of rs227361, but not the A allele, was predicted using the TRANSFAC database to constitute a binding motif for the retinoid X receptor alpha (RXRα) (Figure 9A–C). RXRα forms a multiple-subunit corepressor complex that stabilizes repressive local chromatin structure, and blocks access of the transcriptional machinery to the promoter. Therefore, RXRα is known as a nuclear receptor that binds to response elements in the promoter region of target genes, and inhibits transcription without an agonist ligand.36, 40 RXRα was found to be expressed in both the HepG2 and Jurkat cell lines (Figure 9D), and, hence, a supershift assay using the nuclear extract of HepG2 was performed to investigate the binding of RXRα with the G allele of rs227361. Concordant with the prediction, a supershift was found for the G allele of rs227361 after pre-incubation with an anti-RXRα antibody (Figure 9E).

Figure 9

Identification of transcription factors targeting the G allele of rs227361. (A and B) Prediction of transcription factor binding using the TRANSFAC database. The G allele of rs227361 produces a (A) RXRα binding site, however, this sequence is disrupted in the (B) A allele. (C) The binding motif of RXRα, with the position of rs227361 shown as a rectangle. (D) RXRα expression in cancer cell lines. Both Jurkat and HepG2 cells showed abundant expression of RXRA. Data were obtained using GeneCards (http://www.genecards.org/). (E) Supershift was observed by incubation of HepG2 cell nuclear extracts with an anti-RXRα antibody. Three independent experiments were performed in each assay. Ab, antibody.

The effect of rs227361 genotype on endogenous gene expression levels also was compared using the GTEx portal database as well as rs17032850. Among these 11 genes, NFKB1 (muscle skeletal: P = 4.5 × 10-6), MANBA (breast mammary tissue: P = .000035; small intestine terminal ileum: P = .000079; spleen: P = 2.5 × 10-6), UBE2D3 (esophagus mucosa: P = .000016; lung: P = .000043; pancreas: P = 8.9 × 10-6; stomach: P = 7.8 × 10-6), CISD2 (artery tibial: P = 1.2 × 10-8; nerve tibial: P = .000045; skin not sun exposed: P = 5.3 × 10-7; skin sun exposed: P = 8.0 × 10-11; thyroid: P = 3.8 × 10-12; uterus: P = 6.9 × 10-6), and BDH2 (esophagus muscularis: P = .000034; muscle skeletal: P = .000060) showed changes in expression that were associated significantly with the genotypes of rs227361 (Figure 10). These results indicated that NFKB1, MANBA, UBE2D3, CISD2, and BDH2 are candidate effector genes of rs227361, a PBC-susceptibility allele that has higher binding affinity for the inhibitory transcription factor RXRα.

Figure 10

rs227361 genotypes are associated with differences in endogenous Comparisons between the genotype of rs227361 and endogenous expression levels of (A) NFKB1, (B) MANBA, (C) UBE2D3, (D) CISD2, and (E) BDH2 in all tissues registered in the GTEX database. Statistical significance levels after Bonferroni multiple comparison correction: P = .000108 (0.05/[11 genes × 45 organs]).

Next, the effect of rs227361 on gene expression levels was assessed using the rs227361 genotype knock-in versions of Jurkat cell lines constructed using a CRISPR/Cas9 system. Among the candidate effector genes, only MANBA showed an expression level that differed significantly between rs227361-A/A clones (n = 12) and rs227361-G/G clones (n = 10) (P = .00037; Mann–Whitney U test) (Figure 11A–E). In addition, the difference in MANBA expression level between rs227361-A/A clones and rs227361-G/G clones was eliminated by the knock-down of RXRα (Figure 11F). These results confirmed that rs227361 is a primary functional variant in the MANBA locus, such that the rs227361 SNP alters RXRα binding to, and regulation of, the MANBA gene.

Figure 11

The effects of rs227361 on gene expression in the knock-in clones of the Jurkat cell line generated using a CRISPR/Cas9 system. (A–E) Expression levels of (A) MANBA, (B) NFKB1, (C) UBE2D3, (D) CISD2, and (E) BDH2 comparing between rs227361-A/A knock-in clones and rs227361-G/G knock-in clones. (F) MANBA expression levels in each knock-in clone after the knock-down of RXR-α. *P < .01, **P < .001 (Mann–Whitney U test).

Co-segregation of the Effects of rs17032850 and rs227361 to the Disease Susceptibility of Top-Hit SNPs in This Locus

Finally, to assess the effects of these 2 putative primary functional variants on the disease susceptibility of top-hit SNPs, multiple regression analysis was performed. When rs17032850 and rs227361 were incorporated into a regression model (Homozygote of risk alleles = 1, other genotypes = 0), these SNPs entirely explained the association signal of rs230534 (TT = 2, CT = 1, CC = 0; P = 7.55 × 10-193). In addition, the risk genotypes of rs17032850 and rs227361 co-segregated with the rs230534-T allele (Figure 12). These results indicated that most of the effects of rs230534 were covered by rs17032850 and rs227361.

Figure 12

Co-segregation of the risk alleles of rs17032850 and rs227361 with the rs230534 T allele. The 2 putative primary functional variants co-segregated with the rs230534 T allele, which strongly associates with the susceptibility to PBC in the Japanese population. For each rs230534 genotype, the graph shows the percentages of individuals who are heterozygous or homozygous for at least 1 PBC-protective allele. homo, homozygote; hetero, heterozygote.

Discussion

In the present study, rs17032850, located in intron 11 of NFKB1, and rs227361, located in intron 10 of MANBA, were identified as the putative primary functional variants for disease susceptibility to PBC. The disease susceptibility allele of rs17032850 was shown to disrupt a LEF-1 binding site and was associated with the endogenous expression levels of NFKB1. The disease-protective allele of rs227361 on the other hand, disrupted a RXRα binding site and was associated with significantly higher MANBA endogenous expression than in individuals without this allele.

There were some limitations to this study. First, although the transcription binding sites are believed to be critical components underlying disease susceptibility, important regulatory variants in autoimmune disease are only within transcription binding sites at 10% to 20%. Second, the RegulomeDB is likely out of date, given that this database was last updated in 2014. Third, the PICS algorithm, which is an LD-based method for predicting SNPs that are most likely to be the causal variants responsible for associations with gene expression levels, indicated that rs227361 showed the highest PICS score in the MANBA locus, but did not indicate this for rs17032850 in the NFKB1 locus. Fourth, in the logistic regression analysis used in the present study, most of the effects of rs17032850 and rs227361 were covered by rs230534, but this coverage was not complete. Therefore, some other primary functional variants presumably exist. However, results in this study indicate that the identified PBC susceptibility alleles of NFKB1 and MANBA could be the functional variants that alter the transcription of target genes by the reduction of the effect by these limitations using in vitro functional evaluation analysis.

In the present study, individuals with the PBC-susceptible allele (G allele) of rs227361, in the intron of MANBA, showed significantly lower endogenous levels of MANBA expression. To date, the contribution of β-mannosidase to the pathogenesis of digestive disorders has not been reported, and the relative level of MANBA expression in PBC patients and healthy controls has not been assessed. However, the phenotype of human β-mannosidosis, which is caused by the deficiency of β-mannosidase in lysosomes, may provide us with clues to the potential role of β-mannosidase in the onset of digestive disorders. Lysosomes in β-mannosidosis patients reportedly suffer from lysosomal membrane permeabilization, and this type of lysosomal damage can induce defects in autophagy. Dysfunction of autophagy is well known to play a role in the pathogenesis of digestive disorders. Additional studies are warranted to improve our understanding of the relationship between PBC pathogenesis and β-mannosidase.

To date, the following pathogenic pathway related to adaptive immune responses has been proposed in PBC: (1) release of intact PDC-E2, a well-defined, specific PBC autoantigen, from biliary epithelial cells during apoptosis; (2) endocytosis of PDC-E2 by antigen-presenting cells; (3) deficits in immune tolerance as the result of genetic and environmental factors; (4) presentation of degraded PDC-E2 peptide to autoreactive CD4+ T cells via HLA class II, and autoreactive CD8+ T cells via HLA class I; (5) interaction between autoreactive CD4+ T cells and autoreactive B cells that produce antimitochondrial autoantibody, and (6) immune-mediated destruction of biliary epithelial cells in the small bile duct. Innate immune responses mediated by natural killer cells also contribute to the pathogenesis of PBC. Indeed, most of the susceptibility genes to PBC identified to date have functions related to this disease-development pathway. For example, the protein products of HLA-DP, HLA-DQ, and HLA-DR, reported as the strongest susceptibility genes for PBC,44, 45 are known as HLA class II molecules expressed in antigen-presenting cells. The protein product of TNFSF15 has been reported to be involved in polarization to Th1 and Th17 cells.46, 47 The product of IL12RB2 is a subunit of interleukin 12, which activates both innate and acquired immune responses. Although the exact disease pathway in which β-mannosidase operates remains to be elucidated, the present study provides important clues in regard to their role in the pathogenesis of PBC. Crohn’s disease and PBC in Asian populations have been reported to share a genetic background related to tumor necrosis factor signaling,49, 50 and NFKB1/MANBA has been reported as the susceptibility region for both of these diseases.16, 19 Both human genetic studies and mouse studies51, 52, 53, 54 have reported that defects in autophagy are involved in the pathogenesis of Crohn’s disease, and, therefore, rs227361, which was identified as the putative primary functional variant in MANBA in the present study, also may contribute to the dysfunction of autophagy via lysosomal damage in Crohn’s disease. Although the top-hit SNP in the NFKB1/MANBA region in Crohn’s disease (rs13126505) is not found in Asian populations, the top-hit SNP in the same region in ulcerative colitis (rs3774959) shows strong LD (r2 = 0.99, D’ = 0.99) with rs17032850, one of the putative primary functional variants in the present study. Thus, rs17032850 may operate as a functional variant in both PBC and ulcerative colitis.

GWASs have identified thousands of genetic variants associated with complex human diseases, including digestive disorders. However, these genetic variants often explain only a relatively small proportion of the heritability of the disease. To overcome the missing heritability problem, SNP imputation analysis using a large-scale, whole-genome sequence reference panel has been used as an efficient approach for the supplementation of GWAS data. Particularly for Japanese studies, ethnicity-matched reference data have been required for appropriate SNP imputation analysis because of the higher proportion of variants that are specific to the Japanese population compared with other East Asian populations. In this regard, we have constructed a large Japanese population reference panel based on whole-genome sequences (1KJPN). By using this reference panel, we previously identified primary functional variants in PBC susceptibility gene loci by high-density association mapping.16, 59, 60 Through the application of high-density association mapping based on SNP imputation, using ethnicity-matched reference panels, additional primary functional variants in other digestive disease susceptibility gene loci are likely to be found in future studies.

In conclusion, high-density association mapping based on SNP imputation analysis, together with in silico/in vitro functional analyses, showed the molecular mechanisms by which rs17032850 and rs227361 may contribute to the development of PBC. Altered binding of the transcription factors LEF-1 and RXRα to the enhancer regions, which includes rs17032850 and rs227361, respectively, may affect the transcription of NFKB1 and MANBA. These results indicate that NFkB/LEF-1 and β-mannosidase/RXRα independently have important roles in the development of PBC. In regard to PBC susceptibility genes in the Japanese population, we previously identified the primary functional variants of TNFSF15, PRKCB, and the locus on chromosome 17q12-21, and described their molecular roles.16, 59, 60 Similar approaches for other susceptibility genes are needed to further clarify the molecular mechanisms of disease development.

Materials and Methods

Subjects, Guidelines, and Regulations

Information regarding the participants of this study has been described previously. We obtained written informed consent from all participants, and the study was approved by the committee on research ethics and genetically modified organisms of the Graduate School of Medicine at the University of Tokyo.

Genotype Imputation

The methodology used for SNP filtering and genotype imputation using a phased reference panel of 1KJPN was described in our previous study.

Notably, we applied the variable interval size (depending on the density of SNPs) to the reference panel. Starting from a minimal interval size of 5 Mb, the size was extended until 50 SNPs were included in the interval. The length of the buffer region was set to 1 Mb regardless of interval size.

Databases

The RegulomeDB database, the UCSC genome browser, LD link (https://ldlink.nci.nih.gov/), and NCBI Variation viewer (https://www.ncbi.nlm.nih.gov/variation/view/) were used to evaluate the probability that candidate functional variants might influence transcription regulation, and to obtain LD information. PICS online was used to obtain the probability score of SNPs that are likely to be the causal variants responsible for the association. TRANSFAC Professional was used for the prediction of transcription factor binding. SWISS model, from the University of Basel, and Polyphen2, were used for the prediction of protein structure and the functional consequences of nonsynonymous substitutions.

Gene expression levels in each cell line were examined using data available from GeneCards at the Weizmann Institute of Science (http://www.genecards.org/).

The correlation between rs17032850 or rs227361 genotypes and gene expression was examined using data from the GTEx portal database version 7. P values less than .05, adjusted for multiple testing (Bonferroni correction), were regarded as statistically significant.

EMSA

EMSA was performed using the LightShift Chemiluminescent EMSA Kit (Thermo-Fisher Scientific, Waltham, MA) and biotin-labeled, double-stranded, oligonucleotide probes corresponding to each major and minor allele (Supplementary Table 4), according to the manufacturer’s instructions. These oligonucleotide probes (10 fmol/μL) were incubated with a nuclear extract (2.5 μg/mL) of HepG2 or Jurkat cells (Nuclear Extract Kit; Active Region, Carlsbad, CA) for 30 minutes at 25°C.

The supershift assay was performed by incubating anti–LEF-1 antibody ab53293 (Abcam, Cambridge, UK) and anti-RXRα antibody sc-46659 (Santa Cruz Biotechnology, Santa Cruz, CA) for 24 hours at 4°C with the nuclear extracts, before subsequently incubating these complexes with the biotin-labeled probe for 30 minutes at 25°C. Each assay was performed independently 3 times.

Luciferase Reporter Assay

Specific polymerase chain reaction (PCR) primers (Supplementary Table 5) were used to amplify the 11th intron of NFKB1 and the 7th intron of MANBA from human genomic DNA. Amplicons then were subcloned into the luciferase reporter pGL4.23 (luc2/minP) vector (Promega, Madison, WI). pGL4.23 constructs of each allele (500 ng) and the pGL4.74 (hRluc/TK) vector (50 ng), which was used as an internal control, were transfected into Jurkat and HepG2 cells using Lipofectamine 3000 (Thermo-Fisher Scientific). For the measurement of luciferase activity, the Dual-Luciferase Reporter Assay system (Promega) was used. Differences in relative luciferase activity were compared between major and minor alleles of each SNP using the Student t test. P values less than .05 were regarded as statistically significant. Each figure shows representative data of experiments that were performed independently 3 times. The data in the figures represent the means ± SD of triplicate assays in a single experiment.

Gene Editing Using CRISPR/Cas9

The sequences of the guide RNAs (gRNAs) were designed using CRISPR search and design tools (https://apps.thermofisher.com/apps/crispr/index.html#/search), and gRNAs were synthesized using the GeneArt Precision gRNA Synthesis Kit (Thermo-Fisher Scientific) (Supplementary Table 6). gRNAs, Protein GeneArt Platinum Cas9 Nuclease (Thermo-Fisher Scientific), and donor DNA for each allele were transfected into HepG2 cells using Lipofectamine CRISPR-MAX (Thermo-Fisher Scientific).

After the single-cell cloning and the extraction of genomic DNA using a PureLink Genomic DNA Mini Kit (Thermo-Fisher Scientific), gene editing of the target site and of off-target sites (site with homologies of >75% compared with the target site) were checked by Sanger sequencing (ABI prism 3730-XL, Thermo-Fisher Scientific) using the primer sets shown in Supplementary Table 6.

Quantitative Reverse-Transcription PCR

For Jurkat clones treated either with or without transfection of RXRα small interfering RNA sc-36447 (Santa Cruz Biotechnology) and LEF-1 small interfering RNA sc-35804 (Santa Cruz Biotechnology), total RNA was extracted using a RNeasy kit (Qiagen) and first-strand complementary DNA was synthesized using a High-Capacity Complementary DNA Reverse-Transcription Kit (Thermo-Fisher Scientific). Reverse-transcription PCR to detect each transcript was performed using primers as shown in Supplementary Table 7 and FastStart Taq polymerase (Sigma-Aldrich, St. Louis, MO). To achieve linear amplification, 21 and 27 cycles were found to be optimal in preliminary experiments for GAPDH and other genes, respectively. Quantitation of each transcript was performed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA). These experiments were repeated 3 times with essentially identical results.