Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Stochastic variation of transcript abundance in C57BL/6J mice

BMC Genomics volume 12, Article number: 167 (2011) Cite this article

Abstract

Background

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.

Results

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.

Conclusions

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.

Background

Variation in transcript abundance between individuals has important implications for microarray experimental design and significance testing [1]. 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) [2] 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) [3] 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) [4] 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.

Figure 1
figure 1

Experimental design. Twelve C57BL/6J male mice were co-housed as pairs in six cages from weaning until 10 weeks of age (A). Liver, heart, kidney, and adipose tissues were collected from each mouse and split into two samples per tissue per individual and processed separately to generate RNA (B). In the case of kidney, the samples consist of the entire left or right kidney, respectively. Other tissues were chopped into small pieces, which were separated before placing them into collection tubes.

Full size image

Results

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).

Table 1 Variability of transcript abundance
Full size table

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.

Figure 2
figure 2

Distribution of within- and between-mouse variance components. The total variance (A) and estimated between-mouse (B) and within-mouse (C) variance components are shown as smoothed density histograms after square root transformation. The vertical bars (left panel of B) represent the percentage of probes with estimated variance components less than 1×10-4. These probes are not included in the density plot (right panel of B). Tissues are indicated by colour (Adipose: blue; Heart: red; Kidney: green; Liver: black).

Full size image

We decomposed total variance for each gene into within-mouse (sw 2 ) and between-mouse (sb 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 [5] (as in Pritchard et al. (2001) [3]) and found few genes with significant between-mouse variation (Table 1C). We applied a false discovery rate (FDR) adjustment [6] (as in Pritchard et al. (2006) [4]) 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 [7] 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 [10] 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).

Figure 3
figure 3

Co-expression module clustering of variable genes. The topographical overlap similarity [9] of gene pairs is shown as a heatmap for adipose (A), heart (B), kidney (C), and liver (D). Gene pairs with low similarity are shown in yellow, and those with high similarity are shown in red. Coexpression modules define highly connected sets of correlated expression profiles as indicated by colour codes in the margins of the maps. The colour names are given in Table 2. Genes are ordered within module from left to right by increasing topographical overlap with the module eigengene. Genes that are not assigned to modules are not included in this plot.

Full size image
Table 2 Module highlights
Full size table
Figure 4
figure 4

Transcript abundance profiles of co-expression modules. Module eigengenes representing average transcript abundance across genes are shown for the modules listed in Table 2. Mice are indexed from 1 to 12 and vertical partitions indicate cage pairings. Eigengene profiles are zero-centered and represent an average log2 fold-change (fc) across all genes in the module. All panels are shown with the same y-axis scaling for ease of comparison. Horizontal lines connect mean values for each mouse. Vertical line segments connect the two within-mouse samples. Sample 1 is indicated by the upright triangle and sample 2 by the inverted triangle. For adipose, heart, and liver, samples should not be compared across mice or across tissues, as there is no correspondence between sample labels at this level of the experiment. In the kidney, sample 1 is the left kidney and sample 2 is the right kidney. Modules with 25 or fewer genes are not shown.

Full size image

For each gene, we computed the intraclass correlation coefficient, c ≡ sb 2/(sw 2 + sb 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 [11] 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, rb , 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, rb , 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 (rb = 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, rb , 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 (rb ≥ 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 (rb = 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 (rb = - 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, rw . 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 (|rw | ≥ 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]).

Discussion

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 (rb > 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 [12]. 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 [3].

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 [2]. 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 [27]. Vascularised adipose tissue tends to be more metabolically active [28]. 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 [2] [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.

Figure 5
figure 5

Within-tissue correlation of eigengene transcript abundance profiles. The eigengenes of the adipose-green and adipose-magenta modules exhibit negative within-mouse correlation and positive between-mouse correlation (rw = -0.65, rb = 0.69) Profiles are shown for (A) Ccl6 (l) and Ccl9 (r) of the positively correlated adipose-green module and KEGG cytokine-cytokine receptor interaction pathway, (B) Cd8b1 (l) and Elf1 (r) of the positively correlated adipose-magenta module and GO immune response biological process, and (C) Lep (l) and Trp53inp2 (r), of the negatively correlated adipose-magenta module. The graphical features of the eigengene plots of Figure 4 apply to these gene plots. Profiles are coloured by module membership. Summary statistics for these genes are available [Additional file 10: Supplemental Table S3].

Full size image

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 [31] indicates that other fat depots are at least as variable as the inguinal depot. The Koza et al. study [18] 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 (rb = 0.86; rw = 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 [34].

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[35]. 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 (rw = 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) [36]) 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 [37] suggesting that androgenic effects in individuals with higher Srd5a2 expression could be more pronounced. Hsd11b1 facilitates the conversion of testosterone to adrenosterone [38] and has been shown to be androgen-responsive in mice [39]. 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 [40]. Gusb has significant between-mouse variation that is correlated with the kidney-green module eigengene (rb = 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 [41]. 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[45]. 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 [46]. The effect of growth hormone signalling appears to extend beyond these modules, however. Among 71 genes previously identified as growth hormone responders [47], 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 [48]. This gene has been reported to be variable between strains in liver [4], 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.

Conclusions

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 [57].

Methods

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.

Microarray processing

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 [58]. 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 [59] 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 [60]. We selected 45905 probes which are mapped to 22869 genes based on the R/illuminaMousev1p1BeadID.db package [61]. 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 [59]. The probe sets were mapped to genes using the R/mogene10sttranscriptcluster.db package [62]. 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, yikg, , satisfies yikg ~ 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 [63]. 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 χ223 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 yikg is written as the sum of the average transcript abundance for that gene, μg , a mouse-specific effect, big , and a within-mouse term, wikg .

(1)

The within-mouse term absorbs variation from the mean not accounted for by other terms on the right side of (1). The terms big , and wikg are assumed to satisfy big ~ N(0, σbg2 ) and wikg ~ 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, sbg 2 and swg 2 , for these components were obtained by residual maximum likelihood (REML) estimation from R/lme4 [64]. A modified F statistic [65] 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 [5], and the Benjamini and Hochberg method [6]. The qvalue software package [7] 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 yikg is written as the sum of the average transcript abundance for that gene, μg , a cage-specific effect, cig , a mouse-within-cage term, dj(i)g , and a within-mouse term, wikg .

(2)

The Pritchard et al. (2001) [3] data were revised to correct a processing error as previously reported [56]. 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, rs2, 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 [66]. 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 [10]. 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 ≡ sb 2/(sw 2 + sb 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, rb and rw , for between-mouse (rb ) and within-mouse (rw ) 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 (rb ) and n = 12 (rw ). For significance α < 0.05, this required |rb | > 0.66 and |rw | > 0.64.

Gene set enrichment

Each module of the coexpression networks was tested for enrichment within the Gene Ontology (GO) gene sets [67] 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.

Module overlap

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.

Across-experiment comparison

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 [70]. 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:

http://cgd.jax.org/individualvariation.shtml.

References

  1. Churchill G: Fundamentals of experimental design for cDNA microarrays. Nature Genetics Supplement. 2002, 32: 490-495. 10.1038/ng1031.

    CAS  Google Scholar 

  2. 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.

    PubMed  PubMed Central  Google Scholar 

  3. 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.

    CAS  Google Scholar 

  4. 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.

    PubMed  PubMed Central  Google Scholar 

  5. Westfall PH, Young SS: Resampling-based Multiple Testing. 1993, New York.: Wiley

    Google Scholar 

  6. 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.

    Google Scholar 

  7. 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.

    CAS  Google Scholar 

  8. 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.

    Google Scholar 

  9. 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.

    Google Scholar 

  10. Horvath S, Dong J: Geometric Interpretation of Gene Coexpression Network Analysis. PLoS Comput Biol. 2008, 4 (8): 10.1371/journal.pcbi.1000117.

    PubMed  PubMed Central  Google Scholar 

  11. 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.

    Google Scholar 

  12. 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.

    CAS  PubMed  PubMed Central  Google Scholar 

  13. Raser JM, O'Shea EK: Control of stochasticity in eukaryotic gene expression. Science. 2004, 304 (5678): 1811-1814. 10.1126/science.1098641.

    CAS  PubMed  PubMed Central  Google Scholar 

  14. 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.

    CAS  PubMed  Google Scholar 

  15. 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.

    CAS  PubMed  Google Scholar 

  16. 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.

    CAS  PubMed  PubMed Central  Google Scholar 

  17. Raser JM, O'Shea EK: Noise in gene expression: origins, consequences, and control. Science. 2005, 309 (5743): 2010-2013. 10.1126/science.1105891.

    CAS  PubMed  PubMed Central  Google Scholar 

  18. Ramanathan S, Swain PS: Tracing the sources of cellular variation. Dev Cell. 2005, 9 (5): 576-578. 10.1016/j.devcel.2005.10.004.

    CAS  PubMed  Google Scholar 

  19. 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.

    CAS  PubMed  Google Scholar 

  20. Longo D, Hasty J: Dynamics of single-cell gene expression. Mol Syst Biol. 2006, 2: 64-10.1038/msb4100110.

    PubMed  PubMed Central  Google Scholar 

  21. 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.

    CAS  PubMed  PubMed Central  Google Scholar 

  22. Mar JC, Quackenbush J: Decomposition of gene expression state space trajectories. PLoS Comput Biol. 2009, 5 (12): e1000626-10.1371/journal.pcbi.1000626.

    PubMed  PubMed Central  Google Scholar 

  23. 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.

    CAS  PubMed  Google Scholar 

  24. 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.

    PubMed  Google Scholar 

  25. 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.

    CAS  PubMed  PubMed Central  Google Scholar 

  26. 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.

    Google Scholar 

  27. 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.

    PubMed  PubMed Central  Google Scholar 

  28. 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.

    CAS  PubMed  Google Scholar 

  29. 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.

    CAS  PubMed  Google Scholar 

  30. 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.

    CAS  PubMed  Google Scholar 

  31. 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.

    CAS  PubMed  PubMed Central  Google Scholar 

  32. Christian M, Parker MG: The Engineering of Brown Fat. J Mol Cell Biol. 2009

    Google Scholar 

  33. 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.

    CAS  PubMed  PubMed Central  Google Scholar 

  34. Cinti S: The adipose organ. Prostaglandins, Leukotrienes, and Essential Fatty Acids. 2005, 73 (9-15):

    CAS  Google Scholar 

  35. 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.

    CAS  PubMed  Google Scholar 

  36. 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.

    CAS  PubMed  Google Scholar 

  37. 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.

    Google Scholar 

  38. Kanehisa M, Goto S: KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Research. 2000, 28: 27-30. 10.1093/nar/28.1.27.

    CAS  PubMed  PubMed Central  Google Scholar 

  39. 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.

    PubMed  PubMed Central  Google Scholar 

  40. Paigen K: Mammalian ß-glucuronidase: genetics, molecular biology, and cell biology. Progress in Nucleic Acid Research and Molecular Biology. 1989, 37: 155-205. full_text.

    CAS  PubMed  Google Scholar 

  41. 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.

    CAS  PubMed  Google Scholar 

  42. 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.

    CAS  PubMed  Google Scholar 

  43. 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

    Google Scholar 

  44. 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.

    CAS  PubMed  PubMed Central  Google Scholar 

  45. 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.

    CAS  PubMed  PubMed Central  Google Scholar 

  46. 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.

    CAS  PubMed  PubMed Central  Google Scholar 

  47. 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.

    CAS  PubMed  Google Scholar 

  48. 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.

    CAS  PubMed  PubMed Central  Google Scholar 

  49. Lockwood J, Turney T: Social dominance and stress-induced hypertension: strain differences in inbred mice. Physiology & Behavior. 1981, 26: 547-549.

    CAS  Google Scholar 

  50. 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.

    PubMed  Google Scholar 

  51. 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.

    CAS  PubMed  Google Scholar 

  52. 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.

    PubMed  PubMed Central  Google Scholar 

  53. 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.

    PubMed  Google Scholar 

  54. 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.

    CAS  PubMed  Google Scholar 

  55. Vermeulen A, Goemaere S, Kaufman JM: Testosterone, body composition and aging. J Endocrinol Invest. 1999, 22 (5 Suppl): 110-116.

    CAS  PubMed  Google Scholar 

  56. 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

    Google Scholar 

  57. Efron B: Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction: Cambridge University Press. 2010

    Google Scholar 

  58. 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.

    Google Scholar 

  59. 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.

    CAS  PubMed  Google Scholar 

  60. 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

    Google Scholar 

  61. Dunning M, Barbosa-Morais N, Ritchie M: Illumina Mousev1p1 annotation data (chip illuminaMousev1p1BeadID). illuminaMousev1p1BeadIDdb. R package version 1.2.0. edn. 2009

    Google Scholar 

  62. Li A: Affymetrix Mouse Gene 1.0-ST Array Transcriptcluster Revision 5 annotation data (chip mogene10sttranscriptcluster) assembled using data from public repositories. 2011, Bioconductor

    Google Scholar 

  63. 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.

    Google Scholar 

  64. Bates D: Linear mixed model implementation in lme4: University of Wisconsin-Madison, Department of Statistics. 2008

    Google Scholar 

  65. 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.

    PubMed  Google Scholar 

  66. Langfelder P, Horvath S: WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008, 9: 559-10.1186/1471-2105-9-559.

    PubMed  PubMed Central  Google Scholar 

  67. 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.

    CAS  PubMed  PubMed Central  Google Scholar 

  68. 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.

    CAS  PubMed  Google Scholar 

  69. 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.

    CAS  PubMed  Google Scholar 

  70. Human and Mouse Prostrate Information. [http://pedb.org]

Download references

Acknowledgements

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.

Author information

Authors and Affiliations

  1. The Jackson Laboratory, 600 Main Street, Bar Harbor, Maine, 04609, USA

    Peter T Vedell, Karen L Svenson & Gary A Churchill

Authors
  1. Peter T Vedell
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Karen L Svenson
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Gary A Churchill
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Gary A Churchill.

Additional information

Authors' contributions

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.

Electronic supplementary material

Additional file 1:Supplemental Figure S1. P-value histograms for between-mouse significance tests. (PDF 66 KB)

Additional file 2:Supplemental Figure Captions. Captions for supplemental figures. (DOC 230 KB)

Additional file 3:Supplemental Table S1. Enrichments scores of the gene coexpression modules. (XLSX 6 MB)

Additional file 4:Supplemental Figure S2. Graphical model showing relationships between modules. (PDF 76 KB)

Additional file 5:Supplemental Table S2. Pairwise Pearson correlations for module eigengenes. (XLSX 39 KB)

Additional file 6:Supplemental Figure S3. Cross-platform comparison of liver eigengene profiles. (PDF 67 KB)

12864_2010_10022_MOESM7_ESM.PDF

Additional file 7:Supplemental Figure S4. Transcript abundance profiles for fatty acid metabolism genes in the liver. (PDF 84 KB)

Additional file 8:Supplemental Figure S5. Transcript abundance profiles for circadian rhythm genes. (PDF 142 KB)

12864_2010_10022_MOESM9_ESM.PDF

Additional file 9:Supplemental Figure S6. Transcript abundance profiles for growth-hormone regulated genes in kidney and liver. (PDF 86 KB)

12864_2010_10022_MOESM10_ESM.DOC

Additional file 10:Supplemental Table S3. Within-mouse correlation statistics for selected genes in adipose. (DOC 160 KB)

12864_2010_10022_MOESM11_ESM.PDF

Additional file 11:Supplemental Figure S7. Transcript abundance profiles for variable genes reported in adipose tissue. (PDF 89 KB)

12864_2010_10022_MOESM12_ESM.PDF

Additional file 12:Supplemental Figure S8. Transcript abundance profiles for variable brown fat signature genes in white fat tissue. (PDF 82 KB)

12864_2010_10022_MOESM13_ESM.PDF

Additional file 13:Supplemental Figure S9. Transcript abundance profiles showing region-specific variation of gene expression in heart. (PDF 73 KB)

12864_2010_10022_MOESM14_ESM.PDF

Additional file 14:Supplemental Figure S10. Transcript abundance profiles for androgen-regulated variable genes in the kidney. (PDF 47 KB)

12864_2010_10022_MOESM15_ESM.DOC

Additional file 15:Supplemental Table S4. Transcript abundance variation statistics for Pritchard et al (2001) dataset. (DOC 150 KB)

Additional file 16:Supplemental Figure S11. Transcript abundance profile for Cfd gene in kidney. (PDF 41 KB)

Authors’ original submitted files for images

Rights and permissions

This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Cite this article

Vedell, P.T., Svenson, K.L. & Churchill, G.A. Stochastic variation of transcript abundance in C57BL/6J mice. BMC Genomics 12, 167 (2011). https://doi.org/10.1186/1471-2164-12-167

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/1471-2164-12-167

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us