Age and sex dependent changes in liver gene expression during the life cycle of the rat

The liver's ability to process exogenous chemicals and protect itself from injury is influenced by many factors, including an individual's age [1], sex [2], ethnicity [3], genotype [4], and diet [5]. Age- and sex-dependent factors are decidedly innate and comprise two of the most important variables in gauging liver function and responsiveness to chemical insult. Because the liver is the first and primary line of defense against the potentially toxic effects of oral xenobiotics [6, 7], inquiry into its gene expression profiles to assess potential susceptibilities was performed. Not all potential causes of susceptibility were examined, only those gene expression changes associated with sex and age. Some important gene regulatory events may not be evident at the RNA level, in the tissue examined (liver), or in the absence of chemical challenge. Mechanisms of susceptibility may also be dependent on non-transcriptional events such as protein expression and modifications, as well as polymorphisms in specific genes and proteins, which lie beyond the scope of the current study. Nonetheless, capturing baseline gene expression profile measurements is a necessary first step in characterizing the suite of genes differentially expressed at various ages.

Age-dependency in susceptibility to drug induced liver injury (DILI) has been established clinically. For example, age has been shown as a predisposing factor for adverse effects of the antiepileptic drugs carbamazepine and valproate [8]. In the case of carbamazepine, elderly patients appear to be at increased risk of developing blood dyscrasias and liver reactions whereas risk for adverse liver reactions to valproate is associated with young patients. Sex-differences in DILI have also been reported: for example, acetaminophen seems to exhibit female-specific susceptibilities in clinical and epidemiology studies [9, 10]. However, it remains to be seen if higher rates of female DILI are due to intrinsic (biology) versus extrinsic (increased use) factors. Whether susceptibilities to DILI are generally more common in males versus females have not been as clearly demonstrated clinically when considering all drugs [11–13]. This suggests that when sex-specific susceptibilities are present, the sex-specificity varies from drug to drug.

It is well understood that the liver's drug metabolizing capacity is a function of its gene expression levels [14]. Furthermore, transcriptome profiling continues to play an increasing role in estimating toxicity and efficacy of various drugs in the liver in preclinical and clinical evaluations [15, 16]. Several prior studies have evaluated hepatic gene expression differences between young and old animals or evaluated sex differences on a relatively limited scale [17–22]. However, a more comprehensive evaluation of the whole genome during time points spanning immature, pubertal, early adult and aging adult animals, in both sexes, would provide greater interpretive power to assess life cycle liver gene expression and its impact on drug safety and disease. This study provides a comprehensive look at liver gene expression in untreated male and female rats throughout the entire life cycle in order to assess the basal level profiles of liver genes that may underlie patterns of susceptibility due to sex and age.

A comprehensive time-course study, comprising nine age groups spanning the entire life cycle (between 2 weeks and 2 years of age) of Fischer 344 rats was performed in both male and female animals (Figure 1A). Terminal body weights were recorded for each animal (Figure 1B) and the corresponding growth curves are consistent with historical findings [23]. An original number of sixteen animals of both sexes was included in the 104 week group. However, consistent with previous reports [24], 10 male animals were found to be moribund or died prior to 104 weeks and were thus removed from the study. No female animals were removed early. Pathological examinations were performed on animals at 52, 78 and 104 weeks of age. Neoplastic histopathological findings in the liver which showed sex-related trends included a hepatocellular adenoma in only one 104 week male. Non-neoplastic findings included basophilic foci in females (3 and 12 animals at 78 and 104 weeks, respectively) but not in males. Also, male predominant findings of bile duct hyperplasia with increasing incidence with age were observed. In summary, the incidences of neoplastic and non-neoplastic lesions were typical of those observed in historical control Fischer 344 rats of comparable age [25].

Whole genome profiling of rat liver gene expression was performed (n = 5) in each sex and age group using Agilent microarrays. These rat whole genome microarrays contain 43,379 features or spots of which 23,233 have Entrez Gene IDs, 16,327 of which are unique. Gene expression data were entered into ArrayTrack™ [26], the Food and Drug Administration's database for microarray data storage, processing, analysis and visualization that was created at the National Center for Toxicological Research (NCTR). The composite liver gene expression profiles for all genes (unfiltered) displayed good reproducibility in general between animals of the same sex and age group with average Pearson's correlations of R = 0.986 between biological replicates of the same sex and age. Normalized intensity values were analyzed using a two-way ANOVA (sex, age) to calculate a p-value for statistical evaluation. Data for each spot across all arrays were mean-scaled to calculate relative fold-change differences in expression for each sex and age. The initial filtering criteria of p < 0.05 for sex or age difference with an absolute fold-change greater than or equal to 1.5 in relative expression at any age was used to define a set of 7,951 differentially expressed features comprising 3,770 unique Entrez Gene IDs or differentially expressed genes (DEGs). Using these 7,951 differentially expressed features, global expression pattern analysis was performed using three-dimensional principal component analyses (3D-PCA) within ArrayTrack, wherein the three principal components capturing the highest amount of variability in the DEG data were plotted (Figure 2). These three principal components captured approximately 58% of the total variability present in the data (31.4%, 16.2% and 10.2% for PC1, PC2 and PC3, respectively). Three-dimensional PCA facilitates the visualization of individual animal expression data within each sex and age group and the global relationships between one another. Male and female animals showed clear separation from each other for most of the life cycle. Furthermore, most, if not all age groups displayed temporal continuity with neighboring age groups in a contiguous manner (i.e., expression profiles from 6 week old animals of each sex lay in-between those of 5 week and 8 week animals in 3D-space, etc.). Lastly, the largest sex-differences in expression occurred between 15 and 52 weeks of age, whereas the biggest age-related differences occurred, as expected, between the youngest animals (2 weeks) and adult animals (52 weeks).

A second clustering method, k-means clustering (JMP Genomics 4.0, SAS. Cary, NC), was performed to show individual gene profiles and to find groups of genes with similar patterns of expression associated with specific stages of development or aging (Figure 3). Male and female data for each array feature were allowed to cluster independently such that male and female data for a given gene may group into separate clusters. The number of clusters was determined empirically according to the fewest clusters required to achieve a minimum correlation radius of 0.7 between any individual profile and its cognate cluster members, resulting in 30 clusters. K-means cluster analysis allowed for grouping of genes based upon similarity of expression profiles across all ages. The number of features included in each cluster ranged from 16 to 2,252. Many clusters reflected temporal patterns that could be readily interpreted in the physiological context, such as high and low perinatal expression (clusters 2 and 10, respectively) or up- (clusters 6, 9, 14, 23) and down- (clusters 13, 18, 25, 26) regulation in the aging animals. Clusters also displayed both sustained (clusters 2, 10, 11, 21) and transient (clusters 3, 4, 15, 22, 24) patterns of regulation. Although not all of the clusters are readily interpretable, knowledge of rat growth patterns, sexual maturation, and basic liver function assist in interpreting time-course gene expression data. For example, cluster 2 and cluster 10 (Figure 3) are negatively correlated with each other yet both exhibit early and substantial changes (between 2 weeks and 5 weeks) followed by negligible change for the rest of the life cycle. Clusters 14 and 25 also negatively correlate with each other and show clear and consistent upward and downward trends, respectively, in the late age groups (between 52 weeks and 104 weeks). Thus, temporal patterns can be clearly distinguished for genes potentially involved in early development vs. late aging based upon these life cycle clusters.

Sets of genes for which expression was at least 2-fold greater in the opposite sex were determined for each age group (p-value < 0.05; fold-change > 2 between sexes) from the set of 3,770 DEGs described above. The numbers of these sex-predominant genes at each age are presented in Figure 4. The resulting plots indicate a notable period from 6 to 21 weeks of age in which more than 70% of the genes showing sex-predominance are expressed at a higher level in females than males. From the ages of 21 to 78 weeks, sex-predominant gene expression was higher in males than females.

An unsupervised, pathway-analysis approach was performed to identify functional categories of gene response. Analysis software Ingenuity Pathway Analysis (IPA) was used to assist in biological interpretation of the expression data from the initial list of 3,770 DEGs with unique EntrezGene IDs. Pathways which exhibited large numbers of DEGs included metabolism of xenobiotics by cytochrome P450, fatty acid metabolism, tryptophan metabolism, and acute phase response. Genes involved with liver cholestasis, liver damage and liver necrosis/cell death were highly represented in toxicity-related functional analysis of the DEGs. Furthermore, gene networks represented by the DEGs included seven high-scoring IPA categories, including cell mediated immune response, cell growth and proliferation, lipid metabolism, organ morphogenesis, cell movement, and hepatic system disease.

Since data from all age groups were included in the pathway analysis, the results are generic to the entire life cycle and do not provide meaningful insight into specific periods of the life cycle which may feature unique expression patterns that would underlie age-related susceptibilities. Thus, a more focused functional analysis of targeted age groups was implemented. Adult aging-related clusters (up-regulated with age: clusters 6, 14, 23; down-regulated with age: clusters 13, 18, 25, 26; Figure 3) from the k-means cluster analysis were selected and their corresponding gene lists were queried for pathways. Genes associated with organ injury/abnormality (Crp (C-reactive protein), Il6st (interleukin 6 signal transducer), Stat3 (signal transducer and activator of transcription 3)) and lipid metabolism (Pon1 (paraoxonase 1), Apoa1 (apolipoprotein A1), Cyp27a1, Abcg1 (ATP-binding cassette G1)), as well as toxicity-related functions of liver cholestasis, necrosis and cell death, were all up-regulated during the second year of the rat life cycle. These gene family results overlap with the top pathways represented in the up-late clusters including acute phase response signaling, complement system genes and xenobiotic metabolism. The clustering of up-regulated gene families related to immune response and cell death signaling pathways match well with a presumed accumulation of aging hepatocytes [27], the potential hepatic accumulation of fatty acids, and increased demand for metabolic capacity for lipids and xenobiotics which accumulate in the fat of aging animals [28]. Genes that clustered in a pattern of down-regulation in the adult aging rats are more difficult to interpret as loss of expression may itself be a product of the aging process [29]. However, a number of represented gene networks in the aging-related down-regulation clusters included cell signaling genes particularly associated with the nervous system. A major hub of many of the top represented gene networks in both the up- and down-regulated, aging-associated clusters was NF-κB (nuclear factor κB).

As a result of the global cluster analysis methods, various gene relationships in the expression data were revealed from the analysis of correlated and anti-correlated expression profiles (Figure 5A and 5B). Examples of genes showing sex-divergent expression for both male and female specific regulation were identified (Figure 5C). Patterns of expression associated with specific developmental stages or periods of adult aging were also evident and suggest age-specific regulation in the liver (Figure 5D and 5E). Array based expression profiles were verified using quantitative real time PCR for four genes, using beta actin as a housekeeping control (Figure 6). Cyp2d4, the rat ortholog of human gene CYP2D6, and Por (P450 oxidoreductase) are important genes involved in Phase I xenobiotic metabolism. The other two genes (Dbp and Arntl) are involved in circadian rhythm regulation of gene expression in the liver. An average correlation of R = 0.914 between qPCR and array data was calculated for these four genes, demonstrating high agreement between the two methods. Expression profiles for specific examples of genes found in xenobiotic metabolism pathways (Cyp2c11, Cyp2e1, Cyp3a23/3a1, Gstm1, glutathione-s-transferase mu 1; Slc22a8, solute carrier family 22 member 8) and energy metabolism (Phgdh, 3-phosphoglycerate dehydrogenase), which exhibit notable sex and age related differences, are shown in Figure 7.

The assessment of gene expression changes at key developmental periods of the life cycle, as well as at the ages most commonly used in pre-clinical drug safety testing, was the focus of this study. Fischer 344 rats are frequently used in the National Toxicology Program and one of the most commonly used inbred strains for preclinical pharmaceutical in-life testing [25, 33]. The two early time points, 2 and 5 weeks, span the time of weaning (3 weeks). The National Toxicology Program (NTP) 2-year bioassays are typically initiated in 6-8 week old rats [34], thus 6 and 8 week age groups were included to provide baseline for future NTP bioassays. The 15-week age group captures expression data after sexual and physical maturity [35]. Furthermore, a 13-week NTP sub-chronic study, initiated in 8 week old animals would be completed at 21 weeks. Reproductive senescence and aging phenotypes begin to be apparent during the second year [35], thus 52, 78, and 104 week age groups were assessed.

Data analysis began with an unsupervised, or data-driven, approach to interpreting the expression data. Clustering methods were used to describe the global relationships within the data on a per-group and per-gene basis. Group relationships were illustrated in the 3D-PCA (Figure 2) which succinctly captures two major trends in the expression data. Namely, the clear sex- and age-related differences in expression shown by the discrete spatial separation of groups visualized in these parameters. Data from females follow a generally flat linear trajectory while results from males show conspicuous divergence from females during adulthood followed by convergence at late age. This pattern of late male convergence with female data is consistent with previous reports of increasingly "feminized" liver expression patterns in older rats [18]. Such a large divergence between males and females in the composite data suggests a substantial proportion of the 3,770 individual genes influencing this clear, sex-specific trend may be related to sex steroid biosynthesis pathways. In fact, upon analyzing the component or factor loadings for the first three principal components (Figure 2), genes from each principal component exhibited a unique characteristic. Principal component 1 (PC1), when ranked by factor loadings, consisted of genes exhibiting high female expression and low male expression (female-specific top genes: Cyp2c12 (cytochrome P450, family 2, subfamily c, polypeptide 12); Akr1b8 (aldo-keto reductase family 1 member b8); Sult2a1 (sulfotransferase family 2a member 1); Prlr (prolactin receptor)). PC2-ranked genes identified male-specific expression with concurrent low female expression (male-specific top genes: Mup5 (major urinary protein 5); Cyp2c13; Pgcl1/Pgcl3/Pgcl5 (alpha-2u globulins 1, 3, and 5); Cyp2c11). Whereas PC3-ranked genes included genes that exhibited large-fold temporal differences, yet were positively correlated between sexes (sex-conserved top genes: Cyp17a1; A2m (alpha-2-macroglobulin); Col1a1 (collagen, type I, alpha 1); Phgdh (3-phosphoglycerate dehydrogenase). The top 100 genes ranked by PC1 were analyzed for gene ontology representation (ArrayTrack's Gene Ontology For Functional Analysis, GOFFA) with "regulation of hormone levels", "steroid metabolic process", and "hormone metabolic process" as the top three GO terms which would support sex steroid hormone biosynthesis pathways explaining the biggest sex differences in PC1.

Linear continuity between neighboring age groups within each sex was maintained in the 3D-PCA, (e.g., 6 week animal data spatially appear between 5 week and 8 week animals for both males and females). Additionally, 3D-PCA permits the viewing of individual animal data as they contribute to the average or composite expression for one sex and age group. The relatively tight clustering observed among animals (biological replicates) of the same sex and age group illustrates the degree of reproducibility in sex- and age-related expression profiles. In fact, correlation coefficients (R) varied from 0.962 to 0.993 within a group when using DEGs, demonstrating high biological similarity. Subsequently, a k-means clustering method was used to bin individual genes into groups according to their temporal profiles (Figure 3). Such sorting algorithms establish the presence and prevalence of various patterns, thus suggesting global mechanisms of transcriptional control that can be assessed in the context of biological necessity or functionality in the liver. These clustered temporal profiles gain more interpretive potential when correlated with gross biological trends (such as the growth curves of animal body or liver weight) or developmental windows (sexual and physical maturation). The k-means clusters can then be associated with various developmental stages such as acute growth, puberty, early adulthood, and adult aging. For example, the transition from perinatal to pubertal control of growth and differentiation is believed to be directed in part by the switch from, Igf2 (insulin-like growth factor 2) expression to predominantly Igf-1 expression in the liver [36] during early development. These gene expression changes are illustrated in the different patterns shown in cluster 2 (Figure 3) for Igf-2 versus clusters 8 and 12 for Igf-1. Igf-2 exhibits high expression at 2 weeks (100- to 500-fold change above average for female and male, respectively) followed by relatively low and constant expression at the remaining age groups. Igf-1, however, shows low expression at 2 weeks followed by peak expression between 8 and 21 weeks in both sexes (Figure 5E). Another example, Phgdh, ranked high in the principal component 3 (PC3) factor loadings which exhibited sex-conserved, age-specific expression profiles. In this case, Phgdh fell into clusters 22 and 29 (male and female, respectively), demonstrating transient up-regulation between 5 and 6 weeks with declining expression in subsequent age groups. This expression pattern is consistent with the known association of Phgdh with highly replicating cells as the first and rate-limiting step in de novo biosynthesis of serine [37], which also plays a role as nucleotide precursor in regenerating liver tissue and hepatomas [38]. Thus Phgdh's transient up-regulation during early liver growth phase (5-6 weeks) is consistent with this physiological role. These patterns illustrate the temporal relationships between their expression and roles in liver development.

Several k-means clusters exhibit patterns of up-regulation with age. Clusters 6, 14 and 23 (Figure 3) in particular show an upward trend between 52 and 104 weeks, suggesting a group of genes important in the adult aging process. Examination of these genes using Ingenuity Pathway Analysis suggests many may be involved in compensatory mechanisms adapting to higher levels of liver injury, DNA damage or loss of function. For example, Mgmt (O6-methylguanine DNA methyltransferase) is a DNA damage repair enzyme which acts to directly reverse alkylation adducts in DNA [39]. Mgmt was found in clusters (clusters 11 and 23, male and female, respectively) which exhibit a trend toward of increasing expression in late adulthood (between 21 weeks and 104 weeks) (Figure 5D) which is consistent with previous reports associating Mgmt with aging phenotype [40]. Likewise, Por (P450 oxidoreductase) is a member of the "up-late" clusters. Por's microarray expression patterns were further verified by Taqman quantitative real time PCR (qRTPCR) (Figure 6). Interestingly, Por, along with cytochrome b5 [41], is a key electron donor for microsomal cytochrome P450 monooxygenases (CYP enzymes) [42], many of which function as Phase 1 metabolizing enzymes during xenobiotic metabolism. Por essentially acts to reload CYP enzymes in their monooxygenase activities. Por's reductase activities also target a number of other oxygenase enzymes involved in sterol and cholesterol biosynthesis, heme degradation to bilirubin, and fatty acid metabolism [42]. Thus, Por's consistent up-regulation from 52 weeks to 104 weeks could be explained by increased demand for CYP enzyme reloading due to accumulated xenobiotic exposure or fatty acid accumulation as a function of the adult aging process.

Other genes associated with "up-late" clusters (k-means clusters 1, 6, 9, 14, 23) included several oxidative stress related genes including Ephx1 (epoxide hydrolase 1), Fmo1, Fmo5 (flavin containing monooxygenase 1 and 5) and Hmox1 (heme oxygenase 1). Hmox1 is an important antioxidant protein expressed in response to oxidative stress [43], and thus its continuous increase in expression in older animals suggests a compensatory mechanism specifically to address increased levels of reactive oxygen species in the liver.

Expression level is not always an indicator of protein activity. For example, a previous report showed changes in male rat hepatic Cyp2e1 activity between 3 and 18 months of age despite the relative stasis of mRNA and protein levels [44]. However, examination of Cyp2e1 gene expression on a finer scale, as was done in this study, revealed sex-conserved spikes in young animals (between 2 weeks and 8 weeks) followed by greater than 2-fold differences in expression level between sexes at 21 weeks (Figure 7). Therefore, care must be taken in interpreting changes in mRNA expression levels on their own. This broader examination using more time points during the life cycle may prevent over-interpretation of expression data captured in studies of smaller scope.

Male rodents are preferred in many studies evaluating toxicological or pharmacologic response in vivo due to adult female variability in response during estrous cycle [45]. However, the current study using unsynchronized females is a crucial starting point to evaluating sex differences in liver expression over the larger scale of the life cycle. The number of replicates per age group (5 per age group) and relative consistency across multiple adult age groups increases our confidence in the sex-specific differences. Whereas Mori et al. [17] were limited to comparisons of relatively young (32 weeks) and old males (84 weeks), a full life-cycle analysis in both sexes provides comprehensive data on the temporal dynamics of gene expression and evaluates how well correlated male and female data are throughout the life cycle. If and when male and female data do not comport with each other, these sex-differences may inform our assessments of clinical (human) response which may influence the sex and age of rodent models used for pre-clinical testing. Wauthier et al [22] provide comprehensive lists of male- and female-specific liver gene expression in approximately 13 week old rats. A comparison of their top 100 male- and female-specific genes with the top 100 male/female-specific genes from the 15 week age group of the current data set was performed. An overlap of 73% (46/63, matched by Agilent probe ID) and 61% (43/71) for male- and female-specific expressed genes, respectively, were identified, suggesting a good degree of overlap of the most sex-specific genes. Furthermore, Lee et al. [14] complemented their male data with a meta-analysis of sex-specific genes from a separate study [46] and showed age-related changes in xenobiotic metabolism genes in the course of 6, 11, 18 and 24 month old rats. However, because the male and female data were not from the same study and the time and sex analyses were evaluated separately, inter-study confounding factors can not be adequately controlled for. Thus, analysis of sex- and age-related changes between animals within the same in-life study, as reported here, removes any uncertainties with regard to such confounding factors and allows for more robust comparisons.

Efforts have been made using differing metrics (comparing times of weaning, sexual maturity, skeletal maturity, etc.) to align the life cycle of normal human development and maturation with that of the rat [35], illustrating the compressed timescale that exists in rodent models. Properly extrapolating information about specific windows of development or aging from rodents to humans is therefore only approximate; however, the comprehensive gene expression profiles throughout the life cycle of the rat reported here may help further refine the life-cycle alignment of the rat model to human.

The FDA has assembled a listing of valid genomic biomarkers in the context of approved drug labels that are used in a variety of specific uses, including evaluating clinical response and differentiation; risk identification; dose selection guidance; susceptibility, resistance and differential disease diagnosis; and polymorphic drug targets [47]. Of the 20 gene-specific biomarkers on the FDA list, 13 of the genes were present on the Agilent arrays and 8 of them were differentially expressed by age or sex. Of these 8 genes, CYP2C9 (r Cyp2c11), CYP2D6 (r Cyp2d4), and EGFR (epidermal growth factor receptor), exhibited notable sex-differences in gene expression patterns in this study. Cyp2c11 showed much greater (200-250 fold difference) male-specific expression from 8 to 52 weeks of age. Both Cyp2d4 and Egfr showed ~1.5 and 2-fold higher expression, respectively, in males from 15 to 52 weeks. The relative expression of cytochrome P450 enzymes may have sex- and age-specific implications for the liver; however, Egfr was validated as a biomarker for cancer in non-hepatic tissues. Two other of the 9 DEGs, Dpyd (dihydropyrimidine dehydrogenase) and Nat1 (N-acetyltransferase 1), exhibited sex-conserved repression (~-2-fold) at early ages (2 and 5 weeks) followed by minimal changes at later ages. The remainder of the validated biomarkers that were differentially expressed during the rat life cycle showed few or minimal changes in expression that were very close to the screening criteria threshold. Thus, the expression levels of these biomarkers vary through the life-cycle of the rat model and this may have implications for their appropriate use.

Circadian rhythms have been shown to influence the expression of xenobiotic metabolism genes during the 24-hr cycle [48]. There are also reports that changes in circadian control during the aging process may, in addition, impact xenobiotic metabolism capacity during the life cycle [49]. Master circadian regulators Clock and Arntl/Bmal1 (aryl hydrocarbon receptor nuclear translocator-like) control the expression of the PARbZip family of transcription factors which is comprised of three members: Dbp (albumin site D-binding protein), Tef (thyrotroph embryonic factor), and Hlf (hepatocyte leukemia factor). In turn, these PARbZip transcription factors have been shown to regulate several enzymes and regulators involved in detoxification and drug metabolism, such as cytochrome P450 enzymes, carboxylesterases and constitutive androstane receptor [50]. The genes for two of these transcription factors, Dbp and Tef, were found to be differentially expressed throughout the life cycle in an age- and sex-specific manner (Figure 5A) and also to be highly correlated between the sexes (Pearson, R = 0.94). Furthermore, humans exhibit sexually dimorphic patterns of suprachiasmatic nucleus (SCN) neuron deterioration, with women showing no change throughout adulthood while men exhibit complex and fluctuating patterns with increasing age [51]. It is interesting to note that the pattern of SCN-controlled clock genes, Dbp, Tef, and Hlf (Hlf's 1.3-fold change did not meet criteria for inclusion in list of DEGs) display a conspicuous correlation with these reports (Figure 5A). Namely, while female expression remains consistent between 52 weeks and 104 weeks, males exhibit comparable levels of expression to female at 52 and 104 weeks but a noticeable decrease at 78 weeks. The microarray expression profiles of Dbp and Arntl/Bmal1 were further verified by qRTPCR (Figure 6). Based upon the known common pathway (Dbp and Tef controlled by the SCN) and life cycle expression data exhibiting correlated and sex-specific patterns, these data suggest a shared regulatory mechanism between the PARzip family members' expression and the sex-specific pattern of SCN neuronal deterioration. This further supports the concept that the circadian clock may influence life cycle "clocks" and thus influence downstream xenobiotic metabolism capacity.

The major mechanisms controlling sex differences in liver gene expression have been well-characterized [52] and are dictated by the distinct patterns of continuous vs. pulsatile plasma growth hormone (GH) secretion in females and males, respectively [53]. These temporal differences in plasma GH levels influence the activity of transcription factors Stat5b (signal transducer and activator of transcription 5b), Hnf4α (hepatic nuclear factor 4, alpha), and Stat5a (signal transducer and activator of transcription 5a), which then differentially activate sex-specific cytochrome P450s (CYPs) such as Cyp2c family members. Several Cyp2c family members analyzed in this study exhibited clear and conspicuous differences in expression, including Cyp2c7, Cyp2c11, Cyp2c13, Cyp2c12, Cyp2c79, and Cyp2c22. Previous studies of rat liver gene expression have also shown age-related changes in xenobiotic metabolism genes, although only male rats were examined, [17] and included Cyp2c11 and Cyp3a2. Both of these genes were represented on the Agilent array used in this study and showed comparable changes in gene expression at similar age groups along with related family members (Cyp2c13 and Cyp3a23/3a1, Cyp3a9 and Cyp3a18). A number of xenobiotic metabolism related receptors of interest were present on the array but were either not differentially expressed (pregnane X receptor, Nr1i2; constitutive androstane receptor, Nr1i3; aryl hydrocarbon receptor, Ahr; peroxisome proliferator-activated receptor alpha and gamma, Ppara/Pparg; retinoid x receptor alpha and beta, Rxra/Rxrb) or minimally passed threshold criteria for differential expression (retinoid X receptor, Rxrg). Cyp1a2, a notable xenobiotic and drug metabolizer and known target of Ahr, exhibited roughly parallel expression between sexes with maximum relative expression (~1.5 fold) peaking around 6 to 8 weeks followed by declining expression with age. Of the 34 phase I and II metabolizing genes found to be changed with age in male animals by Mori et al. [17], 27 were present in this study's list of DEGs and showed comparable expression changes across corresponding time scales. However, it is important to note that the female patterns differed remarkably from males in over 50% of those genes, raising uncertainties regarding how well data from male animals reflect the female response. Data for Cyp2c11 and Cyp3a23/3a1 are shown in Figure 7.

It has also been shown that GH regulates Cux2 (cut-like homeobox 2) which has been suggested [54] to act as a suppressor of male-specific gene expression. Laz et al. [54] employed in silico methods to predict binding sites of Cux2 suggesting target proximal genes for repression. In the current data set, when life cycle profiles for Cux2 were compared to the 16 suggested target genes, a clear sex-specific repressive relationship was observed with one of the genes, Cyp4a8 (R = -0.699 in males; Figure 5B). Cux2 showed approximately 12-fold expression ratio difference between males and females from puberty through early adulthood, while Cyp4a8 showed approximately 11-fold sex difference (male/female) in expression level. These comparative life cycle expression data provide further evidence of Cux2's role in sex-specific repression of Cyp4a8. Cyp4a8 is also a target of PPARα (peroxisome proliferator activated receptor alpha) signaling [55] and its expression has been suggested to explain sex and strain differences in the susceptibility to hypertension and target organ damage [56]. Thus these data suggest integrated sex-specific control of genes like Cyp4a8 that are also involved in xenobiotic receptor mediated responses.

To discover further large-scale patterns of sex-specific life cycle gene expression, the data were analyzed on a per-age group basis assessing the number of instances where the expression difference of a given gene between the sexes was at least 2-fold, with the condition that the data at that individual age group first met the screening criteria of p < 0.05 and great than or equal to 1.5 relative fold-change (Figure 4). The periods of the life cycle in which one sex predominates over the other in expression level appear to be non-random and appear at discrete times of the life cycle. For example, female predominant expression prevails during or soon after sexual maturation and into early adulthood (8 and 15 weeks) whereas male predominant expression appears to peak between 21 and 78 weeks of age. Interpreting this result in light of the 3D-PCA analysis, the data suggest that differences between males and females are observed immediately at 2 weeks, with more noticeable sex-divergence at approximately 6 weeks (Figure 2) which extended through sexual maturation. This sex-divergence seems to be due at first to genes being expressed at higher levels in females than males (Figure 4), followed by a transition where male expression is higher than females throughout the majority of adulthood. For example, Stat3 exhibits transient, female-specific predominance at 15 weeks while Cdh17 (cadherin 17) shows clear male-specific expression between 8 and 104 weeks (Figure 5C). Only one clear example, Cyp3a23/3a1 (Figure 7), was found among the 122 xenobiotic metabolism genes that exhibited both the female predominant expression at 15 weeks and male predominant expression at 52 weeks, suggesting that the sex-specific predominance of expression observed in Figure 4 primarily occurs in different sets of genes. The large-scale patterns observed in the global analyses (Figures 2 and 4) are also observed in individual gene profiles (Figure 7), including several cytochrome P450 family members in addition to other xenobiotic metabolizing proteins (Table 1).

This study comprises one of the most comprehensive data sets to date for rodent liver gene expression with regards to sex, life span, biological replication, and genome coverage. A broad description of the dataset in the context of the normal physiology of the rat life cycle suggests important implications for toxicological susceptibility as related to xenobiotic metabolism, DNA damage repair and oxidative stress response. Namely, principal component analysis identified and ranked genes which exhibited the greatest sex-specific and sex-conserved differential expression. Furthermore, specific periods of the life cycle were identified that display a higher number of genes with sex-specific expression which included putative toxicological susceptibility related genes (e.g., Cyp2d4, Cyp3a23/3a1, Gstm1, Slc22a8). Cluster analyses, in combination with functional annotation analysis, identified gene clusters that positively correlated with adult aging which also included several susceptibility related genes (i.e., Por, Mgmt). Lastly, life cycle expression data provided experimental support for gene regulatory relationships previously based upon in silico predictions alone (i.e., Cux2, Cyp4a8).

The ability to make informed regulatory decisions about the impact of drugs and other chemicals on the liver is dependent upon a growing understanding of liver biology and more specifically the expression profiles of the molecular components which play a role in the pharmacokinetics and pharmacodynamics in liver cells. This large-scale study captures comprehensive baseline hepatic gene expression profiles during the life cycle of male and female rats which should prove to be an increasingly valuable resource for toxicogenomics-related research and ongoing evaluations of pre-clinical and clinical biomarkers of toxicity and efficacy.