Multi-Analytic Approach Elucidates Significant Role of Hormonal and Hepatocanalicular Transporter Genetic Variants in Gallstone Disease in North Indian Population

Objective Cholesterol gallstone disease (CGD) is a multifactorial and multistep disease. Apart from female gender and increasing age being the documented non-modifiable risk factor for gallstones the pathobiological mechanisms underlying the phenotypic expression of CGD appear to be rather complex, and one or more variations in genes could play critical roles in the diverse pathways further progressing to cholesterol crystal formation. In the present study we performed genotyping score, Multifactor dimensionality reduction (MDR) and Classification and Regression Tree analysis (CART) to identify combinations of alleles among the hormonal, hepatocanalicular transporter and adipogenesis differentiation pathway genes in modifying the risk for CGD. Design The present case-control study recruited total of 450 subjects, including 230 CGD patients and 220 controls. We analyzed common ESR1, ESR2, PGR, ADRB3, ADRA2A, ABCG8, SLCO1B1, PPARγ2, and SREBP2 gene polymorphisms to find out combinations of genetic variants contributing to CGD risk, using multi-analytical approaches (G-score, MDR, and CART). Results Single locus analysis by logistic regression showed association of ESR1 IVS1-397C>T (rs2234693), IVS1-351A>G (rs9340799) PGR ins/del (rs1042838) ADRB3-190 T>C (rs4994) ABCG8 D19H (rs11887534), SLCO1B1 Exon4 C>A (rs11045819) and SREBP2 1784G>C (rs2228314) with CGD risk. However, the MDR and CART analysis revealed ESR1 IVS1-397C>T (rs2234693) ADRB3-190 T>C (rs4994) and ABCG8 D19H (rs11887534) polymorphisms as the best polymorphic signature for discriminating between cases and controls. The overall odds ratio for the applied multi-analytical approaches ranged from 4.33 to 10.05 showing an incremental risk for cholesterol crystal formation. In conclusion, our muti-analytical approach suggests that, ESR1, ADRB3, in addition to ABCG8 genetic variants confer significant risk for cholesterol gallstone disease.


Introduction
Cholesterol Gallstone disease (CGD) corresponds to one of the most recurrent and costly gastroenterological disorder. It is worldwide health problem representing 10% to 15% of the adult population in industrialised countries [1,2] whereas a prevalence of 6% have been reported from North India [3]. The female gender and increasing age are the documented non-modifiable risk factors for gallstones [4], the pathobiological mechanisms underlying the phenotypic expression of CGD appear to be rather complex, and one or more defects could occur in genes that play critical roles in the diverse pathways leading to cholesterol gallstone formation. The genetic determinants of gallstone formation have only recently been dissected in humans [5] Compelling evidence for familial clustering and an increased concordance of the trait in monozygotic twins as compared to dizygotic twins [6] further confirms the heritability of gallstones. Thus 'Gallstone genes' are continuously being corroborated in genome-wide association studies (GWAS), in case-control cohorts, and in family studies [7,8].
In human physiology, the gender disparity commences with puberty and continues through the childbearing years [9,10,11,12] which suggests that female sex hormones, could be an important risk factor for the formation of cholesterol gallstones [13]. The actions of these hormones such as estrogen, progesterone and catecholamines are executed through one or more of their respective receptors as estrogen receptors (ESRs), progesterone receptor (PGR) and adrenergic receptor (ADR) Allelic variants of ESR, PGR and ADR genes have been shown to be associated with susceptibility or progression with various disorders such as myocardial infarction [14,15], cholesterol gallstones and biliary tract diseases [16].
Another area of interest is hepatocanalicular transporters namely ATP-binding cassette transporters (ABC transporters) and organic anion transporters both encoded by ABC and SLCO1B1 genes respectively. Mutations in genes encoding these transporters have been implicated in cholesterol gallstones formation owing to their ability to influence bile composition and causing retention of substances normally secreted in bile.
Peroxisome proliferator-activated receptor c 2 (PPARc2) orchestrate the adipocyte differentiation process whereas sterol regulatory element binding protein 2 (SREBP-2) is involved in adipocyte differentiation followed by cholesterol homeostasis. Series of previous observations have suggested that regulatory interactions between the SREBPs and PPARc2 can coordinate cholesterol and fatty acid metabolism. Therefore, sequence variation in these genes may further disrupt the cholesterol homeostasis which in turn may nurture the development of CGD.
Previously, we have studied the role of some individual genetic variants with CGD susceptibility in a North Indian population [17,18,19]. Individual SNPs have little predictive value because of their modest effect on risk, but combinations of gene variants may improve the predictive ability and could be used to model susceptibility to CGD. Therefore, the current study aimed to search for gene-gene interactions in the selected pathways (hormonal, hepatocanalicular and adipogenesis differentiation) as a key contributory factor in the disease outcome.

Population Characteristics
The demographic profile of gallstone patients with respect to their age and gender matched controls are presented in Table 1.

Overall Frequency Distribution of Selected Hormonal, Hepatocanalicular Transporter and Adipogenesis Differentiation Gene Polymorphisms in GSD Patients and Healthy Subjects
Association of hormonal pathway gene polymorphisms with gallstone patients. Table 2 shows the risk of gallstones in relation to each of the SNPs of ESR1, ESR2, PGR, and ADR in hormonal pathway. On comparing the genotype frequency distribution of our study groups i.e. gallstone patients with that of healthy subjects (HS), the homozygous variant genotypes of  Table 2) when compared with homozygous wild-type DD genotype. Furthermore, on subdividing the study groups on the basis of gender we observed that ESR1 IVS1-397C.T and ADRB3 -190 T.C conferred increased risk for gallstones in female gender (Table S5).
Association of hepatocanalicular transporter pathway gene polymorphisms with gallstone patients. Table 3 shows the risk of gallstones in relation to each of the SNPs of ABCG8 and SLCO1B1 in hepatocanalicular transporter pathway. We found that in single locus analysis, the variant genotypes (GC+CC) of ABCG8 145 G.C and (CA+AA) of SLCO1B1 463 C.A were significantly associated and conferred increased risk of gallstone disease (p = ,0.019; [OR], 2.4: p = 0.007; [OR], 2.6). On the contrary, no significant difference were observed in the distribution of SLCO1B1 521 T.C (rs4149056) (ptrend = 0.416; MCS = 0.298) polymorphism, both at genotypic and allelic levels and therefore conferred no risk for developing gallstones.
Association of adipogenesis differentiation pathway gene polymorphisms with gallstone patients. Table 4 shows the genotype and allele frequency distribution of sequence variants in SREBP2 1784 G.C and PPAR c2 C.G. A borderline statistical significance was observed when the homozygous variant genotypes of SREBP2 1784 G.C (rs2228314) was compared i.e gallstone patients with that of healthy subjects (HS) (p = 0.045; [OR], 4.8). Furthermore, no significant difference was observed in the distribution of PPAR c2 C.G (rs1801282) (ptrend = 0.256; MCS = 0.218).

Haplotype Analysis
Linkage disequilibrium and haplotypes analysis of ESR1 and ESR2 in case and control groups. On LD analysis, ESR1 rs2234693 and rs9340799 were found to be in strong linkage disequilibrium (D' = 0.575). Haplotypes were constructed for the three polymorphisms in ESR1 gene including IVS1-397C.T, IVS1-351A.G and Ex4-122C.G. The haplotypes comprising   the homozygous wild alleles were taken as reference and the difference in the frequencies of haplotypes between patients and controls were tested using chi-square test.
The results of the studied three polymorphisms of ESR1 revealed that distribution of T,G,C haplotypes was significantly higher in gallstone patients (25.1% v/s 13.7) in comparison to controls and was conferring high risk for gallstone disease (p = 0.0012; [OR], 2.2). Global haplotypes analysis indicated a statistically significant difference between cases and controls based on the distribution pattern of the ESR1 haplotypes (p = ,0.001). Furthermore, none of the ESR2 haplotypes conferred risk for gallstones presenting the global haplotypes association p-value = 0.65 (Table S2; S3).
Linkage disequilibrium and haplotypes analysis of SLCO1B1 in case and control groups. SLCO1B1 Exon4 C.A and Ex6+40T.C were found to be in strong linkage disequilibrium (D' = 0.8916). Haplotypes analysis of these two polymorphisms gave rise to three haplotypes, of which C, T was the most common haplotypes in control population. On comparing the haplotypes frequencies in controls and gallstone cases, A, T haplotypes was more commonly distributed in gallstone patients and was imposing risk for the disease (p = 0.017; OR = 2.21) (Table S4).
G-score. For each individual, we counted the number of riskincreasing alleles. The number of risk alleles ranged from 1 to 11 in overall 450 subjects (Figure 1). The mean (6SD) G-score was 5.4361.96 in gallstones subjects and 4.6361.95 in controls (pvalue = ,0.001) ( Table 5). At the more extreme ends of the risk distribution, CIs around risk estimates became very wide because of small numbers. The number of risk alleles ranged from 1 to 11 with a median of 4 among control subjects and 6 among cases ( Figure 1). The risk for gallstone disease was estimated for each number of risk alleles, relative to the median number of risk alleles of 4, and ranged from an OR of 2.27 (95% confidence interval [CI], 1.0-4.6) for 5 risk alleles to an OR of 8.27 (95% CI, 0.90-75.2) for 11 risk alleles. The average relative risk increase per risk allele, when treated as an ordinal variable, however, could be estimated with a high level of precision, and was 2.7 (95% CI, 1.12-1.16). This corresponded to several fold difference in risk between the lowest and the highest number of risk alleles in our population.
Association of High-Order Interactions with GSD Risk by MDR Analysis Table 6 shows the best interaction model by MDR analysis. The one-factor model for predicting GS risk was ABCG8 145 G.C SNP (testing accuracy = 0.515, CVC = 7/10, permutation p = 0.027). The two-factor model of ESR1 IVS1 351A.G and ADRB3 -190 T.C had an improved testing accuracy of 0.578 (permutation p = ,0.001) however, the CVC was increased (10/ 10). The three factor model was the three-factor model including  Table 7 shows the CART, which included all investigated genetic variants of the selected pathways. The final tree structure contained seven terminal nodes as defined by single-nucleotide polymorphisms of the hormonal, hepatocanalicular transporter and adipogenesis differentiation pathway genes. Consistent with the MDR best one-factor model, the initial split of the root node on the decision tree was ESR1 IVS1-397C.T, suggesting that this SNP is the strongest risk factor for GSD among the polymorphisms examined. Individuals carrying ESR1 IVS1 -397CC, ADRB3-190 TT, ABCG8 145 GG and ADRA2A GG genotypes had the lowest case rate of 17.2%, considered as reference. Further inspection of the tree structure revealed distinct interaction patterns between individuals carrying the ESR1 IVS1-397 variant and those with the ADRB3 variant and SLCO1B1 463 C.A wild genotypes. Using the terminal node with lowest case rate as reference, individuals carrying the combination of ESR1 IVS1-397TT, SLCO1B1 Exon4CC, ESR2 1082GG, ESR1 IVS1-351AA and ESR1 Ex4-122GG exhibited a significantly higher risk for

Discussion
In order to achieve a more comprehensive evaluation of CGD risk, present analysis was performed in order to identify high and low intrinsic risk sets of sequence variants. Of the included 13 polymorphisms, some of them were found to be significantly associated with CGD risk in our previous studies [17,18,19] while others showed little or no influence on the risk for CGD development. Moreover, accumulating evidence supports the importance of adipogenesis differentiation and adrenergic receptor pathways in cholesterol associated diseases [22,23,24]. Therefore, we further extended our work by incorporating these two pathways.
Based on the candidate SNPs in genes involved in the gallstone pathway. We created a consolidated Genotype Score (G-score) from the number of risk alleles as previously reported for risk assessment of cardiovascular events and diabetes [25,26]. Our assumption was that individuals with a high G-score might have a higher probability of gallstone development as compared to those with the low G-score. The overall G-scores for the three selected pathways obtained were highly significant and conferred increased risk for gallstone development. Further calculating the G-score individually in respective pathways we found both hormonal and hepatocanalicular transporter pathway conferred increase risk. These results suggest significant role of hormonal receptor and hepatocanalicular transporters in gallstone disease.
For the higher order gene-gene interaction analysis, we employed statistical approaches namely MDR and CART analysis to find out the particular combinations of genetic variants contributing to CGD risk. In MDR analysis, we observed the best four-factor interaction model consisting of ESR1 IVS1-397C.T, ESR1 IVS1 351A.G, ADRB3 -190T.C and ABCG8 145 G.C with highest testing accuracy compared with the onefactor model.
In CART analysis, which is a non-parametric statistical approach for conducting regression and classification analyses by recursive partitioning. [17], study subjects were grouped according to different risk levels on the basis of the different gene polymorphisms. From this analysis, we found that development  of CGD involves complex genetic interactions among the hormonal and hepatocanalicular transporter genetic variants. As our results from CART analyses consistently suggested that ESR1IVS1-397TT, ABCG8GC+CC, ESR1IVS1-351GG and ADRB3 TC+CC polymorphisms are the most important single susceptibility factor for CGD development.
The association between hormonal receptor gene polymorphisms and risk of gallstones are biologically convincing. It has been assumed that the gallbladder is a female sex hormone responsive organ, and these hormones might be involved in the pathogenesis of gallbladder diseases. Elaborating on estrogen receptor the animal studies have shown that ESRs are present in the hepato-pancreatic-biliary tree [27,28,29] including bile duct epithelial cells and gallbladder. In addition, immunohistochemical and quantitative RT PCR studies have also revealed that the expression level of ESR1 gene is approximately 50 fold higher compared to ESR2. In animal models, 17beta estradiol promoted gallstone formation which further involves the upregulation of hepatic expression of ERalpha but not ERbeta. These studies show that ESR-1 is key player and findings may offer a new approach to treat gallstones by inhibiting hepatic ER activity with a liver-specific, ERalpha-selective antagonists.
The literature regarding the ADRB3 confirms that it is localized in the smooth muscles of the vasculature and the muscularis propria of the gallbladder [30] where it is thought to mediate relaxation and increase mucosal blood flow. The T.C polymorphism results in lowered responsiveness to potent agonists including endogenous catecholamines. [31] The mutated receptor had less ability to stimulate adenylyl cyclase and therefore less accumulation of cAMP. [31] Activation of ADRB3 also results in smooth muscle relaxation in the guinea-pig common bile duct, [32] and since the ductal smooth muscle appear to be more sensitive to activation of the ß3-adrenoceptor, there is the possibility that these receptors may be involved in the regulation of tone in the ductal smooth muscle and hence the outflow of bile. Thus the inhibiting variant C in ADRB3 might result in gallstone formation by impairing the relaxation of the gallbladder and probably the biliary tree too, setting the stage for crystal formation.
In the selected hepatocanalicular transporters ABCG8 145G.C conferred increased risk both individually and in combination to hormonal receptors. A genome wide scan carried out by Buch et al., [33] identified a variant D19H in the hepatic cholesterol transporter (ABCG8) as major susceptibility factor for human gallstone disease. Subsequently, this association has been replicated in various populations [9,19,34,35].
The phenomenon that a combination of polymorphisms within genes of unrelated pathways may elevate the risk for CGD could be explained by two hypotheses. One possibility is that some connection between these genes or proteins exists but still remains to be discovered. Another hypothesis, more credible in our opinion, is that the genes influencing risk for CGD may as well comprise a set of alterations located within genes not related to each other.
Our multi-analytic approach revealed that the combination of genotypes of respective polymorphisms as ESR1 IVS1-397 variant, ABCG8 145 variant, ESR1 IVS1-351 variant and ADRB3 190 variant pose a significant risk for developing gallstone. Comparing to the results of single locus analysis the role of SLCO1B1 Exon4 C.A SREBP2 1784 G.C and PGR ins/del was diminished when the overall analysis of 13 selected polymorphisms was performed. It also suggests that ESR1, ADRB3 and ABCG8 have significant incremental risk factors for gallstone disease. Thus, the application of these multi analytical approaches allowed creating a decision that has more sensitivity or specificity and was more accurate with reasonable power as compared to single strategy employed in calculating risk allele for disease prediction.
Also a prominent significant role of hormonal pathway was elucidated when the means of genotyping scores of selected pathways was calculated separately or all together. Therefore, exhaustive analysis of multi-analytic approaches as MDR, CART and G-scores are well recognized methods in understanding complex traits, such as disease susceptibility and also the etiology of complex diseases.
In summary, this is the first comprehensive study to use a multigenic analysis for cholesterol gallstone disease, and the data suggest that individuals with a higher number of genetic variations in hormonal and hepatocanalicular transporter pathway genes are at an increased risk for cholesterol gallstone disease, confirming the importance of taking a multigenic pathway based approach to risk assessment. The finding also indicates that the development of gallstone involves complex genetic interactions and follows different pathways depending on the specific genetic background of the subjects. The present study provided evidence supporting the cholesterol supersaturation contribution of hormonal, hepatocanalicular transporter and adipogenesis differentiation pathway genes, of which interaction between ESR1, ADRB3 and ABCG8 genes were the most important.
Thus, our results support the concept that genetic polymorphisms can be used as cholesterol gallstone risk predictors and multiple polymorphisms allow more precise delineation of risk groups and suggest the future direction of association studies. However, the present study included only North Indian individuals, therefore the results need to be replicated in other ethnic groups.

Ethics Statement
The institutional ethical committee of Sanjay Gandhi Post Graduate Institute of Medical Sciences (SGPGIMS) approved the present study protocol and the authors followed the norms of World's Association Declaration of Helsinki. All the participants provided written informed consent.

Study Population
The case control study recruited a total of 450 subjects, including 230 cholesterol gallstone patients (GS) and 220 healthy subjects. From the year June 2006 to September 2011 symptomatic cholesterol GS patients attending the Department of Gastrosurgery, Sanjay Gandhi Post Graduate Institute of Medical Sciences and Department of Surgical Oncology Lucknow India, were approached for participation in the present study. All subjects were unrelated and confirmed to North Indian ethnicity.
Phenotype data. For each individual, ultrasound examinations were conducted at the Department of Radio-diagnosis and Imaging SGPGIMS, Lucknow. Participants were considered as having gallstones when one of the subsequent diagnostic criteria was satisfied: (1) Gallbladder lumen with mobile nodular or dependent layering echoes that exhibited posterior acoustic shadowing, or (2) Gallbladder with hyperechoic shadowing material filling the gallbladder lumen with an appearance of the WES triad (i.e., the gallbladder wall, the echo of the stone, and the acoustic shadow-a specific ultrasonographic sign of gallstones used to make a reliable diagnosis of cholelithiasis [36]. The healthy controls were randomly selected from a pool of healthy volunteers that visited the general health check-up center at SGPGIMS Lucknow, during the same period. In addition to a self-reported gallstone history, transabdominal ultrasound was performed to validate gallstone status and to identify silent gallstones. Inclusion criteria for controls also included absence of asthma, coronary artery disease, diabetes mellitus determined through maternal and paternal family history. At recruitment, informed consent was obtained from each subject and the information on demographic characteristics, such as sex and age was collected by questionnaire. Both patients and controls had similar ethnicity. The blood sample and the clinical details were collected from each participant at recruitment.

DNA Samples and Genotyping
Genomic DNA was isolated from peripheral blood leukocytes using salting out method [37]. The polymorphisms were genotyped using the PCR or PCR restriction fragment length polymorphism method. The details of genotyping for studied polymorphisms are shown in Table S1. Ten percent of masked, random sample of cases and controls were tested twice by different laboratory personnel and the reproducibility was 100%.
Genotype Score Calculation (G-score) A Genotype score (G-score) was defined as the cumulative number that counts the total number of risk-increasing alleles in individuals. Genotyping of 13 selected SNPs in candidate genes involved in hormonal, hepatocanalicular and adipogenesis differentiation pathways was performed and G-score was computed from the number of variant alleles. A value of 2, 1 and 0 was allotted to homozygous variant, heterozygous and homozygous wild type genotypes respectively. Variant genotype was considered as risk conferring. Using these 13 SNPs a Genotype Score (Gscore) was constructed ranging from 0 to 26 on the basis of the number of risk alleles. For each sample a consolidated G-score was calculated by adding the values from all 13 SNPs together.

Statistical Analysis
Descriptive statistics were presented as mean and standard deviation [SD] for continuous measures while absolute value and percentages were used for categorical measures. Differences in genotype and allele frequencies between study groups were estimated by chi-square test. Unconditional logistic regression was used to estimate odds ratios [ORs] and their 95% confidence intervals [CIs] adjusting for age and sex. A two-tailed p-value of less than 0.05 was considered a statistical significant result. All statistical analyses were performed using SPSS software version 16.0 (SPSS, Chicago, IL, USA). Ptrend and Monte Carlo Simulation (MCS) were calculated through Cochrane Armitage trend using XLstats whereas haplotype analysis was performed using SNPstats (http://bioinfo.iconcologia.net/SNPstats).
Furthermore, higher-order gene-gene interactions associated with CGD risk were determined through multifactor dimensionality reduction (MDR) using software version 2.0 beta8 and classification regression tree analysis (CRT) using SPSS software version 16.0.