Bilirubin is the major metabolite of heme, the iron-binding tetrapyrrole ring found in hemoglobin, myoglobin, and cytochromes.^[@R1]^ The straight-chain compound biliverdin is produced through the oxidation of heme porphyrin ring by microsomal heme oxygenase. Biliverdine reductase, a nicotinamide adenine dinucleotide phosphate (NADPH)-dependent enzyme, reduces biliverdin to produce bilirubin.^[@R1]^ After being captured by the hepatocyte through its membrane surface in contact with the sinusoids, bilirubin is transported to the smooth endoplasmic reticulum and becomes the substrate of the UDP-glucuronosyltransferase 1 family, polypeptide A1 enzyme (UGT1A1), which catalyzes the esterification of the propionic acid side chains of bilirubin with glucuronic acid (present as uridine diphosphoglucuronic acid) to form mainly the diglucuronide conjugate, a water-soluble conjugated molecule.^[@R1]^ Conjugated bilirubin is then actively secreted into the bile canaliculi by a membrane ATP-dependent transporter, designated as multidrug resistance-associated protein 2 (MRP2).^[@R2]^

It is now well established that unconjugated hyperbilirubinemia is a risk factor for gallstones.^[@R3],[@R4]^ In patients with sickle cell disease, it has been shown that genetic variation in the promoter of *UGT1A1* may be a risk factor for symptomatic gallstones in older people.^[@R5]^ However, no studies are available in the healthy population, especially among older subjects.

Genome-wide association studies (GWASs) systematically evaluate common genetic variants, typically with a minor allele frequency (MAF) \>5% and have been extensively used to dissect the genetic architecture of complex diseases and quantitative traits.^[@R6]^ The high number of genetic variations identified in GWASs contrasts with their low effect on disease risk or quantitative trait variation. Thus, GWASs generally fail to translate into functional understanding or clinical practice. The "missing heritability" observed in GWASs could be explained by the fact that they do not assess low-frequency (MAF, 0.5%--5%) and rare (MAF \< 0.5%) genetic variants that play a major role in human pathology.^[@R6]^ Recent evidence suggest that low-frequency and rare variants are associated with complex diseases.^[@R6]--[@R8]^

It is estimated that the protein-coding regions of the human genome constitute about 85% of the disease-causing mutations.^[@R9]^ The Illumina Human Exome BeadChip provides coverage of functional exonic variants. Around 250,00 markers on this BeadChip represent SNPs in *RefSeq* genes, non-synonymous SNPs, and SNPs in coding regions (including untranslated regions, UTRs). It has been demonstrated that exome-wide genotyping identified additional medically actionable variant calls and helped resolve ambiguous single-nucleotide variants in comparison with genome-wide approaches.^[@R10]^

To date, there are no available data regarding the exome-wide association study approach that evaluated potential-associated variants with unconjugated, conjugated, and total serum bilirubin levels in physiological conditions. To determine the role of rare and low-frequency coding variants in traits reflecting unconjugated, conjugated, and total serum bilirubin level, we evaluated putative functional coding variants using the Illumina HumanExome BeadChip on a well-characterized cohort of ambulatory elderly subjects from Italy.^[@R11]--[@R13]^ Furthermore, we assessed the potential influence of bilirubin-related variants on hyperbilirubinemia risk and the risk of gallstone-related cholecystectomy.

METHODS
=======

Study Populations
-----------------

Between January 1, 1999 and December 31, 2006, ambulatory subjects were included from 3 communes from Sicily, Italy, in the same Province of Enna: Troina, Cerami, and Gagliano Castelferrato (the Oasi cohort Sicily) (Supplemental Digital Content: Figure 1, <http://links.lww.com/MD/A287>).^[@R13]^ Institutional review board approval was obtained from the ethical committees of the Medical Centre of Troina (*Instituto di Ricovero e Cura a Carattere Scientifico*, the Oasi Maria Institute, Institute for Research on Mental Retardation, Troina, Sicily, Italy) and the University Hospital of Nancy (Vandoeuvre-lès-Nancy, France).^[@R13]^ All subjects involved in the study gave valid informed consent in accordance with approval from the relevant local ethical committees or institutional review boards and received a systematic clinical evaluation to rule out cancer, cardiovascular, renal, hepatic, or genetic disease, and any vitamin supplementation before inclusion in the study.

Blood Sampling
--------------

Phlebotomy was performed and fasting venous blood was collected in ethylenediaminetetraacetic acid-containing tubes. Samples were centrifuged immediately at 2500*g* for 15 minutes at room temperature. Aliquots were stored at −70°C until analysis.

Bilirubin Assays
----------------

All specimens were tested in a single laboratory (Vandoeuvre-lès-Nancy, France). Quantitative estimation of total and conjugated fraction of serum bilirubin was assayed by the Synermed (Rosny-Sous-Bois, France) reagent using the Olympus AU2700 chemistry-immuno analyzer (Rungis, France) according to manufacturer\'s protocol. Direct bilirubin reagent measures the absorbance of an oxidation product of bilirubin at 700 nm. The absorbance of the cholecyanin oxidation product is directly proportional to bilirubin concentration. Conjugated bilirubin is oxidized in the presence of the electron transfer agent, KHC3, to a blue-green infrared light absorbing chromophore. Unconjugated bilirubin is not significantly detected due to the fact that it is not soluble under the conditions of the direct reaction mixture and remains bound to albumin. Serum total bilirubin is reported in mg/dL. The intra-assay coefficients of variation for total bilirubin were 0.00% and 0.99% for level 1 and level 2, respectively. The corresponding figures for conjugated bilirubin were 0.00% and 1.07%, respectively.

Genotyping and Quality Controls
-------------------------------

Genomic DNA was extracted from the buffy coat and genotyping was performed at the National Institute of Health and Medical Research Unit-954 (Vandoeuvre-lès-Nancy, Lorraine, France) using the Illumina Infinium HumanExome BeadChip v1.2 (Illumina Inc, San Diego, CA) according to the manufacturer\'s recommendations. The Illumina HumanExome Beadchip includes 247,870 markers focused on protein-altering variants selected from \>12,000 exome and genome sequences representing multiple ethnicities and complex traits.^[@R14]^ Nonsynonymous genetic variants had to be observed ≥3 times in at least 2 studies, and splicing and stop-altering variants had to be observed ≥2 times in at least 2 studies.^[@R14]^ Additional, array content includes variants associated with complex traits in previous GWAS, HLA tags, ancestry-informative markers, markers for identity-by-descent estimation, and random synonymous SNPs.^[@R14]^ Data were analyzed using Illumina\'s GenomeStudio Genotyping Module.

Bead intensity data were processed and normalized for each sample in GenomeStudio. The Genome Reference Consortium GRCh37 (hg19; Feb. 2009) map was used (Illumina manifest file Immuno_BeadChip_11419691_B.bpm), and normalized probe intensities were extracted for all samples that passed standard laboratory quality-control thresholds. Genotype calling was performed with the GenomeStudio GenTrain 2.0 algorithm (Illumina\'s GenomeStudio data analysis software; Illumina Inc.).

Sample quality control measures included sample call rate, overall heterozygosity, relatedness testing. All patients included in this study also had a genotyping success rate \>90%. Duplicated samples and those showing cryptic relatedness were assessed by calculating identity by descent (IBD). We estimated IBD sharing and carried out PCA on an linkage disequilibrium (LD)-pruned set of 28,509 autosomal variants obtained by removing large-scale high-LD regions,^[@R15],[@R16]^ variants with a call rate \<0.99, variants with a minor allele frequency \<5% or variants with Hardy-Weinberg equilibrium (HWE) *P*-value \<10^−4^. By estimating the probability of sharing zero, 1, or 2 alleles for any 2 individuals as P(Z = 0), P(Z = 1), and P(Z = 2), respectively, a proportion of IBD was calculated using the following formula: PI-HAT = PI = P(Z = 1)/2 + P(Z = 2). Eleven pairs of individuals with a relatedness measure (PI-HAT) value \>0.99 were considered to be a duplicate sample. From these pairs, one of the individuals was randomly removed from the data. A PI-HAT threshold \>0.2 was used to suggest cryptic relatedness. Three samples had high positive inbreeding coefficients *f* suggesting truly inbred samples and were removed before analysis. Principal-component analysis (PCA) was performed on the samples merged with 7 HapMap phase 3 populations (CEU, TSI, YRI, MEX, JPT, CHD and CHB) as reference populations to identify ancestry outliers. Inspecting the first 10 principal components, we identified 3 population outliers that were excluded from subsequent analyses (Supplemental Digital Content: Figure 2, <http://links.lww.com/MD/A287>). All genotyping data sets underwent the same rigorous quality checks before cases and controls were compared. Genetic variants were removed from the primary analysis if they had a call rate \<95%, a significant differential missingness in cases and controls (*P* \< 0.05), a significant departure from Hardy--Weinberg equilibrium (exact HWE *P* \< 10^−4^), or a minor allele frequency (MAF) \<0.05%.

Replication Study
-----------------

The replication study was performed in an independent population of 227 Italian elderly subjects originating from Gagliano Castelferrato commune (Sicily, Enna). Genomic DNA was isolated from peripheral leukocytes using Qiagen kit (Qiagen-France, Courtaboeuf cedex, France). Determination of the *UGT1A6* gene polymorphisms rs6759892 (p.Ser7Ala), rs2070959 (p.Thr181Ala), and rs1105879 (p.Arg184Ser) were performed by LightCycler 480 instrument (Roche Molecular Biochemicals, Lyon, France) using commercial LightSNiP assays (SimpleProbe, TIB MolBiol, Berlin, Germany) according to the manufacturer\'s recommendations. Samples were set up in a final volume of 6 μL, containing 10 ng of DNA, 0.6 μL of LightCycler FastStart DNA MasterPLUS HybProbe Mix (Roche Molecular Biochemicals, Lyon, France), 3 mmol/L of MgCl~2~, and 0.3 μL of LightSNiP reagent mix (Tib-MolBiol, Berlin, Germany). The cycling conditions were as follows: initial denaturation at 95°C for 10 min, followed by 45 cycles of denaturation at 95°C for 10 s, annealing at 60°C for 10 s, and extension at 72°C for 15 s. Fluorescence was monitored at the end of each annealing phase at 60°C. After completion of the PCR, a melting curve of the amplification products was plotted by denaturation at 95°C for 30 s, 40°C for 2 min, and then slowly heating the sample to 75°C with a ramp rate of 2 °C/s and continuous fluorescence acquisition.

Statistical Analysis
--------------------

All quantitative variables are described as medians and percentiles (Interquartile range \[IQR\], 25--75th percentile). All proportions are expressed as percentages with 95% confidence intervals (95% CI). In the exome-wide stage, the association between Log transformed total, unconjugated, and conjugated bilirubin and genetic variants was determined by applying linear regression analysis and correlation/trend test using an allelic model with Bonferroni adjustment for multiple tests. The correlation/trend test tests the significance of any correlation between 2 numeric variables (or 2 variables that have been encoded as numeric variables). This test may also be thought of as any "trend," which either one of the numeric variables may have taken against the other one. The association between binary outcomes and genetic variants was determined by applying logistic regression analysis. We visually inspected the cluster plots for the most associated genetic variants to confirm genotyping quality. Linkage disequilibrium pairwise analysis was performed using a matrix output for both the expectation-maximization (EM) algorithm and composite haplotype method and both R2 and D' values. The comparison of bilirubin values across genotypes subgroups was performed using the Kruskal--Wallis test. All quality controls and genotypic analyses were performed using SNP & Variation Suite (SVS) 7.7.8 (Golden Helix Inc, Bozeman, MT).

Gene-Based Analysis
-------------------

We performed gene-based tests to further investigate the influence of rare and low-frequency variants on serum bilirubin level. Gene-based tests offer an alternative to single-variant tests, which are often underpowered to detect association with rare variants.^[@R14],[@R17]^ Accordingly, exome-wide data were analyzed using the combined multivariate and collapsing (CMC) method with regression analysis as described by Li and Leal.^[@R17]^ The CMC approach uses regions such as genes. The difference is that within each region, CMC first bins variants according to a criterion such as minor allele frequency, then collapses the variants within each bin, and finally performs multivariate testing on the counts across the various bins.^[@R17]^ We collapsed variants with MAF \<0.01 and variants with larger MAFs were investigated separately for each gene. In total, we tested 14,87 genes using the CMC method.

Model Building and Structural-Based Analysis
--------------------------------------------

Three-dimensional (3D) modeling of the human wild-type and variant proteins was performed using SWISS-MODEL, an automated homology modeling program.^[@R18],[@R19]^ The homology-modeling server SWISS-MODEL was accessed at <http://swissmodel.expasy.org/>.^[@R18],[@R19]^ The 3D models were visualized by the DeepView/Swiss-PdbViewer 4.1.0 and OpenAstexViewer 3.0 softwares.^[@R20],[@R21]^ For each 3D model of the native and variant proteins, narrow grooves, hydrophobic patches, and electrostatic potential on charged residues using the Coulomb method were studied using the DeepView/Swiss-PdbViewer 4.1.0 software.^[@R21]^

RESULTS
=======

In the initial study, we successfully genotyped 400 subjects originating from the Troina commune (Sicily, Enna) in the population-based OASI study for 242,901 variants on the Illumina HumanExome Beadchip and performed in silico replication study on 373 subjects originating from the Cerami commune (Sicily, Enna). Clinical characteristics of the 773 participants are summarized in Supplemental Digital Content: Table 1, <http://links.lww.com/MD/A287>. A total of 239,059 variants passed quality controls and were tested for association with unconjugated, conjugated and total serum bilirubin level, assuming additive allelic effects.

Genetic Variants Associated With Total Serum Bilirubin Level
------------------------------------------------------------

### **Exome-Wide Analysis in the Initial Study (n** **=** **400, Troina, Sicily, Enna)**

The Exome-wide association study established a significant association between 6 genetic variants on the chromosome 2 locus 2q37.1 and total serum bilirubin level. Among these 6 variants, 3 are located on the *UGT1A1* gene (UDP glucuronosyltransferase 1 family, polypeptide A1; Entrez ID: 54658) and 3 are located in the *UGT1A6* gene (UDP glucuronosyltransferase 1 family, polypeptide A6; Entrez ID: 54578) (Table [1](#T1){ref-type="table"} and Figures [1](#F1){ref-type="fig"} and [2](#F2){ref-type="fig"}). Functional annotation and quality controls of the *UGT1A1* and *UGT1A6* variants associated with serum total bilirubin level are reported in Supplemental Digital Content: Table S2, <http://links.lww.com/MD/A287>. The 3 variants in the *UGT1A1* gene (rs887829, rs4148325, rs6742078) are part of the same haplotype block (Supplemental Digital Content: Table 3, <http://links.lww.com/MD/A287> and Figure 3). The 3 variants on *UGT1A6* are non-synonymous (rs6759892, c.19T\>G, p.Ser7Ala, *P* = 8.34 × 10^−11^; rs2070959, c.541A\>G, p.Thr181Ala, *P* = 1.36 × 10^−11^; and rs1105879, c.552A \> C, p.Arg184Ser, *P* = 1.54 × 10^−11^) and are part of the same haplotype block on the first exon of *UGT1A6* (Supplemental Digital Content: Table and Figures 4 and 5, <http://links.lww.com/MD/A287>).

###### 

Genetic Variants Associated at the Exome-wide Level With Serum Total Bilirubin Level in the Troina Cohort Subjects

![](medi-94-e925-g001)

![(A) Manhattan plot for total serum bilirubin level. Association results of the single-variant analysis (−log10 P) are plotted against genomic position (NCBI build 37). The blue horizontal line indicates the significance threshold for genome-wide association.](medi-94-e925-g002){#F1}

![Bilirubin metabolism and *UGT1A1* and *UGT1A6* variants associated with unconjugated, conjugated, and total serum bilirubin level in physiological conditions. Coding variants are denoted in red font. *UGT1A1*: UDP-glucuronosyltransferase 1 family, polypeptide A1; *UGT1A6*: UDP-glucuronosyltransferase 1 family, polypeptide A6; OATP: organic anion transporter; MRP2: hepatic multidrug resistance protein.](medi-94-e925-g003){#F2}

### **In Silico Replication Study (n** **=** **373, Cerami, Sicily, Enna)**

In the in silico replication study, the 6variants that reached exome-wide significance for the association with serum total bilirubin level in the initial study were significantly associated with total serum bilirubin level (Table [1](#T1){ref-type="table"}). The pooled analysis of initial and in silico replication studies reported the same figures for the 6 variants on *UGT1A1* and *UGT1A6* genes with the following *P* values for the 3 non-synonymous variants on *UGT1A6*: (rs6759892, *P* = 1.98 × 10^−26^; rs2070959, *P* = 2.87 × 10^−27^; and rs1105879, *P* = 3.27 × 10^−29^) (Table [1](#T1){ref-type="table"}). Results were similar after adjusting regression analysis for age and sex (Supplemental Digital Content: Table 4, <http://links.lww.com/MD/A287>). Using the pooled data from both initial and in silico replication studies, for each of the 6 variants, the median total serum bilirubin values differed significantly between the 3 subgroups of genotypes (homozygous major, heterozygous, and homozygous minor) on the Kruskal--Wallis test after Bonferroni correction with higher total serum bilirubin level observed in subjects with homozygous minor genotype (Supplemental Digital Content: Table 5 and Figure 6, <http://links.lww.com/MD/A287>).

Genetic Variants Associated With Unconjugated Serum Bilirubin Level
-------------------------------------------------------------------

### **Exome-Wide Analysis in the Initial Study (n** **=** **400, Troina, Sicily, Enna)**

Consistently with the exome-wide analysis on total serum bilirubin level, the same top variants on *UGT1A1* and *UGT1A6* genes were significantly associated with serum unconjugated bilirubin level (Supplemental Digital Content: Table 6 and Figure 7, <http://links.lww.com/MD/A287>) with the following *P* values for the non-synonymous variants on *UGTA16*: (rs6759892, *P* = 4.75 × 10^−11^; rs2070959, *P* = 9.17 × 10^−12^; and rs1105879, *P* = 9.02 × 10^−12^).

### **In Silico Replication Study (n** **=** **373, Cerami, Sicily, Enna)**

In the in silico replication study, the same 6 variants that reached exome-wide significance for the association with serum unconjugated bilirubin level in the initial study were also significantly associated with unconjugated serum bilirubin level (Supplemental Digital Content: Table 6, <http://links.lww.com/MD/A287>). The pooled analysis of initial and in silico replication studies reported the same figures for the 6 variants with the following *P* values: (rs6759892, *P* = 4.19 × 10^−23^; rs2070959, *P* = 1.26 × 10^−24^; rs1105879, *P* = 1.41 × 10^−26^) (Supplemental Digital Content: Table 6, <http://links.lww.com/MD/A287>). Results were similar after adjusting regression analysis for age and sex (Supplemental Digital Content: Table 4, <http://links.lww.com/MD/A287>). As reported for total bilirubin level, using the pooled data from both initial and *1* replication studies, for each 1 of the 6 variants, the median unconjugated serum bilirubin values differed significantly between the 3 subgroups of genotypes (homozygous major, heterozygous, and homozygous minor) on the Kruskal--Wallis test after Bonferroni correction with higher unconjugated serum bilirubin level observed in subjects with homozygous minor genotype (Supplemental Digital Content: Table 5 and Figure 6, <http://links.lww.com/MD/A287>).

Genetic Variants Associated With Conjugated Serum Bilirubin Level
-----------------------------------------------------------------

### **Exome-Wide Analysis in the Initial Study (n** **=** **400, Troina, Sicily, Enna)**

Consistently with the exome-wide analysis on both total and unconjugated serum bilirubin level, the same top variants on *UGT1A1* and *UGT1A6* genes were significantly associated with serum conjugated bilirubin level (Supplemental Digital Content: Table 7 and Figure 8, <http://links.lww.com/MD/A287>) with the following *P* values for the non-synonymous *UGT1A6* variants: (rs2070959, *P* = 4.37 × 10^−8^, rs1105879, *P* = 4.49 × 10^−8^, rs6759892, *P* = 1.33 × 10^−7^).

### **In Silico Replication Study (n** **=** **373, Cerami, Sicily, Enna)**

In the in silico replication study, the 6 variants were significantly associated with conjugated serum bilirubin level (Supplemental Digital Content: Table 8, <http://links.lww.com/MD/A287>). The pooled analysis of initial and in silico replication studies reported the same figures for the 6 *UGT1A1* and *UGT1A6* variants with the following *P* values: (rs6759892, *P* = 6.51 × 10^−21^; rs2070959, *P* = 1.11 × 10^−20^; rs1105879, *P* = 1.66 × 10^−21^) (Supplemental Digital Content: Table 7, <http://links.lww.com/MD/A287>). Results were similar after adjusting regression analysis for age and sex (Supplemental Digital Content: Table 4, <http://links.lww.com/MD/A287>). As reported for both total and unconjugated serum bilirubin level, using the pooled data from both initial and in silico replication studies, for each of the 6 variants, the median conjugated serum bilirubin values differed significantly between the 3 subgroups of genotypes (homozygous major, heterozygous, and homozygous minor) on the Kruskal--Wallis test after Bonferroni correction with higher conjugated serum bilirubin level observed in subjects with homozygous minor genotype (Supplemental Digital Content: Table 5 and Figure 6, <http://links.lww.com/MD/A287>).

Combined Multivariate and Collapsing Method
-------------------------------------------

Exome-wide data analysis using the CMC analysis approach on the pooled data from the 773 subjects retrieved 6 loci in the 2q37.1 region (position from 2:234590584 to 2:234681945), namely: *UGT1A4* (*P* *=* 9.99 × 10^−19^), *UGT1A3* (*P* *=* 1.27 × 10^−18^), *UGT1A1* (*P* *=* 1.28 × 10^−18^), *UGT1A5* (*P* *=* 1.02 × 10^−13^), *UGT1A6* (*P* *=* 6.35 × 10^−10^), *UGT1A7* (*P* *=* 7.67 × 10^−10^). All the 6 loci *P* values were Bonferroni significant (Supplemental Digital Content: Table 8 and Figure 9, <http://links.lww.com/MD/A287>). From the Illumina Infinium HumanExome BeadChip, 54 variants corresponded to the 6 loci spanning from position 2:234590584 to 2:234681945. Among them, only the 6 above-mentioned *UGT1A1* and *UGT1A6* variants were significantly associated with serum total bilirubin level (rs887829, rs4148325, rs6742078 on the *UGT1A1* locus and rs1105879, rs2070959, and rs6759892 on the *UGT1A6* locus) and 37 variants were monomorphic.

The LD analysis revealed that the 3 top intronic SNPs rs887829, rs4148325, and rs6742078 on *UGT1A1* locus were in strong LD with 4 missense rare and low-frequency coding variants located in the first exon of the *UGT1A1* gene, namely: rs4148323 (MAF = 0.06%; p.Gly71Arg), rs144398951 (MAF = 0.06%; p.Ile215Val), rs35003977 (MAF = 0.78%; p.Val225Gly), and rs57307513 (MAF = 0.06%; p.Ser250Pro) (Supplemental Digital Content: Table 9 and Figure 10, <http://links.lww.com/MD/A287>). Due to the low MAF of these 4 *UGT1A1* exonic variants, only rs35003977 (p.Val225Gly) was studied in exploratory post-hoc analysis in the initial study and was associated with an elevated risk of having a serum bilirubin level beyond the upper normal limit (total bilirubin\>1.2 mg/dL \[odds ratio, OR 4.82; 95% CI, 1.29--18.11; *P* = 4.31 × 10^−2^\]; unconjugated bilirubin \>0.7 mg/dL \[OR = 4.52, 95% CI, 1.42--14.38; *P* = 1.78 × 10^−2^\]; and conjugated bilirubin \>0.3 mg/dL \[OR = 7.45, 95% CI, 1.58--35.11; *P* = 3.99 × 10^−2^\] Supplemental Digital Content: Table 10, <http://links.lww.com/MD/A287>).

Model Building and Structure-Based Analysis of UGT1A1 and UGT1A6 Proteins Variants
----------------------------------------------------------------------------------

The 3D models of the UGT1A1 proteins harboring the rare variants amino acid substitutions (p.Gly71Arg, p.Ile215Val, p.Val225Gly, and p.Ser250Pro) show alterations in the 3D protein structure, narrow grooves, hydrophobic patches, and electrostatic potentials (Figure [3](#F3){ref-type="fig"}). Similarly, the 3D models of the UGT1A6 proteins harboring the amino acid substitutions (p.Ser7Ala, p.Thr181Ala, and p.Arg184Ser) show alterations in the 3D protein structure, narrow grooves, hydrophobic patches, and electrostatic potentials (Figure [4](#F4){ref-type="fig"}).

![(A) Three-dimensional modeling of the human wild-type UGT1A1 protein and the low-frequency variants (p.Gly71Arg \[MAF = 0.06%\], p.Ile215Val \[MAF = 0.06%\], p.Val225Gly \[MAF = 0.78%\], and p.Ser250Pro \[MAF = 0.06%\]) performed using OpenAstexViewer 3.0. Secondary structure succession was colored using the 'rainbow' method. Elements of secondary structure are colored from the N-terminal to the C-terminal end in the order violet, blue, green, yellow, orange, and red. (B) Three-dimensional modeling of the human wild-type UGT1A1 protein and the low-frequency variants (p.Gly71Arg \[MAF = 0.06%\], p.Ile215Val \[MAF = 0.06%\], p.Val225Gly \[MAF = 0.78%\], and p.Ser250Pro \[MAF = 0.06%\]) performed using DeepView/Swiss-PdbViewer 4.1.0. Secondary structure succession was colored using the "rainbow" method. Narrow grooves in each model are depicted as a 3-dimensional network in a different color than the underlying secondary structure of the protein. (C) The negative electrostatic potential map is displayed as a red-colored 3-dimensional network surrounding the protein (D). The negative electrostatic potential map is displayed as a red-colored three-dimensional network surrounding the protein.](medi-94-e925-g004){#F3}

![(A) Three-dimensional modeling of the human wild-type UGT1A6 protein and the protein variants p.Ser7Ala (MAF = 36.3%), p.Thr181Ala (MAF = 32.7%), and p.Arg184Ser (MAF = 35.9%), and p.Ser250Pro (MAF = 0.06%) performed using OpenAstexViewer 3.0. Secondary structure succession was colored using the "rainbow" method. Elements of secondary structure are colored from the N-terminal to the C-terminal end in the order violet, blue, green, yellow, orange, and red. (B) Three-dimensional modeling of the human wild-type UGT1A6 protein and the protein variants p.Ser7Ala (MAF = 36.3%), p.Thr181Ala (MAF = 32.7%), and p.Arg184Ser (MAF = 35.9%) performed using DeepView/Swiss-PdbViewer 4.1.0. Secondary structure succession was colored using the "rainbow" method. Narrow grooves in each model are depicted as a 3-dimensional network in a different color than the underlying secondary structure of the protein. (C) The negative electrostatic potential map is displayed as a red-colored three-dimensional network surrounding the protein. (D) The negative electrostatic potential map is displayed as a red-colored three-dimensional network surrounding the protein.](medi-94-e925-g005){#F4}

**Replication Analysis on the 3 Top Genetic Variants on *UGT1A6* Gene (n** **=** **227, Gagliano Castelferrato, Sicily, Enna)**
-------------------------------------------------------------------------------------------------------------------------------

The independent replication study was performed on 227 Italian subjects from the Gagliano Castelferrato commune (Sicily, Enna). As the 3 intronic variants of *UGT1A1* have been previously reported and validated for their association with serum bilirubin level,^[@R22]--[@R29]^ we report here replication data on the 3 exonic missense variants on *UGT1A6* (rs6759892, p.Ser7Ala; rs2070959, p.Thr181Ala; and rs1105879, p.Arg184Ser). These variants were significantly associated with total serum bilirubin level, as well as unconjugated and conjugated serum bilirubin levels (Table [2](#T2){ref-type="table"}).

###### 

Association Between UGT1A6 Missense Variants and Serum Bilirubin (Total, Unconjugated, and Conjugated) in the Replication Study (Gagliano Castelferrato, Sicily, Enna)

![](medi-94-e925-g006)

Cholecystectomy Risk According to *UGT1A1* and *UGT1A6* Variants
----------------------------------------------------------------

Among the 773 OASI cohort participants from the Troina (Sicily, Enna) and Cerami (Sicily, Enna) communes, standardized collected data regarding gallstone-related cholecystectomy was available on 669 participants. Among them, 8 underwent cholecystectomy for cholelithiasis (1.20%; 95% CI, 0.37--2.02). All *UGT1A1* intronic variants (rs887829, rs6742078, and rs4148325) and *UGT1A6* exonic variants (rs6759892, rs2070959, and rs1105879) were significantly associated with gallstone-related cholecystectomy risk (Table [3](#T3){ref-type="table"}). The *UGT1A6* variant rs2070959 (p.Thr181Ala) was associated with the highest risk of gallstone--related cholecystectomy (OR, 4.58; 95% CI, 1.58--13.28; *P* = 3.21 × 10^−3^) (Figure [5](#F5){ref-type="fig"}). None of the 4 low-frequency and rare exonic variants of *UGT1A1* (rs4148323, p.Gly71Arg; rs144398951, p.Ile215Val; rs35003977, p.Val225Gly; and rs57307513, p.Ser250Pro) has been replicated due to their low MAF.

###### 

Effect-sizes of Gallstone-related Cholecystectomy Risk Among the OASI Cohort Participants According to UGT1A1 and UGT1A6 Genetic Variants

![](medi-94-e925-g007)

![Gallstone-related cholecystectomy risk, according to *UGT1A6* variants: (A) (rs6759892, p. Ser7Ala), (B) (rs2070959, p. Thr181Ala), and (C) (rs1105879, p. Arg184Ser); (D) gallstone-related cholecystectomy risk effect sizes according to *UGT1A1* and *UGT1A6* variants.](medi-94-e925-g008){#F5}

DISCUSSION
==========

The UGT1A enzymes family is derived from a single gene locus (*UGT1A*) spanning about 210 kb on chromosome 2 (2q37), which is composed of 17 exons including 13 unique alternate first exons followed by 4 common exons.^[@R30]^ Because of alternative splicing, only one of 13 different exon-1 sequences on the locus is associated with the 4 downstream exons, which are common to all *UGT1A* gene isoforms.^[@R30]^ Among the 13 exon-1 sequences, 9 code for functional proteins (UGT1A1, UGT1A3 to UGT1A10) and 4 correspond to pseudogenes (UGT1A2p, UGT1A11p, UGT1A12p, and UGT1A13p).^[@R30],[@R31]^ UGT1A proteins are composed of 265 (UGT1A6) to 534 amino-acid residues (UGT1A3, 4, and 5) with a molecular weight of 50 to 57 kDa.^[@R30]^ The exon-1 sequence of UGTs codes for the substrate-binding domain (N-terminal half of the protein), whereas the 4 common exons code for the cosubstrate-binding domain (C-terminal half of the protein).^[@R30]^

Several SNPs in the *UGT1A1* locus have been reported in previous studies to be significantly associated with serum bilirubin.^[@R22]--[@R29],[@R32]--[@R34]^ The reported *UGT1A1* SNPs are intronic (rs6742078, rs887829, rs4148324, and rs4148325) and are in complete linkage disequilibrium within the same haplotype block.^[@R24]^ To date, outside of published polymorphisms on *UGT1A1* gene isoform, no data have been reported regarding potential variants on other isoforms of *UGT1A* gene complex in physiological conditions.

We identified 4 low-frequency and rare missense variants located in the first exon substrate-binding domain coding region of the *UGT1A1* gene (rs4148323, p.Gly71Arg; rs144398951, p.Ile215Val; rs35003977, p.Val225Gly; and rs57307513, p.Ser250Pro). These coding variants are in strong linkage disequilibrium with 3 intronic *UGT1A1* SNPs, which are highly associated with serum bilirubin concentration. The *UGT1A1* intronic SNPs (rs6742078, rs887829, rs4148324) serve as proxy markers for the low-frequency and rare *UGT1A1*-coding variants, thereby providing mechanistic explanation to the relationship between *UGT1A1* intronic SNPs and the UGT1A1 enzyme activity.

Among the 4 low-frequency and rare variants that we found on the exon-1 of the *UGT1A1* gene, only 1 (rs4148323, p.Gly71Arg) has been reported to be associated with serum bilirubin level in Asian populations.^[@R22],[@R29],[@R33],[@R35],[@R36]^ Using the 1000 genomes database, the MAF of the rs4148323 (p.Gly71Arg) variant was 0.9% in Europeans versus 17.3% in East Asians. Regarding the 3 other missense variants that we found on *UGT1A1* first exon, the MAFs in Europeans for rs144398951 (p.Ile215Val), rs35003977 (p.Val225Gly), and rs57307513 (p.Ser250Pro) were 0% (Exome Sequencing Project, ESP6500), 0.1% (ESP6500), and 0% (1000 genomes), respectively. None of these variants has been reported to date as potentially associated with serum bilirubin level.

It has been demonstrated that 2 alternative orientations for bilirubin binding to UGT1A1, either toward the loop 5 and the membrane or toward the loop 3 and the endoplasmic reticulum lumen were possible.^[@R37]^ In the loop 5-binding configuration, bilirubin is in contact with residues Gln84, Asp87, Lys115, and Ile116, whereas in the loop 3-binding configuration, bilirubin is in contact with Lys115, Ile116, Ala174, Val193, and Asp224.^[@R37]^ It should be mentioned, however, that the substrates in the known glucuronosyltransferase structures are found closer to loop 5.^[@R37]^ It is noteworthy that 3 of the 4 variants on exon-1 of *UGT1A1* gene (rs144398951, p.Ile215Val; rs35003977, p.Val225Gly; and rs57307513, p.Ser250Pro) are close to the loop 5 on the corresponding UGT1A1 protein and therefore can potentially alter the binding site of the UGT1A1 protein to bilirubin.

In our study, we report for the first time a relationship between 3 nonsynonymous exonic SNPs of the *UGT1A6* gene (rs6759892, p.Ser7Ala; rs2070959, p.Thr181Ala; rs11058792, p.Arg184Ser) with unconjugated, conjugated, and total serum bilirubin level. To date, no study has reported association between *UGT1A6* gene variations and serum bilirubin in the general population. It has been demonstrated in whites that the 2 common variants, rs2070959 (p.Thr181Ala) and rs11058792 (p.Arg184Ser), were associated with 30% to 50% lower UGT1A6 enzyme activity.^[@R38]^ In a recent study on subjects with Gilbert syndrome---characterized by intermittent unconjugated hyperbilirubinemia without structural liver damage, affecting about 10% of the white population---UGT1A6 protein variants (p.Ser7Ala, p.Thr181Ala, and p.Arg184Ser) were observed in 80%, 79%, and 80% of subjects, respectively versus 9% in healthy blood donors.^[@R39]^ A humanized transgenic *UGT1A* mouse model containing a gene fragment from human chromosome 2 spanning *UGT1A1* up to a portion of *UGT1A10*, which encodes the most prevalent wild-type sequence (*htgUGT1A wild-type*) and a fragment of equal length containing *UGT1A6* SNPs (*htgUGT1A SNP*), showed lower *UGT1A6* messenger RNA (mRNA) expression and UGT1A6 protein synthesis in comparison with the corresponding wild-type mouse model.^[@R39]^

To date, 11 GWASs have reported the association with unconjugated, conjugated, and/or total bilirubin; however, no study reported significant association between *UGT1A6* locus and serum bilirubin the general population (Supplemental Digital Content: Table 11, <http://links.lww.com/MD/A287>).^[@R22]--[@R29],[@R32]--[@R34]^ In our study, using an exome-array approach, we demonstrated for the first time that the *UGT1A6* locus significantly influence total serum bilirubin as well as unconjugated and conjugated serum bilirubin in physiological conditions. Milton et al^[@R28]^ performed a GWAS on total bilirubin levels and cholelithiasis risk in patients with sickle cell anemia. In this study, total bilirubin level and cholelithiasis risk were associated with *UGT1A6* locus, confirming its potential direct implication in serum bilirubin determinism.^[@R28]^ However, results from sickle cell anemia patients' cohorts could not be inferred to general population and bilirubin metabolism in physiological conditions because the *UGT1A* alleles and haplotypes undergo protozoa-driven selective pressure.^[@R40]^

A major strength of the OASI cohort is that it is focused on the elderly and has assessed the phenotypic expression of genetic variants that require decades for their clinical translation. We reported, for the first time, an association between *UGT1A1* intronic variants and *UGT1A6* exonic variants (rs6759892, rs2070959, and rs1105879) and gallstone-related cholecystectomy risk. Indeed, using standardized collected data on a well-characterized cohort of ambulatory elderly subjects with a median age of 72 years, we showed for the first time high effect sizes for the significant association between *UGT1A1* and *UGT1A6* variants and the risk of gallstone-related cholecystectomy. In a meta-analysis on 2 GWASs, using hospital discharge diagnosis records Johnson et al^[@R24]^ reported no association between *UGT1A1* intronic SNP (rs6742078) and cholelithiasis-related admission, history of gallbladder disease, or gallbladder removal. Milton et al^[@R28]^ reported a small effect size significant association between *UGT1A1* SNPs and cholelithiasis risk in a young African American sickle cell anemia patients (median age around 33 years). In that study, cholelithiasis risk was defined by the presence of gallstones, non-functional gallbladder, or cholecystectomy and no data were reported separately on the risk of cholecystectomy as a hard criteria.^[@R28]^ Moreover, results on sickle cell anemia patients' cohorts could not be inferred to general population and bilirubin metabolism in physiological conditions.

Our study was limited in its ability to look at very rare variants because of the content of the exome array. Although sequencing will still be required to completely assess variants associated with serum bilirubin, this study provides proof of concept that exome array genotyping is a powerful approach to identify low-frequency functional variants.^[@R14]^

In conclusion, this study demonstrates for the first time that exome-array genotyping is a valuable approach to identify protein-coding variants that contribute to bilirubin metabolism and sheds new light on the genetic architecture of bilirubin metabolism, thus providing new evidence on the relationship between bilirubin metabolism and the risk of gallstone disease.

The authors thank all the data managers and hospitals for participating in the study and for their great efforts in entering data.

Abbreviations: 3D = three-dimensional, 95% CI = 95% confidence interval, EM = expectation-maximization, GWAS = genome-wide association study, HWE = Hardy--Weinberg equilibrium, IBD = identity by descent, LD = linkage disequilibrium, MAF = minor allele frequency, MRP2 = multidrug resistance-associated protein 2, OR = odds ratio, PCA = principal-component analysis, SNP = single nucleotide polymorphism, UGT1A1 = UDP-glucuronosyltransferase 1 family, polypeptide A1, UGT1A6 = UDP-glucuronosyltransferase 1 family, polypeptide A6.

The study was supported by grants from the Region of Sicily, the French National Institute for Health and Medical Research (INSERM) and the French National Agency for Research (ANR, Nutrivigène program)

The authors declare no conflicts of interest.

Supplemental digital content is available for this article. Direct URL citations appear in the printed text and are provided in the HTML and PDF versions of this article on the journal\'s Website ([www.md-journal.com](www.md-journal.com)).
