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
]. 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 ) 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 ), 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 ). 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 ) 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 ). 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 ) 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 ) 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 ). 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 ). 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 ) 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 ). 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 ). 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 .
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 ). 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 ). 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 ) 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 ), 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 ). Only one clear example, Cyp3a23/3a1 (Figure ), 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 primarily occurs in different sets of genes. The large-scale patterns observed in the global analyses (Figures and ) are also observed in individual gene profiles (Figure ), including several cytochrome P450 family members in addition to other xenobiotic metabolizing proteins (Table ).