Definition, conservation and epigenetics of housekeeping and tissue-enriched genes

HKGs and TEGs

We identified 1,522 HKGs from a total of 18,149 genes in 42 normal human tissues monitored on the microarray (see Methods). This list of HKGs was used for analysis of genomic and epigenetic features (see Additional file 1). We also identified 975 TEGs from a subset of 29 representative tissues. These TEGs are expressed at much higher level in one single tissue than any other tissues (see Methods. Additional file 2). TEGs were found in 26 tissues, while no TEGs meeting our criteria were identified in spleen, colon and CD4+ T-cells.

The 20 HKGs with the highest and most consistent expression (See Methods) were selected from this list as the best candidates to serve as reference HKGs in laboratory experiments (Figure 1A). Three of the 20 highest quality HKGs, GAPDH, ACTB and UBC, are commonly used HKGs for experimental controls. Expression data for several other commonly used or commercially available HKGs ([11] and http://invitrogen.com) not included in the top 20 HKG lists (Figure 1B) illustrate that genes commonly used as controls do not necessarily show high expression with low variance across diverse tissues, and in fact, the expression of some of these genes varies by more than an order of magnitude across tissues.

There have been other recent efforts to identify HKGs in human tissues to serve as experimental references or internal controls [12–16], but these studies have significant shortcomings. One study surveyed a much smaller set of 7012 potential candidate genes [12], while another one used a much smaller set of 15 tissues [16]. Others were based on microarray data from heterogeneous sources lacking systematic experimental controls [15] or were based on in silico predictions [14]. The HKGs identified here are based on high quality microarray expression data systematically gathered from a large and diverse set of tissues. This is a systematically and experimentally defined set of HKG which have both high expression and low fluctuation across all major organ/tissues.

The human and mouse transcriptome in multiple tissues have been surveyed in microarray studies [17–19] that built foundations for studies on housekeeping and tissue-specific genes [13, 19, 20]. Large collections of EST and SAGE data have also been used to identify HKG [21] and tissue-specific genes [16]. Comparing HKGs of this study and other studies based on microarray [13] or EST datasets [21], significant portions of genes in the three HKG sets overlap, while the HKG list described here has the fewest genes unique in a single study, suggesting our fluctuation-controlled microarray approach is more conservative than the other methods that either depend on sampling or representation in an EST dataset [21] or lack control of variation across different tissues [13] (see Additional file 3). We also compared our TEGs in testis, prostate, liver and skin with tissue-specific genes from another study [19] (see Additional file 4). A significant portion of our TEGs overlap with this human tissue-specific gene study, particularly for liver, in which 70% of the defined TEGs are identical. Other tissues are more discrepant, as a result of different tissue selection for the surveys and different criteria used to identify these genes.

We next investigated the functions of HKGs by testing for enrichment of Panther Biological Process annotations [22]. As expected, HKGs are enriched for biological processes related to basic maintenance of the cell, including protein biosynthesis, pre-mRNA processing, cell cycle and intracellular protein trafficking (Figure 2). In contrast, TEGs were observed to be enriched largely for tissue-specific biological processes, as expected (see Additional file 5). For instance, TEGs of bone marrow, brain, kidney, liver and skeletal muscle are enriched in immunity and defense, synaptic transmission, ion transport, lipid metabolism and muscle contraction, respectively. These biological processes are presumably less likely to be essential for cell maintenance and survival.

Conservation of functions and expression patterns in HKGs across species

Distribution of segmental duplication, copy number variation sites and ultraconserved elements in HKGs and TEGs

We also calculated the distribution of ultra conserved elements (UCEs) which are sequences that are absolutely conserved (100% identical) between orthologous regions of the human, rat, and mouse genomes [32]. Consistent with the slower evolution of HKGs, UCEs are significantly enriched in HKGs and not changed in TEG (Table 3). HKGs have been found to evolve more slowly than TEGs at the sequence level point mutation [5]. The distribution of SD, CNV and UCE demonstrated that HKGs are also more conserved than TEGs with respect to genomic structural changes.

Enriched CpG islands at HKG transcription start sites

Chromatin structure and epigenetic modifications in HKGs

We next examined differences in chromatin structure and epigenetic modifications, including nucleosome occupancy, histone modifications, and DNA methylation between HKG and TEGs as possible mechanisms contributing to the differential expression patterns of these two groups of genes.

DNase I hypersensitive (HS) sites, formed by nucleosome-free chromatin regions, have been used to identify locations of many different types of regulatory regions including enhancers, silencers, promoters, insulators, and locus control regions [35, 36]. Figure 3A shows the distribution of HS sites in CD4+ T cells [35] at TSS for HKGs, TEGs and RefSeq genes. The density of HS sites peaks at TSS, as expected, and rapidly drops to background levels beyond 1 kb from TSS. HKGs show an elevated density of HS sites relative to RefSeq genes, indicating that TSS of HKGs are less likely to be packaged into nucleosomes and are more exposed to transcription and regulatory factors. HS site density in TEGs is low, consistent with the low levels of expression of these genes in CD4+ T cells, the cell type from which the HS sites data were generated.

The HKGs identified in this study are highly expressed in CD4+ T cells (and all other tissues), leading to the possibility that the differences observed in HS site density seen for HKGs and TEGs may reflect only the overall expression level of these genes in CD4+ T cells, rather than the difference in expression patterns across tissues. To address this question, we partitioned both the HKGs and RefSeq genes into subgroups based on their expression level in CD4+ T cells: low, intermediate and high (Figure 3B). While a correlation between HS site density and expression level is still observed across the subgroups for either HKGs or RefSeq genes, the HKG-low expression subgroup (average expression intensity: 2.55) has a higher HS site density than the RefSeq-high expression subgroup (average intensity: 3.48), clearly demonstrating that the HS site density is not simply a function of gene expression level in CD4+ T cells, but also correlates with the high levels of expression across different tissues. A recent study showed the positive association of CpG density with the distribution of HS sites across different tissues [37], suggesting that the increase in HS sites in HKGs may be related to high CpG density. Another possible explanation for this observation is that HKGs may contain sequence elements at their TSS that inhibit formation of nucleosomes, leading to high promoter accessibility and higher expression levels of these genes across different tissues. Further investigation of TSS sequences and more HS site mapping in other tissues would be necessary to test this hypothesis.

Histone acetylation regulates gene expression by allowing transcription factors access to promoters in the chromatin by neutralizing the positive charge of histone tails and weakening their contact with negatively charged DNA [38]. Using histone acetylation data of chromosome 21 and 22 derived from liver Hep G2 cells[39], we show that the acetylation ratio of genes peaks at the TSS (Figure 4). The histone acetylation density observed for HKGs is higher than RefSeq genes on average, while the density observed for TEGs, which are not expressed in liver Hep G2 cells, is almost at the background level (Figure 4A). This observation is consistent with recent genome-wide studies in which histone acetylation is positively correlated with transcription factor binding or gene expression [40–42]. To study the correlation of gene expression with histone acetylation, HKGs and RefSeq genes were partitioned into subgroups based on their expression level in Hep G2 cells. 74% of all RefSeq genes have low expression (intensity < 1) in the Hep G2 cells and the higher levels of histone acetylation levels for HKGs relative to RefSeq genes was maintained. The difference in histone acetylation between RefSeq genes and HKGs with medium (1 < intensity < 3) and high (intensity > 3) expression levels, however, virtually disappears. It appears that histone acetylation density depends only on gene expression level, indicating that histone acetylation of TSS does not contribute to the expression pattern of HKGs across tissues (Figure 4B). Interestingly, while histone acetylation density is low for genes with low expression levels, it is essentially identical for all subgroups with expression intensity >1. This observation of a threshold-like relationship suggests that histone acetylation of TSS may be a necessary condition for gene expression and may serve as a "transcriptional switch" that opens the chromatin structure and allows other transcription factors to regulate gene expression. Significant differences between HKGs and TEGs in both nucleosome occupancy and histone acetylation suggest that the regulation of gene expression for these different groups is affected by multiple epigenetic factors.

We also compared the transcription factor binding and histone methylation density between HKGs and TEGs using data collected in CD4+ T cells [43] (Figure 5). As expected, binding of the transcription factor pol II and regulatory factor CTCF reaches a peak level at the exact position of TSS. CTCF is a zinc finger protein with diverse regulatory functions, including a role in mediating chromatin interactions to form the genomic three-dimensional structure[44]. Recent studies on HKG clusters around the α-Globin and β-Globin loci suggest chromatin loop/hub forms during the transcription of the gene clusters [45, 46]. Our results suggest CTCF plays a role in HKG transcription. Transcription factor binding is highest for HKGs among all gene groups, while transcription factor binding of TEGs, which are not expressed in CD4+ T cells, is near background levels (Figure 5A, B). The patterns observed for different histone methylation sites of HKGs and TEGs are complex, likely reflecting the complex relationship between histone methylation and transcription (Figure 5C–I). Differences in the shape of the distribution are observed between HKGs, TEGs and RefSeq genes for some histone methylations. When similar distribution shapes are observed, histone methylations may be either positively or negatively correlated with the expression level. These features suggest that histone methylation likely contributes to the differential gene expression patterns of these genes in a complex fashion.

DNA methylation in mammals occurs on cytosine residues of CpG dinucleotides, which may lead to formation of heterochromatin, imprinting and transcriptional repression [47]. The distribution of genome-wide DNA methylation [48] in HKGs, RefSeq genes and TEGs (see Additional file 6) shows that DNA methylation peaks at TSS for all gene groups and that there is no significant difference of DNA methylation levels between HKGs, TEGs and RefSeq genes in either sperm cells or fibroblast cells. Additionally, comparison of the list of methylated genes from another recent study [49] with our HKGs and TEGs did not yield any significant overlapping genes (data not shown). Based on these two pieces of evidence, HKGs do not appear to be enriched for DNA methylation, despite enrichment for CpG islands. This observation is consistent with previous reports that CpG islands in normal tissues are protected from methylation and that methylation of CpG islands is one of the mechanisms of tumorigenesis [50–52].

Tissues

mRNA from human tissues was purchased from commercial vendors, including Clontech, Ambion, and Biochain. Most samples were pooled from multiple donors, typically twelve.42 normal tissues were tested, in cluding adipose, adrenal gland, bladder, activated CD4-positive T-lymphocyte, activated CD8-positive T-lymphocyte, bone marrow, brain, fetal brain, cerebellum, cerebral cortex, hippocampus, thalamus, pituitary gland, cervix uteri, colon, epididymis, heart, kidney, fetal kidney, liver, fetal liver, lung, fetal lung, trachea, lymph node, mammary gland, skeletal muscle, ovary, placenta, prostate, retina, salivary gland, skin, duodenum, ileum, jejunum, spinal cord, spleen, stomach, testis, thymus, and thyroid gland. These selected tissues cover most major organs and normal tissue types. Four fetal tissues of brain, kidney, liver and lung were included.

Microarray expression profiling

Human tissue microarray expression profiling was performed as described previously [53]. In brief, purchased mRNA pooled from multiple normal individuals was amplified and labeled using a full-length amplification protocol and hybridized in duplicate against a common reference pool in a two-color dye swap experiment [54]. Each gene is represented by 3 microarray probes placed at exon-exon junctions or in exons. Gene expression was calculated as the median probe intensity, after normalization by the pool of all data. The dataset is available at National Center for Biotechnology Information's Gene Expression Omnibus database [GEO accession: GSE16546].

Selection of HKGs and TEGs

We used fairly conservative criteria to identify HKGs: the intensity of the gene must be greater than the median intensity of all genes in the microarray in at least 41 out of 42 tissues and the coefficient of variance (CV, standard deviation/average) of the gene intensity across tissues must be less than 1. The intensity and CV of the 18,149 genes monitored in the microarray are distributed over a wide range, with average intensity of all genes 1.04 ± 1.94 (SD) and average CV of all genes 0.83 ± 0.77 (SD). A recent study shows that genes' breadth of expression in tissues is positively correlated with the expression level of the genes [24]. Therefore it is reasonable to select HKGs from among those genes with higher intensity. While the CVs of most genes (76% of all genes) are below 1, some genes' expression is very volatile across tissues, with CV as high as 6. Our criteria guarantee the HKGs are highly expressed in vast majority of tissues with limited fluctuation in intensity level across tissues.

More stringent criteria were used to identify a reference HKG list for laboratory experimental controls. We required that the intensity of each HKG be greater than the median of all genes in each of the 42 tissues and CV of intensity less than 0.35. A total of 362 HKGs meet these criteria. The top 20 genes ranked by their average intensity across all 42 tissues were selected as the experimental housekeeping genes reference.

To identify TEGs, we selected 29 representative tissue types, removing fetal and redundant tissues from the set of 42 tissues described above. The resulting set was as follows: adipose, adrenal gland, bladder, bone marrow, brain, cervix uteri, colon, heart, kidney, liver, lung, trachea, mammary gland, ovary, skeletal muscle, lymph node, placenta, prostate, retina, salivary gland, skin, spinal cord, spleen, stomach, testis, thymus, thyroid gland, jejunum, and CD4-positive T-lymphocyte. To be identified as a TEG, the intensity of the gene in the relevant tissue was required to meet three criteria: 1) among the top 25% percentile of all genes in that particular tissue; 2) greater than 50% of the sum of intensities for that gene in all other tissues in the set of 29; and 3) greater than three times of intensity of the gene of interest in any other tissue.

Conservation of functions

We used the number of orthologs of human genes in other eukaryotic species as identified by NCBI HomoloGene [23] as an indication of functional conservation across species. We mapped human HKGs, TEGs and all genes represented in the microarray to orthologs in mouse, rat, dog, fly (D. melanogaster), worm (C. elegans) and budding yeast (S. cerevisiae). The numbers of human genes that map to genes of other species through HomoloGene are counted. Student's T-tests were applied between orthologs of HKG and all genes and between orthologs of TEG and all genes.

Distribution of SD and CNV in genes

We required at least a quarter of the total genomic length of a gene to overlap the SD or CNV region (Table 3). The p-values, indicating the statistical significance of the overlap for HKGs and TEGs relative to all RefSeq genes, were calculated according to the hypergeometric distribution with a Bonferroni correction.

CpG islands

CpG islands coordinates were obtained from UCSC genome browser http://genome.ucsc.edu human CpG island track. The number and length of CpG islands located within 500 bp upstream and downstream of transcription start sites and end sites are calculated for HKGs, TEGs and RefSeq genes. CpG density is indicated by the fraction of base pairs occupied by CpG islands. The hypergeometric distribution with Bonferroni correction is applied to determine the significance of the enrichment or depletion of CpG islands relative to the density seen for RefSeq genes.

Chromatin structure and epigenetics modifications

Data of DNase I hypersensitive (HS) sites, histone acetylation, methylation, transcription binding sites and DNA methylation were obtained from recent publications [35, 39, 43, 48]. The density of each feature is calculated in a 500 bp sliding window advancing 100 bp each time near transcription start sites for HKGs, RefSeq genes, TEGs. The average intensity of all genes in each group is plotted as a function of the distance to transcription start site.