GWAS meta-analysis of intrahepatic cholestasis of pregnancy implicates multiple hepatic genes and regulatory elements

Abstract

Intrahepatic cholestasis of pregnancy (ICP) is a pregnancy-specific liver disorder affecting 0.5-2% of pregnancies. The majority of cases present in the third trimester with pruritus, elevated serum bile acids and abnormal serum liver tests. ICP is associated with an increased risk of adverse outcomes, including spontaneous preterm birth and stillbirth. Whilst rare mutations affecting hepatobiliary transporters contribute to the aetiology of ICP, the role of common genetic variation in ICP has not been systematically characterised to date. Here, we perform genome-wide association studies (GWAS) and meta-analyses for ICP across three studies including 1138 cases and 153,642 controls. Eleven loci achieve genome-wide significance and have been further investigated and fine-mapped using functional genomics approaches. Our results pinpoint common sequence variation in liver-enriched genes and liver-specific cis-regulatory elements as contributing mechanisms to ICP susceptibility.

Fig. 1. Manhattan plot for the genome-wide association study meta-analysis of intrahepatic cholestasis of pregnancy (ICP).

Data shown are from the combined meta-analysis including 1138 cases and 153,642 controls, all of European ancestry. The chromosomes are ordered on the x axis; the y axis shows the −log₁₀(P) values for the association tests. Eleven loci achieved genome-wide significance (P < 5 × 10⁻⁸, indicated by the red line). The prioritized gene at each locus is shown. The effector gene functional prioritization strategy followed in this study is presented in Fig. 2. Association testing was performed using a generalized logistic mixed model to account for population stratification and saddle-point approximation to control for type 1 error rates due to unbalanced case–control ratios, followed by meta-analysis weighting the effect size estimates using the inverse of the standard errors.

Fig. 2. Flowchart describing the variant prioritization strategy followed for the eleven loci reaching genome-wide significance in the GWAS meta-analysis of ICP.

Sequential steps used to prioritize the variants are shown with the different approaches for nonsynonymous, synonymous, and noncoding variants indicated.

Fig. 3. Analysis of intrahepatic cholestasis of pregnancy (ICP) association signals uncovers three coding risk variants affecting the SERPINA1, GCKR, and HNF4A genes.

a–c Regional association and linkage disequilibrium (LD) plots for the ICP signals at SERPINA1, GCKR, and HNF4A, respectively. The y axis represents the −log₁₀(P) values for the association test from the GWAS meta-analysis of ICP (1138 cases and 153,642 controls) and the x axis the chromosomal positions (GRCh38/hg38). The SNP with the lowest p-value in the locus is indicated by a purple diamond. The remaining SNPs in the region are colored according to their r² with the lead SNP (r² calculated using 1000 Genomes European (EUR) reference). The dotted gray line represents the genome-wide significance threshold of 5 × 10⁻⁸. Plots were generated using LocusZoom. Association testing was performed as described in Fig. 1. Zoomed insets reveal the prioritized coding variant at each locus, with functional posterior probability (PP_func) above 0.99. PP_func were calculated using PAINTOR. All three variants showed Combined Annotation Dependent Depletion (CADD) scores^, above 10, which is usually considered the threshold for pathogenic variants. Further details on coding variant analysis are provided in Supplementary Data 3. 95% PP ICP variants, variants within the 95% genetic credible set in the ICP meta-analysis. LD ICP variants, variants with EUR r² > 0.8 with the meta-analysis lead variant.

Fig. 4. Three ABC transporter gene loci are associated with intrahepatic cholestasis of pregnancy (ICP).

a, c, e LocusZoom plots showing the association signals at ABCG8, ABCB1/4, and DHRS9, respectively. The LocusZoom plots are as per Fig. 3 and association testing was performed as described in Fig. 1. At each locus, a single noncoding variant was prioritized as likely causal following the strategy outlined in Fig. 2. All epigenomic and transcriptomic datasets shown were retrieved from ENCODE. Adult liver ATAC-seq peaks track (green regions) corresponds to regions accessible in at least two out of four ENCODE adult liver ATAC-seq (see Methods). 95% PP ICP variants, variants within the 95% genetic credible set in the ICP meta-analysis. LD ICP variants, variants with EUR r² > 0.8 with the meta-analysis lead variant. PP_func, functional posterior probability calculated using PAINTOR. CADD, combined annotation-dependent depletion. a Zoomed inset shows a coding variant at ABCG8 overlapping a hepatic transcriptional enhancer, as revealed by strong H3K27ac enrichment and accessible chromatin in adult female liver tissue. Metadome analysis did not suggest that this variant affects protein function (see Supplementary Fig. 7). b Allele-specific effect of the prioritized risk variant rs4148211, assessed by Survey of Regulatory Elements (SuRE) in HepG2 cells. The results suggest that the risk variant associates with decreased transcriptional activity. Mean signal and two-sided Wilcoxon rank-sum test. Previously determined 5% false-discovery rate threshold: P < 0.00173121. Source data are provided as a Source Data file. c, e Zoomed insets show two prioritized regulatory risk variants for ICP: one affecting a hepatic enhancer between ABCB1 and ABCB4 (c), and one affecting a hepatic enhancer in an intron of DHRS9 (d). We note however that DHRS9 is not expressed in liver tissue, as shown in the strand-specific RNA-seq tracks. In contrast, ABCB11 is highly expressed in this tissue making it a likely effector transcript in this locus.

Fig. 5. GWAS meta-analysis of ICP reveals association signals affecting the bile acid homeostasis genes CYP7A1 and SULT2A1.

a, d LocusZoom plots showing the association signals at CYP7A1 and SULT2A1, respectively. The LocusZoom plots are as per Fig. 3 and association testing was performed as described in Fig. 1. All epigenomic and transcriptomic datasets shown were retrieved from ENCODE. Adult liver ATAC-seq peaks track (green regions) corresponds to regions accessible in at least two out of four ENCODE adult liver ATAC-seq (see “Methods”). 95% PP ICP variants, variants within the 95% genetic credible set in the ICP meta-analysis. LD ICP variants, variants with EUR r² > 0.8 with the meta-analysis lead variant. PP_func, functional posterior probability calculated using PAINTOR. CADD, combined annotation-dependent depletion. a Functional fine-mapping with PAINTOR identified a single variant (rs10504255) with a posterior probability over 0.99 overlapping an active liver enhancer located between the CYP7A1 gene, which encodes the rate-limiting enzyme for bile acid synthesis, and UBXN2B. b, e GTEx [GTEx Analysis Release V8 (dbGaP Accession phs000424.v8.p2)] human liver eQTLs identify the top prioritized variant in the CYP7A1 locus (rs10504255) as an eQTL associated with the eGene UBXN2B (b) and show that the SULT2A1 promoter ICP variant (rs296361) is an eQTL for both SULT2A1 and its downstream lncRNA LINC01595 (e). Violin plots represent the density distribution of the samples in each genotype (n for each genotype is indicated below in blue). Box plots show normalized gene expression in median (white line), first and third quartiles. Gene-level adjusted P values calculated with FastQTL are shown. c The risk allele of the ICP-associated variant rs10504255 (G) disrupts a DMRTA1-binding motif. We note that this allele also associates with weaker transactivation in liver tissue, as observed by eQTL analysis (b). d Overlap of ICP risk variants with adult liver accessible chromatin sites revealed one variant in the promoter of SULT2A1 (rs296361). f Transcription factor motif analysis revealed that the ICP variant rs296361 affects a motif for the SOX-D transcriptional repressor family.

Fig. 6. Analysis human liver eQTLs uncovers TMEM147 and ENPP7 as effector transcripts of ICP susceptibility.

a, d LocusZoom plots showing the association signals at TMEM147 and ENPP7, respectively. The LocusZoom plots are as per Fig. 3 and association testing was performed as described in Fig. 1. All epigenomic and transcriptomic datasets shown were retrieved from ENCODE. Adult liver ATAC-seq peaks track (green regions) corresponds to regions accessible in at least two out of four ENCODE adult liver ATAC-seq (see “Methods”). 95% PP ICP variants, variants within the 95% genetic credible set in the ICP meta-analysis. LD ICP variants, variants with EUR r² > 0.8 with the meta-analysis lead variant. PP_func, functional posterior probability calculated using PAINTOR. CADD, combined annotation-dependent depletion. Overlap of ICP-associated variants in these two loci with adult liver accessible chromatin sites enabled the prioritization of two ICP regulatory variants in intronic regions of GAPDHS (rs4806173) (a) and ENPP7 (rs9916601) (d). b, e The ICP variants prioritized with functional liver chromatin annotations are eQTLs in liver tissue and pinpoint TMEM147 and ENPP7 as effector transcripts in these loci. Violin plots represent the density distribution of the samples in each genotype (n for each genotype is indicated below in blue). Box plots show normalized gene expression in median (white line), first and third quartiles. Gene-level adjusted P values calculated with FastQTL are shown. Liver eQTL plots were retrieved from GTEx [GTEx Analysis Release V8 (dbGaP Accession phs000424.v8.p2)]. We note that GAPDHS is a testis-specific gene and its interrogation as an eGene for rs4806173 was therefore not possible (b). c The ICP risk allele of rs4806173 (G), which associates with less transcription of TMEM147, is predicted to create a binding site for glucocorticoid receptors, a family of context-dependent transcriptional regulators that can act as either activators or repressors.

Fig. 7. Tissue-specificity and 3D chromatin conformation analyses prioritize SCARB2 as an ICP susceptibility gene.

a LocusZoom plot showing the association signal at SHROOM3. The LocusZoom plots are as per Fig. 3 and association testing was performed as described in Fig. 1. All epigenomic and transcriptomic datasets shown were retrieved from ENCODE. Adult liver ATAC-seq peaks track (green regions) corresponds to regions accessible in at least two out of four ENCODE adult liver ATAC-seq (see “Methods”). 95% PP ICP variants, variants within the 95% genetic credible set in the ICP meta-analysis. LD ICP variants, variants with EUR r² > 0.8 with the meta-analysis lead variant. CADD, combined annotation-dependent depletion. Overlap of ICP risk variants in this locus with adult liver accessible chromatin sites enabled the prioritization of a single ICP-regulatory variant in an intronic region of the SHROOM3 gene (rs4859682). None of the other ICP variants (EUR r² > 0.8 with lead variant or in the 95% PP credible set) overlapped accessible chromatin sites. b Expression of all genes within the TAD containing the ICP lead variant rs13146355 in the human liver. Expression data was retrieved from GTEx [GTEx Analysis Release V8 (dbGaP Accession phs000424.v8.p2)]. Violin plots represent the density distribution of expression of the indicated genes in human liver tissue (n = 226 donors). Box plots show gene expression in median (white line), first and third quartiles. TPM, transcripts per million. Note that the gene harboring the ICP risk variant, SHROOM3, is lowly expressed in the liver. c Scatter plot showing the expression of all genes within the TAD containing the ICP lead variant rs13146355 in the human liver versus their liver-specificity. Tissue-specificity scores (SPECS) were precalculated by Everaert et al.. Source data are provided as a Source Data file. This plot highlights two highly expressed and liver-specific genes, STBD1, and SCARB2. d, Analysis of long-range chromatin interactions in human liver detected by promoter capture Hi-C, showing that the genomic region containing the ICP risk variant rs4859682 interacts with SCARB2, but not SHROOM3. Long-range chromatin interactions stemming from the promoter of STBD1 were not detected in liver tissue.