A Genome-Wide Association Meta-Analysis of Circulating Sex Hormone–Binding Globulin Reveals Multiple Loci Implicated in Sex Steroid Hormone Regulation

1Sex hormone-binding globulin (SHBG) is a protein secreted mainly by the liver that binds to the sex steroids, testosterone, dihydrotestosterone, and estradiol, transports them in the circulation, and influences their action in target tissues by regulating their bioavailability. SHBG thereby influences the expression of sex hormone sensitive phenotypes including sexual characteristics and reproductive function in men and women. In addition to regulating sex steroid hormone effects, SHBG may exert independent effects through its own receptor. Variation in SHBG concentration has also been associated with various chronic diseases including cancers, polycystic ovary syndrome (PCOS), and type 2 diabetes (T2D),. Although SHBG is estimated to have a heritable component (∼50%), little is known about the genetic regulation of SHBG. Polymorphisms at the SHBG gene locus have been associated with SHBG concentrations,, but much remains unknown about specific genetic variants that may determine circulating SHBG concentrations. Identifying genetic factors that influence SHBG may provide insights into the biology of sex steroid hormone regulation, metabolism and tissue effects that underlie their relationship with chronic diseases such as T2D as well as hormone-sensitive cancers such as breast and prostate cancer.

2We identified nine loci associated with SHBG concentrations at the genome-wide significance threshold of p = 5×10−8 (Table 1 and Figure 1) in a genome-wide association study (GWAS) meta-analysis of circulating SHBG concentrations in 21,791 men and women from 10 studies (Table S1). All nine lead SNPs at these loci had effects in the same direction (seven with p<0.05) in the validation dataset of 7,046 men and women from six additional studies (Table S2). The strongest association was within the SHBG locus (rs12150660, p = 2×10−106). Together, these nine lead SNPs explained 7.2% of the genetic variance (assuming 50% heritability) in SHBG concentrations.

3We next performed a series of additional analyses to explain more of the phenotypic variance (Figure 2). First, we hypothesized that genetic effects may be different in men and women, as SHBG concentrations are >50% higher in females than males, and may be differentially regulated between sexes. In a sex stratified analysis, three of the nine loci showed evidence of sex-differentiated effects at p<0.02 when we would not expect any signals to have reached this level of significance by chance. The associations at the 17p13.1-SHBG and Xq22.3 loci were stronger in males whereas the association at the 8q21.13 locus was stronger in females. To investigate the apparent differential sex effect for the X chromosome further we ran a recessive regression model for the X chromosome SNP rs1573036 in women in the Framingham Heart Study and found no association with SHBG suggesting the sex-differentiated effect is not the result of a recessive inheritance pattern. Sex stratified GWAS identified one novel signal in men, which showed no association in women (4q13.2: men p = 2.5×10−8, women p = 0.66, heterogeneity p = 0.003).

4A series of conditional analyses were performed to identify statistically independent signals. At the SHBG locus, three apparently independent additional signals separate from the main index SNP were observed, based on low (r2<0.05) pairwise correlations in HapMap (rs6258 p = 2.7×10−46, rs1625895 p = 1.2×10−14 and rs3853894 p = 2.5×10−11). A series of iterative conditional analyses (Table 2) involving SNPs at the SHBG locus generated a final regression model including six statistically independent SHBG SNPs. Four of these SNPs (#1–4 Table 2) retained GWS when conditioned against the other five, and two were nominally associated (SNP#5 p = 0.0001, SNP#6 p = 0.01). Re-running the GWAS meta-analysis adjusting for these six SNPs revealed evidence for three additional statistically independent (pairwise HapMap r2<0.01) signals at the SHBG locus (SNP#7 p = 1.5×10−7, SNP#8 p = 4.6×10−5, SNP#9 p = 9.9×10−6) (Figure 3). There were also two additional trans signals located at 2p16.3 and 7q21.3 (Table 1). Although the 2p16.3 signal dropped below GWS when combined with follow-up samples (p = 1×10−7), the index SNP at 2p16.3 is ∼300 kb away from a strong candidate gene, the luteinizing hormone receptor gene (LHCGR).

5The majority of pair-wise correlations for the nine SHBG locus SNPs highlighted by our conditional analyses showed very low HapMap r2 values. However, the pairwise D′ values are often high (Table S3) indicating that no or few recombination events have occurred between some SNPs, and that combinations of SNPs may be tagging un-typed variants on a common haplotype. To investigate this possibility, we performed more extensive analyses in a single study (NFBC1966, n = 4467). We used a denser set of SNPs imputed from the June 2011 version of the 1000 Genomes data and performed model selection analyses. Model selection identifies a set of SNPs that best explain phenotypic variation, while simultaneously penalizing each SNP included in this set, and therefore correlated SNPs tend to be excluded from the final model. These analyses consistently included at least seven SNPs in the model, although it is hard to estimate the false-negative rate of using the reduced sample size. While we are underpowered to accurately pinpoint the exact number of independent signals, these analyses support the results of the conditional analysis and suggest that multiple variants at the SHBG locus are independently associated with SHBG concentrations.

6Data from an independent study, the InCHIANTI study, was used to calculate the proportion of genetic variance in SHBG concentrations explained when accounting for sex specific effects, the multiple signals of association at the SHBG locus, and the additional trans signals identified post conditional analysis. In men and women we explained ∼15.6% and ∼8.4% of the heritable component respectively. The SHBG locus accounted for ∼10% and ∼6.6% of the genetic variation in men and women respectively with the lead SNP in isolation accounting for ∼7.8% and ∼3.3% of the variation in men and women, respectively.

7We identified genes near the associated SNPs and explored their biologic relevance to SHBG. The genes associated with identified SNPs included the SHBG locus (rs12150660, 17p13.1, p = 1.8×10−106), PRMT6 (rs17496332, 1p13.3, p = 1.4×10−11), GCKR (rs780093, 2p23.3, p = 2.2×10−16), ZBTB10 (rs440837, 8q21.13, p = 3.4×10−09), JMJD1C (rs7910927, 10q21.3, p = 6.1×10−35), SLCO1B1 (rs4149056, 12p12.1, p = 1.9×10−08), NR2F2 (rs8023580, 15q26.2, p = 8.3×10−12), ZNF652 (rs2411984, 17q21.32, p = 3.5×10−14), TDGF3 (rs1573036, Xq22.3, p = 4.1×10−14), LHCGR (rs10454142, 2p16.3, p = 1.3×10−07), BAIAP2L1 (rs3779195, 7q21.3, p = 2.7×10−08), and UGT2B15 (rs293428, 4q13.2, p = 5.5×10−06) (Figure 1).

8We used the online tool STRING (www.string-db.org) to perform pathway analyses to explore possible interactions between the SHBG gene and the proteins encoded by the 11 most plausible genes nearest the 11 SNPs listed above. There was an interaction noted between GCKR and JMJD1C which were associated with the lipoprotein fractions VLDL and HDL, respectively. In an expanded analysis, we assessed protein interactions among SHBG and 67 genes within 500 kb of our 11 identified SNPs and uncovered additional protein interaction pathways. An interaction between two proteins encoded by GTF2A1L and STON1 was found; these proteins are co-expressed in testicular germ cells in the mouse. An interaction between LHCGR and BRI3 encoded proteins that are associated with the G-protein coupled receptor complex in the human luteinizing hormone receptor was also identified. Finally, an interaction between LHCGR and IAPP (amylin) proteins which are components of a ligand/G-protein receptor/G-protein alpha subunit complex was found (database: www.reactome.com).

9Targeted analysis of two strong candidate genes, hepatocyte nuclear factor-4α (HNF4α) and peroxisome-proliferating receptor γ (PPARγ) did not identify any SNPs at HNF4α but did identify one SNP, rs2920502, at PPARγ that reached statistical significance (p = 9.9×10−5) and a second SNP at PPARγ, rs13081389, that reached nominal significance (p = 0.01).

10In total, we identified 12 genomic regions associated with circulating SHBG concentrations, including extensive allelic heterogeneity at the SHBG locus itself. Conditional meta-analyses carried out at the SHBG locus, identified nine genome-wide significant SNPs with low correlation (r2<0.01) between them. Two of these signals (rs6258 and rs6259) are missense variants and two are low frequency variants (MAF ∼2%). Furthermore, rs12150660 is highly correlated (r2>0.95) with a pentanucleotide repeat, which affects SHBG expression in-vitro. To our knowledge, the magnitude of secondary signals observed at this locus are the largest seen for any complex trait.

11The proportion of genetic variance in SHBG serum concentrations explained when accounting for sex specific effects, the multiple signals of association at the SHBG locus, and the additional trans signals identified post conditional analysis was ∼15.6% in men and ∼8.4% in women. The SHBG locus accounted for ∼10% and ∼6.6% of the genetic variance in men and women, respectively, with the lead SNP explaining most of the genetic variation at ∼7.8% for men and ∼3.3% for women. Thus additional signals at the SHBG locus identified through conditional analyses approximately doubled the variance of the trait explained. While we provide evidence for multiple variants associated with SHBG concentrations, further studies are needed to pinpoint the causal loci and functional variants. For the 11 regions outside the SHBG locus, most have biologically plausible related genes within 300 kb.

12Several genes near the identified SNPs regulate sex steroid production and function. The NR2F2 locus (15q26.2) encodes a nuclear receptor important in testicular Leydig cell function, the primary source of gonadal testosterone production, and has been linked to male infertility. NR2F2 has also been associated with estrogen receptor alpha (ERα) signaling and may influence hormone responsivity in breast cancer. PRMT6 (1p13.3) also encodes a nuclear receptor regulatory protein that mediates estrogen signaling as a co-activator of the estrogen receptor. LHCGR (2p16.3) encodes the luteinizing hormone receptor which was associated with polycystic ovary syndrome (PCOS) in a recent GWAS,. PCOS is both a reproductive and metabolic disorder characterized by higher testosterone serum concentrations as well as an increased prevalence of obesity, insulin resistance, and T2D in women. Inappropriate secretion of luteinizing hormone leads to increased ovarian production of testosterone. Coincident lower SHBG concentrations contribute to increased bioavailable testosterone concentrations and the expression of both reproductive and metabolic phenotypes in PCOS,,.

13The SLCO1B1 locus encodes a liver-specific transporter of thyroid hormone as well as estrogens which impact liver production of SHBG. JMJD1C (10q21.3), also known as TRIP 8 (thyroid hormone receptor interactor protein 8 ), may impact SHBG concentrations via thyroid hormone effects on liver protein production. Thyroid hormone may alter SHBG production through effects on HNF4α which is known to regulate SHBG transcription,.

14Many of the genes identified are involved in carbohydrate and lipid metabolism and liver function. The GCKR locus (2p23.3) encodes a protein that regulates glucokinase activity and has been associated with T2D in several ethnic populations,,,. GCKR has been associated with metabolic and inflammatory traits including triglyceride concentrations and other lipid fractions,, fasting plasma glucose,, insulin concentrations, uric acid, c-reactive protein (CRP), and non-alcoholic fatty liver disease which are all characteristic of the metabolic syndrome and T2D,,,,,,,,. The SLCO1B1 locus (12p12.1) codes for a protein, hepatocyte protein anion-transporting polypeptide 1B1, involved in liver metabolism of both endogenous and exogenous compounds. Consistent with SLCO1B1's role in liver metabolism, the same SNP (rs4149056) has been associated with circulating bilirubin concentrations in previous GWAS. BAIAP2L1 (7q21.3) encodes a protein important in cytoskeleton organization that has been associated with the inflammatory marker CRP in patients with arthritis. BAIAP2L1 is also known as IRTKS (insulin receptor tyrosine kinase substrate) which is involved in insulin receptor signaling and may relate to insulin resistant states including obesity and T2D,,,,,,. We conducted a targeted analysis of PPARγ, a gene that influences SHBG gene expression in the liver, and is associated with T2D,. Our analysis identified one significant SNP (rs2920502, p = 9.9×10−5) and a second nominally significant SNP (rs13081389, p = 0.01) at PPARγ. Some of the identified genes involved in hepatic metabolism of lipids and carbohydrates may be affect SHBG concentrations indirectly through effects on the SHBG transcription regulator HNF4α although HNF4α itself was not identified in this meta-analyses,,,.

15The UGT2B15 locus (4q13.2) was significantly associated with SHBG concentrations in men but not women in this meta-analysis. UGT2B15 belongs to a family of genes (the UGT2B gene family) that code for enzymes involved in the metabolism of sex hormones through glucuronidation which allows for excretion of sex steroids through the kidney and the gut via bile excretion,, primary clearance mechanisms for sex steroids. UGT2B15 is involved in the conjugation and inactivation of testosterone. An association between rs293428 in the UGT2B15 locus and circulating SHBG concentrations in men is supported by a previous study demonstrating that a non-synonymous SNP in UGT2B15 (rs1902023; D85Y) is associated with serum SHBG concentrations in younger adult men. UGT2B15 is thought to play a significant role in local tissue inactivation of androgens in androgen dependent prostate cancer,. The mechanism behind the influence of genetic variants in UGT2B15 on SHBG concentrations is unknown, but one may speculate that UGT2B15 affects the local androgenic environment in selected tissues, which in turn results in regulation of SHBG concentrations.

16In addition to UGT2B15, three other genes near the identified SNPs are associated with carcinogenesis, particularly in the prostate and breast. ZBTB10 (8q21.13), has been linked to breast cancer. In breast cancer cell lines ZBTB10 is suppressed by ROS-microRNA27a thereby enhancing ERα alpha expression and mediating estrogen effects. The ZNF652 (17q21.32) locus codes for a DNA binding protein thought to act as a tumor suppressor gene in breast cancer,, that is also co-expressed with the androgen receptor in prostate cancer. TDGF3, teratocarcinoma derived growth factor 3, is the only significant region identified on the X chromosome ((Xq22.3). TDGF3 is a pseudogene of TDGF1 located on chromosome 3p23-p21 that has been associated with testicular germ cell tumors.

17This GWAS meta-analysis incorporated data from approximately 22,000 men and women from 16 epidemiologic cohorts. The overall size of the study yields power but the meta-analysis of data from different epidemiologic studies requires the inclusion of different laboratory methods. The different studies used a variety of assay methodologies to measure serum SHBG concentrations although the vast majority were immunoassays (Tables S1 and S2, Text S1) with similar methodologies. Variation introduced by the use of different SHBG assays would result in loss of statistical power and likely bias toward the null. Additionally, the majority of women were post-menopausal as ascertained by self-report in all studies (Table S1). SHBG concentrations, like testosterone, decline only slightly across the menopause so adjustment for menopause status is not necessary. SHBG may also increase with ovulation and be slightly higher in the luteal versus the follicular phase of the menstrual cycle in premenopausal women, but most studies did not collect data on menstrual phase at the time of SHBG measurement so adjustment for menstrual phase was not possible. Finally, individuals were not excluded based on health status, therefore some individuals with chronic conditions that may affect hepatic production of or clearance of proteins including SHBG such as liver disease, renal disease, or severe malnutrition, may have been included in this analysis.

18SHBG synthesis in the liver is known to be affected directly or indirectly by estrogens, androgens and thyroid hormones and has been observed to be inversely associated with the higher insulin concentrations characteristic of insulin resistant states such as T2D,. In summary, the results of this GWAS reflect these influences. Three regions map to proteins related to hepatic function (12p12.1-SLCO1B1, 2p23.3-GCKR and 10q21.3-JMJD1C ). In addition, 2p23.3-GCKR and 7q21.3-BAIAP2L1 [alias insulin receptor tyrosine kinase substrate (IRTKS)] are involved in susceptibility to T2D and insulin signaling, respectively. Two signals also mapped to loci involved in thyroid hormone regulation (10q21.3-JMJD1C and 12p12.1-SLCO1B1). One signal mapped to the receptor for luteinizing hormone 2p16.3-LHCGR, the hormone that stimulates testosterone production. Five regions mapped to genes previously implicated in androgen and estrogen signaling (1p13.3-PRMT6, 8q21.13-ZBTB10, 12p12.1-SLCO1B1, 15q26.2-NR2F2, 4q13.2-UGT2B15 ).

19We have combined a conventional GWAS approach with detailed additional analyses, including sex stratification, conditional analysis and imputation from 1000 Genomes. Our results demonstrate that these approaches can lead to an appreciable gain in heritable variance explained. It does however highlight the complexity of elucidating individual variant causality through statistical approaches. In addition to the extensive allelic heterogeneity at the SHBG locus, our data identify loci with a role in sex steroid hormone metabolism, which may help elucidate the role of sex steroid hormones in disease, particularly T2D and hormone-sensitive cancers.

20We performed a genome wide association study (GWAS) meta-analysis of 21,791 individuals (Table S1: 9,390 women, 12,401 men) from ten observational studies. Data from an additional six studies totaling 7,046 individuals (Table S2: 4,509 women; 2,537 men) were used for validation. The proportion of variance explained was estimated in an independent study (InCHIANTI, n = 1,129). The individual study protocols were approved by their respective institution's ethics committee/institutional review board and all participants provided informed consent prior to participation. Individuals known to be taking hormonal contraceptives or hormone replacement therapy at time of SHBG measurement were excluded from analysis. Age, sex and body mass index (BMI) were included as covariates. After applying standard quality control measures, imputed genotypes were available for approximately 2.5 M SNPs. See Figure 2 for an overview of the analytic plan and the Text S1 for further information for individual studies included in this meta-analysis.

21The initial starting point for the conditional analysis was the four SHBG locus SNPs that all showed low Hapmap LD (r2<0.05) with each other: rs12150660 (lead SNP Table 1), rs6258 p = 2.7×10−46, rs1625895 p = 1.2×10−14 and rs3853894 p = 2.5×10−11. Each cohort fitted a single regression model, fitting SHBG concentrations against these four genome-wide significant SHBG locus SNPs (rs12150660, rs6258, rs1625895 and rs3853894), in addition to age, sex and BMI. After meta-analyzing the results from all cohorts, three of the SNPs retained genome wide significance when regressed against each other, with the fourth SNP narrowly missing that threshold (rs3853894, p = 4.1×10−6).

22We next performed a conditional GWAS meta-analysis, where each study included, as additional covariates to the original analysis plan, the ten genome-wide significant autosomal SNPs (the eight ‘Main’ signals from Table 1 and the two unique SHBG locus signals described above in addition to the lead SNP rs12150660: rs6258 and rs1625895). Three additional signals (independence based on HapMap r2<0.05) at the SHBG locus reached genome-wide significance (rs1641537 p = 7.8×10−32, rs6259 p = 1.5×10−12 and rs10432029 p = 3×10−8), giving a total of six independent signals in this gene region. In addition, two novel signals reached genome-wide significance in the conditional analysis, at 7q21.3 (rs3779195 p = 1×10−8) and 2p16.3 (rs10454142 p = 3×10−8). After replication, only rs3779195 at the BAIAP2L1 locus retained genome-wide significance.

23Given the six signals observed at the SHBG locus (three through conditional analysis #1 rs12150660, rs6258, rs1625895, three through LD estimates from conditional GWAS #1: rs1641537, rs6259, rs10432029), we sought to confirm which of these six were truly independent by a second round of conditional analysis. All discovery and replication cohorts fitted a single regression model of the six SNPs (SNPs # 1–6, Table 2) against SHBG concentrations, using the same parameters and covariates as conditional analysis #1. Four of the six SNPs (#1–4: rs12150660, rs6258, rs1641537, and rs1625895) retained genome-wide significance when conditioned against each other, with two showing nominal evidence of association (SNP #5 rs6259, p = 0.0001; SNP #6 rs10432029, p = 0.01).

24Finally, we performed a second conditional GWAS analysis, adjusting for the six SHBG locus SNPs which had evidence of association from conditional analysis #2. All the discovery cohorts were used in this analysis, in addition to three replication cohorts (total sample size 24,354). This analysis revealed evidence for a further three independent signals at the SHBG locus (based on HapMap r2<0.01), SNP #7 rs9901675 p = 1.5×10−7, SNP #8 rs8077824 p = 4.6×10−5, and SNP #9 rs9393218 p = 9.9×10−6.

25We performed a sensitivity analysis using samples from the 1966 Northern Finland Birth Cohort (NFBC1966) study to further investigate allelic heterogeneity at the SHBG locus (Text S1). The conditional meta-analysis showed evidence for up to nine signals at the SHBG locus, but it is possible that these signals could be explaining a much smaller number of causal variants in the region. Since 1000 Genomes imputation allows us to assess the genetic variation associated with a phenotype across a much denser set of markers, it increases our power to detect allelic heterogeneity within a region. Therefore, 1000 Genomes imputation was carried out on all the samples in the NFBC1966 study and forward selection was used to identify the set of SNPs that best explain the variation in the SHBG phenotype. 1000 Genomes imputation was carried out using IMPUTE2. The mean genotype probabilities for each SNP were calculated and used in the model selection step. Only SNPs 250 kb upstream and 250 kb downstream from the SHBG locus (7283453–7786700 bp) were used in the analysis. All SNPs with MAF <0.1% or an imputation quality score less than 0.4 were excluded from the analysis. In total, 1978 SHBG region SNPs measured or imputed in 4467 samples from the NFBC1966 study were used in the sensitivity analysis. Forward selection was implemented in R (version 2.13.0) using the stepAIC package to estimate the Akaikie Information Criterion (AIC), an inclusion parameter. Given the high degree of correlation between the SNPs in this region, we increased the penalty (k) on the number of terms included in the model to 12 (where it is usually two), to minimize possible over fitting. The final model included seven SNPs, adjusted for sex and BMI.

26We examined potential interactions among the proteins encoded by the SHBG locus and the proteins encoded by the 11 genes (ZBT10, TDGF1, ZNF652, PRMT6, JMJD1C, GCKR, BAIAP2L1, LHCGR, SLCO1B1, UGT2B15, NR2F2) closest to the 11 identified SNPs using pathway analysis with Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) Pathways Analysis (www.string-db.org). The interactions explored by STRING include direct (physical) and indirect (functional) associations. We then expanded the analysis to examine protein interactions among the SHBG gene and the proteins encoded by 67 genes within 500 kb of the 11 identified SNPs.

27We conducted targeted analysis of two strong candidate genes, hepatocyte nuclear factor-4α (HNF4α) and peroxisome-proliferating receptor γ (PPARγ). Statistical significance thresholds were set correcting for the number of SNPs tested in each gene region (±100 kb).