Evidence for large domains of similarly expressed genes in the Drosophila genome
© Spellman and Rubin, licensee BioMed Central Ltd 2002
Received: 28 March 2002
Accepted: 17 May 2002
Published: 18 June 2002
Transcriptional regulation in eukaryotes generally operates at the level of individual genes. Regulation of sets of adjacent genes by mechanisms operating at the level of chromosomal domains has been demonstrated in a number of cases, but the fraction of genes in the genome subject to regulation at this level is unknown.
Drosophila gene-expression profiles that were determined from over 80 experimental conditions using high-density oligonucleotide microarrays were searched for groups of adjacent genes that show similar expression profiles. We found about 200 groups of adjacent and similarly expressed genes, each having between 10 and 30 members; together these groups account for over 20% of assayed genes. Each group covers between 20 and 200 kilobase pairs of genomic sequence, with a mean group size of about 100 kilobase pairs. Groups do not appear to show any correlation with polytene banding patterns or other known chromosomal structures, nor were genes within groups functionally related to one another.
Groups of adjacent and co-regulated genes that are not otherwise functionally related in any obvious way can be identified by expression profiling in Drosophila. The mechanism underlying this phenomenon is not yet known.
The regulation of gene expression is a fundamental process within every cell that often allows exquisite control over a gene's activity (for review see ). Altering transcription rates is an effective strategy for regulating gene activity. It is well established that transcription of a given gene is dependent upon a promoter sequence located within a few hundred base pairs of the transcriptional start site. Promoter activity is modulated by sequence-specific transcription factors that physically interact either with the protein complexes that make up the core transcriptional machinery or with the promoter sequence itself.
In eukaryotes, the activity of a promoter can be modified by transcription factors binding to DNA sequences (frequently termed cis-regulatory modules or enhancers) that are located from hundreds to hundreds of thousands of base pairs away from the promoter. These regulatory modules can either increase or decrease the rate of transcription for a target gene, depending on the cellular state and the activities of the bound transcription factors. There are several mechanisms by which transcription factors bound to regulatory modules exert their effects. First, many transcription factors interact directly with the core transcriptional machinery by recruiting the latter's protein complexes to the promoter. Second, transcription factors may bend or twist the DNA, altering the way in which other transcription factors interact with the DNA. Finally, transcription factors can alter local chromatin structure by modifying histones (typically through methylation, acetylation, and substitution of histone subunits) to permit or restrict access to the DNA. Modifications of chromosome structure also occur at much larger scales. Most eukaryotes exhibit distinct chromosomal regions that are usually either transcriptionally active (euchromatin) or inactive (heterochromatin). In animals, heterochromatin is typically found near centromeres and other regions of low sequence complexity.
Less clear are the mechanisms by which the regulation provided by a cis-regulatory module is restricted to specific target genes. Several examples of insulators - sequences that prevent neighboring modules from affecting transcription - have been identified (reviewed in ). Insulators seem to function not by deactivating cis-regulatory modules but by preventing their influence from being propagated along the chromosome. It is not known how common insulators are in the Drosophila (or any other) genome. Some insulator-binding proteins localize to a few hundred chromosomal positions, and these positions coincide with genomic sequences that are not heavily compacted by chromatin structure (the 'interbands' of polytene chromosomes) . There is substantial evidence that, although gene expression can be tightly controlled, neighboring genes or chromatin regions are important for the expression of individual genes. For example, otherwise identical transgenes inserted into different chromosomal sites show varying levels of expression .
Two recent observations lend credence to the idea that genomes may be divided into domains important for controlling the expression of groups of adjacent genes. First, there is evidence from budding yeast that some genes are found in pairs or triplets of adjacent genes that display similar expression patterns . Second, about 50 much larger regions of the human genome show a strong clustering of highly expressed genes , which is caused by clustering of genes that are expressed in nearly all tissues . We have examined the fraction of genes in the Drosophila genome that are subject to regulation that reflects large domains, using data from high-density oligonucleotide microarrays that reflect over 80 experimental conditions, and have found more than 20% of the genes clustered into co-regulated groups of 10-30 genes.
Many neighboring genes show similar expression patterns
We collected relative gene-expression profiles covering 88 distinct experimental conditions from 267 Affymetrix GeneChip Drosophila Genome Arrays (see Materials and methods section). When the genes in this dataset were organized according to their positions along the chromosome, we observed numerous groups of physically adjacent genes that shared strikingly similar expression profiles. We sought to measure the magnitude of this effect by identifying all groups of physically adjacent genes that showed pair-wise correlations between their expression profiles that were higher than expected by chance.
The number of ten-gene groups of adjacent, similarly expressed genes that are found in ordered and randomized datasets, or are expected to be found in a randomized dataset
Significance (p value)
The number of groups of genes, and total numbers of genes in groups, that are identified at various levels of significance (p values)
Significance (p value)
Gene groups are not explained by gene function or homology
The number of groups of genes, and total numbers of genes in groups, from a dataset containing no physically close homologs
Significance (p value)
We considered an extreme model to account for our observations - that evolutionary selection has organized gene groups according to the biological processes the genes are involved in, so that their expression can be coordinately regulated. We sought to test this model using the Gene Ontology (GO) database [9,10] as a source of annotations of biological processes. We first used the hypergeometric distribution to calculate the probability of observing each GO term as enriched in each group, on the basis of the number of genes in the group, the number of genes in that group that are annotated with that GO term, and the number of genes in the genome that are annotated with that GO term. We then selected all GO 'process' terms associated with a group at p < 0.05 where at least two genes had the selected GO term. Of the 211 groups identified in our full dataset and the 176 groups from the 'homologs-removed' dataset, 43 and 11 GO terms, respectively, have associations to groups that meet the above criteria. These numbers are modestly higher than would be expected by placing a random selection of genes into groups, where we would expect 7 ± 2 from the full dataset and 4 ± 2 from the homologs-removed dataset. The observed enrichment is clearly dependent on homologs, however, given the nearly four-fold decrease in observed associations when homologs are excluded from the analysis. Thus, with the present level of functional annotation, the vast majority of gene groups we observe are not composed of genes with similar biological processes, and the extreme model is not supported.
Similarly expressed gene groups can be identified from smaller datasets
The number of genes within groups identified in either 'adult' or 'embryo' experiments
Significance (p value)
The correlation between sets of genes identified in the adult, embryo and combined datasets
Significance (p value)
Correlations with known chromosome structures
We attempted to determine whether the locations of similarly expressed gene groups correlate with known chromosome structures. Polytene chromosomes show a distinct, reproducible pattern of extended and compacted regions. The compacted regions contain the vast majority of the DNA, although the amount of DNA in each band can vary by more than one order of magnitude. The mean DNA content of each band is approximately 25 kbp [11,12] as compared with approximately 125 kbp for each group of co-expressed genes. We calculated the number of bands that overlap (or are contained in) each group and compared this with the number of bands that overlap (or are contained in) a randomly placed group matched for size. There was very little difference in the average number of bands overlapping each co-expressed group or each randomly placed group (5.9 versus 6.6).
It has been proposed that Drosophila chromosomes are attached to a nuclear scaffold at precise locations , but there is very limited mapping data on the position of these attachments. Mirkovitch et al.  mapped four attachment sites in a 320 kbp region near the rosy gene on chromosome 3R, dividing the region into a number of discrete domains of average size 50 kbp, each containing many genes. We wished to determine whether the groups we identified might correspond to distinct regions between attachment sites, as several of our groups fall in the region studied by Mirkovitch et al. . We attempted to align these regions but there are no clear overlaps; the sizes and positions of the domains identified between attachment sites did not correspond to the groups we found.
We have found that over 20% of the genes in the Drosophila genome appear to fall into groups of 10-30 genes such that the genes within each group are expressed similarly across a wide range of experimental conditions. Our data do not reveal the mechanism(s) responsible for the observed similarities in expression of adjacent genes but we believe the findings are most consistent with regulation at the level of chromatin structure, for the following reasons. First, the regions showing similarities in expression are quite large, containing on average 15 genes, with each gene presumably having its own core promoter. Second, it is frequently the case that one or two genes in a group display a high level of differential expression (see Figure 2c). If the chromatin in a region of the chromosome that contained many genes was 'opened' so that a single target gene could be expressed, it might increase the accessibility of the promoters and enhancers of other genes to the transcriptional machinery, leading to modest parallel increases in their expression. Such an effect could account for the observations we have made.
Discussions of transcriptional regulation often emphasize the belief that the process is tightly controlled and essentially error-free. We believe that the degree of precision, at least at a quantitative level, may be less than is generally assumed. For example, only a few genes show an obvious phenotype when heterozygous, and heterozygosity generally results in a two-fold reduction in expression level . Moreover, there are numerous examples in the literature of genes that, when misexpressed either temporally or spatially, do not generate a phenotype. Although it is difficult to prove that individuals carrying such traits are as fit as their normal relatives, it is likely that the precise regulation of many genes is allowed to vary considerably. If we presume that the groups we have observed arise because of selection on the regulation of a small subset of genes in each group, then the vast majority of genes are in effect being 'carried along for a ride'. The regulation of transcription may be precise when it is needed and sloppy when it is not important.
If coordinated gene expression is unimportant, there should be no selection that drives the groups of co-regulated genes we observed to be evolutionarily conserved. It will be possible to test this when the D. pseudoobscura sequence is completed. If the groups of genes we identify here are found to be more syntenic in the D. melanogaster and D. pseudoobscura genomes than expected, that would support the idea that the observed coordinated expression is advantageous.
Although we have assayed a relatively large number of biological samples, we cannot infer the profiles of unique cellular states. As further experiments are carried out it may be that our observation of similarly regulated groups will grow to include all genes - that is, the entire euchromatic genome may be structured in such domains.
Materials and methods
We collected a dataset composed of 88 experimental conditions hybridized to a total of 267 GeneChip Drosophila Genome Arrays (Affymetrix, Santa Clara, CA, USA) . This dataset came from six independent investigations that will be described in detail elsewhere (A. Bailey, personal communication; M. Brodsky, personal communication; ; E. De Gregorio personal communication; A. Tang, personal communication; and P. Tomancak, personal communication), which study five different experimental questions - aging, DNA-damage response, immune response, resistance to DDT, and embryonic development. Supplemental data including software used in this study and the underlying expression dataset is available at our website  and from the ArrayExpress database  with the accession id E-RUBN-1.
Genes are represented on the GeneChip Drosophila Genome Array by one or more transcripts, which in turn are represented by a probe set. Each probe set has 14 pairs of perfect match (PM) and mismatch (MM) oligonucleotides. Data were collected at the level of the transcript, but for ease in the text, the data are referred to by gene. Intensity data for each feature on the array were calculated from the images generated by the GeneChip scanner, using the GeneChip Microarray Suite. These intensity data were loaded into a MySQL database where information on each of the features was also stored. The difference between the PM and MM oligonucleotides (probe pair) was calculated, and the mean PM-MM intensity for each array was set to a constant value by linearly scaling array values. The mean intensity of individual probe pairs was calculated across all arrays, and the log2 ratio of each value to this mean was stored. Next, all log ratios for each probe pair set (transcript) were averaged, creating one measurement for each transcript on each array. The final dataset was generated by averaging data for each transcript on replicate arrays and subtracting the average log ratio of each gene in the dataset.
Definition of homologs
BLAST scores based on predicted protein sequence were obtained from Gadfly (Release 2) . We used these scores to define a homolog pair as those gene pairs for which BLAST E values are less than 10-7.
Identification of adjacent similarly regulated genes
We calculated the average pair-wise correlation of gene-expression profiles for all genes that were within n genes (an n-gene window) of one another using the Pearson correlation. Significance (p values) was estimated by sampling random sets of n genes 1 million times to determine the likelihood distribution for the dataset. We also calculated the average pair-wise correlation for a random dataset in which the associations between genes and expression profiles were shuffled. We have calculated the number of genes in groups at each of the three p values, namely 10-2, 10-3, and 10-4, for window sizes ranging from 2 to 25 genes.
Next, we set out to show that homologs did not account for the increase in the number of gene groups with higher than expected average correlations. We searched for cases in which homologs (as defined above) were near each other in the genome by scanning the set of genes for each chromosome from one end to another. If a gene showed homology to another gene that appeared less than 10 genes ahead, it was removed from the dataset, although no break in gene order was created. For example, in a set of 11 genes where the third and fourth were homologs, gene 3 would be removed, and a ten-gene group would consist of genes 1, 2 and 4 through 11. In total, 1,369 genes were removed from the dataset. This 'homologs removed' dataset was subjected to the average pair-wise algorithm, as was a randomized version of it.
We also constructed two non-overlapping subsets of the total data matrix. All hybridizations were divided into either the 'embryo' or 'adult' dataset on the basis of the source of the RNA used in that hybridization. In total, 35 experiments remained in the embryo dataset and 53 experiments remained in the adult dataset. The random pair-wise correlation algorithm was applied independently to each of these datasets as well as to randomized versions of each dataset.
Significantly enriched GO terms among gene groups
GO terms for all genes were obtained from the GO database . Using the hypergeometric distribution, the probability of observing each GO term with each group was calculated. Briefly, the probability p that a GO term is significantly enriched among a specified set of genes can be calculated with the following formula:
where k is the number of genes in the group, G is the total number of genes, n is the number of genes in the group with a given annotation and A is the total number of genes with a given annotation. Because many sets of GO terms (> 1,000) were tested on many groups of genes (> 200), there is a problem of multiple testing. All GO terms significantly associated with a group of similarly expressed genes at a p value of less than 5 × 10-4 were recorded.
Correlation of groups with known chromosomal structures
We determined the number of polytene bands present in each group of similarly expressed genes. The coordinates of each group were determined by using the transcription start sites (from GadFly Release 2)  of the genes at each end of a group. We then determined how many bands overlapped each group based on the positions reported [11,12]. We also calculated the number of bands that overlap randomly placed groups (with the same sizes as the real groups).
Additional data files
The following are provided as supplemental materials; a tab-delimited text file of the underlying expression data; the perl scripts used to process the data; and a text file used to generate Figure 2. All expression data are reported as log base 2 and are mean centered (the mean expression value for each gene in all experiments is zero). The first column of each expression data file is the CT identifier of each transcript. The second column is a description field, which includes the CT identifier, CG identifier, gene name, and brief Gene Ontology annotations. The remainder of the columns contain expression data, classified by the column header (either adult or embryo). The data used to generate Figure 2 can be loaded into the TreeView software  to visualize individual groups (null data rows indicate boundaries between groups). The software and underlying expression dataset are also available at our website  and from the ArrayExpress database  with the accession ID E-RUBN-1.
We thank Adina Bailey, Michael Brodsky, Amy Tang, and Pavel Tomancak for sharing data prior to publication. P.T.S. was a recipient of an NSF Biocomputing postdoctoral fellowship. G.M.R. is an investigator of the Howard Hughes Medical Institute.
- Emerson BM: Specifity of gene regulation. Cell. 2002, 109: 267-270.View ArticlePubMedGoogle Scholar
- Bell AC, West AG, Felsenfeld G: Insulators and boundaries: versatile regulatory elements in the eukaryotic genome. Science. 2001, 291: 447-450. 10.1126/science.291.5503.447.View ArticlePubMedGoogle Scholar
- Zhao K, Hart CM, Laemmli UK: Visualization of chromosomal domains with boundary element-associated factor BEAF-32. Cell. 1995, 81: 879-889.View ArticlePubMedGoogle Scholar
- Spradling AC, Rubin GM: The effect of chromosomal position on the expression of the Drosophila xanthine dehydrogenase gene. Cell. 1983, 34: 47-57.View ArticlePubMedGoogle Scholar
- Cohen BA, Mitra RD, Hughes JD, Church GM: A computational analysis of whole-genome expression data reveals chromosomal domains of gene expression. Nat Genet. 2000, 26: 183-186. 10.1038/79896.View ArticlePubMedGoogle Scholar
- Caron H, van Schaik B, van der Mee M, Baas F, Riggins G, van Sluis P, Hermus M-C, van Asperen R, Boon K, Vouöte PA, et al: The human transcriptome map: clustering of highly expressed genes in chromosomal domains. Science. 2001, 291: 1289-1292. 10.1126/science.1056794.View ArticlePubMedGoogle Scholar
- Lercher MJ, Urrutia AO, Hurst LD: Clustering of housekeeping genes provides a unified model of gene order in the human genome. Nat Genet. 2002, 31: 180-183. 10.1038/ng887.View ArticlePubMedGoogle Scholar
- Eisen MB, Spellman PT, Brown PO, Botstein D: Cluster analysis and display of genome-wide expression patterns. Proc Natl Acad Sci USA. 1998, 95: 14863-14868. 10.1073/pnas.95.25.14863.PubMed CentralView ArticlePubMedGoogle Scholar
- 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: 25-29. 10.1038/75556.PubMed CentralView ArticlePubMedGoogle Scholar
- The Gene Ontology Consortium. [http://www.geneontology.org/]
- Ashburner M, de Grey A: Cytological table used to infer a genetic map position from a published cytogenetic map position. [http://fly.ebi.ac.uk:7081/maps/lk/cytotable.txt]
- Saura AO, Saura AJ, Sorsa V: Electron micrograph maps of Drosophila melanogaster polytene chromosomes. [http://www.helsinki.fi/~saura/EM/index.html]
- Mirkovitch J, Spierer P, Laemmli UK: Genes and loops in 320,000 base-pairs of the Drosophila melanogaster chromosome. J Mol Biol. 1986, 190: 255-258.View ArticlePubMedGoogle Scholar
- Lindsley DL, Sandler L, Baker BS, Carpenter AT, Denell RE, Hall JC, Jacobs PA, Miklos GL, Davis BK, Gethmann RC, et al: Segmental aneuploidy and the genetic gross structure of the Drosophila genome. Genetics. 1972, 71: 157-184.PubMed CentralPubMedGoogle Scholar
- Affymetrix GeneChip Drosophila Genome Array. [http://www.affymetrix.com/products/arrays/specific/fly.affx]
- De Gregorio E, Spellman PT, Rubin GM, Lemaitre B: Genome-wide analysis of the Drosophila immune response by using oligonucleotide microarrays. Proc Natl Acad Sci USA. 2001, 98: 12590-12595. 10.1073/pnas.221458698.PubMed CentralView ArticlePubMedGoogle Scholar
- Spellman PT, Rubin GM: Web supplement to "Identification of adjacent gene groups showing similar expression". [http://www.fruitfly.org/expression/dse/]
- ArrayExpress. [http://www.ebi.ac.uk/microarray/ArrayExpress/arrayexpress.html]
- Gadfly. [http://www.fruitfly.org/annot/index.html]