Deciphering genetic causes for sex differences in human health through drug metabolism and transporter genes

Sex differences in genetic architecture and genetic effects of DMET genes

Sex differences in human phenotypes can be driven by sex-differentiated genetic effects^{4, 18}. However, such effects can be masked in GWASs because individuals of both sexes have typically been analyzed together⁷. In this study, we evaluated the sex differences in the genetic architecture of human complex traits in each sex separately using sex-stratified GWAS summary statistics from the UKBB (http://www.nealelab.is/uk-biobank, Fig. 1a). We analyzed 564 traits (421 binary/categorical traits, 143 continuous traits, Supplementary Data 2), for which at least one DMET region SNP shows significant trait association based on the GWAS catalog. DMET genes, which encode 222 metabolism enzymes and 150 transporters, were retrieved from a previous publication (Supplementary Data 1)¹⁰.We first estimated the narrow-sense heritability of each trait separately in each sex to quantify the proportion of phenotypic variance explained by the common genetic variation¹⁹. Sex differences in trait heritability suggest a different molecular mechanism and/or a different degree of environmental effect between sexes on the trait. Of 564 traits, 83 traits (14.72%) showed significant sex differences in heritability (FDR < 0.05, Fig. 1b and Supplementary Data 3). Among them, 56 traits (67.4%) had higher heritability in males, including neuropsychiatric traits, such as ever attempted suicide and anxiety, and other diseases, such as gout, cardiovascular complications (coronary atherosclerosis) and diabetes. In comparison, traits such as hypothyroidism, osteoporosis and gallbladder disorders have higher heritability in females. These findings are robust across different methods employed to estimate heritability (Supplementary Fig. 1). When expanding the genome-wide heritability analysis to additional 1222 traits that have sufficient samples for both sexes in UKBB and are not known to related to DMET genetic regions, we found 13.7% (167/1222) of them showing sex differences in global heritability (Supplementary Fig. 2, Supplementary Data 4). Interestingly, similar 13.40% (71/530) traits showing significant differences in their heritability between two sexes have been reported by an independent study⁷.

We also estimated genetic correlation between sexes for each trait to evaluate potential sex differences in genetic architecture. Genetic correlation estimates the proportion of variance that two groups share due to genetic causes²⁰. Low genetic correlation between males and females in the same trait suggests sex differences in the genetic architecture. We focused on only those heritable traits (defined by greater than median estimated heritability across all traits through sex-combined GWAS analysis) to generate robust estimates of the male-female genetic correlation. We identified 253 traits for which male-female genetic correlation differed significantly from 1 (Supplementary Fig. 3, Supplementary Data 3, FDR < 0.05), suggesting a global difference in genetic architecture between the sexes. These traits with sex-differentiated architecture include gout (r_g = 0.459; FDR = 1.26 × 10⁻⁹), heart attack (r_g = 0.573; FDR = 6.39 × 10⁻¹⁰), and cholelithiasis (r_g = 0.634; FDR = 1.18 × 10⁻⁵).

For variants located in DMET gene regions, to quantify sex differences in SNP-trait associations, we calculated the z-score and P values for each SNP that mapped into the cis-DMET regions to each trait. We defined SNPs with sex-differentiated effects (SDEs) as those trait-associated SNPs with a significant sex difference in genetic effects (based on z-score FDR < 0.05). In total, we identified 25 traits harboring at least one SDEs (Fig. 1c, Table 1 and Supplementary Data 5) and mapped 954 SDEs onto 109 tagging loci (r² < 0.2) and 125 genes. The two disease traits with the largest numbers of SDEs are self-reported gout and hypothyroidism/myxedema (Table 1). We indeed found SDEs mapped in genes previously implicated in disease, such as ABCG2 in gout¹⁵.

15 traits were found to have sex differences both in heritability (global) and SDEs in DMET gene regions (Fig. 1d), affirming the presence of different genetic effects between males and females. However, whether these differential genetic effects are functionally related to human health is unknown. To test this, we employed colocalization analysis, which assess the probability of two GWAS traits sharing a common causal variant in a genomic region²¹. In another word, if for the same trait, the sex-stratified GWAS results were not colocalized, the probability is high that males and females have a different causal mechanism for that trait. Among the aforementioned 15 traits, we excluded three of them (Treatment code: levothyroxine sodium, treatment code: thyroxine product, major coronary heart disease event excludes revascularization) as they represented information redundant to other traits in the list. All 12 traits exhibit genetic correlation coefficients differing from 1 (Supplementary Fig. 4). We detected 22 male-specific causal loci in traits such as gout and major coronary heart disease events, and 13 female-specific causal loci in traits such as hypothyroidism (Fig. 1e, Supplementary Data 6) after colocalization analysis. All support the presence of sex-differentiated genetic effects of traits in cis-DMET gene regions.

Gout and hypothyroidism were found to have higher number of sex-differentiated causal SNPs in the DMET gene regions. Previous findings indicate that accumulation of crystal form of urate at the joint causes gout²². We identified male-specific SNP-trait associations of gout in the SLC22A12 region (Fig. 1f), specifically an upstream variant rs2360872. SLC22A12, also known as urate transporter 1, is involved in regulating urate levels in the blood and its variants are associated with serum uric acid level and gout development²³. Our results suggest that rs2360872 affects gout uniquely in males. Additionally, we detected female-specific SNP-trait associations of hypothyroidism at the 5’UTR of SLC66A1 (Fig. 1f). SLC66A1 is a lysosomal amino acid transporter that mediates the export of cationic amino acids from lysosomes²⁴. The relationship between SLC66A1 and thyroid function remains unknown.

Sex differences in the genetic regulation of gene expression may lead to observed sex differences in human complex traits

Sex-differentiated genetic regulation of gene expression can lead to different manifestations in downstream biological pathways. Therefore, we characterized sex differences of genetic regulation on gene expression in human liver samples using sex-stratified cis-eQTLs produced by the GTEx Consortium²⁵.

Sex-differentiated eQTLs were defined as expression associated variants who have significantly different effect size between males and females (quantified through a z-score testing). We identified 31 sex-differentiated eQTLs (FDR < 0.1) in and near 2 DMET genes (Supplementary Fig. 5, Supplementary Data 7). Among them, rs34109652, is associated with UGT2B17 gene expression only in males (β_male = 0.84; P_male = 5.44 × 10⁻⁷; β_female = 0.27; P_female = 0.078; Z-score = −2.95; Fig. 2a). This gene encodes an enzyme that transforms steroid hormones such as testosterone²⁶. A previous study reported higher UGT2B17 enzyme activity in males²⁷. This male-specific genetic regulation could in part explain the observation of higher UGT2B17 activity in male, though functional work is required for validation.

In addition to sex-differentiated eQTLs, we also identified 40 potential sex-specific eQTLs in the DMET gene regions (FDR < 0.1 in one sex and FDR > 0.1 in the other sex of eQTL, Fig. 2b, Supplementary Data 7). For example, variant rs854572 on chromosome 7, is associated with PON1 gene expression in males only (P_male = 4.13 × 10⁻⁸; P_female = 0.01; Z-score = 0.33 Fig. 2a). This gene encodes serum paraoxonase and arylesterase 1 enzyme, an enzyme exerting protective effects in Major Adverse Cardiovascular Events (MACE = death, myocardial infarction, stroke). In a sex-combined GWAS, this variant is associated with the activity of paraoxonase and arylesterase 1²⁸. Moreover, higher arylesterase activity has been reported in females²⁹. Although sex-specific association of this variant with serum enzyme activity has not been tested, this sex-specific eQTL provides a putative explanation for the observed higher PON1 serum activity in females. Together, our results support the presence of sex-differentiated genetic regulation of DMET gene expression.

To examine the potential biological consequences of these sex-differentiated/specific eQTLs, we performed colocalization²¹ between them and the 564 complex traits described above. We identified 11 traits with a sex-differentiated/specific eQTL colocalization (PPH4 > 0.5 in only a single sex, Fig. 2c, Supplementary Data 8 Binary Traits, Supplementary Data 9 continuous traits). Among them, multiple traits related to alcohol intake are colocalized with a male-specific eQTL for ADH1C (PPH4_male = 0.73, PPH4_female = 0.073, Fig. 2d). This gene encodes enzymes that metabolize a wide variety of substrates, including alcohol³⁰. Alcohol metabolism determines blood alcohol level over time, and the extent of organ exposure to alcohol. Therefore, ADH1C has been implicated in alcohol dependency. Sex differences in alcohol consumption have been widely reported, where males are more likely to drink alcohol and to consume more than females³¹. Genetic variants in ADH1C are associated with alcohol metabolism capacity³². In males, genetic variants of ADH1C showed significant association and a larger effective size with heavy/excessive alcohol drinking habits than females³³. Our results provided a plausible explanation for this phenomenon where sex-differentiated genetic regulation of ADH1C may be the culprit. Importantly, all traits that showed significant sex-dependent colocalization with ADH1C eQTLs were related to alcohol drinking frequency, but not alcoholism. This is concordant with previous reports that ADH1C genetic variants were not associated with alcoholism³³. Of note, the observed sex dimorphism in alcohol consumption can be affected by factors, such as body size, sociocultural behavior. The identification of genetic contributor to this trait should not be interpreted without these other factors. Our results provided a plausible explanation for this phenomenon where sex-differentiated genetic regulation of ADH1C may play a part in this complex issue. Overall, our colocalization results linked the genetic basis of a small number of complex traits to sex-dependent genetic effects on expression levels of specific DMET genes.

Sex differences in the DMET region genetic regulation of human serum biomarkers and their impact on human health outcomes

Clinical laboratory tests are frequently used to diagnose diseases and monitor human health status. Many of these serum biomarkers are transformed in the liver by DMET enzymes¹⁰. Importantly, sex differences have been known for the level of many serum biomarkers³⁴. While the genetic basis of serum biomarkers has been studied in large cohorts in a sex-combined fashion³⁵, sex differences in genetic basis of those biomarkers and putative causal relationships with human diseases have not been extensively studied. Given the relevance of serum biomarkers in human health and the potential masking effect introduced by the sex-combined model, we hypothesized that sex-differentiated causal relationships exist, and can be revealed by sex-aware characterization of the genetic regulation of serum biomarkers and human disease risks. To test this hypothesis, we performed Mendelian Randomization (MR) in both sex-combined and sex-stratified models. We pre-selected traits that harbor at least 1 significant SNP (P < 5 × 10⁻⁸) that are mapped to a DMET gene region. In total, we have 29 serum biomarker traits as exposures and 186 outcomes selected from the 564 human complex traits described above. A sex-specific likelihood of causal relationship is defined if it is significant in only a single sex but not in the sex-combined population (P < 0.05/(186*29) threshold, Fig. 3a).

We first examined the proportion of significant serum biomarker level-associated SNPs (P < 5 × 10⁻⁸) that are located in the DMET regions relative to the whole genome (Fig. 3b, Supplementary Data 10); and found the DMET region variants were associated with levels of 29 serum biomarkers. The DMET region serum biomarker associations account for 20% (in the case of c-reactive protein) to 100% (in case of oestradiol) of overall genomic association with these serum biomarkers. The proportion of serum biomarker-associated variants that mapped to cis-DMET gene regions is greater than random selection of variants from genome sequence of the same length (Supplementary Fig. 6), which suggests the important impact of DMET genes on serum biomarker levels.

Using the same analytical pipeline described earlier on human complex traits (not including any of the serum biomarker traits), we characterized sex differences in the genetic basis of serum biomarker traits. As expected, we observed that hormone-related traits, such as testosterone, harbor the largest sex differences in estimated heritability (Supplementary Data 11), genetic correlations (Supplementary Fig. 8, Supplementary Data 11), and SDEs (Supplementary Fig. 9, Supplementary Data 12)³⁴.

By conducting MR analysis using sex-combined and sex-stratified models, we identified a total of 141 female-specific and 167 male-specific putative causal relationships (Fig. 3c, Supplementary Data 13), which were not seen in the sex-combined model. Serum biomarkers, such as bilirubin and aspartate aminotransferase, have a large number of female-significant causal relationships, whereas testosterone, cholesterol and glucose have a large number of male significant causal relationships (Supplementary Fig. 10). We found sex-specific causal relationships between blood cell counts (e.g., monocyte count, platelet count) and several anthropometric traits (e.g., BMI, waist-to-hip ratio, hip circumference). While sex differences have been known for either the blood cell counts or the anthropometric traits^{36, 37}, these putative sex-specific causal relationships have not been reported. Among diseases, we identified 24 traits that exhibit at least 1 sex-specific causal relationship (Fig. 3d). We observed that testosterone increases the likelihood of high blood pressure only in females; and apolipoprotein B increases the likelihood of coronary atherosclerosis and major coronary heart disease events only in males. Given that DMET genes are largely involved in the transformation of these endogenous compounds, we hypothesized that genetic variations in DMET gene regions affect levels of these serum biomarkers and result in sex-specific human health outcomes. We further conducted MR analysis only considering SNPs in DMET gene regions in traits with sex-specific causal effect; and identified 61/141 (43.2%) female-specific and 66/167 (39.5%) male-specific putative causal relationships that remained valid (Bonferroni adjusted P < 0.05, Supplementary Fig. 11, Supplementary Data 14). Overall, our results highlight the importance of examination of causal relationships in each sex separately, by which sex-specific causal effect can be identified.

To understand the molecular basis of sex differences in this causal relationship, we examined shared loci that associated with testosterone and high blood pressure in each sex. In females, six loci (P < 5 × 10⁻⁸) were shared between testosterone and high blood pressure, and 3 (CYP11B1, SLC16A1, SLC22A7) of them were in DMET gene regions (Fig. 3e). In males, only 4 loci were shared between testosterone and high blood pressure, none of these loci were located in the DMET gene regions (Supplementary Fig. 12). Variant rs7003319 in CYP11B1 3′UTR is associated with both testosterone level and high blood pressure in females, but not males. Colocalization analysis of these two traits revealed a shared genetic basis only in females (PPH4_females = 0.73, PPH4_males = 0.00034; Fig. 3f). CYP11B1-encoded enzymes catalyze glucocorticoid production and lead to the production of 11-ketotestosterone (11-KT)³⁸. A previous study reported that 11β-hydroxylase deficiency (11βOHD), a rare autosomal recessive disorder caused by mutations in the CYP11B1 gene, is associated with hypertension³⁹. Further investigation of the interplay between CYP11B1, testosterone, and high blood pressure in the context of biological sex is warranted.

Sex differences in DMET gene expression are linked to sex differences in drug response

Elimination of drugs and exogenous toxins is the major function of DMET gene products. Such processes mainly take place in the liver. Differences in the abundance and activity of DMET have been shown to play a major role in drug responses and side effects¹⁰. Although sex differences in the gene expression of DMET genes have been studied in different contexts⁴⁰, the impact of these differentially expressed genes on responses to relevant medications is sparsely reported. Here, we hypothesized that sex differences in DMET gene expression would provide insights into observed sex differences in drug responses and side effects.

Using RNA-seq data from GTEx liver tissue, we identified 20 DMET genes with significant sex differences in expression; 10 expressed higher in males and 10 others expressed higher in females (FDR < 0.05, Fig. 4a, Supplementary Data 15). Among them, several well-characterized pharmacogenes, such as CYP1A2, CYP3A4 and CYP2C19, were found to exhibit sex-differential expression. When assessing the reproducibility of our discovery in an independent dataset, we recapitulated the differential expression for 14 of 19 genes (expression of CYP1A2 was not quantified in the validation dataset) (Supplementary Fig. 13)⁴¹. Two additional smaller datasets^{42, 43} were also evaluated and our top differential expressed genes, such as UGT2B17, UGT2A3, CYP3A4, SLC3A1, SLC16A14 are concordant with previous finding (Supplementary Data 16). Protein abundance of two pharmacogenes, CYP1A2 and CYP3A4, was quantified in pooled human liver microsomes (HLMs) that were collected separately from male and female donors. We confirmed higher protein abundance of CYP1A2 in the male-pooled HLMs and the opposite trend for CYP3A4 (Fig. 4b), which is correlated with a previous report¹¹.

To evaluate the potential clinical consequences of these sex differentially expressed DMET genes with respect to medication use, we annotated these genes with FDA-approved drugs. Specifically, using data from PharmGKB⁴⁴ and DrugBanks⁴⁵, both of which have summarized gene and drug relationships from published literature, we found 1,166 drugs that have been linked to at least one of the sex differentially expressed genes identified in this study (Supplementary Fig. 14, Supplementary Data 17). Not surprisingly, CYP3A4, a cytochrome P450 family member which is responsible for metabolizing approximately 50% of drugs on the market, has been linked with the greatest number of drugs.

To further investigate whether these annotated drugs have been reported to have sex differences in response, we performed web scraping to comprehensively search literature in PubMed for evidence supporting sex differences in response to these annotated drugs. Using “sex difference” and drug names as our key search terms, we identified 2,340 journal articles that contained information for 306 drugs. Given the high intensity of inspecting all articles, we focused on those drugs that were annotated with CYP1A2, as a proof-of-concept. This resulted in 519 journal articles which were then manually inspected to determine whether sex differences in drug responses were reported. Among them, we identified 75 articles that reported sex differences in either efficacy, pharmacokinetics, or toxicity of 49 drugs that are metabolized by CYP1A2 (Fig. 4c, Supplementary Data 18). Fluoxetine, for example, was reported to have greater weight reduction in females than in males at the same dosage⁴⁶. Furthermore, antipsychotic medications, such as clozapine and olanzapine, have been reported to have a slower elimination in females than males⁴⁷, which is in agreement with our findings that males have a higher expression of CYP1A2 (Fig. 4d).

Our analysis identified a number of sex-differential drug responses that were supported by existing literature. For example, flunarizine, another CYP1A2 substrate, which is used in treating epilepsy⁴⁸, was reported to not affect catalepsy in male mice, but attenuated catalepsy in females at the same doses⁴⁹. Another study found that male rats formed two oxidative metabolites of flunarizine at higher rate than female rats⁵⁰. These are in agreement with our findings that higher CYP1A2 expression in male can lead to faster metabolism/breakdown of this drug and therefore less response. Similarly, female mice have been reported to be less susceptible to acetaminophen overdose induced hepatotoxicity than male mice⁵¹. Acetaminophen is also metabolized by CYP1A2 to N-acetyl-p-benzoquinone imine (NAPQI), known to induced hepatotoxicity. In this case, the lower expression of CYP1A2 in female liver would lead to lower production of the toxic metabolite and therefore partially explain the lower toxic response in females.

Given the wide usage of clozapine and the higher rate of clozapine-induced toxicity observed in males⁵², we hypothesized that the higher abundance of CYP1A2 contributes to higher formation of an active metabolite, and subsequently leads to higher rate of clozapine-induced side effects in males (Fig. 4e). We quantified the sex differences in clozapine metabolism using separate pools of HLMs from men and women. Because CYP3A4 also contributes to the metabolism of clozapine, in our HLM experiments we employed ketoconazole, a CYP3A4 inhibitor, to specifically evaluate CYP1A2-mediated metabolite formation. In the presence of ketoconazole, the formation of N-demesthylclozapine, the active metabolite, was significantly higher in the male HLMs (Fig. 4f, P = 0.0309). We observed no sex differences in the formation of an inactive metabolite, clozapine N-oxide (Fig. 4f) and parent drug clozapine (Supplementary Fig. 15). Higher plasma concentration of N-desmethylclozapine has been associated with clozapine side effects, such as agranulocytosis⁵³. Our observation supports the notion that the observed higher frequency of clozapine-induced toxicity in males can be due to elevated formation of N-desmethylclozapine resulting from higher CYP1A2 levels in males (Fig. 4e).