We used COLOC (Giambartolomei et al., 2014) and expression data from The Genotype Tissue Expression (GTEx) project (Carithers and Moore, 2015) to identify cis-eQTL that colocalize with SUA1 and SUA2. eQTL for the genes MAFTRR and LINC01229 (Figure 2, gene body depicted in black) colocalize with the SUA1 urate signal (Table 2) although the lead eQTL variants are not the lead urate variants at SUA1 (Supplementary Figure S2). SUA1 urate-raising alleles associate with increased expression of MAFTRR but with lowered expression of LINC01229 (e.g., rs7188445_G Supplementary Figure S3 and Supplementary Table S1). An eQTL for LINC01229 in the GTEx tissue Skin – Exposed (Lower leg) (Supplementary Figure S2) colocalises with the SUA2 urate signal (Table 2). Urate-raising variants within SUA2 increase LINC01229 expression (Supplementary Table S2). Colocalization of the SUA1 and SUA2 urate signals with MAFTRR and LINC01229 eQTL suggests that the causal variant(s) at these regions for serum urate and lincRNA expression are shared. Collectively these results are consistent with expression of MAFTRR and LINC01229 being functionally important in serum urate control.

MAFTRR eQTL variants (all GTEx tissues) that have the greatest effect on MAFTRR expression lie immediately 3′, or within the 3′ region of the MAFTRR transcript (Supplementary Figure S2). By contrast, the locations of eQTL for LINC01229 differ in genomic location between tissues (Supplementary Figure S2). eQTL for MAF are not found within SUA1 or SUA2 and instead lie within the WWOX transcript (+375 kb of MAF), downstream and within the MAF transcript, and within the promoter of DYNLRB2 (-1 Mb of MAF) (Supplementary Figure S4).

GTEx currently has no renal eQTL data, therefore we queried the NepheQTL database (Gillies et al., 2018) to assess whether the SUA1 and SUA2 MAFTRR and LINC01229 eQTL exist in the kidney tubule. MAFTRR has a very strong eQTL signal (FDR adjusted p = 8.12 × 10^-18) in the tubulointerstitium (Supplementary Figure S2) that coincides with SUA1. Four maximally associated variants (R²> 0.98) lie immediately 3′ to the MAFTRR transcript consistent with the GTEx data. COLOC indicates that the eQTL signal for MAFTRR and the SUA1 serum urate signal are different (H3 posterior probability = 0.976). Nonetheless, SUA1 urate-raising variants associate with increased expression of MAFTRR in kidney tubule (Supplementary Table S1). No SUA2 variants were eQTL for MAFTRR in the kidney tubules within the NepheQTL database (Gillies et al., 2018) (Supplementary Table S2). Furthermore, there were no significant gene eQTL signals for MAF or LINC01229 within NepheQTL.

An additional whole kidney-derived eQTL dataset (Ko et al., 2017) and Japanese blood eQTL sample set¹ [monocytes, CD4+, CD8+, natural killer and B cells (Kanai et al., 2018)] were queried. eQTL for MAF were located downstream of MAF (consistent with the GTEx data), however, there were no eQTL data for MAFTRR and LINC01229 and no SUA (SUA1-4) variants were significant eQTL in these tissues (Supplementary Figure S5).

Intergenic disease-associated variants could influence disease traits by regulating gene expression. We used PAINTOR (Kichaev et al., 2014) in conjunction with ENCODE (Encyclopedia of DNA Elements) (ENCODE Project Consortium, 2012) to identify functional variants with regulatory activity at the SUA regions upstream of MAF. Fine mapping with PAINTOR using kidney cell type annotations (Finucane et al., 2015) as a prior identified the most likely causal variants within SUA1 and SUA2 (Table 3). We interrogated ENCODE to identify putative regulatory elements that overlap the SNPs identified by PAINTOR (Figure 2, Supplementary Tables S1, S2, and Supplementary Figures S6, S7). Functional annotations of the human genome were scanned to identify DNaseI hypersensitivity regions [DHS; characteristic of enhanced chromatin accessibility (Thurman et al., 2012)], the coexistence of histone H3 monomethylated on lysine 4 and histone H3 acetylated on lysine 27 (H3K4me1, H3K27ac; hallmarks of active cis-regulatory elements) (Heintzman et al., 2009) and transcription factor binding sites (Wang et al., 2012, 2013). We also looked for hyper-accelerated regions (HARs), which represent conserved genomic loci with elevated divergence in humans (Doan et al., 2016) and can therefore indicate potential human-specific enhancers.

Our ENCODE analyses indicate that rs58722186 within SUA1 identified by PAINTOR (Table 3, Figures 2A,B, PAINTOR posterior probability = 0.31) is in close proximity to regulatory elements. rs58722186 is adjacent to five predicted transcription factor binding sites based on chromatin immunoprecipitation (ChIP) studies from HEPG2 and K652 cells, falls within a GATA2 and GATA3 binding site and is bound by GATA2 in SH-SY5Y cells (ENCODE; Figure 2B and Supplementary Figure S6). rs58722186 is also flanked by DHS peaks and enriched for histone 3 monomethylation on lysine 4 (H3K4me1) (Figure 2B and Supplementary Figure S6). rs58722186 is one of the maximally associated serum urate variants (p = 7.10 × 10^-13) in the Kanai et al. (2018) dataset. rs58722186 is an eQTL for MAFTRR in multiple tissues with a similar effect on expression as the lead SUA1 urate SNP rs7188445 (p = 1.3 × 10^-36 and p = 9.6 × 10^-37, respectively). Chromatin state segmentation by hidden Markov model (ChromHMM) (Ernst and Kellis, 2017) indicates that rs58722186 is within heterochromatin (Figure 2B and Supplementary Figure S6) but is in proximity to a weak enhancer in K652 cells. We find minimal evidence for regulatory elements that overlap urate-associated SNPs identified in the European dataset in SUA1 (Supplementary Figure S2B and Supplementary Table S1) consistent with our PAINTOR analyses (Supplementary Table S3).

At SUA2 the lead PAINTOR SNP rs192158533 (Table 3, Figures 2A,C, PAINTOR posterior probability = 0.89) does not mark any significant regulatory elements. rs192158533 is monomorphic in East Asian populations yet there is a strong urate association at SUA2 signal in Kanai et al. (2018) data. The second most likely SUA2 PAINTOR SNP (Table 3, Figures 2A,C, PAINTOR posterior probability = 0.27) is rs4077450 and the maximally associated urate SNP in the Kanai et al. (2018) Japanese data (p = 3.9 × 10^-14).

ENCODE identifies that rs4077450 along with SUA2 variant rs4077451 encompass a region that is enriched for H3K4me1 and H3K27ac (Figure 2C and Supplementary Figure S7) (Heintzman et al., 2009). DHS peaks overlap this region in 33 cell lines, including HRCEpiC renal cortex cells (Figure 2C, Supplementary Figure S7, and Supplementary Table S4) (Thurman et al., 2012). The ENCODE data identifies ∼30 predicted transcription factor binding sites within the rs4077450_rs4077451 SNP region in chromatin immunoprecipitation (ChIP) studies in HEPG2 and K652 cells (Figure 2C and Supplementary Figure S7) (Gerstein et al., 2012; Wang et al., 2012, 2013).

ChromHMM (Ernst and Kellis, 2017) predicts that the rs4077450_rs4077451 SNP region marks a strong enhancer in liver carcinoma HEPG2 and leukemia K562 cells (Figure 2C and Supplementary Figure S7). In these cells there is also evidence for strong active transcription signatures at the MAFTRR/LINC01229 TSS but the MAF TSS is marked by repressed chromatin signatures. Our eQTL analyses indicate that SUA2 SNPs rs4077450 and rs4077451 are within the colocalized LINC01229 eQTL (p = 1.2 × 10^-5 and p = 2.0 × 10^-5, respectively, Supplementary Figures S8A,B) and the urate lowering alleles associate with lowered LINC01229 expression (Table 2 and Supplementary Table S2).

Conserved non-coding regions of the genome might indicate developmental enhancers (Polychronopoulos et al., 2017). LINSIGHT (Huang et al., 2017) indicates that the SUA2 SNP region is under sequence constraint (Supplementary Figure S9). SUA2 also encompasses a HAR (Bird et al., 2007) which is located ∼2 kb upstream of the rs4077450 SNP region at Chr16: 79934248–79934397 (Supplementary Figure S9). In combination, our in silico analyses indicate that the SUA2 SNP region is a strong candidate for enhancer function and follow up in functional assays.

We performed transient enhancer assays in zebrafish embryos to investigate the ability of the SUA2 region to act as an enhancer in vivo. The region encompassing the DHS and TFBS hotspot and spanning rs4077450 and rs4077451 was amplified from 1000 Genome DNA samples heterozygous for each of the allele variants (Figure 2C). Allele fragments containing rs4077450_T-rs4077451_A and rs4077450_G-rs4077451_T are herein named major and minor allele fragments, respectively. The minor alleles associate with increased urate levels and increased LINC01229 expression.

The major and minor allele fragments of the SUA2 region are both capable of driving kidney-specific GFP expression in zebrafish larvae (Figure 3). Cells that are destined to form the zebrafish pronephric duct arise from the ventral mesoderm and then populate the intermediate mesoderm where the differentiation of the pronephric duct into the glomerulus, proximal tubule and distal tubule begins (Horsfield et al., 2002; Drummond, 2003). GFP was initially observed in the zebrafish kidney proximal tubule at 24 h post-fertilization (hpf) (Supplementary Figure S10). By 3 days post-fertilization (dpf), the pronephric duct can be clearly divided into four sections; the glomerulus, the proximal convoluted tubule (PCT), the proximal straight tubule (PST) and the distal tubule (DT) (Figure 3A) (Wingert et al., 2007). At 3 dpf, the presence of mosaic GFP showed that the enhancer is active in the PCT, PST and DT (Figure 3B). By 5 dpf GFP kidney expression was confined to the PCT (Figure 3C).

The extent of GFP localization in the PCT, PST, and DT was quantified in 3 dpf embryos. The minor allele fragment drove GFP expression in the cells of the PCT in 34% of embryos, the PST in 26% of embryos and the DT in 13% of embryos at 3 dpf. The major allele fragment drove GFP expression in the PCT in 34% of embryos, PST in 22% of embryos and DT in 5% of embryos at 3 dpf (Figure 3D). The percentage of embryos with GFP in the DT between the major and minor allele fragments was significantly different (Fishers exact test, p = 0.0025) indicating that minor (urate-raising) allele fragment has subtly greater spatial enhancer activity during development.

The SUA2 major and minor allele fragments were cloned into luciferase reporter constructs and evaluated for differential enhancer activity in the human embryonic kidney cell line HEK293. The major and minor allele fragments both exhibited significant (p = 0.0216 and p = 0.0002, respectively) enhancer activity in HEK293 cells when compared to the vector only control (Figure 3E). Notably, the minor urate-raising allele fragment had elevated enhancer activity when compared to the major allele (p = 0.0032).

Collectively, our results indicate that the genomic region containing the rs4077450_rs4077451 SNPs is a kidney-specific regulatory element (henceforth called the SUA2 enhancer) that shows allele-specific enhancer activity. The minor (urate-raising) alleles increased the range of tissue expression in zebrafish, and also increased the amplitude of expression in HEK293 cells, providing the first indication of a mechanism for influence of SUA2-specific genetic variation on serum urate levels.

CoDeS3D which leverages genome connectivity (Hi-C datasets) and gene expression associations [eQTL data from the GTEx catalog (Carithers and Moore, 2015)] provides evidence that variants within the SUA2 enhancer regulate the LINC01229 locus via long-range interactions. Genome connectivity data from Hi-C datasets for K562, NHEK, IMR90, HMEC, GM12878, and KBM7 cells (Rao et al., 2014) shows that an interaction fragment containing rs4077450 and rs4077451 within SUA2 physically interacts with 45 restriction fragments that fall within a ∼130 kb interaction region [Figure 3F, Chr16: 79712175–79841647 (build 37.7) FDR < 0.02, Supplementary Table S5]. This interaction region includes the transcripts for the lincRNA genes LINC01229, LINC01228, and MAFTRR and a large part of SUA1 (Figure 3F). Hi-C data from renal cells (Caki-2) (Yang et al., 2018) provides further confirmation that a restriction fragment containing rs4077450 interacts with restriction fragments containing rs7188445 (SUA1) and rs889472 (SUA3) (Supplementary Figure S11A). Hi-C data also shows that topologically associated domains (TADs) at the MAF locus shift in location depending on cell/tissue type. In NHEK cells MAF is found on the TAD boundary (Supplementary Figure S11B) and in HUVEC, IMR90 and KBM7 cells MAF shares a TAD with SUA1, MAFTRR, LINC01229, and SUA2 (Supplementary Figure S11B).

We used in situ hybridization to examine the spatial location of Danio rerio maf transcripts in zebrafish embryos and assess whether maf expression overlaps with the activity of the kidney specific enhancer. We found that maf is expressed in the developing kidney and proximal tubules at 12 hpf, 24 hpf, and 48 hpf (Figures 4A–C). GFP expression driven by the SUA2 enhancer in the proximal tubules coincides with the developmental expression of zebrafish maf at 48 hpf (Figure 4D).

The HNF4α transcription factor is robustly expressed in the human kidney proximal tubule² (Ponten et al., 2008) and regulates the expression of genes encoding proteins that control serum urate (Prestin et al., 2014), including PDZK1 (Ketharnathan et al., 2018). Furthermore, variation at the HNF4A locus is associated with serum urate levels (P∼10^-6) (Köttgen et al., 2013). SUA2 variants rs4077450 and rs4077451 flank a core HNF4α consensus motif (Supplementary Figure S7) (Wang et al., 2012; Huang et al., 2017) and fall inside an HNF4α ChIP-seq peak in HEPG2 cells (Supplementary Figure S7) (Gerstein et al., 2012; Wang et al., 2012, 2013). We conducted ChIP-qPCR using an anti-HNF4α antibody and confirmed that the SUA2 enhancer is physically bound by HNF4α in HEPG2 cells (Figure 4E). HNF4α binding was not assessed in HEK293 cells because HNF4A is not expressed in this cell line (Ketharnathan et al., 2018). GFP expression driven by the SUA2 enhancer also coincides with hnf4a expression in the zebrafish pronephric duct at these stages (Wingert and Davidson, 2011). Combined, our analyses of the SUA2 kidney enhancer indicate that it functionally and physically connects with the MAFTRR/LINC01229 lincRNA region, is likely responsive to HNF4α, and could contribute to MAF expression. However, several other transcription factor consensus motifs, including those for TFAP2A and HNF4G, are also in close proximity to SUA2 variants rs4077450 and rs4077451, therefore it is possible that the variants could influence the binding of one or more transcription factors.

CoDeS3D and Hi-C data for K562, NHEK, IMR90, HMEC, GM12878 and KBM7 cells provides spatial data supporting the rs7188445 MAFTRR cis-eQTL (Supplementary Table S6). Hi-C data also provide evidence that the lead MAFTRR eQTL SNPs which lie within the MAFTRR and LINC01229 transcripts physically connect with MAF, however, there are no significant MAF eQTLs associated with these connections (Supplementary Table S7).

We also tested whether genes in trans have functional genomic regulatory connections at SUA1 (Rinn and Chang, 2012). CoDeS3D identified 12 spatially supported trans-eQTL with SUA1 variants located within the MAFTRR and LINC01229 transcripts (FDR < 0.05, Supplementary Table S8) and also with the lead SUA1 SNP rs7188445 (Supplementary Table S6). The lead SNPs for these trans-eQTL at SUA1 are different than the lead cis-eQTL SNPs for MAFTRR (Supplementary Figures S2, S12 and Figure 5). Importantly COLOC analysis (Giambartolomei et al., 2014) indicates that these trans-eQTL signals colocalize with the SUA1 serum urate signal (Table 4) indicating that SUA1 could control serum urate levels via trans-eQTL connections. Of the 12 trans-eQTL, only 2 did not colocalize (H4 < 0.5) with the SUA1 signal (CCDC79 and TBC1D10B). Of the 10 remaining trans-eQTL, four genes (i.e., SLC5A8, CCDC6, EHHADH, and DLGAP1) are potentially involved in urate metabolism based on our Gene Ontology (GO) analyses³ and previous urate associations (Köttgen et al., 2013; Kanai et al., 2018) (Figure 5).

CCDC6 is in very close proximity to a urate GWAS signal that sits in the intergenic region between CCDC6 and SLC16A9 (Supplementary Figure S13A). CCDC6 and SLC16A9 are strongly co-regulated (p = 6.9 × 10^-9, GeneNetwork.nl) and the lead Köttgen SNP at CCDC6/SLC16A9 (rs1171614) is an eQTL for both SLC16A9 and CCDC6 [GTEx (Carithers and Moore, 2015)]. SUA1 lead urate-raising variant rs7188445_G associates with lowered expression of CCDC6 (Supplementary Figure S13B). At SLC5A8 there is a weak urate GWAS signal in the Japanese population (Kanai et al., 2018) (Supplementary Figure S13A). SLC5A8 is a monocarboxylate transporter that has been strongly implicated in urate transport in the kidney (Merriman, 2015) and associates with the GO biological function term urate transport (p = 6.3 × 10^-3) (Supplementary Table S9). The urate-raising allele of rs7188445 associates with increased expression of SLC5A8 (Supplementary Figure S13B). The only other genes with a trans-eQTL at SUA1 with evidence for a urate signal was DLGAP1 in the Japanese population (Supplementary Figure S13A). The serum urate-raising and gout risk allele rs7188445_G associates with increased expression of DLGAP1 (Supplementary Figure S13B). The proteins encoded by both DLGAP1 and SLC5A8 contain PDZ3 domains and interact with PDZK1 (Hu et al., 2009), a scaffold protein involved in urate reabsorption within the kidney tubule (Merriman, 2015).

Of the 10 genes with colocalized SUA1 trans-eQTL, EHHADH has the strongest signal in kidney tubule (gene level FDR = 0.078) (Gillies et al., 2018). Expression of EHHADH is enriched in the proximal tubule (Habuka et al., 2014) and associates with the GO biological function term urate transport (p = 7.7 × 10^-3) in the GeneNetwork database (Supplementary Table S9). The serum urate-raising and gout risk allele rs7188445_G associates with increased MAFTRR expression and lowered EHHADH expression (Supplementary Figure S13B). These analyses collectively indicate that SUA1 variants associate with the expression of serum urate and kidney relevant genes in trans.

Our analyses indicate that LINC01229 and MAFTRR, along with genes in trans, likely have functions in serum urate regulation. In addition, the SUA2 enhancer recruits HNF4α and coincides with the expression of zebrafish hnf4a and maf, implicating MAF in serum urate regulation. Supporting this co-expression analysis the GeneNetwork database (Kumar et al., 2013) indicates that MAFTRR and LINC01229 are strongly co-expressed with MAF (p = 5.0 × 10^-27 and p = 2.4 × 10^-26, respectively). However, MAFTRR and LINC01229 are not co-expressed (Supplementary Tables S10, S11). Interestingly, HNF4A is co-expressed with LINC01229 (p = 0.011) but not with MAFTRR (Supplementary Tables S10, S11) consistent with the SUA2 variant effect on LINC01229 expression. We also assessed whether the SUA1 trans-eQTL genes are co-expressed with MAFTRR, MAF and/or LINC01229. SLCO3A1, TNKS and ENTPD3 (the likely target of ENTPD3-AS1) are co-regulated with MAF (p = 0.0088, p = 0.011, and p = 0.022, respectively) (Supplementary Table S12) and rs7188445_G associates with increased SLCO3A1, TNKS, and ENTPD3 expression (Supplementary Figure S13B).

To test if MAFTRR and LINC01229 can contribute to the regulation of MAF transcription, we depleted MAFTRR and LINC01229 using siRNA in HEK293 cells. Reduction of MAFTRR RNA levels (Figure 6A) resulted in upregulation of endogenous MAF mRNA (Figure 6B). Reduction of LINC01229 RNA levels (Figure 6D) also led to upregulation of endogenous MAF mRNA (Figure 6E). Our results indicate that LINC01229 and MAFTRR normally repress MAF expression in kidney HEK293 cells.

Our analyses indicate that LINC01229 and MAFTRR have the ability to regulate MAF in kidney cells. We assigned putative phenotypes and biological processes to LINC01229 and MAFTRR using GeneNetwork (Kumar et al., 2013), which infers functional enrichment from co-expressed genes⁴. The top enriched human phenotype ontology terms for MAFTRR and LINC01229 are ‘tubulointerstitial abnormality’ (p = 3.9 × 10^-4) and ‘accelerated skeletal maturation’ (p = 4.3 × 10^-5), respectively. The top GO biological process term for MAFTRR is ‘chondroitin sulfate catabolic process’ (p = 8.4 × 10^-3). The top GO biological process term for LINC01229 ‘cell development’ (p = 7.9 × 10^-5) implicates LINC01229 in a pathway with MAF and three other serum urate-associated loci, INHBB, INHBC and INHBE, detected by GWAS (Köttgen et al., 2013) (Supplementary Figures S14A,B). Additionally, HNF4A and LINC01229 share GO enrichment in lipid-related pathways including ‘chylomicron assembly,’ ‘phospholipid efflux,’ and ‘high-density lipoprotein assembly’ (Supplementary Figure S14A). LINC01229 also associates with inflammation- and gout-relevant GO terms including ‘interleukin β secretion,’ ‘T-cell differentiation,’ and ‘interleukin 2 biosynthetic process’ (Supplementary Figure S14A).

The intergenic region upstream of MAF encompassing the lincRNAs is conserved in mammals. Two non-coding RNA loci upstream of MAF in mouse (i.e., Gm30925 and Gm31037) are positionally conserved with respect to MAF, expressed in kidney and homologous with the human lincRNA region (BLAST E-value = 3 × 10^-59 and 2 × 10^-16, respectively) (Supplementary Figure S15). Moreover, sequence homologous (E-value = 6 × 10^-20) to the SUA2 enhancer in mouse is located 276,670 bp upstream of Maf. Our analyses suggest that the ∼300 kb MAF upstream intergenic region marked by SUA1 and SUA2 represents an evolutionarily conserved regulatory network active in kidney and serum urate regulation.

Publicly available summary level statistics from the Köttgen et al. (2013) GWAS for serum urate were used. Conditional association analysis on rs7188445 was performed using the Genome-Wide Complex Trait Analysis (GCTA) software (Yang et al., 2011). 6,654 HapMap2 imputed genotypes from European participants of the Atherosclerosis Risk in Communities study were used as a reference. LD R² values were calculated using LDlink (Machiela and Chanock, 2015) with the 1000 Genomes European and East Asian reference panels.

ImpG (v1.0) (Pasaniuc et al., 2014) was used to impute Z-scores into the European urate GWAS (Köttgen et al., 2013) summary statistics. For the reference haplotypes, the latest release of the 1000 Genomes project was used, and only bi-allelic SNP markers having a minor allele frequency greater than 0.01 in the relevant population were included. All imputed markers with a predicted LD R² of less than 0.8 were removed. PAINTOR leverages functional genomic annotation data, in addition to genetic association strength, to prioritize functional variants. The LD matrix used for fine mapping of the Köttgen GWAS data in PAINTOR was calculated based on European data from the 1000 Genomes Project Phase 1. The annotation matrix used in PAINTOR was kidney cell-type DHS (Finucane et al., 2015). PAINTOR (v3.0) (Kichaev et al., 2014) was carried out using a LD matrix matched locus file for SUA1 and SUA2.

We used COLOC (Giambartolomei et al., 2014) to assess the similarity between the urate-associated loci SUA1 and SUA2 and publicly available eQTL data from the Genotype Tissue Expression Project (GTEx). COLOC is a Bayesian method that compares four different statistical models at a locus and gives a posterior probability for each test. The four tests are: H0: no causal variant in the GWAS or an eQTL region H1: a casual variant in either the GWAS or the eQTL region but not both, H3: a causal variant in the GWAS and the eQTL region that is different, H4: a causal variant in the GWAS and the eQTL region that is shared.

The Contextualize Developmental SNPs using 3D Information (CoDeS3D) algorithm (GitHub⁵) (Fadason et al., 2017) was used to identify long-distance regulatory relationships for rs4077450 and rs4077451 (SUA2), rs7188445 (SUA1) and SUA1 SNPs that overlap the MAFTRR transcript. This analysis leverages known spatial associations from Hi-C databases (Rao et al., 2014) and gene expression associations [eQTL data from the GTEx catalog (GTEx 2013)] to assess regulatory connections. Briefly, SNPs were mapped onto Hi-C restriction fragments, the genes that physically interact with these restriction fragments identified and collated (SNP-gene spatial pairs). SNP-gene pairs were screened through GTEx to identify eQTLs. The false discovery rate (FDR) was calculated using a stepwise Benjamini-Hochberg correction procedure and incorporated the number of tests and eQTL value list. A FDR value of <0.05 was accepted as statistically significant (Fadason et al., 2017).

A DNA fragment (562 bp) harboring rs4077450_rs4077451 (Chr16: 79931421–79931980) was amplified from the genomic DNA of individuals who were heterozygous for the major alleles (rs4077450_T-rs4077451_A) and the minor alleles (rs4077450_G-rs4077451_T), respectively. Primers F: 5′-CCTCCATACAGTGTCCAGCA-3′ and R: 5′-TGGACCGTTTTGGCTTTTAC-3′ were used for amplification. Amplicons were cloned into the pCR^®8/GW/TOPO (Invitrogen) entry vector and gateway-cloned into the destination vectors pGL4.23 (addgene 60323) and ZED (Bessa et al., 2009) using LR Clonase^® enzyme (Invitrogen) for the enhancer assays in human cell lines and zebrafish, respectively.

Human embryonic kidney HEK293 cells (ATCC) and human hepatocellular carcinoma HEPG2 cells (ATCC) were grown (37°C in a 5% CO₂) and maintained in Dulbecco’s Modified Eagle’s Medium (DMEM) and Eagle’s Minimum Essential Media (EMEM), respectively, (Life Technologies), supplemented with 10% fetal bovine serum.

For luciferase assays, HEK293 cells were transiently transfected with pGL4.23 constructs and Renilla, using Lipofectamine 3000 (Invitrogen). Luminescence was measured 48 h post-transfection using the Dual Glo Luciferase Assay (Promega), on the Perkin Elmer Victor X4 plate reader. Values were normalized to the expression of Renilla luciferase.

For the siRNA knockdown experiments, HEK293 cells were reverse-transfected with ON-TARGETplus Non-targeting Control Pool siRNA (5 nM) (Dharmacon Cat# D-001810-10), Lincode SMARTpool LINC01229 siRNA (5 nM) (Dharmacon Cat# R-191862-00-0005), and Lincode SMARTpool MAFTRR (5 nM) (Dharmacon Cat# R-192817-00-0005) using RNAiMax (Invitrogen) for gene expression analysis.

Total RNA was isolated from control and siRNA-treated HEK293 cells at 48 h post-treatment using the NucleoSpin RNA kit (Macherey-Nagel). cDNA was synthesized with qScript cDNA SuperMix (Quanta Biosciences). MAF, MAFTRR and LINC01229 expression were measured using TaKaRa SYBR Premix Ex Taq^TM (Clontech) on a LightCycler400 (Roche Diagnostics). Primers MAF F: 5′-AGCAAGTCGACCTCAAG-3′ and R: 5′-CGAGTGGGCTCAGTTATGAAA-3′, MAFTRR F: 5′-CCTGGACAATGCTGGTTTTT-3′ and R: 5′-CACGTCCTTCCATTTTGCTT-3′ and LINC01229 F: 5′-ATGGGAGCTCCACACAGGT-3′ and R: 5′-TGGGTGCCTTTAAACAAGAGA-3′ were used for amplification. Gene expression analyses were carried out on qBase Plus (Biogazelle) and were normalized relative to the mean of reference genes encoding glyceraldehyde 3-phosphate dehydrogenase (GAPDH) and ribosomal protein L13a (RPL13A).

Chromatin was extracted from 10^-7 HEPG2 cells, sonicated (Vibra Cell VCX130 Sonicator, Sonics) to approximately 500 bp (determined experimentally) and diluted 10 times with immunoprecipitation (IP) buffer as described previously (Veto et al., 2017). Chromatin immunoprecipitation (ChIP) was carried out according to Vaisanen et al. (2005) with minor adjustments. Equal amounts of diluted chromatin in 2 ml IP buffer were pre-cleared (to reduce background) overnight with Dynabeads protein G (Thermo Fisher). Immunoprecipitations were performed overnight at 4°C with 10 μl of a ChIP grade HNF4α (ab41898, Abcam) antibody conjugated to 50 μl of Dynabeads G. The immunocomplexes were then pelleted by centrifugation at 4°C for 1 min at 100 × g and washed sequentially for 5 min by rotation with 1 ml of the following buffers: low-salt wash buffer [20 mM Tris–HCl (pH 8.1), 0.1% SDS, 1% Triton X-100, 2 mM EDTA, 150 mM NaCl], high-salt wash buffer [20 mM Tris–HCl (pH 8.1), 0.1% SDS, 1% Triton X-100, 2 mM EDTA, 500 mM NaCl] and LiCl wash buffer [10 mM Tris–HCl (pH 8.1), 0.25 mM LiCl, 1% (v/v) Nonidet P-40, 1% (w/v) sodium deoxycholate, 1 mM EDTA]. Finally, the beads were washed twice with 1 ml of TE buffer [10 mM Tris–HCl (pH 8.0), 1 mM EDTA] and eluted in ChIP elution buffer (50 mM Tris, pH 8, 10 mM EDTA, 1% SDS) at room temperature for 30 min with rotation. Crosslinking was reversed with Proteinase K (final concentration 40 μg/ml) at 50°C for 2 h, DNA was recovered using phenol:chloroform:isoamyl alcohol (Invitrogen) and precipitated with 0.1 volume of NaCl and 2.5 volumes of ethanol using Ambion^® Linear Acrylamide (5 mg/mL) as a carrier. For qPCR analyses, 1 μl of pre-cleared or immunoprecipitated chromatin was used for each reaction. HNF4α binding at the SNP region, negative binding site B-globin and positive binding site ABCC6 (Veto et al., 2017) was expressed relative to the pre-cleared input chromatin after subtraction from no antibody control. Primers SNP region F: 5′-CCTCCATACAGTGTCCAGCA-3′ and R: 5′-TGGACCGTTTTGGCTTTTAC-3′, B-globin F: 5′-AGGACAGGTACGGCTGTCATC-3′ and R: 5′-TTTATGCCCAGCCCTGGCTC-3′ and ABCC6 F: 5′-AGCCCATTGCATAATCTTCTAAGT-3′ and R: 5′-ATGGAGACCGCGTCACAG-3′ were used for amplification.

Wild type (WIK) fish lines were maintained according to established protocols (Westerfield, 2000). The University of Otago Animal Ethics Committee approved all zebrafish research.

An enhancer test vector ZED (Bessa et al., 2009) was used to investigate the enhancer capacity of the genomic region (Chr16: 79931421–79931980) marked by rs4077450 and rs4077451. Single-cell WIK embryos were injected with 1 nl of 30 ng/μl ZED (Bessa et al., 2009) DNA construct containing either of the minor and major fragments for the rs4077450_rs4077451 SNP region and 90 ng/μl Tol2 transposase mRNA (Kawakami, 2005). Injected embryos were screened for Green Fluorescent Protein (GFP) expression 24–120 h post-fertilization (hpf) using a Leica M205 FA fluorescence microscope.

Full length maf riboprobe synthesis and whole mount in situ hybridization was performed as described previously (Monnich et al., 2009). The cDNA of MAF on a pCMV SPORT 6.1 plasmid was acquired from Life Biosciences as the Danio rerio maf cDNA DR.81288. The plasmid was linearized with ApaI and the antisense DIG-labeled riboprobe was generated by in vitro transcription using SP6 RNA polymerase (Roche).

GraphPad PRISM 7 (GraphPad software, San Diego, CA, United States) was used for performing all statistical analysis. One-way ANOVA, (Tukey’s multiple comparisons tests), Fisher exact tests and unpaired t-tests were used for estimating the statistical significance of the luciferase, enhancer assays and quantitative PCR data, respectively. All data are presented as mean and error bars represent standard error of the mean (SEM).