Reverse transcriptional profiling: non-correspondence of transcript level variation and proximal promoter polymorphism

Transcriptional profiling via whole genome arrays has revealed both impressive diversity and unexpected commonalities in gene expression. Some fraction of this variation results in the various functional, developmental, and reproductive phenotypes that are downstream of gene expression [1, 2]. The source of this variation, however, is in the nucleotides that regulate gene expression (as reviewed in [3]). Given the relationship between these nucleotides and gene expression, it should in principle be possible to reverse this relationship and retrodict the nucleotide sequences that have given rise to variation in gene expression. Indeed, this is the basis for numerous bioinformatic efforts to identify novel cis-regulatory motifs if not to elucidate the entire cis-regulatory code [4, 5] from sequence that is conserved among genes with similar expression patterns. Here we ask a somewhat different but equally fundamental question of sequence that varies among genes with dissimilar expression: Are regulatory regions different for genes differing in expression, and similar for genes not differing in expression?

We ask this question of cis-regulatory sequence immediately upstream of the transcription start site and of similarities or differences in gene expression between members of the same species. Candidate cis-regulatory variation includes polymorphisms in transcription factor binding sites, multiallelic tandem repeat variation, such as microsatellites, single nucleotide polymorphisms (SNPs), and transposable element insertions. These polymorphisms can alter transcription rates and affect mRNA expression levels in vivo [6–10]. Functional cis-regulatory variation is widespread in the human genome [11] and these sites have a high degree of variability (about 0.6% of these sites are polymorphic [12]). Clearly the source of variation in gene expression need not lie in the proximal promoter and could be either in cis but far up- or downstream, or in trans. Lack of correlated variation in proximal promoter sequence and gene expression can be viewed as implicating these alternative sources.

We focus on variation in gene expression among members of the same species because this is the raw material of evolution. For evolutionary change to ensue by ordinary Darwinian mechanisms, variation within a species must exist, be heritable, and be consequential for fitness. Both classical [13, 14] and recent [15–17] and others reviewed in [18] works implicate variation in gene expression among members of a population, controlled by regulatory as opposed to coding sequence, as a principal source of such variation. Although simple sequence analysis of regulatory regions from diverse and related members of a species can readily ascertain whether sequence variation is present and heritable, ascertainment of impact on gene expression from sequence alone is less precise. Although differences in gene expression are not necessarily consequential for fitness, if variation in regulatory sequence affects gene expression (or that similarity in sequence reliably corresponds to no variation in gene expression), then intraspecific variability in regulatory sequence could be used as a preliminary screen for evolvability. Indeed, sequence analysis of regulatory regions may potentially assess the evolutionary forces that shape them [15].

To test this possibility, we mined the whole-genome transcriptional profiles of two near-isogenic strains (Oregon R and Russian 2b) of Drosophila melanogaster, here representing variation within species. The original data [1] are for genes whose expression does or does not vary between non-virgin adults of the Oregon R and Russian 2b strains. Of the 12017 genes examined, 527 do not differ between strains and 483 differ between strains regardless of sex. With no a priori knowledge of these genes' sequence variation or function, we have selected 15 and 19 representatives of each class, respectively. We chose genes whose regulation is yet uncharacterized to avoid bias due to over- or under-sampling of genes with known transcriptional regulation, and included both known transcription response elements/conserved motifs and uncharacterized sequence in our analysis. We now report whether or not their proximal promoter sequence varies and, if so, whether or not this variation corresponds to variation in gene expression.

Many genes in the first geneset are false positives. The expected number of candidate genes in the first geneset, 325, is less than 483, the observed number resulting in a False Discovery Rate (FDR) of 67.4%. Thus, of the 483 genes chosen in the first gene set, 67.4% or ~325 genes are called as significantly different when they are not. In other words, the expression of ~158 of the 483 genes in the first geneset presumably differs significantly between the two strains.

Polymorphisms are numerous in both genesets. Of the 34 promoters analyzed, six do not differ between strains; however, they belong to genes whose transcripts differ between strains. The remaining 28 proximal promoters contain at least one polymorphism, with one promoter containing as many as 37 polymorphisms (Figures 1 and 2). The mean, median, and mode of polymorphisms per gene are 8.5, 6, and 0, respectively. Of the 288 total polymorphisms detected in at least 1 kb of the proximal promoters, 258 (89.6%) were SNPs and 30 (10.4%) were indels. Over half of the SNPs were transitions within a nucleotide class, while 43.8% were transversions between nucleotide classes (Figure 3). 239 putative binding sites for transcription factors were created or removed by these proximal promoter polymorphisms; thus 258/288 or 90% of proximal promoter polymorphisms fall within putative transcription factor binding sites (Figures 1 and 2).

Although few, indels varied between the strains by kind and from 1 to 43 nt in length [see Additional File 1]. Indels were classified as direct repeats (dr), homopolymer repeats (hpr), microsatellites (mcs), or non-repetitive (nr) according to designations of Schaeffer [19], page 165. In most cases, indels in one strain were the same as in the Celera strain used as a reference. However, six repeats in five different genes resulted in differences in sequence among the Russian 2b, Oregon R, and Celera strains. In bin (CG18647), five additional Ts are found in the Oregon R strain in comparison with the poly-T₄ tract in the Russian 2b strain. This caused a transition and insertion from T₄CT₃ in the Celera strain. Distal to this in bin, the ATACCCGTACCCGTACCCAT sequence in the Russian 2b strain was shortened to ATACCCGTACCCAT in the Celera strain but absent altogether in the Oregon R strain. In tacc (CG9765), the poly-A tract varies from A_14–17 nt among individuals in the Russian 2b strain, to A_23–26 nt among individuals in the Oregon R strain, to A₂₂ in the Celera strain. In qkr58E-3 (CG3584), a SNP and variation in the length of a poly-T tract resulted in T₃GT₄ in the Oregon R strain, T₉ in the Celera strain and T₁₁ in the Russian 2b strain. In KP78a (CG6715), the dinucleotide microsatellite AC was repeated 9X in the Russian 2b strain, 10X in the Celera strain, and 11X in the Oregon R strain. In stan (CG11895), the homopolymer repeat was T₁₀ nt long in the Celera strain, T₁₁ in the Russian 2b strain, and T₉ in the Oregon R strain.

Indels in Scab, Cry, Ih, and bin (CG8095, CG16963, CG8585, CG18647) contained small regions (12–15 nt long) that shared sequence similarity with known transposable elements. Despite these matches, the sequences are not long enough to discriminate confidently between a TE footprint and chance occurrence of the same sequence. Flanking the indel in Cry listed in the Supplementary Table [see additional file 1], a 152 bp sequence from -889 to -1041 relative to the translational start site in the Oregon R strain matched a DNA LINE retroelement. In the same location in the Russian 2b strain a 152 bp sequence matched two overlapping DNAREP1_DM LINE elements [20].

The variation in proximal promoter sequence among the Oregon R, Russian 2b, and Celera strains, while extensive, typifies that among individuals and populations of eukaryotes. In a similar survey of 107 transcriptionally active genes in the human genome, Rockman and Wray [11] identified 140 experimentally validated cis-regulatory polymorphisms resulting in two-fold or greater variation in transcription rate and subsequent gene expression. In another survey of regulatory variation, Cowles et al. [21] found a relatively high frequency (6%) of cis-regulatory polymorphism (including SNPs, complex nucleotide repeats, insertions/deletions, and microsatellites) in the 1 kb region 5' of the predicted transcription start site of genes in four inbred mouse strains with allelic differences in expression by 1.5-fold or greater.

The frequencies of SNPs and indels were not independent. In fact, indels, including microsatellites, homopolymer repeats, and possibly TE footprints (see next paragraph), co-occurred with SNPs but SNPs could occur in the absence of indels. This suggests that indels occur infrequently and only after a promoter is poised for mutation via the presence of SNPs [22, 23].

Some promoters have diverged remarkably between the two strains. For example, the Cry promoter contains 28 SNPs, 7 indels, and a 152 bp sequence including and flanking one indel that matches two different classes of transposable elements a DNA LINE element in the Oregon R strain and 2 overlapping DNAREP1_DM LINE elements in the Russian 2b strain [20]. This variation suggests that the Cry promoter is susceptible to mutation.

The polymorphisms are in regions within the proximal promoter that ought to affect gene expression. 90% of them fall within putative transcription factor binding sites (Figures 1 and 2). In addition, all regulatory information necessary for transcription in Drosophila is generally present within 1 kb of the basal promoter [24].

Surprisingly, genes differing in expression and genes not differing in expression had the same diversity and frequency of polymorphisms in their proximal promoters. Indeed, the proximal promoters from six genes whose expression differed were identical in sequence between the two strains. Because this region includes core regulatory sequences, polymorphisms between the two genesets (or their absence) ought reasonably be correlated with gene expression. In support of this expectation, 74% of transcription factor binding sites (as proxies for proximal promoter polymorphisms) in yeast genes are between 100 to 500 bases upstream of the translational start site and fewer than expected lie outside this region [25]. Thus, the non-concordance between proximal promoter sequence and gene expression is unexpected.

One possible explanation for this outcome is that the 34 genes examined are anomalous or unrepresentative of Drosophila genes in general. This explanation is unlikely. As discussed above, these genes' proximal promoters resemble other genes' in variation among individuals. Choice of these genes was also unbiased with respect to prior knowledge of transcriptional regulation, known transcriptional response elements, and mode of gene effect (e.g., dominance, underdominance, or additivity). Indeed, 5/34 (15%) of the blindly chosen genes were dominant, underdominant, or additive in their impact on gene expression, which compares favorably with the proportions of such effects for the Drosophila genome in general [1].

A second explanation for non-concordance between proximal promoter sequence and gene expression is that the expression of the genes in the first geneset actually does not differ between strains. To exclude this possibility, the first geneset included only genes whose expression profiles were different between strains (P < 0.01), and not between sexes or by the sex*strain interaction (P > 0.05 for both) [1]. These selection criteria confer a rather high False Discovery Rate of 67.4%. The actual statistical significance of variation between the two strains averages <0.0009 for the 19 genes in the first geneset. An additional 333 genes also meet this criterion. For these 352 genes (total), the FDR becomes ((0.0009)*(0.95)*(0.95)*(12017)*(3))/352 = 29.28/352 = 0.0832 or 8.3%. This low rate indicates 91.7% certainty that the 19 genes in this geneset truly differ in expression.

A third explanation is that some of these promoter variants may affect gene expression but that neutral sites dilute their regulatory activities. Functionally important regulatory sites may make up only a small fraction of the many promoter variants detected in the 1 kb region and the vast majority may be neutral. On the other hand, promoter variants whose expression levels are similar may be functionally equivalent (or neutral) as a result of stabilizing selection or epigenetic forces may even out promoter variation resulting in similar expression levels. Functional assays should be able to distinguish between these outcomes.

Another possibility is that the regulatory elements responsible for the observed variation in transcript abundance occur outside the proximal promoter, including in introns and 5' and 3' untranslated regions, and/or in trans. Several recent localizations of regulatory sequence also support this possibility. In a recent assignment of 142 expression phenotypes in humans to both cis- and trans-acting loci, 27 (19%) had a single cis-acting regulator (defined as regulators that map within 5 megabases of the target gene), 110 (77.5%) had a single trans-acting regulator (those that map elsewhere), and 5 (3.5%) had two or more regulators acting in cis and/or trans [26]. In budding yeast, trans- acting regulators were linked to 365 of 570 (64%) genes whose expression diverged between a laboratory and a wild strain of Saccharomyces cerevisiae whereas only 36% of divergent gene expression was linked to individual cis-acting loci [27]. Of 2294 genes whose expression differs between laboratory and wild yeast strains, 1716 (75%) map to 100–200 distinct trans-acting loci with widespread genomic effects and molecular functions but that do not encode transcription factors [28]. Of 3546 highly heritable transcripts in the same two strains of yeast, 3% map to a single locus, 17–18% to 1–2 loci, and half to more than 5 loci [29]. Using gene expression differences as quantitative traits, Schadt et al. [30] determined that genes involved in determining patterns of obesity in mice and humans map to QTLs other than the genes themselves. Many genes map to the same QTL, suggesting these are trans-regulatory hotspots.

This study leads us to suggest that the relationship of coding sequence conservation and functional similarity may not be true for cis-regulatory sequences. Indeed, recent work as shown this to be the case for enhancers [31]. The problem, however, is that there is not a precise quantitative framework for the interpretation of cis-regulatory variation.

In summary, we have characterized the frequency and diversity of proximal promoter polymorphisms on a genome-wide scale and shown that they do not differ between D. melanogaster strains with divergent gene expression. We conclude that sequences elsewhere, such as in trans, may cause differential gene expression. Indeed, in the very strains examined in the present study, extensive nonadditivity of gene expression among these strains and their reciprocal F₁ hybrids indicates that transcription is controlled predominantly by trans-regulatory factors [1]. Linkage analysis might allow for mapping of regulatory regions outside the proximal promoter that cause expression variation within the D. melanogaster transcriptome.