physiological differences between the sexes have their roots in sex differences in gene expression (37, 58). Although the sexes possess identical autosomes and differ by a single Y chromosome, genes exhibiting expression differences are broadly distributed throughout the genome. In the liver, many of these differences are driven primarily by the pituitary growth hormone (GH) (59), which acts through signal transducer and activator of transcription (STAT)5b (13, 29) and hepatocyte-enriched nuclear factor (HNF)4α (30). These transcription factors, through direct and indirect mechanisms, are thought to regulate expression of many steroid and xenobiotic metabolizing enzymes, leading to observed patterns of sexually dimorphic expression in a number of important biological pathways (58).

Several investigators have examined sex differences in gene expression in the mouse liver with microarrays. Rinn et al. (41) surveyed the expression of 13,977 genes in male and female Swiss-Webster mice in the hypothalamus, kidney, liver, and reproductive tissues. Using criteria of fold change (FC) between sexes > 3 and P < 0.001, they found six genes to be sexually dimorphic in the liver, Cyp2b13, Cyp3a16, Cyp4a12, Hsd3b5m, Elovl3, and Keg1. Clodfelter et al. (12) studied liver gene expression changes in Stat5b knockout and wild-type mice, finding 1,603 differentially regulated genes, with 850 being male- and 753 female biased (P < 0.05 and FC > 1.5). A large study consisting of 344 mice comprising an F₂ cross between C57B/6J.apoE^−/− and C3H/HeJ.apoE^−/− strains (∼50% from each sex) produced two reports (57, 61) that examined sexually dimorphic gene expression in adipose tissue, brain, liver, and muscle. It was reported that 9,250 genes are dimorphic in the liver (P < 0.01 and FC > 1). In an effort to uncover regulatory pathways driven by genetic variation, the authors also performed gene expression quantitative trait locus (eQTL) mapping to find genomic loci responsible for transcriptional regulation of clusters of genes. The authors reported that one locus on medial chromosome (Chr) 5 regulated the expression of ∼80 genes, and the genes in this cluster were enriched for cell cycle, cell proliferation, and DNA replication. In addition, van Nas et al. (55) also illustrated that there are differences in both expression levels and gene coexpression between males and females. While these and other previous investigations into sex differences in gene expression in the liver made crucial contributions to the knowledge of sex-specific gene expression patterns in the liver, only a few studies examined the genetic basis of these effects (55, 57).

Even though sex-independent polymorphic local and distant QTLs, including several loci that control expression of large numbers of genes, were identified with the eQTL approach in the BXD panel (19) and confirmed in an independent cohort of inbred strains (20), relatively little is known about the extent of genetic control of sex differences in liver gene expression networks in the mouse, or in other species. We hypothesize that sex differences in liver xenobiotic metabolizing functions are largely due to variation in the expression of steroid metabolizing enzymes. We examined differences in constitutive liver gene expression between males and females across a genetically diverse population of BXD recombinant inbred (RI) mice derived from the C57BL/6J and DBA/2J strains. We examined both genes and pathways that are differentially regulated or differentially coexpressed between sexes in the mouse. eQTL mapping was performed separately for each sex, and while sex-specific eQTLs were observed, no sex-specific eQTL hotspots of transcriptional regulation were identified. Finally, we compared sexually dimorphic genes between murine and human liver by drawing a comparison to the data from Schadt et al. (44), who reported on liver gene expression in over 400 human subjects.

METHODS

Gene Expression Data

The details of mouse breeding, housing, RNA isolation, and gene expression are described elsewhere (19). Tissue collection (synchronized with respect to the time of day) was conducted at the University of Tennessee at Memphis and approved by the Institutional Animal Care and Use Committee. Briefly, gene expression in the livers of 37 strains of male and female C57BL/6J × DBA/2J (BXD) RI mice, the C57BL/6J and DBA/2J parentals, and the C57BL/6J × DB2/2J F1 generation (B6D2F1) mice was measured with the Agilent (Santa Clara, CA) G4121A microarray (20,868 transcripts). Additional details regarding animals, microarrays, data acquisition, processing, and analyses can be found at WebQTL.org (http://webqtl.org/dbdoc/LV_G_0106_B.html). The data have been deposited in the Gene Expression Omnibus (GSE17522).

Gene Expression Data Analysis

The raw data files from the Agilent scanner were uploaded into the University of North Carolina Microarray Database (http://genome.unc.edu), and a custom normalization was applied to the log₂ of the Lowess normalized Red/Green ratio of background-subtracted intensities. A nested, mixed analysis of variance (ANOVA) model (where batch, sex, and strain were considered as main effects and sex by strain was an interaction term) was fit to the data, and least-squares means were produced for each sex and strain along with P values for batch, strain, sex, and strain by sex interaction. The Significance Analysis of Function and Expression (SAFE) tool (5) was used to find significant differentially regulated Gene Ontology (GO) categories, as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Rather than employing an arbitrary significance cutoff, SAFE uses the entire data set and array permutation to produce empirical P values for each category or pathway. Significant categories were selected based on q values of ≤0.25, a 25% false discovery rate (FDR) for the entire list.

Sex-specific correlated gene clusters.

The Spearman correlation of each transcript with all 20,867 other transcripts on the array was calculated for males and females independently. Transcripts with correlations ≥0.7 in either males or females were retained separately for each transcript as highly correlated clusters. This threshold corresponds to a P value of 3.44 × 10⁻⁷, which is equivalent to a Bonferroni-corrected P value (for 20,868 transcripts) of 0.0072. Gene sets with absolute mean differences in correlation between sexes ≥0.35 were selected as sex-specific coexpression clusters. Both cutoffs were selected to produce reasonably short lists of interpretable results. The highly correlated sex-specific gene clusters from the above analysis were searched for GO and KEGG pathway enrichment with Fisher's exact test. Significant categories were selected with q values ≤0.1.

Transcription factor binding site enrichment.

oPOSSUM (28) was used to search for enrichment in transcription factor binding sites (TFBSs) among sets of genes. The options used, in the order that they appear on the website, were vertebrate taxonomic supergroup, using the top 10% of conserved regions, matrix threshold = 80%, upstream/downstream sequence = 5,000/5,000, Z-score cutoff ≥ 10. oPOSSUM does not provide P values for results.

Gene Expression Quantitative Trait Locus Mapping

eQTL mapping was performed by using single marker mapping in FastMap (21), which uses sample permutation to produce per-gene significance thresholds. Heritability was estimated with the strain intraclass correlation. This calculation was performed on variance components obtained with restricted maximum likelihood estimation in a random-effects ANOVA model implemented with SAS v9.1.2. In RI lines, this analysis can be performed by using the replicates within strain and sex as detailed in Reference 10. Within each strain, these replicates should be similar if the gene expression trait is heritable, or if it has a high proportion of variance explained by genetic factors relative to technical and environmental noise. The top quartile (5,217) of transcripts with the highest heritabilities (h² > 0.136) were kept. Mapping was performed on these 5,217 transcripts by using a set of 3,795 markers downloaded from WebQTL (http://www.genenetwork.org) for the BXD mice. Permutations (1,000) of the expression data were employed to establish significance cutoffs for each gene. Mapping was performed with the strain means, the male data, the female data, and the strain differences between males and females.

Mapping with the strain differences allowed us to test sex by marker interactions as follows. For a given transcript, we observe the expression G_ij, where i = 1 … n indexes the mouse strains and j = 0, 1 denotes the sex. Let a (fixed) genotype be denoted as S_i = 0,1, i = 1 … n, as it does not depend on the sex j.

The full model is given by:

For each sex separately (j = 0 and j = 1) the model takes the form:

The model (Eq. 2) contains S_i as the regressor, estimated on the sex differences (G_i1 − G_i0). This model produces the same estimate of b₃ as in the full model (Eq. 1). The permutation scheme for the full model (Eq. 1) involves permuting S_i relative to G_ij, keeping the strain pairs (G_i0, G_i1) intact in order to correctly estimate the empirical P value (11). To select a significance threshold that accounts for multiple testing across transcripts, we applied the q value algorithm (50) to the maximum P value for each transcript and selected all transcripts with q < 0.1. We also investigated whether polymorphisms that lie within the probe sequence produce false local eQTLs (14). Since C57BL/6J is the reference strain used to create probe sequences, probes querying the DBA/2J allele should have lower probe intensities among local eQTLs (2), creating false local eQTLs. We tested for this by performing Fisher's exact test, using a 2 × 2 table with the high-expression allele (C57BL/6J or DBA/2J) on one axis and SNP-in-probe/no-SNP-in-probe eQTL on the other, and found no evidence of an effect (P = 0.359). To assess a stringency of eQTL band analysis, eQTL bands were derived from the above analysis repeated 100 times, permuting the markers each time and counting the number of eQTL hotspots of each size. The probability of seeing an eQTL band of a given size was assessed from the empirical distribution of the number of eQTL bands of each size produced in all 100 permutations.

RESULTS

Characterization of Sexually Dimorphic Liver Genes in BXD Mouse Panel

Sex differences in gene expression were identified with a two-way (sex and strain of each sample) ANOVA model (Supplemental Table S1).¹P values for sex and strain were filtered separately, using q values (50) of 0.01, 0.05, and 0.1, which represent the FDR of the resulting gene lists, corresponding to FDRs of 1%, 5%, and 10%. At the smallest q value, we find a similar number of differentially expressed genes as that in a previous report (61) and replicate the finding that there are more female-biased genes than male-biased genes (Table 1). To further explore the differences, we plotted the mean expression value for each gene in males versus that in females (Fig. 1) and observed that, while there are more female-biased genes, there is also a large cluster of male-biased genes (Fig. 1, circle). The genes in this cluster (q ≤ 0.01, 0 ≤ male ≤ 2.6, −0.5 ≤ female ≤ 1.5) are significantly enriched for olfactory receptor activity (GO: 0004984, P = 9.6 × 10⁻²²), which consists of G protein-coupled receptors involved in signaling.

As expected, some of the most sexually dimorphic genes (Table 2) are located on the X and Y chromosomes, including inactive X specific transcripts (Xist), DEAD box polypeptide 3, Y-linked (Ddx3y), and eukaryotic translation initiation factor 2, subunit 3, structural gene Y-linked (Eif2s3y). However, most genes that exhibit considerable fold difference in expression (>2) between sexes reside on autosomes. Although most sexually dimorphic genes were found on the autosomes, there was no statistically significant enrichment for sexually dimorphic genes on autosomal chromosomes compared with the sex chromosomes. We have compared our results with two independent microarray studies that examined sex differences in liver gene expression in mice. The study of Clodfelter et al. (12) was conducted on Stat5b-null and wild-type mice. Using a threshold of P ≤ 0.01 for sex bias in both data sets, we found that 51.5% (410 of 796 genes from that study) were sexually dimorphic in both studies. A similar comparison with a study that used a much larger number of mice in an F₂ cross (61) yields a 27.2% overlap (9,250 genes were reported as sexually dimorphic with P ≤ 0.01, of which 8,087 mapped to Agilent 4121A microarray and 2,196 matched our results). There was no significant GO category or KEGG pathway enrichment (q ≤ 0.1) among the genes that replicate between the three studies.

We also examined the differential expression of several cytochrome P-450 (CYP) enzymes previously reported to be sexually dimorphic. Consistent with previous reports (42), we found that liver expression of Cyp2c55, Cyp2d9, Cyp4a12, Cyp7b1, and Cyp8b1 is male biased and expression of Cyp2a4, Cyp2b10, Cyp2b13, Cyp2c40, Cyp3a16, and Cyp3a41 is female biased. However, in contrast to some previous observations, we find that Cyp2j13 is male biased and Cyp2e1 is female biased in the BXD mouse panel.

To find transcription factors responsible for sexually dimorphic gene expression, we submitted lists of male- and female-biased genes to oPOSSUM (28), a web-based tool that searches for TFBS enrichment in lists of genes. In females, Hnf1α was significant (Z-score = 16.98) and was found in 395 transcripts. Hnf1α was more highly expressed in females (P = 7.67 × 10⁻⁴). In males, Pax5 was significant (Z-score = 11.54) and was found in 165 transcripts. Pax5 did not show differential expression between sexes (P = 0.6).

While the examination of expression differences between individual genes increases our understanding of physiological differences between the sexes, it is often informative to search for functionally related gene sets that are differentially expressed. SAFE (5) is a permutation-based software tool that assesses the significance of differentially expressed gene sets, in this case between males and females, using the entire data set rather than a predefined “significant” gene list. We found 89 significant KEGG pathways (q ≤ 0.10, Supplemental Table S2a) of which many are reflective of xenobiotic metabolism pathways such as Drug Metabolism-Cytochrome P-450 (KEGG 00982, q = 0.01; Fig. 2A, Supplemental Table S3a), Retinol Metabolism (KEGG 00830, q = 0.01), Glutathione Metabolism (KEGG 00480, q = 0.05), and Fatty Acid Metabolism (KEGG 00071, q = 0.08). Similarly, of 45 significant GO Molecular Function categories (q ≤ 0.25, Supplemental Table S2b) many are also highly relevant to xenobiotic metabolism in the liver, such as Oxidoreductase Activity (GO: 0016491, q = 0.03; Fig. 2B, Supplemental Table S3b), Cysteine-type Endopeptidase Activity (GO: 0004197, q = 0.04), Iron Ion Binding (GO: 0005506, q = 0.05), Heme Binding (GO: 0020037, q = 0.08), Sulfotransferase Activity (GO: 0008146, q = 0.13), and Glutathione Transferase Activity (GO: 004364, q = 0.14). There were 30 significant GO Biological Process categories (q ≤ 0.25; Supplemental Table S2c). The most significant (q = 0.11) were GO: 0006766 (Vitamin Metabolic Process; Fig. 2C, Supplemental Table S3c) and GO: 0051186 (Cofactor Metabolic Process). There were eight significant GO Cellular Component categories (q ≤ 0.25; Supplemental Table S2d). The most significant categories were Vesicular Fraction (GO: 042598, q = 0.09; Fig. 2D, Supplemental Table S3d) and Microsome (GO: 0005792, q = 0.09), both of which include a large number of sexually dimorphic CYPs that are located on the endoplasmic reticulum. Also included in the list was Peroxisome (GO: 0005777, q = 0.1), which contains the female-biased hydroxyacid oxidase (Hao)3 and ATP-binding cassette D2 (Abcd2), as well as the male-biased Hao1.

Comparison of Sex Differences in Liver Gene Expression Between Mouse and Human

Sex differences in human liver gene expression have been largely unexplored (58). We compared the sexually dimorphic genes in the BXD liver data set with a recently published large human liver gene expression data set (44) obtained from the Gene Expression Omnibus. Expression measurements for over 39,000 transcripts in over 400 human samples (both male and female) were taken with a custom microarray. We determined that 7,342 transcripts were sexually dimorphic in this data set (Student's t-test between sexes for each transcript at q ≤ 0.10 and p ≤ 0.03). Murine orthologs could be mapped with high certainty only for 2,850 of these genes (due to rather insufficient annotation of the custom array used in the human study), and 1,258 (44.1%) of these were also sexually dimorphic in the mouse. We compared the direction of expression difference between sexes in this subset of 1,258 transcripts (Supplemental Table S4) and found that 747 (59.4%) genes were sex biased in the same direction in mouse and human and 511 (40.6%) were biased in the opposite direction. The genes that exhibited unidirectional sex bias between species were enriched in GO Biological Process categories Oxidation Reduction (GO: 0055514, P = 3.05 × 10⁻⁹, q = 5.60 × 10⁻⁶) and Glutathione Metabolic Process (GO: 0006749, P = 1.22 × 10⁻⁴, q = 0.045) and GO Molecular Function categories Monooxygenase Activity (GO: 0016491, P = 1.68 × 10⁻⁷, q = 1.48 × 10⁻⁴) and Catalytic Activity (GO: 0003824 P = 2.47 × 10⁻⁵, q = 7.24 × 10⁻³).

Gene Expression Quantitative Trait Locus Analysis

Quantitative trait locus (QTL) mapping is a statistical technique that finds associations between phenotype and genotype in a genetically segregating population (32). Here, we performed eQTL mapping on the male and female data separately to determine whether genetic polymorphisms are responsible for variation in transcript abundance. First, we filtered the transcript expression data by heritability, retaining the top 5,217 transcripts with heritabilities greater than the 75th percentile (h² > 0.136). We performed eQTL mapping on the male and female data separately and also on the mean and difference of the male and female data. There were 1,638 significant eQTLs (q ≤ 0.1, P ≤ 0.079) found for liver gene expression in males (Supplemental Table S4a) and 2,076 (q ≤ 0.1, P ≤ 0.105) in females (Supplemental Table S4b). Performing eQTL mapping on the mean of the male and female data uses all of the data and increases the power to detect additive effects. Using the full data set, we found 2,925 significant eQTL (q ≤ 0.1, P ≤ 0.0005) (Supplemental Table S4c). The difference between the male and female data produced only five significant eQTLs (q ≤ 0.1, P ≤ 0.199), which suggests either that genotype by sex interactions are rare in the BXD strains or that this population is not of sufficient size.

First, we explored differences in patterns of eQTL locations between sexes by plotting the genomic locations of each eQTL versus the transcript location (Fig. 3, A and B). The strong diagonal band represents local eQTLs (43), for which the eQTL is located within 1 Mb of the transcript. Points off the diagonal represent distant eQTLs (43) for which the eQTL is located farther than 1 Mb from the transcript or on another chromosome. A histogram of the number of significant eQTLs associated with each genomic marker is also shown for each sex (Fig. 3, D and E). Markers associated with a larger number of transcripts than expected by chance, termed eQTL bands, may lie near genes that regulate the expression of clusters of genes. Legitimate concern has been expressed over the validity of eQTL bands, and a permutation approach to assess their significance has been proposed (7). We applied this approach (see methods) and colored the eQTL bands by their empirical P value. In males, there are four significant (P ≤ 0.05) eQTL bands on Chr 2, 7, 12, and 17 and in females four eQTL bands on Chr 1, 2, 12, and 17. The genomic regions covered by the eQTL bands on Chr 2, 12, and 17 overlap in males and females and contain 25, 33, and 20 transcripts in common, respectively.

Second, to further explore the similarities in genetic control of liver gene expression between sexes, we compared the eQTL overlap between males, females, and all of the data as well as the overlap of local and distant eQTLs (Fig. 4, A–C). It is evident that a larger number of local eQTLs are shared between the two sexes than distant eQTLs. The strongest distant eQTL band on distal Chr 12 is shared between sexes and has been previously reported (19). There are 87 transcripts regulated by this locus in the male data and 45 transcripts in the female data; however, the two sets of genes are not identical. The union consists of 101 genes, and the intersection of these two sets consists of only 33 genes. Several GO categories are significantly overrepresented in the union of these genes: Response to Cytokine Stimulus (GO: 0034097, P = 0.0050), Calcium-Mediated Signaling (GO: 0019722, P = 0.0070), Serine-Type Endopeptidase Inhibitor Activity (GO: 0004867, P = 0.0064), and Cyclin-Dependent Protein Kinase Activity (GO: 0004693, P = 0.0098). A clade of serine protease inhibitors contains a putative regulator for these genes, as discussed in Reference 20. GO and KEGG pathway enrichment analysis for other sex-specific eQTL bands on Chr 1, 2, 7, and 17 did not produce significant categories or pathways (q ≤ 0.1). No TFBS enrichment (P ≤ 0.05) exists for genes in theses eQTL bands.

To examine whether the transcripts that have significant eQTLs in both sexes are identically regulated in males and females, we plotted the male versus female eQTL locations for 1,030 transcripts (of 1,152) with known genomic locations (Fig. 4D). For 757 (73.5%) of the transcripts, the eQTL is in the same location (diagonal band in Fig. 4D). For 273 (26.5%) transcripts, one sex had a single significant eQTL while the other had two, one at the same location and another elsewhere (off-diagonal points in Fig. 4D), suggesting that sex is a factor in the regulation of these transcripts. Still, for all transcripts that had a single significant eQTL in one sex and two in the other, we found that either sex had a less significant (P ≥ 0.05) second eQTL. We examined the direction of the eQTL effect for each allele in the 757 eQTLs that replicate between males and females. We found that transcript expression was higher when the C57BL/6J allele was present for 394 (52%) eQTL and that expression was higher when the DBA/2J allele was present for 363 (48%) eQTLs. In all cases, the effect size was in the same direction in both sexes (Supplemental Fig. S1).

Third, we compared the sex-specific eQTL locations obtained from BXD RI mice with those from an F₂ cross study of Wang et al. (57). While the comparisons of eQTL results are challenging because of the varying genetic backgrounds between populations, differences in microarray probe sequences, experimental conditions, and statistical techniques, and likely environmental differences inherent in different animal facilities, it is likely that eQTLs that overlap may serve as strong candidate loci for further investigation (39). Wang et al. (57) reported mouse liver eQTLs for 8,349 transcripts (P < 0.001), and we subselected the 3,855 most significant (P < 10⁻⁶) eQTLs for 2,667 transcripts. An eQTL was deemed reproducible between two studies if the maximal eQTL interval intersected, and we found that ∼25% of local and 8.5% of distant eQTLs satisfied these criteria (Table 3). This result is consistent with previous observations that local eQTLs are more reproducible than distant eQTLs (39).

It has been reported that genetic polymorphisms that are located within transcript probe sequences may alter apparent transcript intensity in an allele-specific manner and hence create false local eQTLs (2, 14). To investigate this in our data, we intersected the microarray probe sequences with a high-density single nucleotide polymorphism (SNP) data set (52). We did not find evidence for this in our data (see methods). Transcripts with SNPs in the probe sequences are identified in Supplemental Tables S5a, b, and c.

Gene Correlation Exhibits Distinct Patterns from Gene Expression Between Sexes

While it is clear that many genes differ in expression between the sexes, there is also a more subtle type of variation that involves variability in the correlation among genes within each sex, often without pronounced difference in expression. The sex-specific clusters of such differentially correlated genes may represent complex networks of transcriptional regulation that are important in physiology and pathophysiology.

On a genomewide level, we examined the extent of correlation in gene expression between sexes. We selected genes with Spearman correlation ≥0.7, an arbitrary cutoff that reflects a balance between lower values, which quickly produced so many genes that analysis proved computationally intractable, and higher values, which yield only a small set of genes. We searched for genes exhibiting three types of patterns. The first type contains genes with stronger correlations in males than in females (i.e., male-biased correlations). Member of RAS oncogene family Rsb30 (Fig. 5A) is a representative gene that does not differ in expression between sexes yet exhibits strong correlation in males. The genes that are highly correlated (>0.7) with Rsb30 in males do not show significant enrichment for any GO category or KEGG pathway. However, the GO categories with the lowest P values (P < 0.0023) relate to male sexual development, suggesting that there may be some relevance to liver metabolism in males. The second type of pattern, female-biased correlations, is exemplified by neutrophil cytosolic factor (Ncf)2 gene (Fig. 5B). Ncf2 is involved in superoxide production, and the genes highly correlated (>0.7) with Ncf2 in females are enriched for Immune Response (GO: 0006955, P = 8.24 × 10⁻¹¹). Finally, there are genes that have an equal level of correlation in both sexes regardless of expression levels, exemplified by potassium channel K5 (Kcnk5) gene (Fig. 5C).

To establish whether genes that have the highest difference in expression between sexes (Table 2) also would exhibit strong correlation, we selected genes that had the largest number of correlations with other genes in one of the sexes (Table 4). It is noteworthy that there is only one gene (Elovl3) that is both differentially expressed between sexes and exhibits strong correlation in only one sex (males).

There are also differentially correlated gene clusters in which a set of genes is highly correlated (≥0.7) with each other in one sex but not in the other (Fig. 6). We focused on the genes that had strong correlation in one sex but considerably weaker correlation in the other (mean absolute correlation difference ≥0.35 between the sexes). There were 33 differentially correlated clusters in males consisting of 254 unique genes (Fig. 6, A–C) and 28 clusters in females containing 332 unique genes (Fig. 6, D–F). The male-specific clusters are related through a complex network of correlations (Fig. 6A) that are neither primarily driven by expression differences between sexes (Fig. 6B) nor correlated in females (Fig. 6C). The female clusters exhibit a similar pattern in the opposite direction (Fig. 6, D–F). Finally, we searched for clusters of genes that are strongly correlated (≥0.7) in both sexes. There were 2,639 genes with high correlation in both sexes (Fig. 6, G and I), and most of these genes were also differentially expressed between sexes (Fig. 6H).

DISCUSSION

Sex-Dependent Differences in Gene Expression in Liver: Physiological and Toxicological Significance

In this study, we compared liver gene expression of males and females in a panel of BXD RI mice in an effort to understand whether genetic background affects sexually dimorphic gene expression. Our results show that sex differences in liver gene expression are largely determined by the differences in endogenous steroid metabolism between sexes, and to a lesser extent by the differences in genetic polymorphisms. Indeed, the role for liver steroid metabolism in sex differentiation is well established, and it is known that the neonatal exposure to androgens is a crucial imprinting mechanism (24, 25). The physiological process of sexual differentiation requires differentiation of hepatic processing of steroids in fetal liver, and this differentiation results in an epiphenomenon of sex differences in xenobiotic processing.

It has been long known that females exhibit a higher incidence of drug-induced liver injury (34). There are significant pharmacokinetic differences between the sexes that are likely due to significant differential expression in CYPs, glutathione S-transferases (Gst), sulfotransferases (Sult), UDP glucuronosyltransferases (Ugt), and ATP-binding cassette (ABC) transporters (49). Consistent with previous work, we found a large number of differentially expressed CYPs, including the previously reported female-biased genes Cyp2a4 (54), Cyp2b9, and Cyp3a16 (58) and the male-biased Cyp1a2 (58) and Cyp2d9 (54). Overall, 39 of the 69 CYPs on the array (56%) were differentially expressed, with 26 CYPs biased toward high female expression and 13 male biased. However, we also found Cyp2j13 to be male biased and Cyp2e1 female biased, suggesting that sex differences of important xenobiotic enzymes may also be affected by a complex interaction with genetics.

Higher expression of most Sult genes was found in females, including Sult3a1 (3), Sult1a1, Sult1c1, Sult1c2, Sult1d1, and Sult1e1. In agreement with previous reports, we find that Gstp1 is male biased and Gsta3 is female biased (36). We also find that Gstm7 is male biased and that most differentially expressed Gsts are female biased (11 vs. 23). We found that seven of nine (77.8%) Ugt genes on the array were differentially expressed at q values ≤0.01. Ugt genes were found to be exclusively male biased, which is consistent with the finding that Ugt expression is regulated by male GH production in the liver (8). ABC transporters affect the disposition of many endo- and xenobiotics. About 44% of these genes on the array (out of 50) were more highly expressed in females than in males, whereas only three, including Abcg2 (35), were more highly expressed in males. Similarly, male bias was observed in expression of epoxide hydrolases 1 and 2 (Ephx1, Ephx2) and critical cofactor supply genes hexose-6-phosphate dehydrogenase (H6pd) and nicotinamide nucleotide adenyltransferase 1 (Nmnat1).

The retinol binding protein (Rbp1), alcohol dehydrogenase 1 (Adh1), and retinol dehydrogenases 6 and 9 (Rdh6, Rdh9) exhibit strong female bias in liver gene expression. Estrogen is known to control retinoic acid metabolism in many tissues (33), and there is considerable overlap between alcohol and retinol metabolism (27). Animal models for the study of alcoholism suggest that, in general, females drink more alcohol than males and there are sex-related differences in drinking patterns (15). Interestingly, Rbp1 has been proposed as a locally regulated candidate gene for alcohol preference in the mouse based on the meta-analysis of the brain gene expression data (38). Furthermore, Rbp1 expression in the liver shows a strong negative correlation (−0.783) with alcohol preference in BXD strains (WebQTL Phenotype no. 10158; Ref. 23), further supporting the hypothesis that Rbp1 may play a role in alcohol preference through its involvement in alcohol metabolism in multiple tissues. However, Rbp1 did not have a significant liver eQTL in either sex and did not show significant differential correlation. Likewise, previous QTL studies examining alcohol preference in BXD and AXB mice have not revealed a QTL for alcohol preference in either sex that maps to the genomic location of Rbp1 (15). Thus, while additional studies that address the role of Rbp1 in alcohol preference and metabolism may yield additional information, it is also possible that the genetic predisposition for elevated alcohol consumption extends to thousands of transcripts forming many functional groups (38), which may render a candidate-gene approach futile.

It should be noted that not all sex-specific gene expression differences found in the mouse translate to humans. While the mapping of mouse microarray probes to human array probes limited the comparison, we found that many genes are biased in the same direction in each sex, suggesting that the mouse is an appropriate model for studying human liver disease and function. Still, even though both human and mouse sex-specific liver gene expression differences are prominently enriched with xenobiotic metabolism, important gene-directed distinctions exist. One example is flavin-containing monooxygenase 3 (Fmo3), a key player in metabolism of numerous substrates, including nicotine, tertiary amines, drugs, carbamates, and organophosphates (9). It is a prominent female-biased gene in the mouse liver, but in humans Fmo3 is known to be expressed in livers of both sexes (31). Thus extrapolations between studies and species should be exercised with caution.

Transcriptional Control of Sex Differences in Liver Gene Expression

The sex difference in expression of many drug metabolism and other liver-expressed genes has been reported to be largely regulated by the temporal pattern of plasma GH release by the pituitary gland (13, 58). These differences are most pronounced in rodents, where plasma GH profiles are highly pulsatile in males but are nearly continuous in females. Hnf4α and Stat5b are key transcription factors that respond to sexually dimorphic GH secretion pattern responsible for upregulation of male-biased genes (60). Our study did not identify enrichment in binding sites for GR, Stat5b, or HNF4α among sexually dimorphic genes in the liver of the BXD mouse panel. While Stat5b and Hnf4α undoubtedly play a major role, and can be identified in comparisons of male and female mice of the same strain, our result suggests that the regulation of sex-dependent gene expression on the population level may be complex, and both direct and indirect linkages in the network of regulators that maintain sex-dependent liver homeostasis are yet to be determined. Instead, we found that HNF1α was enriched among genes more highly expressed in females and PAX5 was enriched among genes with high expression in males. Hnf1α knockouts survive only until weaning, suggesting that Hnf1α is critical for normal liver development and function (40). In our data, Hnf1α was more highly expressed in females (P = 7.67 × 10⁻⁴). Hnf1α upregulates the expression of several acute-phase proteins (4), and we found several, such as albumin (P = 5.55 × 10⁻⁰⁹), α-1-antitrypsin (P = 4.29 × 10⁻¹⁰), angiotensin (P = 3.49 × 10⁻¹⁵), and insulin-like growth factor I (P = 7.45 × 10⁻⁰⁸) that are more highly expressed in females. However, the cholesterol oxidizing genes that are regulated by Hnf1α, Cyp7a1 (P = 0.0129), and Cyp27a (P = 0.0828) are not differentially expressed between the sexes. Pax5 did not show differential expression between sexes (P = 0.633), and its function has primarily been studied in the central nervous system and B cell development (1). Pax5 has been shown to regulate gene involved in cell signaling, adhesion, and migration in B cells (46), but further research is needed to validate whether it is involved in the same functions in the murine liver.

Genetic Component of Sex-Specific Gene Regulation

It is likely that mechanisms other than growth and sex hormones may also contribute to regulation of sex-specific hepatic gene expression. QTL mapping has been adapted to gene expression data to find expression QTLs that represent genomic loci responsible for differential transcriptional regulation (45). The sex-specific eQTL hotspots discovered in our study provide some evidence for the sex-specific differences in genetic control that may be also important. We have assessed the existence of eQTLs in males and females by considering gene expression data from each sex separately. While three loci identified as putative regulatory hotspots are shared between males and females, including the strongest regulatory locus on Chr 12 (19, 20), there were few additional significant sex-specific putative regulatory loci.

Ghazalpour et al. (22) also performed eQTL mapping using a C57BL/6J × DBA/2J intercross but did not include all eQTL results in their publication, so we were not able to perform a detailed comparison with their study. It was reported by the authors that significantly enriched eQTL bands on Chr 3, 6, 16, and 19 were identified in females. These bands did not reproduce in our study, and there are several possible explanations for the discord. The microarrays used, while both based on Agilent technology, were different in probes, transcripts, as well as regions of overlapping transcript. Most importantly, whereas the BXD mice in our study were fed a normal chow diet and were killed at a mean of 70 days, a high-fat diet and a 16-mo study length were employed by Ghazalpour et al. (22).

In addition, Yang et al. (61) examined the interaction of sex with genetic background in liver gene expression. It is of relevance that our study was based on a cross between C57BL/6J and DBA/2J strains, whereas Yang et al. (61) used an intercross between C57BL/6J and C3H/HeJ. The genomes of these three strains differ at 37.2% of known SNPs [http://compgen.unc.edu/DisplayIntervals/DisplayIntervals.html], which may at least partially account for the lack of overlap. This also means that only 62.8% of the mouse genome has been explored in these two studies and that much work remains to be done in studying the effect of genetics on sex differences. Collectively, since the number of genes regulated by the QTLs in a sex-specific manner is relatively small and is poorly reproducible between studies, we posit that it is unlikely that genetic polymorphisms play a major role in controlling sex-specific gene regulatory networks.

Our data show that there were several hundred sex-biased eQTLs and slightly more than 2,925 total eQTLs out of 5,217 transcripts that were probed. This raises the question of how transcripts without significant eQTL are regulated. First, it is plausible that most genes may be regulated via signaling and germline effects may be fairly narrow. Second, the genetic diversity of the BXD panel used in this study is rather limited, and other panels should be used to detect additional regulatory loci. Finally, the modest sample size of the BXD panel limits the power and allows detection of eQTLs with large effect sizes. As gene expression is a complex trait, it is expected that genes with heritable expression may also be regulated by multiple loci.

Correlation Analysis Reveals Sex-Specific Gene Expression Networks

While differential expression plays a large role in sexual dimorphism, it was also shown that there is a more subtle, yet important, role played by varying degrees of differential correlation among gene expression networks between the sexes (56). Correlated genes represent functionally related expression networks that are often under common transcriptional regulation (10). Several of the clusters, such as the male immune response cluster and the female ribosomal cluster, show functional coherence that has been annotated in GO or KEGG. However, many of the clusters discovered here have not been previously reported and represent coexpression networks that may help to explain differences between male and female liver metabolism.

One of the male-specific correlated clusters with the largest differential correlation has a mean absolute correlation of 0.68 in males and 0.30 in females. It is significantly enriched in acute-phase response genes (GO: 0006953, q = 4.6 × 10⁻⁸) and contains three serum amyloid A (Saa) genes (Saa1, Saa2, Saa3) along with 16 other highly correlated genes. With the exception of the Saa genes, which are located together on mouse Chr 7, the rest of the genes in this cluster are not colocated in the mouse genome. The Saa genes produce acute-phase proteins in response to inflammatory stimuli (16) and are associated with insulin resistance in mice (47). Several other genes in this cluster have also been implicated in inflammation or the immune system, including lipocalin (Lcn)2 (63), intercellular adhesion molecule (Icam)1 (26), orosomucoid (Orm)1 (48), orosomucoid (Orm)2 (6), LPS-induced TNF-α factor (Litaf) (53), and β-galactoside α2,6 sialyltransferase (St6gal)1 (62). Even though other genes in this cluster have been not implicated in the inflammatory response, their strong correlation may suggest that they may also be involved in immune response in liver.

In females, one of the most differentially correlated clusters (mean absolute male cor. = 0.33, mean absolute female cor. = 0.66) is significantly enriched in ribosomal genes (GO: 0003755, q = 4.9 × 10⁻⁶) and contains 12 genes, including ribosomal proteins L19, L21, S14, and S27. Another cluster of 11 genes (mean absolute male cor. = 0.34, mean absolute female cor. = 0.69) is enriched for antigen processing and presentation of exogenous peptide antigens (GO: 0002478, q = 2.47 × 10⁻⁶), and similar to observation in the male clusters, many genes that highly correlate may have yet undiscovered functional relationships.

The genes that are strongly correlated (≥0.7) in both sexes form two large clusters (Fig. 6, G and I, top right and bottom left). The bottom left cluster, consisting of 620 genes and exhibiting no sex bias in expression, shows no significant GO category or KEGG pathway enrichment (q < 0.05). The top right cluster, consisting of 740 genes with male-biased expression, shows enrichment for ion channel activity (GO: 0005216, q = 6.02 × 10⁻⁴), a category that includes calcium, sodium, and potassium channel transporters. These genes show significant TFBS enrichment for several transcription factors, including myeloid zinc finger 1 (Mzf1) (P = 3.09 × 10⁻⁴), which is known to activate calcium-activated potassium channels (51). The third largest cluster consists of 245 genes with a female-biased expression pattern and shows enrichment for immune response (GO: 0006955, q = 8.32 × 10⁻¹²). These genes show TFBS enrichment for transcription factor SpiB (P = 7.14 × 10⁻⁵) and signal transducer and activator of transcription Stat6 (P = 1.6 × 10⁻³). SpiB is known to be associated with the transcriptional regulation of immunity (17), and Stat6 controls immune response in many organs, including the liver (18).

In conclusion, this study assessed the hepatic gene expression networks between two sexes in a BXD RI mouse population. Unlike studies carried out in only one strain, the genetic diversity in the BXD population, combined with a comparison to human data, sheds light on the global sexually dimorphic gene expression that is highly relevant to the human population. Studies in mouse populations may allow us to generate testable hypothesis with regard to gene regulation networks and modes of action of toxicants. Also, our findings have several implications for toxicological studies. It is not only the genes that show significant sex differences that are important in assessing toxic responses in males and females; a large number of functionally related genes with small differences between sexes could also contribute to many sex-biased phenotypes. Caution should be taken when applying the results from a study that was carried out only in male subjects (which is very common in toxicological studies) and applying the conclusions to the whole population. Finally, sexually dimorphic gene expression, at least partly, may be due to the genetic control. Other genomic information, like genome structure similarities and identical by descent regions, when combined with eQTL mapping, can facilitate the process of finding candidate regulatory genes.

GRANTS

This study was funded, in part, by grants from the National Institutes of Health (R01-ES-015241, R01-AA-016258, U01-AA-014425, and P42-ES-005948). E. J. Chesler and R. Yordanova were supported by the Office of Biological and Environmental Research, U.S. Department of Energy, under Contract DE-AC05-00OR22725 with UT-Battelle, LLC.

DISCLOSURES

The authors declare that they have no competing financial interests.

FOOTNOTES

• ¹Supplemental Material for this article is available online at the Journal website.