- Research article
- Open Access
Stochastic variation of transcript abundance in C57BL/6J mice
BMC Genomics volume 12, Article number: 167 (2011)
Transcripts can exhibit significant variation in tissue samples from inbred laboratory mice. We have designed and carried out a microarray experiment to examine transcript variation across samples from adipose, heart, kidney, and liver tissues of C57BL/6J mice and to partition variation into within-mouse and between-mouse components. Within-mouse variance captures variation due to heterogeneity of gene expression within tissues, RNA-extraction, and array processing. Between-mouse variance reflects differences in transcript abundance between genetically identical mice.
The nature and extent of transcript variation differs across tissues. Adipose has the largest total variance and the largest within-mouse variance. Liver has the smallest total variance, but it has the most between-mouse variance. Genes with high variability can be classified into groups with correlated patterns of expression that are enriched for specific biological functions. Variation between mice is associated with circadian rhythm, growth hormone signaling, immune response, androgen regulation, lipid metabolism, and the extracellular matrix. Genes showing correlated patterns of within-mouse variation are also associated with biological functions that largely reflect heterogeneity of cell types within tissues.
Genetically identical mice can experience different individual outcomes for medically important traits. Variation in gene expression observed between genetically identical mice can identify functional classes of genes that are likely to vary in the absence of experimental perturbations, can inform experimental design decisions, and provides a baseline for the interpretation of gene expression data in interventional studies. The extent of transcript variation among genetically identical mice underscores the importance of stochastic and micro-environmental factors and their phenotypic consequences.
Variation in transcript abundance between individuals has important implications for microarray experimental design and significance testing . Ideally, microarray experiments are designed with samples from multiple individuals in each treatment group. This biological replication provides the variance estimator that is required to establish the statistical significance of between-group differences. In this study, we collected multiple samples of tissues within each of several genetically identical mice. Multiple sampling within individuals is not necessary in an experiment aimed at making between-group comparisons, but it is essential if the aim is to identify significant variation between individuals within the same experimental treatment group. An important procedural detail in this type of study is to determine how to collect and at what stage to divide the tissues to create multiple samples. In this study, we elected to split tissues immediately after dissection and before RNA extraction in order to restrict the possible sources of between-mouse variation to events that occur prior to dissection. With this experimental design, transcript variation can be decomposed into within-mouse and between-mouse variance components. Between-mouse variance reflects differences in whole-tissue transcript abundance between genetically identical mice. Within-mouse variance captures variation due to RNA extraction, array processing, and heterogeneity of gene expression within tissues, which may be amplified by dissection and tissue collection procedures.
Individual variation in gene expression can have important phenotypic consequences. However, only a few studies have previously attempted to characterize gene expression variation in genetically identical mice. Koza et al. (2006)  described gene expression signatures in adipose tissue that are predictive of future adiposity among genetically identical C57BL/6J mice. The use of multiple biopsy samples in this time-course study was essential to establish the link between gene expression variation and late-life adiposity. However, biopsy sampling may be subject to unexpected variation introduced by tissue heterogeneity, as we illustrate below.
Two previous studies have used multiple sampling within individuals to provide a statistical basis for detecting transcript variation between genetically identical mice. Pritchard et al. (2001)  examined 3 tissues in each of 6 C57BL/6J mice and reported that immune function, stress response, and hormone regulation were important sources of biological variation. Pritchard et al. (2006)  examined liver tissue in 3 animals from each of 5 inbred mouse strains and found that genes differentially expressed within strains were enriched for cell growth, cytokine activity, amine metabolism, and ubiquitination. In these experiments, technical replicates were obtained by splitting samples after RNA extraction. This approach confounds variation due to dissection and RNA preparation with variation between mice.
We designed and carried out an experiment to study transcript abundance variation in four tissues among young adult male C57BL/6J mice (Figure 1). Our sampling design enabled us to partition the variance for each gene into within-mouse and between-mouse components, with a division line that corresponds to the step of splitting tissues. We examined within-mouse and between-mouse variation in more than 22,000 protein coding genes and identified groups of genes with shared patterns of variation that are enriched for known biological functions. To facilitate exploration of our data, we have created an on-line resource that includes graphical displays, test statistics, and gene groupings for all transcripts characterized in this study http://cgd.jax.org/individualvariation.shtml.
We performed a microarray experiment using the Illumina Sentrix® Mouse-6 v1.1 BeadChip microarray platform to study transcript variation in 10-week old male C57BL/6J mice (Figure 1). Six pairs of siblings were co-housed from weaning under uniform environmental conditions. From each mouse we obtained duplicate samples of adipose (inguinal fat pad), heart, kidney, and liver tissues by splitting whole organs or tissues prior to homogenization and RNA extraction. Adipose, heart, and liver tissues were coarsely cut into pieces and divided into two samples that were homogenized separately in order to extract RNA. The left and right kidneys were also homogenized separately. We computed a decomposition of variance for each probe on the array (Methods). The within-mouse variance component captures biological variance between two dissected tissue samples as well as technical variance due to sample and microarray processing. The between-mouse variance component reflects differences between individual mice. We repeated gene expression assays on the liver samples, using the Affymetrix Whole-Transcript Mouse Gene 1.0 ST array, to provide validation on a different measurement platform.
Expressed genes and variable genes
We declared a gene to be expressed if the probe intensity was greater than the 95th percentile of the negative control probes for both samples in at least 1 of the 12 mice. A total of 12657 genes, representing 55% of the annotated probes on the array, were expressed in at least one of the four tissues. Across tissues, the number of expressed genes ranged from 8919 (39%) in liver to 11204 (49%) in adipose tissue (Table 1A).
We computed the total variance, s2 , across all samples for each gene in each tissue (Figure 2A). Liver and kidney have relatively few genes of high variability but heart and adipose have many. We tested the hypothesis that the distribution of total variance occurred by chance using a χ2 test (Methods) and found significantly greater variance than expected in each tissue (Table 1B). We applied coexpression network analysis to the top 2500 genes in each tissue, which we refer to as the variable genes.
We decomposed total variance for each gene into within-mouse (s w 2 ) and between-mouse (s b 2 ) components. The distribution of between-mouse variance components was similar across all four tissues (Figure 2B). Adipose tissue showed the greatest number of genes with a large within-mouse component followed by heart, kidney, and liver (Figure 2C). The variable genes include the 313 (adipose), 189 (heart), 405 (kidney), and 990 (liver) genes with the largest between-mouse variance. They also include the 1526 (adipose), 1347 (heart), 593 (kidney), and 221 (liver) genes with the largest within-mouse variance.
Significance of between-mouse variance
Within each tissue, for each gene, we computed a test statistic to assess the significance of the between-mouse variance component relative to the within-mouse variance component. We applied a family-wise error rate correction  (as in Pritchard et al. (2001) ) and found few genes with significant between-mouse variation (Table 1C). We applied a false discovery rate (FDR) adjustment  (as in Pritchard et al. (2006) ) and found no differentially expressed genes in adipose or heart and only modest numbers in kidney (2%) and liver (11%) (Table 1C). We estimated the proportions of differentially expressed genes (1 - π 0 ) using the q-value software  and found similar results (Table 1C; [Additional files 1, 2: Supplemental Figure S1]).
A different picture of the variability in gene expression across tissues emerges when we look at the maximal fold change between mice (Table 1D). In adipose, 2.6% of all genes exhibit maximal fold changes greater than 2, whereas 0.4-0.6% of all genes show fold changes this large in the other three tissues. Although the fold-change is not a statistical criterion, the differences across tissues are dramatic. There are many genes with large maximal fold changes between mice but, particularly in adipose tissue, these same genes also have large within-mouse variance, which reduces their statistical significance.
Variable genes form clusters that are enriched for specific biological functions
We used co-expression network analysis [8, 9] to cluster the variable genes into modules with correlated patterns of expression (Methods) (Figure 3). Module sizes ranged from 34 to 1340 genes with an average module size of 215 genes (Table 2). We identified 8 to 9 modules in each tissue comprising 97% (adipose), 80% (heart), 61% (kidney), and 54% (liver) of the variable genes. For each module, we applied principal components analysis to compute a module eigengene  that represents the dominant pattern of variation (Figure 4). The percentage of variation explained by the module eigengene ranges from 47% to 88%, indicating that the eigengenes are representative of expression profiles of the individual genes in each module. In the following, we refer to modules using a colour code within each tissue (Table 2).
For each gene, we computed the intraclass correlation coefficient, c ≡ s b 2/(s w 2 + s b 2 ), which is the proportion of total variance attributable to the between-mouse component. Median values by module ranged from c = 0.00 (8 modules) to c = 0.79 (liver-pink) (Table 2). Kidney and liver, respectively, have 5 and 8 modules with high intraclass correlation (c ≥ 0.35) indicating substantial between-mouse variance while adipose has two and heart has no modules with high intraclass correlation (c ≥ 0.35). Each tissue also has at least one module with low intraclass correlation (c ≤ 0.02) indicating that differences between samples within mice are greater than differences between mice.
For each module, we computed enrichment scores  for the GO biological process, cellular component, and molecular function terms and for KEGG pathways. The highest scoring enrichment category for each module is listed in Table 2. Each module can be divided into two subsets such that all correlations within a subset are positive (Methods). We also tested for enrichment within each of these subsets [Additional file 3: Supplemental Table S1]. Many of the module enrichment scores are highly significant indicating that correlated groups of variable genes are enriched for specific biological functions.
Most modules in a given tissue share similar features with at least one module in another tissue [Additional files 2, 4: Supplemental Figure S2; Additional file 5: Supplemental Table S2]. Several sets of modules shared similar patterns of between-mouse variation and had significant gene overlap and functional enrichment. Other sets of modules shared similar patterns of within-mouse variation, but with distinct between-mouse variation. Several pairs of modules had significant gene overlap but did not have correlated patterns of variation. Examples of each are described below.
Between-mouse patterns of variation are shared across tissues
Modules from different tissues that are enriched for similar functional categories typically have high intraclass correlation and similar patterns of between-mouse variation. To quantify this similarity, we computed a between-mouse correlation, r b , for all pairs of module eigengenes by averaging the two within-mouse samples before computing the Pearson correlation.
Each of the four tissues has at least one module that is enriched for immune response. The heart-brown (c = 0.03), kidney-gold (c = 0.57), and liver-pink (c = 0.79) modules are enriched for the GO category exogenous antigen presentation (Table 2). The between-mouse correlations, r b , range from 0.53 to 0.80, and the genes in these modules overlap significantly based on a hypergeometric test (Methods). Pairwise overlaps range from 16 to 19 genes and seven genes (Cd274, Cd74, Cxcl9, H2-DMa, H2-Eb1, Igtp and Iigp2) are found in all 3 modules.
The kidney-blue (c = 0.42) and liver-brown (c = 0.74) modules are enriched for GO category extracellular matrix, each containing more than 12 genes of that category. Their between-mouse profiles are correlated (r b = 0.75) and they share 20 genes in common (out of 36) including Adamts2, Col5a1, Col6a1, Col14a1, Ecm1, Igfbp3, Tgfbi and Timp2.
The adipose-red (c = 0.42), heart-blue (c = 0.20), kidney-brown (c = 0.59) and liver-black (c = 0.78) modules are enriched for the GO category apoptosis and have between-mouse correlations, r b , ranging from 0.52 to 0.93. These modules overlap with 16 genes present in at least 3 of the 4 modules including Ccrn4l, Gadd45g, and Map3k6. The liver-blue module (c = 0.77) also has a high between-mouse correlation (r b ≥ 0.64) and significant gene overlap with these adipose, heart and kidney modules including Fkbp5 and Per1.
The kidney-pink (c = 0.69) and liver-magenta (c = 0.74) modules have correlated between-mouse profiles (r b = 0.88), and each contains 18 or more genes of the GO category DNA-dependent regulation of transcription. Their gene overlap (12 out of 47) includes Bcl6, Cish, Rgs3, and Socs2.
The between-mouse profiles of the kidney-green (c = 0.59) and liver-red (c = 0.72) modules are correlated (r b = - 0.73) and each module contains 12 or more genes of the GO category cellular lipid metabolic process. They have 12 genes in common (out of 60) including Acaa2, Acadm, Agxt2l1, Cyp26b1, Cyp4a10, Cyp4a14 and Slc2a2.
Within-mouse patterns are similar across modules of the same tissue
Some modules had similar patterns of within-mouse variation but different patterns of between-mouse variation. To measure similarity of within-mouse variation, we centred the sample values on individual mouse means and then computed a Pearson correlation, r w . This measure is only meaningful for comparisons within the same tissue as there is no correspondence between the duplicate samples from different tissues. Adipose and heart each have multiple highly correlated modules (|r w | ≥ 0.64). The adipose-green, adipose-red, adipose-black, and adipose-magenta modules have distinct patterns of between-mouse variation and different functional enrichment, but they all share high within-mouse correlation (Figure 4A, Table 2). A similar relationship was observed for the heart-green, heart-red, heart-turquoise, heart-blue, heart-brown, and heart-gold modules (Figure 4B, Table 2).
Uncorrelated modules have gene overlap and similar functional enrichment
Some modules share genes and functional enrichment categories but do not have correlated patterns of variation. The adipose-gold (c = 0.12), heart-red (c = 0.00), and kidney-black (c = 0.10) modules have a high gene overlap (adipose-gold & heart-red, 48 out of 72; adipose-gold & kidney-black, 10 out of 29; heart-red & kidney-black, 25 out of 40 and 9 genes in all three including Acaca, Cidea, Cox8b and Ucp1). They are enriched for the GO category fatty acid metabolic process. The adipose-magenta (c = 0.00) and heart-gold modules (c = 0.00) share 118 out of 120 genes including Cd8b1 and Lck and are enriched for the GO category immune system process. The adipose-brown module (c = 0.07) shares 87 out of 182 genes with the heart-green module (c = 0.00) and 31 out of 35 genes with the liver-turquoise module (c = 0.02). These modules are enriched for the GO actin cytoskeleton category and share 8 genes in common including Ckm and Myl1. The adipose-turquoise (c = 0.44) and kidney-turquoise (c = 0.01) modules share 30 out of 40 genes including Apoal, Cyp8b1, and Ugt2b3 and are enriched for the KEGG pathway complementation and coagulation cascades. The adipose-green (c = 0.28) and heart-turquoise (c = 0.17) modules are overlapping in 12 out of 50 genes including Ccl9, Cxcl1, Egr1, Fos, and Hmox1 and are enriched for chemokine activity. The adipose-black (c = 0.00), heart-black (c = 0.00), kidney-magenta (c = 0.00), and liver-green (c = 0.39) modules have pairwise overlaps ranging from 33 to 146 genes. Twenty-two genes are shared among all 4 of these modules and they are enriched for KEGG pathway oxidative phosphorylation and the GO category mitochondrial inner membrane.
Comparison across platform
We repeated the gene expression assays for only the liver samples on a different platform, the Affymetrix Whole-Transcript Mouse Gene 1.0 ST array. To facilitate comparison, we generated a cross-platform probe map based on gene annotation (Methods). Using this map, we computed eigengenes of the previously defined clusters from the Affymetrix data. Correlation of the eigengenes across platforms was very high for 7 of the 9 modules (r > 0.89 for 6 of 9 modules, r = 0.76 for liver-brown; [Additional files 2, 6: Supplemental Figure S3]). Two modules with lower correlation (liver-gold: r = 0.42, liver-green: r = 0.21) had less than 20% of variance explained by the eigengene with Affymetrix data. However, for the liver-gold module, low expression for mouse 6 was a consistent pattern across platforms. The profiles of all 19 genes that are highlighted in the Discussion (below) are highly correlated across platform (r > 0.55 for all 19, r > 0.70 for 16 of the 19, [Additional files 2, 7, 8, 9: Supplemental Figures S4-S6]).
There are several mechanisms that may contribute to between-mouse variation in gene expression in C57BL/6J mice. New mutations that create single nucleotide or copy number variants may result in variable gene expression. We expect such events to be rare. However, we have observed a striking pattern of differential expression (r b > 0.88; p < 0.01) in the insulin degrading enzyme (Ide) with approximately two-fold higher expression in all 4 tissues for the two mice of cage 4. We speculate that these siblings may have inherited a copy number variant at this locus on chromosome 19 for which copy number changes have been observed previously in C57BL/6J mice . Genes that display circadian or other periodic expression patterns can be out of phase in different animals. We attempted to control for cyclical variation by collecting samples in a consistent and narrow time frame for all mice. Variation in feeding behaviour is another possible factor and although we implemented a 4-hour fast prior to tissue collection, some variation in time since last feeding is inevitable. Epigenetic differences may affect the expression of genes as a result of variable access to nutrients in utero, birth order, maternal stress or other pre- or post-partum events. Slight differences in phenotype at birth may be magnified over time. Response to subtle differences in local environment may have an effect on gene expression and finally, the expression of some genes may be sensitive to events just prior to euthanasia .
Within-mouse transcript variation could reflect stochastic variation in gene expression, which has been observed within individual cells and across cell populations [13–20]. However, if it is present, this effect seems to be dominated by other factors in our study. Tissue heterogeneity due, for example, to localization of stem and progenitor cell populations can result in sampling variation [21–24]. This variation may be amplified by dissection, especially in tissues with imprecise boundaries. Even a relatively homogenous and easily isolated tissue such as liver will have internal structure that can influence local gene expression [25, 26].
Phenotypic implications of between- and within-mouse variation in adipose tissue
Adipose tissue is compartmentalized into adipocytes, preadipocytes, and vascular epithelium . The degree of vascularisation can vary significantly across different regions of the same fat pad and is expected to be greater in the portion of the inguinal fat pad that is near the inguinal lymph node . Vascularised adipose tissue tends to be more metabolically active . We found a large number of genes that have within-mouse variation related to vascularisation in the adipose-magenta module (1340 genes, c = 0.00). The positively correlated subset of this module is enriched for GO biological processes immune response, T cell activation, and lymphocyte activation [Additional file 3: Supplemental Table S1] and include genes expressed in lymphocytes such as Lck, Cd8b1, and Elf1 (Figure 5, [Additional file 10: Supplemental Table S3]). Some genes within the adipose-magenta module, which is dominated by within-mouse variation, also have large between-mouse fold changes. These genes, including Bmp3, Sfrp5, Mest, Lep and Trp53inp2, are positively correlated with body weight and were previously found to be predictive for adiposity  [Additional files 2, 11: Supplemental Figure S7]. They are also negatively correlated with the module eigengene, which is consistent with higher expression in the less vascularised region of the inguinal fat pad, suggesting an inverse relationship between vascularisation and adiposity.
We chose to study the inguinal fat pad because it can be efficiently dissected. Gene expression can vary among fat depots [29, 30] and proximity to the inguinal lymph node clearly contributed to heterogeneity in the inguinal fat pad. This limits our ability to generalize our findings. However, our previous experience  indicates that other fat depots are at least as variable as the inguinal depot. The Koza et al. study  identified their adiposity signature, which we have replicated, in epididymal and retroperitoneal fat.
Variable brown fat signature in white fat tissue
Several genes in the adipose-gold module are expressed exclusively in brown fat, including Ucp1, Cidea, and Cox8b [Additional files 2, 12: Supplemental Figure S8A-B]. This module is enriched for fatty acid metabolism and the module eigengene is correlated with Prdm16 (r b = 0.86; r w = 0.74; [Additional files 2, 12: Supplemental Figure S8C]), which is part of a transcriptional complex that promotes brown fat differentiation and suppresses skeletal muscle cell differentiation [32, 33]. The adipose-brown module is enriched with 21 genes of the GO biological process muscle contraction. Genes in this module are expressed in both skeletal muscle and brown fat and many are related to brown fat cell differentiation [32, 33]. We ruled out cross contamination with muscle tissue by inspection of the dissection procedure. The enrichment for muscle contraction appears to be spurious and reflects a potential pitfall of enrichment analysis using GO annotation.
Most of the variation in the adipose-gold (c = 0.12) and adipose-brown (c = 0.07) modules is attributable to the within-mouse component, which suggests a heterogeneous spatial distribution of brown fat within the inguinal fat pad. However, large between-mouse fold changes, including Ckm, with 56-fold change, the largest observed in this study [Additional files 2, 12: Supplemental Figure S8D], suggest that the proportion of brown fat may also vary across mice. Brown fat tissue proportion have previously been shown to vary with age, strain, and environmental conditions .
Region-specific variation of gene expression in heart
The heart is composed primarily of cardiac smooth muscle, but it is differentiated into atrial, ventricular and trabecular regions with a left-right asymmetry. Several genes expressed in atria and trabeculae of the heart are repressed in the ventricles, in part, through activity of the transcription factor, Gata4. The heart-green module (898 genes) is enriched for these genes and shows a pattern of within-mouse variation with little between-mouse variation (c = 0.00). Gata4 is in the heart-red module (c = 0.00), which has a strong within-mouse correlation to the heart-green module (r w = 0.89) [Additional files 2, 13: Supplemental Figure S9]. Gata4 is negatively correlated with the heart-red eigengene such that the within-mouse variation in Gata4 is inversely related to the expression of ventricle-repressed transcripts [Additional files 2, 13: Supplemental Figure S9]. We compared our results with a study of chamber-specific gene expression (Tabibiazar et al. (2003) ) and found that, of the 27 genes previously reported to be more highly expressed in the atria than in the ventricles, 26 are included in the heart-green module. The relatively small magnitude of between-mouse variation in these modules reflects the effect of averaging of the two samples, which together comprise the whole heart. We conclude that much of the within-mouse variation observed for heart tissue is a consequence of variable proportions of anatomical substructures, specifically ventricular tissue, within the samples.
Androgen-regulated genes are variable between mice in the kidney
Many genes are regulated in response to androgens. In mice, Srd5a2 plays a key role in androgen signal amplification  suggesting that androgenic effects in individuals with higher Srd5a2 expression could be more pronounced. Hsd11b1 facilitates the conversion of testosterone to adrenosterone  and has been shown to be androgen-responsive in mice . These genes were found to be variable between mice and cluster together in the kidney-green module (c = 0.59), which is enriched for the KEGG androgen and estrogen metabolism pathway. Other androgen-responsive genes in the kidney-green module include Cyp4a14, Slco1a1, Nudt19, Prlr, Angptl7, Hsd17b11, and Tmco3 [Additional files 2, 14: Supplemental Figure S10].
It is not immediately clear if this variation reflects transient or steady state variation in androgen levels between mice. The expression of a mouse urinary protein, Gusb, is responsive to androgens in the long-term but not in the short-term . Gusb has significant between-mouse variation that is correlated with the kidney-green module eigengene (r b = 0.73) (Figure 4, [Additional files 2, 14: Supplemental Figure S10]). This suggests that other genes in this module also reflect steady state androgen levels, which may have important physiological and behavioural implications.
Between-mouse variation in fatty acid metabolism in the liver
Genes in the liver-red module have either low or high expression in the two mice of cage 3 [Additional files 2, 7: Supplemental Figure S4]. Genes in the low expression subset are involved in oxidation of fatty acids (Acaa2, Acadm, Ces3, Crat, Cyp4a10, Cyp4a14, and Elovl3). Genes in the high expression subset, specifically Tnfrsf1a and Ptgis, are involved in the conversion of the essential fatty acid arachidonic acid to prostaglandins. Thus, we see decreased fatty acid degradation in mice that are actively utilizing fatty acids. The liver-red module also shares genes with the androgen-associated kidney-green module which may reflect the requirement for lipids as precursors in androgen synthesis.
Between-mouse variation in circadian rhythm
The adipose-red, heart-blue, kidney-brown, liver-blue, and liver-black modules are correlated and share multiple genes related to apoptotic activity, which varies following circadian rhythm in mice . Several other genes that are known to vary in a circadian fashion are also found in these modules [Additional files 2, 8: Supplemental Figure S5], including Ccrn4l, Fkbp5, Gadd45g, Map3k6, Per1, Pim3, Mt1, Sgk1, Errfi1, Cdkn1a, Dusp1, and Angptl4. The core circadian gene Per2[42, 43] is found in the adipose-red module. Genes that follow a circadian expression pattern are expected to vary with the time of day and with fasting/feeding cycles. Despite our efforts to control both of these factors, between-mouse variation can be expected to arise if the mice are in slightly different phases of their diurnal cycles.
Angptl4, Cdkn1a, Dusp1, and Fkbp5 vary in circadian fashion [43, 44] and are all located in a 7 Mb region on proximal chromosome 17. This region is the strongest example of coexpression clustering that we found in this study. However, statistical assessment suggests that a cluster of this size could be explained by chance.
Between-mouse variation associated with growth hormone
The genes Socs2, Bcl6, Cish, and Gadd45g have correlated patterns of variation in kidney and liver and are among the genes with the most significant between-mouse variation [Additional files 2, 9: Supplemental Figure S6]. Growth hormone has been shown to elicit a strong transcriptional response in Socs2 (positive), Cish (positive), Bcl6 (negative), and Gadd45g (positive), as well as in the growth hormone responders Igf1 and Il6. Three of these genes (Socs2, Bcl6, and Cish) belong to the kidney-pink and liver-magenta modules, which have 12 overlapping genes and are enriched for genes involved in transcription regulation. Growth hormone signalling affects transcription of genes such as Xbp1 (kidney-pink, liver-magenta), which is critical for the regulation of hepatic lipogenesis . The effect of growth hormone signalling appears to extend beyond these modules, however. Among 71 genes previously identified as growth hormone responders , 49 were classified as variable in our study, indicative of widespread individual variation in growth hormone signalling.
Similarities and differences in transcript abundance for sibling cage mates
Sibling cage mates may be expected to exhibit greater similarity than randomly selected mice of the same strain due to shared developmental or micro-environmental factors. When we further partitioned the between mouse variance into between-cage and within-cage components (Methods), we found more genes with significant between-cage variation (adipose, 318; heart, 294; kidney, 1003; liver, 2066) than within-cage variation (adipose, 91; heart, 77; kidney, 639; liver, 1652). The liver-red module provides a striking example of within-cage similarity [Additional files 2, 7: Supplemental Figure S4]. Enrichment for genes associated with fatty acid oxidation in this module could reflect an extended period of fasting just prior to euthanasia. For example, expression of murine hepatic Cyp4a14 (liver-red module) is known to increase in expression under fasting conditions . This gene has been reported to be variable between strains in liver , but it is not clear whether this is a genuine strain-specific effect or an artefact due to co-housing of mice of the same strain.
Other factors could contribute to greater differences between mice within a cage. Cohabitating outbred male mice form a social structure that includes dominance status even when mice are housed as siblings from birth. Dominance behaviour has been observed within male mice of some inbred strains (e.g. CBA, DBA2) but not C57BL/6J. However, cohabitation has known phenotypic effects on C57BL/6J males including change in body weight, adrenal weight, and aggressiveness [49–53]. The factors that determine the social status of siblings raised together are unclear, but once established, social behaviour can reinforce these minor differences leading to distinct individual phenotypes in adult mice.
In our experiment, we observed within-cage body weight difference of as much as 3g (10% of total body weight). Some of the transcriptional changes that we have observed are likely to be related to these body weight differences. For example, in cage 5 we observed a large body weight difference coincident with a large difference in transcription of signature genes for adiposity [Additional files 2, 11: Supplemental Figure S7], but small differences in signature genes for androgen levels [Additional files 2, 14: Supplemental Figure S10]. In contrast, in cages 3 and 4, body weight differences coincide with a transcriptional signature for androgen response but not for adiposity. This suggests that bodyweight differences may reflect two distinct processes, one that affects adiposity and another that affects androgen levels and lean mass [54, 55]. Moreover, these findings provide evidence for an effect of social context on biological processes that have important consequences for human health.
Comparison to a previous study of transcript variation
We directly compared our results to a previous study of transcriptional variation in C57BL/6J mice [3, 56] by computing variance components and applying the same significance tests to both data sets [Additional file 15: Supplemental Table S4]. We found little correlation in total variation (r < 0.09) which we attribute to the predominance of technical variation, especially in the older study. However, we did find good agreement across these studies when we examined specific genes highlighted in the previous study. Cfd was reported to vary significantly between mice in the kidney for the previous experiment in which effects due to dissection and RNA extraction are included in the between-mouse variance component. We also found it to be a variable gene, but, in contrast, we identified Cfd as a gene with primarily within-mouse variation (c = 0.11) in the kidney-black module [Additional files 2, 16: Supplemental Figure S11]. Both studies identified significant between-mouse variation in several highly variable genes, including Gadd45g, Dusp1, Cish, and Bcl6 [Additional files 2, 8, 9: Supplemental Figures S5-S6]. Our study, with a larger sample size, a more recent array technology, and a different experimental design should provide a more precise and detailed picture of variation in gene expression.
Transcript abundance varies significantly among genetically identical male C57BL/6J mice housed under uniform conditions. Patterns of variation can be tissue specific or shared across multiple tissues and transcripts can vary between tissue samples collected from the same animal. Groups of genes with correlated patterns of between-animal or within-animal variation are often enriched for specific functional annotations. We utilized correlation-based clustering to organize a large number of distinct patterns of variation. Literature search tools and functional annotation aided in the interpretation of our findings. However, annotation of gene function is incomplete and this presented some challenges as exemplified by the finding of a skeletal muscle signature in white adipose tissue, which was due to the presence of a related cell type, brown fat.
This study highlights a number of potential biological and technical sources of variation that practitioners should be aware of for both experimental design and interpretation. Much of the between-animal variation reported here reflects functions that are sensitive to environmental cues, such as androgen response, circadian rhythm, and immune response. External environmental cues tend to elicit similar responses in multiple tissues. Variation of gene expression within tissues reflects their heterogeneous cellular composition, and is also a major factor contributing to variation in gene expression. This underscores the potential for dissection or biopsy procedures to introduce unwanted variation into studies of gene expression. Adipose tissue is especially problematic in this regard as it is a highly dynamic and heterogeneous tissue with few anatomical features to guide consistent dissection.
Our tissue collection procedure involved a coarse separation of tissue fragments which, in retrospect, was useful to reveal within-tissue heterogeneity. An exception to this was our use of the intact left and right kidneys as replicates. This may explain the relatively low within-mouse variation observed for this heterogeneous and highly structured tissue. In future studies, we recommend the use of procedures that more effectively homogenize tissue, such as pulverization and mixing of snap frozen samples. Our finding also raises questions about the potential for introducing systematic variation in the dissection of anatomical substructures. This may be a particular concern for studies of gene expression in the brain, for which we have no data at this time.
The presence of biologically meaningful covariation in a setting with no experimental perturbation underscores the need for replication and careful adherence to statistical design principles in gene expression studies. Seemingly innocuous experimental factors such as co-housing of mice can result in systemic differences that may lead to strong statistical support for incorrect conclusions. Prior knowledge of the categories of genes that are intrinsically variable can help to identify such effects. Our study further demonstrates that the variation used to construct statistical tests (error variance) in microarray experiments can have substantial correlation across large sets of genes. This can have a profound impact on testing procedures, especially those that rely on multiple test adjustment of p-values across many genes .
Animals and RNA isolation
We obtained 12 C57BL/6J male mice from The Jackson Laboratory. Six pairs of littermates were housed together from weaning and put on LabDiet's 5k52 diet (standard chow containing 6% fat) in a facility with a 12 h:12 h light:dark cycle beginning with lights on at 6:00 a.m. Animals had ad libitum access to food and acidified water. At 10 weeks of age, body weight was recorded and the mice were euthanized by cervical dislocation and perfused with RNase-free DEPC-treated PBS. Dissection procedures were started at 11:00 a.m. after a 4-hour period of food deprivation and were completed within a one-hour time window. The Jackson Laboratory Animal Care and Use Committee approved the animal housing and experimental procedures described in this work. Inguinal fat pad, heart, liver, and both kidneys were dissected, cut into pieces not exceeding 0.5 cm in any dimension, divided into two samples and placed in 15 ml conical tubes containing RNAlater solution (Ambion, Austin TX). Each kidney sample consisted of one complete kidney, left or right. Tissues were homogenized in TRIzol™ reagent (Invitrogen, Carlsbad, CA). Total RNA was isolated by standard TRIzol™ methods according to the manufacturer's protocols, and quality was assessed using an Agilent 2100 Bioanalyzer instrument and a RNA 6000 Nano LabChip assay (Agilent Technologies, Santa Clara, CA). The RNA was treated with DNase1 (Qiagen, Valencia, Ca.) according to the manufacturer's methods.
Illumina Sentrix® Mouse-6 v1.1 BeadChip processing
Total RNA was reverse transcribed followed by second strand cDNA synthesis. For each sample, an in-vitro transcription (IVT) reaction was carried out incorporating biotinylated nucleotides according to the manufacturer's protocol for Illumina® Totalprep RNA amplification kit (Ambion). 1.5 μg biotin-labelled cRNA was then hybridized onto Mouse-6 Expression BeadChips (Illumina, San Diego CA) for 16 hours at 55°C. Post-hybridization staining and washing were performed according to manufacturer's protocols (Illumina). Illumina Sentrix® Mouse-6 v1.1 BeadChips were scanned using Illumina's BeadStation 500 scanner. Images were checked for grid alignment and then quantified using the BeadStudio software. Control summary graphs generated by BeadStudio were used as quality assurance tools for hybridization, washing stringency, and background. Integrity of the arrays was investigated using the BeadStudio array images and also using bead level image plots generated using the R/beadarray package. Mean pixel intensities by bead type, were created using BeadStudio v3.1 and processed with the R/beadarray package . We performed the experiment in two blocks of three cages, separated by one month. Within each block, we assayed gene expression in each tissue (12 samples) using two Illumina Sentrix® Mouse-6 v1.1 BeadChips. Samples were randomly assigned to array positions within each chip with the constraint that samples from the same mouse were placed on separate chips. Quantile normalization  was applied within each tissue, and a correction for batch effects was applied separately for each gene using an MM-regression estimator from the R/robustbase software package . We selected 45905 probes which are mapped to 22869 genes based on the R/illuminaMousev1p1BeadID.db package . A transcript was called expressed if the mean intensity across the 2 samples of (at least) 1 mouse was above the 95th percentile of the distribution of the mean intensities for the negative control probes. The data are available in accession series GSE20121 from the Gene Expression Omnibus http://www.ncbi.nlm.nih.gov/geo/.
Affymetrix Mouse Gene 1.0 ST Array processing
Following reverse transcription with random T7 primers (Affymetrix, Santa Clara, CA), double stranded cDNA was synthesized with the GeneChip® WT cDNA Synthesis and Amplification Kit (Affymetrix, Santa Clara, CA). In an in vitro transcription (IVT) reaction with T7 RNA polymerase, the cDNA was linearly amplified to generate cRNA. In the second cycle of cDNA synthesis, random primers are used to generate single stranded DNA in the sense orientation. Incorporation of dUTP in the cDNA synthesis step allows for the fragmentation of the cDNA strand utilizing uracil DNA glycosylase (UDG) and apurinic/apyrimidinic endonuclease 1 (APE 1) that specifically recognizes the dUTP and allows for breakage at these residues. Labeling occurs by terminal deoxynucleotidyl transferase (TdT), where biotin is added by an Affymetrix Labeling Reagent. 2.3 μg of biotin-labeled and fragmented cDNA was then hybridized onto GeneChip® Mouse Gene 1.0 ST Arrays (Affymetrix) for 16 hours at 45°C. Post-hybridization staining and washing were performed according to manufacturer's protocols using the Fluidics Station 450 instrument (Affymetrix). Then, the arrays were scanned with a GeneChip™Scanner 3000 laser confocal slide scanner, quantified, and exported to .CEL file format using the GeneChip® Operating Software. Probes were mapped to 34760 probe sets using the R/mogene10stv1.r3cdf package. The .CEL files were processed using the R/affy package using the Robust Multichip Average normalization method . The probe sets were mapped to genes using the R/mogene10sttranscriptcluster.db package . For this experiment, we used a partially balanced incomplete block design method that accommodated hybridization and washing/staining batch factors. Data are available as part of accession series GSE20121 from the Gene Expression Omnibus http://www.ncbi.nlm.nih.gov/geo/.
Identifying variable genes and estimating variance components
For gene g = 1,...,G; mouse i = 1,2,3,...,12; and sample within mouse k = 1,2, we assumed that the log-transformed transcript abundance, y ikg, , satisfies y ikg ~ N(0, σ g 2 ) and considered the null hypothesis H0: σ g 2 = σ2 where σ2 is a fixed variance common to all genes. The alternative is that some genes, g, have excess variability: Ha: σ g 2 > σ2 . To test this hypothesis, we compared the observed distributions of variance to a χ2 distribution for each tissue. The distributions were scaled by dividing each variance by the robust bias-corrected James-Stein estimate . For each tissue, the frequency of genes in the tail of the scaled distribution was then compared to the frequency of a random sample from a χ2 23 distribution. We identified 2000-4000 genes in each tissue with greater than expected variance (α = 0.05) (Table 1). The 2500 most variable genes in each tissue were designated as variable genes and were used in the coexpression network analysis. We chose this number of genes due, in part, to computational constraints of the coexpression network analysis.
We used random effects ANOVA to decompose total variance into between-mouse and within-mouse variance components. Briefly, each y ikg is written as the sum of the average transcript abundance for that gene, μ g , a mouse-specific effect, b ig , and a within-mouse term, w ikg .
The within-mouse term absorbs variation from the mean not accounted for by other terms on the right side of (1). The terms b ig , and w ikg are assumed to satisfy b ig ~ N(0, σbg2 ) and w ikg ~ N(0, σ wg 2 ), respectively. The terms σ bg 2 and σ wg 2 are the between-mouse and within-mouse variance components in this model. Estimates, s bg 2 and s wg 2 , for these components were obtained by residual maximum likelihood (REML) estimation from R/lme4 . A modified F statistic  was used to identify transcripts with significant between-mouse variation. False positive rates were estimated using p-values that were calculated by permuting model residuals. Two types of multiple test corrections were performed. The p-values were adjusted using the Sidak step-down approach , and the Benjamini and Hochberg method . The qvalue software package  was used to estimate the number of genes that do not have significant between-mouse transcript variation, π 0 . To separately assess significance of between-cage and within-cage variation, the following model was used: Each y ikg is written as the sum of the average transcript abundance for that gene, μ g , a cage-specific effect, c ig , a mouse-within-cage term, d j(i)g , and a within-mouse term, w ikg .
The Pritchard et al. (2001)  data were revised to correct a processing error as previously reported . For comparative purposes, we applied the same tests for significance of between-mouse variation described above to the corrected data.
Coexpression network analysis
Variable genes were analysed separately for each tissue using coexpression networks [8, 9]. Every pair of genes was given a weighted connection, r s 2, equal to the square of their correlation coefficient across all samples. Transcript abundance profiles were hierarchically clustered and modules were obtained by a dynamic dendrogram cutting method and subsequent module merge procedure . We only retained modules with more than 25 members. Modules are referenced by their tissue of origin and by a colour index.
For each module, the first principal component was computed to give a representative profile, referred to as the module eigengene . We determined the sign of the module eigengene to be positively correlated with the majority of genes in the module and refer to this majority as the positively-correlated module genes. The complementary genes are referred to as the negatively-correlated module genes. Module eigengenes were scaled to match the median variance over all genes in the module (Figure 4). For each gene, we computed the intraclass correlation coefficient, c ≡ s b 2/(s w 2 + s b 2 ) as a measure of the relative contribution of the between-mouse variance component. We decomposed each gene profile into a between-mouse profile and a within-mouse profile. The between-mouse profile averages the two samples within each mouse and the within-mouse profile is the difference between sample 1 and the average value for that mouse. To measure similarity of between- and within-mouse profiles, we computed Pearson correlation coefficients, r b and r w , for between-mouse (r b ) and within-mouse (r w ) profiles. When assessing significance of similarity of correlation among eigengenes [Additional file 5: Supplemental Table S2], we applied a Fisher transformation with sample size n = 11 (r b ) and n = 12 (r w ). For significance α < 0.05, this required |r b | > 0.66 and |r w | > 0.64.
Gene set enrichment
Each module of the coexpression networks was tested for enrichment within the Gene Ontology (GO) gene sets  and the Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway gene sets [38, 68, 69]. The universe was defined as the set of variable genes present in the relevant database (either GO or KEGG). Only one probe per gene was included in the set of variable genes. The positively- and negatively-correlated subsets of each module were also tested for enrichment. We considered two modules to have significant overlap of functional enrichment if there were 4 genes in each module from a given category and enrichment p-values were less than p < 0.01 for the category in all modules.
We tested for overlap of modules across tissues on a pairwise basis using the hypergeometric test with a Bonferroni multiple-testing correction (α < 0.05). We also used the hypergeometric test to assess the significance of the overlap between published gene lists and modules in our study. In this case, the universe of genes was defined as those assayed in our experiment.
To compare the results of the replicated liver experiments, a map from Illumina probe to Affymetrix probe-set was created based on gene symbol annotation. Where multiple Affymetrix probe sets had the same gene symbol assignment, we selected the one with highest correlation to the Illumina probe. For Affymetrix module eigengene calculation, we excluded Affymetrix probe sets with average intensity less than 7.
To compare our variance component distributions with those of Pritchard et al. (2001), we generated a map from Illumina probe to gene symbol annotation for the Pritchard et al. experiment . Where multiple probe sets had the same gene symbol assignment, we selected the one with highest intraclass correlation coefficient. For this selection and for our comparison of total variation, we excluded the array component of variation for the Pritchard et al. experiment.
Additional resource: Database
An on-line resource has been created to allow access to the experimental data, graphics, and test statistics for all probes in this study:
Churchill G: Fundamentals of experimental design for cDNA microarrays. Nature Genetics Supplement. 2002, 32: 490-495. 10.1038/ng1031.
Koza RA, Nikonova L, Hogan J, Rim J-S, Mendoza T, Faulk C, Skaf J, Kozak LP: Changes in Gene Expression Foreshadow Diet-Induced Obesity in Genetically Identical Mice. PLoS Genetics. 2006, 2 (5): 10.1371/journal.pgen.0020081.
Pritchard C, Hsu L, Delrow J, Nelson P: Project normal: Defining normal variance in mouse gene expression. Proceedings of the National Academy of Sciences, USA. 2001, 98: 13266-13271. 10.1073/pnas.221465998.
Pritchard C, Coil D, Hawley S, Hsu L, Nelson PS: The contributions of normal variation and genetic background to mammalian gene expression. Genome Biol. 2006, 7 (3): R26-10.1186/gb-2006-7-3-r26.
Westfall PH, Young SS: Resampling-based Multiple Testing. 1993, New York.: Wiley
Benjamini Y, Hochberg Y: Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Roy Statist Soc Ser B (Methodological). 1995, 57: 289-300.
Storey J, Tibshirani R: Statistical significance for genomewide studies. Proceedings of the National Academy of Sciences, USA. 2003, 100 (16): 9440-9445. 10.1073/pnas.1530509100.
Ravasz E, Somera A, Mongru D, Oltvai Z, Barabasi A: Hierarchical organization of modularity in metabolic networks. Science. 2002, 297: 1151-1155. 10.1126/science.1073374.
Zhang B, Horvath S: A general framework for weighted gene co-expression network analysis. Statistical Applications in Genetics and Molecular Biology. 2005, 4 (1): Article 17-10.2202/1544-6115.1128.
Horvath S, Dong J: Geometric Interpretation of Gene Coexpression Network Analysis. PLoS Comput Biol. 2008, 4 (8): 10.1371/journal.pcbi.1000117.
Newton M, Quintana F, Den Boon J, Sengupta S, Ahlquist P: Random-set methods identify distinct aspects of the enrichment signal in gene-set analysis. The Annals of Applied Statistics. 2007, 1 (1): 85-106. 10.1214/07-AOAS104.
Watkins-Chow DE, Pavan WJ: Genomic copy number and expression variation within the C57BL/6J inbred mouse strain. Genome Res. 2008, 18 (1): 60-66. 10.1101/gr.6927808.
Raser JM, O'Shea EK: Control of stochasticity in eukaryotic gene expression. Science. 2004, 304 (5678): 1811-1814. 10.1126/science.1098641.
Laforge B, Guez D, Martinez M, Kupiec JJ: Modeling embryogenesis and cancer: an approach based on an equilibrium between the autostabilization of stochastic gene expression and the interdependence of cells for proliferation. Prog Biophys Mol Biol. 2005, 89 (1): 93-120. 10.1016/j.pbiomolbio.2004.11.004.
Elowitz MB, Levine AJ, Siggia ED, Swain PS: Stochastic gene expression in a single cell. Science. 2002, 297 (5584): 1183-1186. 10.1126/science.1070919.
Swain PS, Elowitz MB, Siggia ED: Intrinsic and extrinsic contributions to stochasticity in gene expression. Proc Natl Acad Sci USA. 2002, 99 (20): 12795-12800. 10.1073/pnas.162041399.
Raser JM, O'Shea EK: Noise in gene expression: origins, consequences, and control. Science. 2005, 309 (5743): 2010-2013. 10.1126/science.1105891.
Ramanathan S, Swain PS: Tracing the sources of cellular variation. Dev Cell. 2005, 9 (5): 576-578. 10.1016/j.devcel.2005.10.004.
Kupiec JJ: On the lack of specificity of proteins and its consequences for a theory of biological organization. Prog Biophys Mol Biol. 102 (1): 45-52. 10.1016/j.pbiomolbio.2009.11.002.
Longo D, Hasty J: Dynamics of single-cell gene expression. Mol Syst Biol. 2006, 2: 64-10.1038/msb4100110.
Lapidus S, Han B, Wang J: Intrinsic noise, dissipation cost, and robustness of cellular networks: the underlying energy landscape of MAPK signal transduction. Proc Natl Acad Sci USA. 2008, 105 (16): 6039-6044. 10.1073/pnas.0708708105.
Mar JC, Quackenbush J: Decomposition of gene expression state space trajectories. PLoS Comput Biol. 2009, 5 (12): e1000626-10.1371/journal.pcbi.1000626.
Enver T, Pera M, Peterson C, Andrews PW: Stem cell states, fates, and the rules of attraction. Cell Stem Cell. 2009, 4 (5): 387-397. 10.1016/j.stem.2009.04.011.
Huang S, Eichler G, Bar-Yam Y, Ingber DE: Cell fates as high-dimensional attractor states of a complex gene regulatory network. Phys Rev Lett. 2005, 94 (12): 128701-10.1103/PhysRevLett.94.128701.
Cox LA, Schlabritz-Loutsevitch N, Hubbard GB, Nijland MJ, McDonald TJ, Nathanielsz PW: Gene expression profile differences in left and right liver lobes from mid-gestation fetal baboons: a cautionary tale. J Physiol. 2006, 572 (Pt 1): 59-66.
Richard ID, Parker JS, Lobenhofer EK, Burka LT, Blackshear PE, Vallant MK, Lebetkin EH, Gerken DF, Boorman GA: Transcriptional Profiling of the Left and Median Liver Lobes of Male F344/N Rats Following Exposure to Acetaminophen. Toxicol Pathol. 2005, 33 (1): 111-117. 10.1080/01926230590522257.
Macqueen HA, Waights V, Pond CM: Vascularisation in adipose depots surrounding immune-stimulated lymph nodes. J Anat. 1999, 194 (Pt 1): 33-38. 10.1046/j.1469-7580.1999.19410033.x.
Gibney MJ, Kearney J: Inter- and intra-fat pad variation in vascularization and the release of 14C-labelled fatty acids in mice. Br J Nutr. 1993, 70 (3): 737-745. 10.1079/BJN19930169.
Vohl MC, Sladek R, Robitaille J, Gurd S, Marceau P, Richard D, Hudson TJ, Tchernof A: A survey of genes differentially expressed in subcutaneous and visceral adipose tissue in men. Obes Res. 2004, 12 (8): 1217-1222. 10.1038/oby.2004.153.
Palou M, Priego T, Sanchez J, Rodriguez AM, Palou A, Pico C: Gene expression patterns in visceral and subcutaneous adipose depots in rats are linked to their morphologic features. Cell Physiol Biochem. 2009, 24 (5-6): 547-556. 10.1159/000257511.
Hageman RS, Wagener A, Hantschel C, Svenson KL, Churchill GA, Brockmann GA: High-fat diet leads to tissue-specific changes reflecting risk factors for diseases in DBA/2J mice. Physiol Genomics. 2010, 42 (1): 55-66. 10.1152/physiolgenomics.00072.2009.
Christian M, Parker MG: The Engineering of Brown Fat. J Mol Cell Biol. 2009
Seale P, Bjork B, Yang W, Kajimura S, Chin S, Kuang S, Scime A, Devarakonda S, Conroe HM, Erdjument-Bromage H, et al: PRDM16 controls a brown fat/skeletal muscle switch. Nature. 2008, 454 (7207): 961-967. 10.1038/nature07182.
Cinti S: The adipose organ. Prostaglandins, Leukotrienes, and Essential Fatty Acids. 2005, 73 (9-15):
Koibuchi N, Chin MT: CHF1/Hey2 Plays a Pivotal Role in Left Ventricular Maturation Through Suppression of Ectopic Atrial Gene Expression. Circulation Research. 2007, 100: 850-855. 10.1161/01.RES.0000261693.13269.bf.
Tabibiazar R, Wagner RA, Liao A, Quertermous T: Transcriptional profiling of the heart reveals chamber-specific gene expression patterns. Circ Res. 2003, 93 (12): 1193-1201. 10.1161/01.RES.0000103171.42654.DD.
Mahendroo MS, Cala KM, Hess DL, Russell DW: Unexpected Virilization in Male Mice Lacking Steroid 5alpha-Reductase Enzymes. Endocrinology. 2001, 142 (11): 1652-1662. 10.1210/en.142.11.4652.
Kanehisa M, Goto S: KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Research. 2000, 28: 27-30. 10.1093/nar/28.1.27.
Llamas B, Verdugo RA, Churchill GA, Deschepper CF: Chromosome Y variants from different inbred mouse strains are linked to differences in the morphologic and molecular responses of cardiac cells to postpubertal testosterone. BMC Genomics. 2009, 10: 150-10.1186/1471-2164-10-150.
Paigen K: Mammalian ß-glucuronidase: genetics, molecular biology, and cell biology. Progress in Nucleic Acid Research and Molecular Biology. 1989, 37: 155-205. full_text.
Ijiri K, Potten CS: The circadian rhythm for the number and sensitivity of radiation-induced apoptosis in the crypts of mouse small intestine. Int J Radiat Biol. 1990, 58 (1): 165-175. 10.1080/09553009014551521.
Kawamoto T, Noshiro M, Furukawa M, Honda KK, Nakashima A, Ueshima T, Usui E, Katsura Y, Fujimoto K, Honma S, et al: Effects of Fasting and Re-Feeding on the Expression of Dec1, Per1, and Other Clock-Related Genes. Journal of Biochemistry. 2006, 140 (3): 401-408. 10.1093/jb/mvj165.
Yan J, Wang H, Liu Y, Shao C: Analysis of Gene Regulatory Networks in the Mammalian Circadian Rhythm. PLoS Comput Biology. 2008, 4 (10): e1000193-10.1371/journal.pcbi.1000193. 1000191-1000113
Almon RR, Yang E, Lai W, Androulakis IP, Ghimbovschi S, Hoffman EP, Jusko WJ, DuBoi DC: Relationships between circadian rhythms and modulation of gene expression by glucocorticoids in skeletal muscle. Am J Physiol Regul Integr Comp Physiol. 2008, 295: R1031-R1047. 10.1152/ajpregu.90399.2008.
Chen Y, Lin G, Huo JS, Barney D, Wang Z, Livshiz T, States DJ, Qin ZS, Schwartz J: Computational and functional analysis of growth hormone (GH)-regulated genes identifies the transcriptional repressor B-cell lymphoma 6 (Bc16) as a participant in GH-regulated transcription. Endocrinology. 2009, 150 (8): 3645-3654. 10.1210/en.2009-0212.
Lee AH, Scapa EF, Cohen DE, Glimcher LH: Regulation of hepatic lipogenesis by the transcription factor XBP1. Science. 2008, 320 (5882): 1492-1496. 10.1126/science.1158042.
Huo JS, McEachin RC, Cui TX, Duggal NK, Hai T, States DJ, Schwartz J: Profiles of growth hormone (GH)-regulated genes reveal time-dependent responses and identify a mechanism for regulation of activating transcription factor 3 by GH. J Biol Chem. 2006, 281 (7): 4132-4141. 10.1074/jbc.M508492200.
Savas U, Machemer DE, Hsu MH, Gaynor P, Lasker JM, Tukey RH, Johnson EF: Opposing roles of peroxisome proliferator-activated receptor alpha and growth hormone in the regulation of CYP4A11 expression in a transgenic mouse model. J Biol Chem. 2009, 284 (24): 16541-16552. 10.1074/jbc.M902074200.
Lockwood J, Turney T: Social dominance and stress-induced hypertension: strain differences in inbred mice. Physiology & Behavior. 1981, 26: 547-549.
Razzoli M, Carboni L, Andreoli M, Ballottari A, Arban R: Different susceptibility to social defeat stress of BalbC and C57BL6/J mice. Behav Brain Res. 2011, 216 (1): 100-108. 10.1016/j.bbr.2010.07.014.
Bartolomucci A, Palanza P, Gaspani L, Limiroli E, Panerai AE, Ceresini G, Poli MD, Parmigiani S: Social status in mice: behavioral, endocrine and immune changes are context dependent. Physiol Behav. 2001, 73 (3): 401-410. 10.1016/S0031-9384(01)00453-X.
Bartolomucci A, Cabassi A, Govoni P, Ceresini G, Cero C, Berra D, Dadomo H, Franceschini P, Dell'Omo G, Parmigiani S, et al: Metabolic consequences and vulnerability to diet-induced obesity in male mice under chronic social stress. PLoS One. 2009, 4 (1): e4331-10.1371/journal.pone.0004331.
Bartolomucci A, Pederzani T, Sacerdote P, Panerai AE, Parmigiani S, Palanza P: Behavioral and physiological characterization of male mice under chronic psychosocial stress. Psychoneuroendocrinology. 2004, 29 (7): 899-910. 10.1016/j.psyneuen.2003.08.003.
Mudali S, Dobs AS: Effects of testosterone on body composition of the aging male. Mech Ageing Dev. 2004, 125 (4): 297-304. 10.1016/j.mad.2004.01.004.
Vermeulen A, Goemaere S, Kaufman JM: Testosterone, body composition and aging. J Endocrinol Invest. 1999, 22 (5 Suppl): 110-116.
Cui X, Churchill G: How many mice and how many arrays? Replication in mouse cDNA microarray experiments. Methods of Microarray Data Analysis III, Papers from CAMDA '02. Edited by: JKaL SM. 2003
Efron B: Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction: Cambridge University Press. 2010
Dunning M, Thorne N, Camilier I, Smith M, Tavare S: Quality control and low-level statistical analysis of Illumina beadarrays. REVSTAT - Statistical Journal. 2006, 4 (1): 1-30.
Bolstad B, Irizarry R, Astrand M, Speed T: A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics. 2003, 19 (2): 185-193. 10.1093/bioinformatics/19.2.185.
Rousseeuw P, Croux C, Todorov V, Ruckstuhl Andreas, Salibian-Barrera M, Verbeke T, Maechler M: robustbase: Basic Robust Statistics. R package version 0.4-5. 2009
Dunning M, Barbosa-Morais N, Ritchie M: Illumina Mousev1p1 annotation data (chip illuminaMousev1p1BeadID). illuminaMousev1p1BeadIDdb. R package version 1.2.0. edn. 2009
Li A: Affymetrix Mouse Gene 1.0-ST Array Transcriptcluster Revision 5 annotation data (chip mogene10sttranscriptcluster) assembled using data from public repositories. 2011, Bioconductor
Tong T, Wang Y: Optimal shrinkage estimation of variances with applications to microarray data analysis. Journal of the American Statistical Association. 2007, 102 (477): 113-122. 10.1198/016214506000001266.
Bates D: Linear mixed model implementation in lme4: University of Wisconsin-Madison, Department of Statistics. 2008
Cui X, Hwang J, Qiu J, Blades N, Churchill G: Improved statistical tests for differential gene expression by shrinking variance components estimates. Biostatistics. 2005, 6 (1): 59-75. 10.1093/biostatistics/kxh018.
Langfelder P, Horvath S: WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008, 9: 559-10.1186/1471-2105-9-559.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, et al: Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000, 25 (1): 25-29. 10.1038/75556.
Kanehisa M, Goto S, Hattori M, Aoki-Kinoshita K, Itoh M, Kawashima S, Katayama T, Araki M, Hirakawa M: From genomics to chemical genomics: new developments in KEGG. Nucleic Acids Research. 2006, 34: D354-357. 10.1093/nar/gkj102.
Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, Katayama T, Kawashima S, Okuda S, Tokimatsu T, et al: KEGG for linking genomes to life and the environment. Nucleic Acids Research. 2008, 36: D480-D484. 10.1093/nar/gkm882.
Human and Mouse Prostrate Information. [http://pedb.org]
Funded by the National Institute of General Medical Sciences (NIGMS) National Centers for Systems Biology: Center for Genome Dynamics [GM076468] to GAC and by the NIGMS Ruth L. Kirschstein National Research Service Award [F32GM087849-01] to PTV. We thank Genevieve Billings and Holly Savage for animal care tissue collection, Sonya Kamdar and The Jackson Laboratory Gene Expression Services for performing the microarray experiments, Ricardo Verdugo for helpful discussions, Jesse Hammer for assistance in preparing the figures, Joel Graber and Imogen Hurley for careful review of the manuscript, and Keith Sheppard for implementing the online data resource.
KLS and GAC conceived the experiment. KLS directed the animal husbandry and tissue collection. GAC and PTV designed and implemented microarray experiments. PTV performed the data analysis. PTV and GAC wrote the manuscript. All authors read and approved the final manuscript.