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

Thermal evolution of gene expression profiles in Drosophila subobscura

Abstract

Background

Despite its pervasiveness, the genetic basis of adaptation resulting in variation directly or indirectly related to temperature (climatic) gradients is poorly understood. By using 3-fold replicated laboratory thermal stocks covering much of the physiologically tolerable temperature range for the temperate (i.e., cold tolerant) species Drosophila subobscura we have assessed whole-genome transcriptional responses after three years of thermal adaptation, when the populations had already diverged for inversion frequencies, pre-adult life history components, and morphological traits. Total mRNA from each population was compared to a reference pool mRNA in a standard, highly replicated two-colour competitive hybridization experiment using cDNA microarrays.

Results

A total of 306 (6.6%) cDNA clones were identified as 'differentially expressed' (following a false discovery rate correction) after contrasting the two furthest apart thermal selection regimes (i.e., 13°C vs . 22°C), also including four previously reported candidate genes for thermotolerance in Drosophila (Hsp26, Hsp68, Fst, and Treh). On the other hand, correlated patterns of gene expression were similar in cold- and warm-adapted populations. Analysis of functional categories defined by the Gene Ontology project point to an overrepresentation of genes involved in carbohydrate metabolism, nucleic acids metabolism and regulation of transcription among other categories. Although the location of differently expressed genes was approximately at random with respect to chromosomes, a physical mapping of 88 probes to the polytene chromosomes of D. subobscura has shown that a larger than expected number mapped inside inverted chromosomal segments.

Conclusion

Our data suggest that a sizeable number of genes appear to be involved in thermal adaptation in Drosophila, with a substantial fraction implicated in metabolism. This apparently illustrates the formidable challenge to understanding the adaptive evolution of complex trait variation. Furthermore, some clustering of genes within inverted chromosomal sections was detected. Disentangling the effects of inversions will be obviously required in any future approach if we want to identify the relevant candidate genes.

Background

Temperature is a fundamental feature that affects all living organisms. Each species, particularly ectotherms, has a non-stressful thermal tolerance range and responds to temperature by physiological, biochemical, and molecular level adjustments that underlie adaptation. For instance, many latitudinal clines exist in Drosophila for allele frequencies at allozyme loci, chromosomal inversions, and microsatellites; as well as for traits like starvation resistance, desiccation resistance, and body size where the differences between populations have a genetic basis and can even persist for many generations under laboratory reared conditions [1–4]. By and large, the empirical evidence suggests that variation in these markers and traits are directly or indirectly related to temperature (climatic) gradients. Perhaps the most pertinent example comes from recent studies on chromosomal inversion polymorphisms showing that the genetic constitution of populations is responding to contemporary rapid global warming [5–8].

The native Palearctic fly Drosophila subobscura spans more than 30° latitude in the Old World: from North Africa to Scandinavia. As a result, its populations experience a strong climatic gradient [9]. In the late 1970s and early 1980s the species invasively spread in North and South America and nowadays spans about 15° latitude on each continent [10, 11]. Remarkably enough, latitudinal clines in the New World for inversion polymorphism and body size parallel to the long standing ones in the native geographic area were evident after a few years since the American colonization. This 'replicated time series experiment of evolution in action' [12, 13] strongly suggests that those traits are indeed subject to selection from temperature-related factors. However, a laboratory natural selection experiment (i.e., an experimental protocol where stocks of organisms are reared under different conditions and allowed to evolve by natural selection over many generations [14]) specifically designed to test the putative role of temperature per se in the evolution of these clines found results conflicting to those expected from the latitudinal gradients for both inversions and body size [3, 15, 16]. Certainly, laboratory experiments are not the best way to mirror what is happening at different latitudes and to reconstruct natural clines. But at present it is unclear whether temperature alone drives the clines. What types of genetic changes are needed for an organism to adapt to new thermal conditions?

A number of authors (e.g. [17, 18]) have argued that changes in the transcriptome constitute a major component of the genetic basis for phenotypic evolution. Gene expression profiling by means of microarrays has become a popular way of finding candidate genes of trait variation and is providing new insights into some old but fundamental questions in evolutionary biology [19–22]. Here we examine global gene expression by measuring the relative abundance of mRNAs in third instar larvae of D. subobscura from 3-fold replicated laboratory thermal selection stocks -derived from the estimated Chilean epicentre (Puerto Montt) of the original New World invasion [15]- that had evolved at three constant temperature regimes during 3 years: cold (13°C), optimum (18°C) and warm (22°C). The connection between the very high dimensional nature of the gene-expression data and the multivariate whole organism phenotype, however, is not straightforward and detailed functional and ecological analyses of candidate genes will obviously be required to understand the genetic basis for thermal adaptation (e.g., [23]).

In holometabolous insects like Drosophila many adaptations to changing environments involve changes in larval behaviour and physiology that may impinge on other phases of the life cycle (e.g. [24–26]). For this reason we used third instar larvae for the microarray experiment. A possible drawback was that both sexes were mixed, so we have overlooked any sex-specific thermal response that might have been present. Total mRNA from each population was compared to a reference pool mRNA independently derived from the optimum (P18) populations (see Methods) in a standard two-colour competitive hybridization using cDNA microarrays with D. melanogaster clones [27]. Heterologous hybridization to study gene expression profiles has been validated between closely related species (divergence time < 10 Mya), and consistent data are also obtained for less closely related taxa (divergence time ~65 Mya; [28]). Drosophila subobscura belongs to the D. obscura group, and the divergence time between the D. melanogaster and D. obscura groups has been estimated to be ~25 Mya [29]. This apparently offers a reasonable warranty to use D. melanogaster arrays in heterologous hybridization with D. subobscura. In addition, we carried out some preliminary tests in order to optimize the experimental conditions. Highly reproducible and consistent gene profiling, comparable to that obtained with homologous hybridization by using some D. subobscura clones added to the arrays (see Methods), was observed.

It is also important to remark here that the lower and upper thermal regimes used in the experiment are not stressful: the temperature range likely covers much of the physiologically tolerable range in this species [9]. Obviously, the constant thermal regimes and light:dark period where the populations have evolved (see Methods) do not mirror the seasonal changes experience by natural populations, but with this experimental protocol we can control that temperature is the only factor differing between the thermal stocks. The thermal stocks had already diverged for inversion frequencies, pre-adult life history components, and morphological traits [3, 16]. The experimental design equated to a four-way factorial analysis of variance (ANOVA) with thermal selection regime and cyanine dyes (Cy3, Cy5) in a flip dye design as fixed effects, replicated populations as a random factor nested in thermal selection regime, and slide (spotted microarray) as a random factor nested in thermal selection, dye, and replicate. The analysis allowed identifying quantitative differences in larval gene expression between cold- (P13) and warm-adapted (P22) populations.

The candidate genes were assigned a biological function and/or biological process when information was available. Also important, a number of genes were mapped by in situ hybridization to the polytene chromosomes of D. subobscura. The karyotype of this species consists of five acrocentric chromosomes and a dot chromosome (see Methods). What is crucial here is that D. subobscura harbours one of the richest inversion polymorphisms in the genus, with a total of 92 chromosomal arrangements (produced from 66 inversions located on all major chromosomes) recorded in the native Palearctic region [9, 30]. This number reduces to 18 arrangements in colonizing populations of the New World, all of them but one segregating in the thermal stocks [16]. Variation in some traits is known to be tied to inversion polymorphisms in Drosophila (e.g. [31, 32]), and quantitative associations between larval gene expression and thermal adaptation could be due to position effects (e.g., the inversion of a chromosomal segment can remove or exchange the regulatory sequences of a gene and alter its expression pattern [33]) or hitchhiking arisen from linkage disequilibrium. In view of the rapidly and consistently evolved latitudinal clines in chromosome inversion polymorphism following the New World invasion by the species [12], and the shifts in inversion frequencies in response to laboratory thermal adaptation [16] and to climate change [7], we expect that a large number of genes will be included inside inverted chromosome segments. Linkage with inversions will highly complicate the identification of chromosome regions that are targets of selection.

Results and discussion

Overall patterns of gene expression in the thermal lines

An important point in the experiment was that the parents of treatment larvae had also been reared at the same temperature of 18°C to control for the possibility of non-genetic parental effects on offspring (see Methods). In order to generate a reliable data set we analyzed mRNA abundance from a highly replicated experiment: four independent batches of 250 optimum-reared larvae each -amounting to 9,000 larvae in total (i.e., 250 larvae × 4 slides per population × 9 experimental populations)- whose mRNAs were competitively hybridized to a reference pooled mRNA from 9,000 control larvae on contiguous duplicated gene spots using a dye-reversal experimental design, thus providing up to 72 gene expression values for each probe. Furthermore, some genes were spotted (in duplicated) several times on the slides, which helped to confirm the quality and consistency of the data as there was a clear correspondence among different spots [see Additional file 1: summary of the microarray results].

The distribution of normalized measures for valid relative gene expression levels in the g = 11,767 cDNA clones is shown in Fig. 1. All clones with less than 57 valid expression observations were excluded from further analyses and, therefore, we will only focus on the right part of Fig. 1 just after the x-axis scale break. As a result, 4,651 cDNA clones (4,310 non-redundant genes) were allowed for differential gene expression analysis and their microarray signature, plotted as expression ratio versus fluorescence intensity, is shown in Fig. 2a. For each clone g = 1, ⋯, 4,651 the dye-reversal experimental design was subjected to a least squares ANOVA model as that shown in Table 1 for (e.g.) gene CG12236. A key premise in the experimental design with 2 degrees of freedom for the main fixed factor temperature was to define a priori the linear contrast between the two furthest apart thermal selection regimes: warm (22°C) vs. cold-adapted (13°C) populations (each comparison or contrast between two means has one degree of freedom).

Figure 1
figure 1

Distribution of valid gene expression levels. For a given probe g the maximum number of valid gene expression values for each thermal selection regime was N = 72. All probes with N < 57 (left part just before the x-axis scale break) were excluded from the statistical analysis.

Figure 2
figure 2

Microarray signature. Gene expression profiles in the three thermal regimes. (a) Expression ratio versus fluorescence intensity of the 316,229 spots (black dots) from the 4,651 clones approved for statistical analyses. The averages for those probes identified as differentially expressed when contrasting the two extreme thermal regimes (i.e.,13°C vs. 22°C) are in blue (P13), green (P18) and red (P22). (b) Box-plots of the average expression ratios for the differentially expressed genes.

Table 1 ANOVA of relative intensity ratios for gene CG12236.

A total of 419 (9%) cDNA clones were identified as 'differentially expressed' when considering a p-value (from 10,000 rounds of permutation) cut-off of 5% for the temperature factor with 2 degrees of freedom (see Table 1), but none of them was labelled as truly significant in terms of the false discovery rate (FDR; [34]) method used in detecting differential gene expression (q-value threshold of 5%; see Methods). On the other hand, from the permutation p-values obtained after the linear contrasts between the two furthest apart thermal selection regimes the number of 'differentially expressed genes' rose up to 950 (20.4%), with 306 (6.6%) remaining significant after a FDR correction (recall that a q-value threshold of 5% means that among all genes considered as significant, 5% of these are truly null on average [35]). Fig. 2 also shows the averages for the expression ratio versus fluorescence intensity of the identified 306 genes differing in gene expression (Fig. 2a), together with the corresponding box-plots (Fig. 2b). The reason why the linear contrasts comparing P22 vs. P13 populations yielded more differentially expressed genes was because the averages of log 2 relative intensity ratios for optimum (P18) populations were normally in between the averages for P13 and P22 populations. Therefore, a substantial proportion of the sum of squares for the temperature factor with 2 degrees of freedom in the ANOVAs was accounted for by the linear contrast (e.g., ~92% in Table 1).

A far more informative plot is shown in Fig. 3, where gene expression values are sorted according to the contrast estimates. It is clear that those populations that have evolved at the optimum thermal regime (18°C) ranged in between, and that the average gene expression difference between warm- and cold-adapted populations was generally low: from 0.6 for CG4183 (heat shock protein 26; Hsp26) to 1.7 for CG4867 (bc10).

Figure 3
figure 3

Transcriptome changes following thermal adaptation. Scatterplots of the average log2 relative intensity ratio for the 306 'differentially expressed genes' detected from the linear contrasts between the two furthest apart thermal selection regimes (i.e., 13°C vs. 22°C) and sorted according to the difference between warm- and cold-adapted (i.e., P22 – P13) populations. The sorted CG names of genes are given at the bottom (1, 3, 5, ⋯) and top (2, 4, 6, ⋯) axes. For each gene, the points in blue (P13), green (P18) and red (P22) give the average log2 relative intensity ratio for the different thermal regimes. The corresponding points have been connected by polynomial fitting to enhance visibility.

Gene expression grouped by cellular function and biological process

Analysis of functional categories defined by the Gene Ontology project [36] using the GOToolBox [37] revealed that our reference dataset (the 4,651 cDNA clones with 4,310 non-redundant genes that were allowed for differential gene expression analysis) includes 989 annotated genes, and only 66 out of 306 genes labelled as 'differentially expressed' were annotated. These genes could be assigned to different cellular or molecular functions: over two-thirds are involved in metabolism processes (41 genes), in transport processes (14 genes), and in regulation of transcription (8 genes). (Note that many genes belong to more than one category.) For each functional category, we compared the actual number of occurrences with the expected one under the null hypothesis that all categories should be equally represented. Namely, the probability of obtaining by chance a number n of annotated genes for a given term among a dataset of size N, knowing that the reference dataset contains m such annotated genes out of G genes, is calculated. This test follows the hypergeometric distribution and the GOToolBox allows for FDR correction, pointing at statistically relevant over- or underrepresented terms within a dataset. The results obtained are shown in Fig. 4 and indicate an overrepresentation of genes involved in carbohydrate metabolism, nucleic acids metabolism and regulation of transcription among other categories. Two categories are apparently underrepresented: organic acid and carboxylic acid metabolism [see Additional file 2: molecular function gene ontology categories of differentially expressed genes].

Figure 4
figure 4

Genes grouped by gene ontology (GO) terms. For all statistical significance categories, the bars in the reference set plot the subset frequencies of annotated genes in the cDNA microarrays belonging to a particular class with valid gene expression data for statistical analysis (Fig. 1), and the bars in the dataset represent those frequencies in the group of genes labeled as 'differentially expressed' when contrasting the two extreme thermal regimes. Two categories are apparently underrepresented (organic acid and carboxylic acid metabolism), whereas the rest are overrepresented in the dataset [see Additional file 2].

It seems reasonable that genetic adjustments to environmental differences may involve changes in metabolism. Thus, when inbred and non-inbred D. melanogaster lines are reared under benign and stressful (high temperature) environmental conditions gene expression patterns of metabolic genes are strongly affected by both inbreeding and temperature stress [38]. Furthermore, previous studies on thermal evolution using the same Drosophila species have shown that differences between cold- and warm-adapted populations can be due to differences in the efficiency of larval growth [24, 39]. In our thermal stocks with D. subobscura we have shown that cold-adapted (P13) populations had longer development times in the whole range of developmental temperatures assayed, and that warm-adapted (P22) populations seem to have evolved faster development [3]. Together with the lack of divergence for adult body size, it seems that cold-adapted D. subobscura stocks achieve the 'target' size by growing more slowly. This apparently agrees with their lower level of gene expression for genes involved in metabolic processes when compared to their warm-adapted counterparts.

Candidate genes for thermotolerance in Drosophila

Expression levels for genes of the heat shock protein group (Hsps) that act as molecular chaperones and are important for cellular housekeeping are known to covary with the thermal regimes experience by populations, species and higher taxa [40]. Hsp70 appears to be the primary protein involved in thermotolerance in D. melanogaster [41] -though apparently not in other Drosophila species [42]-, and Hsp70 allele frequencies show latitudinal clines and change in response to thermal evolution in the laboratory [43]. In addition, Hsp23 and Hsp26 latitudinal variation in the D. melanogaster Australian cline [44], and correlated responses to selection for knockdown resistance at 39°C for Hsp68 [45], have also been found. Besides Hsps, other candidate genes for adaptation to thermal extremes (summarized in [46]) are: Hsrω (heat-shock RNA ω, which produces two RNA products but no known protein product [47]), Hsf (heat-shock transcription factor), Tot (Turandot), mth (methuselah, also a candidate aging gene [48, 49]), Dca (Drosophila cold acclimation gene), Fst (frost; involved in recovering from cold shock [50]), Drs (drosomycin), shark (involved in a signaling pathway for epithelial cell polarity [51]), anon-23Da (encoding a protein with currently unknown function), desat2 (desaturase2; [52]), period (clock gene that determines biological rhythmicity in Drosophila [53]), Ddc (dopa decarboxylase; involved in the catecholamine biosynthesis pathway, which has been implicated in the response to various stressors including temperature [54]), and various metabolic enzymes as Adh (alcoholdehydrogenase; e.g., [55, 56]), Gpdh (Glycerol 3 phosphate dehydrogenase; [56]), Gdh (NAD-dependent glutamate dehydrogenase; [57]) and Treh (trehalase; [57]).

We had valid gene expression values (i.e., N ≥ 57 in Fig. 1) for all but 6 of the formerly listed genes; namely, Hsp70, Hsp23, Hsp26, Hsp68, Hsf, Fst, shark, anon-23Da, period, Ddc, Gpdh (D. subobscura clone added to the microarrays), Gdh and Treh. (We also had valid expression data for desat1 but not for desat2.) Differential gene expression with the q-value cut-off chosen for the linear contrasts between cold- (P13) and warm-adapted (P22) populations was found only for Hsp26, which showed increased expression levels in P13 populations. However, the q-value thresholds for Hsp68 (0.058; P13 > P22), Fst (0.061; P13 < P22), and Treh (0.061; P13 < P22) were low enough as to suggest that these three candidate genes also diverged in gene expression levels in our populations. (q-value thresholds for the other genes were: Hsp70, 0.506; Hsf, 0.505; shark, 0.293; anon-23Da, 0.135; desat1, 0.148; period, 0.505; Dcd, 0.353; Gdh, 0.096 [see Additional file 1].) It seems, therefore, that our survey in D. subobscura apparently links thermal adaptation in a temporally stable environment to some specific candidate genes (Hsp26, Hsp68, Fst, and Treh) previously associated to thermotolerance in D. melanogaster. It is also worth saying that the Hsp60 gene had a q-value threshold of 0.082 (P13 > P22) and might have also been associated with thermal adaptation. On the other hand, Hsp83 is known to be expressed during normal development in D. subobscura but increased transcription occurs when flies are reared at heat-shock temperatures from 26 to 34°C [58]. Consistent with this finding we did not observe a statistically significant differential gene expression between cold- and warm-adapted populations for Hsp83 (q-value threshold 0.467).

Mapping of differentially expressed genes by in situ hybridization and correlated expression patterns

In situ hybridizations to the polytene chromosomes were routinely carried out after crossing wild-type D. subobscura males with virgin females from the ch-cu marker strain. A total of 106 genes out of 306 differing in gene expression as concluded from the a priori contrasts were used as probes, and most of them (83%) yielded hybridization signals [see Additional file 3]. In all cases the negative results were due to failures in the amplifications. In no instance an exchange of genes among the different chromosomal elements of D. melanogaster and D. subobscura has been detected, which agrees with the well supported evidence for the established chromosomal homologies previously proposed for each Muller/Sturtevant/Novitski element [59]. The order of genes within each chromosomal element, however, is known to have widely changed among different species via the fixation of paracentric inversions (e.g.; [60–62]).

From the chromosomal homologies between D. melanogaster and D. subobscura [59], and the distribution of the ~120-megabase euchromatic portion of the D. melanogaster genome on each chromosomal arm [63], we tested the null hypothesis that the location of genes differing in gene expression between warm- and cold-adapted populations on the different D. subobscura chromosomes was at random (Table 2). The G-test for goodness of fit [64] detected a marginally nonsignificant random distribution (p = 0.059. If candidate genes Hsp68 and Fst on chromosome O, Hsp60 on chromosome U, and Treh on chromosome E are included, then p = 0.076), with chromosome J apparently being overrepresented and chromosome O underrepresented. It seems interesting to contrast these results with the previously reported chromosomal inversion shifts in the thermal populations [16]. Inversions on chromosomes J, E and, to a lesser extent, chromosomes A and O showed clear shifts in frequency according to the thermal regime, whereas those on chromosome U showed no trend whatsoever. It is not at all evident how the distribution of the differentially expressed genes on chromosomes matches with these patterns as, for example, chromosome E and U are well represented in Table 2 but their behaviour after two years of thermal evolution was completely different.

Table 2 Chromosomal distribution of genes differing in gene expression

Figs. 5, 6 show a schematic representation of the location of genes along the D. subobscura chromosomes, together with the chromosomal regions covered by the gene arrangements segregating in the thermal stocks. To test whether or not there is any clustering of genes within the inverted segments of the chromosomes we have relied on published tables where inversion lengths of D. subobscura are given as percentages of the total length of all chromosomes (except for the dot chromosome; [65]). Since overlapping inversions in chromosomes U, E and O introduce potential confusion, we have only considered the following heterokaryotypes in the analysis: Ast/A2, Jst/J1, Ust/U1+2, U1+2/U1+8+2(i.e., segment covered by inversion U8), Est/E 1+2+9 , E 1+2+9 /E 1+2+9+3 (segment covered by E3), E 1+2+9 /E 1+2+9+12 (segment covered by E12), Ost/O 3+4 , O 3+4 /O3+4+2(segment covered by O2), and O 3+4 /O3+4+7(segment covered by O7). The G-test for goodness of fit detected a highly significant deviation from a random distribution (G (9) = 61.55; p < 0.001), with all but O7 chromosome segments covered by inversions having a higher than expected number of genes.

Figure 5
figure 5

Physical mapping of differentially expressed genes on chromosomes A, J, U of Drosophila subobscura. Schematic representation of the location of 17 differentially expressed genes mapped along chromosome A (the sex chromosome), 14 mapped along chromosome J (homologous to arm 3L in D. melanogaster), and 20 mapped along chromosome U (homologous to arm 2L in D. melanogaster). The centromere is placed on the left (black circle) and the telomere on the right. The linear order of genes is that in the standard gene arrangements, and the chromosomal regions covered by inversions segregating in the thermal stocks (labelled on the right-hand side next to the segments; overlapping inversions underlined) are indicated. (For further details on the formation of the gene arrangements by overlapping inversions, see [84].)

Figure 6
figure 6

Physical mapping of differentially expressed genes on chromosomes E and O of Drosophila subobscura. Schematic representation of the location of 20 differentially expressed genes mapped along chromosome E (homologous to arm 2R in D. melanogaster), and 16 mapped along chromosome O (homologous to arm 3R in D. melanogaster). The centromere is placed on the left (black circle) and the telomere on the right. The linear order of genes is that in the standard gene arrangements, and the chromosomal regions covered by inversions segregating in the thermal stocks are indicated (labelled on the right-hand side next to the segments; overlapping inversions underlined). Inversions O5 and O7 were also sporadically found but are ignored here: the first is associated to a lethal gene and the second is probably the result of a recombination event in the O3+4+7/Ost heterokaryotype [16]).

We further analyzed the correlated expression levels in cold- (P13) and warm-adapted (P22) populations and compared the correlation matrices of gene expression levels between P13 and P22 populations with a Mantel test [64] by randomly permuting 10,000 times the rows and columns of one of the matrices. For those physically mapped genes (Figs. 5, 6) the Mantel tests showed that correlated patterns of expression were similar for all chromosomes (ranging from Mantel's r = 0.418, p = 6.0 × 10-4 for chromosome U to Mantel's r = 0.891, p = 2.0 × 10-4 for chromosome J). The same conclusion was also obtained when simultaneously comparing the correlation patterns of the 306 differently expressed genes (Mantel's r = 0.783, p = 1.0 × 10-4).

To summarize, the physical mapping has shown that a larger than expected number of differentially expressed genes are located inside inverted chromosomal segments, and that the past thermal selection regime does not seem to have significantly changed the correlated patterns of gene expression. Since inversions in D. subobscura are know to influence temporal patterns of linkage disequilibrium between allozymes [66, 67], disentangling the effects of inversions will be obviously required in any future approach if we want to identify the relevant candidate genes underlying thermal adaptation in D. subobscura.

Conclusion

By looking at changes in the transcriptome this study has identified a large number of genes that appear to be involved in thermal adaptation in Drosophila subobscura, with a substantial fraction implicated in metabolism. Interestingly, expression levels of four previously reported candidate genes for thermotolerance in Drosophila (Hsp26, Hsp68, Fst, and Treh) were found to be correlated with past thermal selection regime. The highest experimental temperature used in our stocks (22°C) was not expected to be stressful. This was apparently confirmed by the lack of induction in the P22 populations of a heat shock response known to occur in D. subobscura [58]. The data also suggest an association with polymorphic inversions as some clustering of genes within inverted chromosomal sections was detected. This result is probably not surprising in view of the rapidly evolved latitudinal clines in inversion frequencies after the introduction of the species into the New World [12], as well as the quick response of inversion polymorphism to laboratory temperature [16]. The challenge now is to elucidate what associations are causal and what are due to correlated responses or hitchhiking arisen from linkage disequilibrium with the inversions.

Methods

Microarray experiment

a) Experimental populations and sampling protocol

All nine laboratory populations used here were initiated from an ancestral population of Drosophila subobscura derived from a large outbred stock collected in November 1999 at the estimated Chilean epicentre of the original New World invasion (Puerto Montt, Chile, 41° 28' S, 73° 00' W [68]). From that ancestral population three sets (P13, P18 and P22) of three replicate populations each (R1, R2 and R3) were set up in May 2001 and have since kept at three experimental temperatures on a discrete generation, controlled breeding under constant larval density (~5 larvae/mL of food) and constant 12:12 light:dark period: cold (13°C), optimum (18°C) and warm (22°C), respectively. The number of breeding adults per population is typically well over 1,500 flies. Complete details of the derivation and maintenance of these populations have been previously described [15, 16].

Because we wanted to compare cDNA microarray gene expression patterns for all populations in a nested ANOVA design (see below), the use of a common reference mRNA to be competitively hybridized to treatment mRNAs coming from the different populations seemed to be more appropriate for statistical analysis than a design involving the correct pairing of all nine samples [69]. Therefore, the three 'optimum' P18 populations were sampled in April 2004 (two generations prior the experimental flies from the same populations were obtained) by placing a large number of eggs (± 4 h) in 120-mL plastic chambers (~100 eggs per chamber) with spoons containing 30 mL of David's killed-yeast Drosophila medium [70] stained with 0.05% bromophenol blue. This dye has no effect on larval growth and allows for accurate staging of third instar larvae just prior to pupation at their maximum size [71]. Larvae with clean guts that stopped feeding and started to wander on the wall of the plastic chambers were gently removed with a spatula, cleaned several times with distilled water, placed in bunches with 25 larvae each inside microcentrifuge tubes containing 500 μ L of TRI Reagent® (Molecular Research Center, Inc.), and immediately homogenized and stored at -80°C. These larvae (hereafter referred to as C18) provided the mRNA used as reference.

Subsequently, samples from all nine populations were obtained in May-June 2004 (25 generations at 13°C, 35 at 18°C, and 46 at 22°C) by placing eggs into twelve 130-mL bottles (~200–250 eggs per bottle). These bottles were cultured at 18°C and emerging adults were dumped into Plexiglas cages for egg collections. Eggs for the experiment were collected over a six days period by placing Petri dishes containing non-nutritive agar with a generous smear of live yeast in the cages. As before, ~100 eggs (± 4 h) were placed at 18°C in 120-mL plastic chambers with stained Drosophila medium to sample the treatment third instar larvae for further mRNA extraction.

b) RNA extraction and gene expression analysis

Total RNA was extracted from the frozen homogenized larvae by using TRI Reagent® (MRC, Inc.), and mRNA was extracted by using Promega PolyATtract® isolation system following the manufacturer's specifications. Three mRNA extractions were performed from 9,000 C18 reference larvae (i.e., 3,000 larvae from each replicated population) to obtain a single reference pooled mRNA, and four independent extractions from 250 larvae each were made for each treatment population. This procedure ensured true replication in the experiment; namely, the reference mRNA was always hybridized with treatment mRNAs coming from independent larvae and extractions.

Relative mRNA levels were determined by parallel two-colour hybridization to cDNA microarrays from the Canadian Drosophila Microarray Centre (CDMC D12Kv1; [72]). The arrayed elements in these slides represent approximately 10,500 D. melanogaster genes to which seven D. subobscura genes were added as positive controls. As previously discussed, we were confident that the use of D. melanogaster arrays in heterologous hybridization with D. subobscura could offer a reasonable warranty to obtaining reliable data. However, the divergence between both species (~25 Mya; [29]) may have been the reason for the relatively high number of genes that failed to hybridize (Fig. 1), even though other factors like lack of gene expression in third instar larvae may have also been important. 1 μg of poly(A)+-enriched RNA was labelled using the SuperScript Indirect cDNA labelling system (Invitrogen Corporation, California, USA). mRNA concentrations and cDNA synthesis were checked with Agilent 2100 bioanalyzer (Agilent Technologies, California, USA). Equal amounts of labelled cDNA were combined with 10 μg of yeast tRNA, speed vac dried and re-suspended in 230 μl Dig Easy Hyb solution (Boehringer-Roche). The solution was incubated at 65°C for 10 min to denature the probes.

Hybridizations and washes were performed using the automatic system Lucidea SlidePro (Amersham, UK). The hybridization was allowed to proceed for 15 h at 25°C, and the slides were sequentially washed three times at 50°C for 10 min with medium stringency buffer (1 × SCC, 0.1% SDS), twice at room temperature for 1 min with high stringency buffer (1 × SCC), post wash buffer (0.1 × SCC) and air dried. Then each slide was scanned using an Axon GenePix 4000B microarray scanner (Axon Instruments, Union City, California, USA). Data were extracted from the scanned images using GenePix® Pro (Axon Instruments) microarray image analysis. Labelling, hybridization and scanning were carried out at the Plataforma de Transcriptòmica from the Parc Cientific de Barcelona and Universitat de Barcelona (PCB-UB; [73]).

The raw data were adjusted using lowess normalization software (TIGR MIDAS ver. 2.19; [74]) with a tri-cube weight function and 0.33 smooth parameter applied to the C18 reference mRNA dye-labelled green (Cy3) or red (Cy5). For each experimental population four microarrays were independently hybridized and scanned, adding to 36 arrays in total. There were 14,440 duplicated spots on each array, and only the spots that passed a quality control of image analysis (i.e., array elements with intensities significantly different from background) were used in the differential expression analysis. The gene spots were further filtered by excluding those with less than 57 out of 72 (i.e., <79%) valid expression observations (Fig. 1), leaving 4,651 probes for differential gene expression analysis. The data acquired from these procedures were relative measures of gene expression independent of larval size differences among the thermal stocks.

c) Experimental design and data analysis

The unit of analysis here is the population, and the three replicated populations (R1, R2 and R3) of each thermal selection stock were treated as a random factor nested within experimental temperature (13, 18 and 22°C), which was a fixed effect [64]. Given any treatment population pop = 13R1, 13R2,⋯, 22R3, and any probe g = 1, ⋯, G for which valid expression levels were obtained, we use notation Z c o n t ( G ) g MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGAbGwdaqhaaWcbaGaem4yamMaem4Ba8MaemOBa4MaemiDaqNaeiikaGccbaGae83raCKaeiykaKcabaGaem4zaCgaaaaa@37C7@ , Z c o n t ( R ) g MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGAbGwdaqhaaWcbaGaem4yamMaem4Ba8MaemOBa4MaemiDaqNaeiikaGccbaGae8NuaiLaeiykaKcabaGaem4zaCgaaaaa@37DD@ to denote normalized and background adjusted gene expression from the reference C18 mRNA sample that was dyed green (G) or red (R); and Z p o p ( G ) g MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGAbGwdaqhaaWcbaGaemiCaaNaem4Ba8MaemiCaaNaeiikaGccbaGae83raCKaeiykaKcabaGaem4zaCgaaaaa@3674@ , Z p o p ( R ) g MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGAbGwdaqhaaWcbaGaemiCaaNaem4Ba8MaemiCaaNaeiikaGccbaGae8NuaiLaeiykaKcabaGaem4zaCgaaaaa@368A@ for gene expression intensities obtained from the treatment mRNAs dyed green or red. For each treatment population Z p o p ( G ) g MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGAbGwdaqhaaWcbaGaemiCaaNaem4Ba8MaemiCaaNaeiikaGccbaGae83raCKaeiykaKcabaGaem4zaCgaaaaa@3674@ / Z c o n t ( R ) g MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGAbGwdaqhaaWcbaGaem4yamMaem4Ba8MaemOBa4MaemiDaqNaeiikaGccbaGae8NuaiLaeiykaKcabaGaem4zaCgaaaaa@37DD@ , Z p o p ( R ) g MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGAbGwdaqhaaWcbaGaemiCaaNaem4Ba8MaemiCaaNaeiikaGccbaGae8NuaiLaeiykaKcabaGaem4zaCgaaaaa@368A@ / Z c o n t ( G ) g MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGAbGwdaqhaaWcbaGaem4yamMaem4Ba8MaemOBa4MaemiDaqNaeiikaGccbaGae83raCKaeiykaKcabaGaem4zaCgaaaaa@37C7@ are the relative intensity ratios measured from the corresponding slides. The fully balanced dye-reversal experimental design can be written as the linear model:

y i j k l m g = μ g + T i g + ℛ j ( i ) g + D k g + T D i k g + A l ( i j k ) g + ε i j k l m g , ( 1 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWG5bqEdaqhaaWcbaGaemyAaKMaemOAaOMaem4AaSMaemiBaWMaemyBa0gabaGaem4zaCgaaOGaeyypa0dcciGae8hVd02aaWbaaSqabeaacqWGNbWzaaGccqGHRaWkcqWGubavdaqhaaWcbaGaemyAaKgabaGaem4zaCgaaOGaey4kaSYenfgDOvwBHrxAJfwnHbqeg0uy0HwzTfgDPnwy1aaceaGae43gHi1aa0baaSqaaiabdQgaQjabcIcaOiabdMgaPjabcMcaPaqaaiabdEgaNbaakiabgUcaRiabdseaenaaDaaaleaacqWGRbWAaeaacqWGNbWzaaGccqGHRaWkcqWGubavcqWGebardaqhaaWcbaGaemyAaKMaem4AaSgabaGaem4zaCgaaOGaey4kaSIae4haXh0aa0baaSqaaiabdYgaSjabcIcaOiabdMgaPjabdQgaQjabdUgaRjabcMcaPaqaaiabdEgaNbaakiabgUcaRiab=v7aLnaaDaaaleaacqWGPbqAcqWGQbGAcqWGRbWAcqWGSbaBcqWGTbqBaeaacqWGNbWzaaGccqGGSaalcaWLjaGaaCzcamaabmaabaGaeGymaedacaGLOaGaayzkaaaaaa@78AC@

where for probe g, μgis the overall grand mean of the log 2 relative intensity ratios; T i g MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGubavdaqhaaWcbaGaemyAaKgabaGaem4zaCgaaaaa@30BC@ is the fixed effect of the i th experimental treatment (P13, P18, P22); ℛ j ( i ) g MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaat0uy0HwzTfgDPnwy1egaryqtHrhAL1wy0L2yHvdaiqaacqWFBeIudaqhaaWcbaGaemOAaOMaeiikaGIaemyAaKMaeiykaKcabaGaem4zaCgaaaaa@3D6A@ is the random effect of the j th replicate population (R1, R2, R3) within treatment i; D k g MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGebardaqhaaWcbaGaem4AaSgabaGaem4zaCgaaaaa@30A0@ is the fixed effect of dye k (Cy3, Cy5); T D i k g MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGubavcqWGebardaqhaaWcbaGaemyAaKMaem4AaSgabaGaem4zaCgaaaaa@332C@ is the interaction term; A l ( i j k ) g MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaat0uy0HwzTfgDPnwy1egaryqtHrhAL1wy0L2yHvdaiqaacqWFaeFqdaqhaaWcbaGaemiBaWMaeiikaGIaemyAaKMaemOAaOMaem4AaSMaeiykaKcabaGaem4zaCgaaaaa@40CE@ is the random effect of slide l = 1, 2 within treatment i, replicate j and dye k; and ε i j k l m g MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaiiGacqWF1oqzdaqhaaWcbaGaemyAaKMaemOAaOMaem4AaSMaemiBaWMaemyBa0gabaGaem4zaCgaaaaa@36B9@ is the residual error associated with the corresponding log 2 relative intensity ratio of the ijklm th spot. This linear model easily allows partitioning all sources of experimental variation: biological (temperature and replicated population effects) and technical (dye and slide effects).

Notice that for the treatment effect we are interested in (i.e., the T i g MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGubavdaqhaaWcbaGaemyAaKgabaGaem4zaCgaaaaa@30BC@ component due to thermal adaptation) the linear model (1) can be conveniently reduced to the following two-level nested ANOVA model:

y i j k g = μ g + T i g + ℛ j ( i ) g + e i j k g , ( 2 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWG5bqEdaqhaaWcbaGaemyAaKMaemOAaOMaem4AaSgabaGaem4zaCgaaOGaeyypa0dcciGae8hVd02aaWbaaSqabeaacqWGNbWzaaGccqGHRaWkcqWGubavdaqhaaWcbaGaemyAaKgabaGaem4zaCgaaOGaey4kaSYenfgDOvwBHrxAJfwnHbqeg0uy0HwzTfgDPnwy1aaceaGae43gHi1aa0baaSqaaiabdQgaQjabcIcaOiabdMgaPjabcMcaPaqaaiabdEgaNbaakiabgUcaRiabdwgaLnaaDaaaleaacqWGPbqAcqWGQbGAcqWGRbWAaeaacqWGNbWzaaGccqGGSaalcaWLjaGaaCzcamaabmaabaGaeGOmaidacaGLOaGaayzkaaaaaa@5B3B@

where the sum of squares for the error term e i j k g MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGLbqzdaqhaaWcbaGaemyAaKMaemOAaOMaem4AaSgabaGaem4zaCgaaaaa@339A@ is simply the sum of the sum of squares for the remainder terms in (1). The usefulness of this model reduction is obvious to efficiently perform randomization tests to test the null hypothesis about treatment effects in a randomized (i.e., random assignment) experiment [75]. Permutation tests are far less sensitive to the presence of outliers and are particularly necessary with unequal sample sizes; i.e., when some data points are missing as is usually the case with microarray experiments. The null hypothesis of no treatment or evolutionary thermal regime effect was tested here after performing random permutations among replicate and selection temperature for the among selection temperature F-statistics. Each test used 10,000 random permutations of the log 2 relative intensity ratios (recall that when N = 72 there are 72 ! ( 24 ! 24 ! 24 ! ) ≈ 2.56 × 10 32 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaadaWccaqaaiabiEda3iabikdaYiabcgcaHaqaaiabcIcaOiabikdaYiabisda0iabcgcaHiabikdaYiabisda0iabcgcaHiabikdaYiabisda0iabcgcaHiabcMcaPaaacqGHijYUcqaIYaGmcqGGUaGlcqaI1aqncqaI2aGncqGHxdaTcqaIXaqmcqaIWaamdaahaaWcbeqaaiabiodaZiabikdaYaaaaaa@44BF@ possible assignments of observations).

A planned comparison between the two treatment means from the stocks at the two extreme thermal regimes (i.e., P13 vs. P22) that had already diverged for 71 generations was also performed for each probe g. The permutation tests were performed following [76]; namely, we first calculated the F g 1 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGgbGrdaqhaaWcbaGaem4zaCgabaGaeGymaedaaaaa@3035@ statistic for the observed data and next the residuals of the log 2 relative intensity ratios from the populations at P13 and P22 were randomly allocated to both treatment temperatures. From B = 10,000 random permutations we got a set of null statistics F g 0 b MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGgbGrdaqhaaWcbaGaem4zaCgabaGaeGimaaJaemOyaigaaaaa@3180@ , b = 1,2, ..., B; and the p-value was computed as:

p g = ∑ b = 1 B # { [ F g 1 , F g 0 b ] ≥ F g 1 } B + 1 . ( 3 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGWbaCdaWgaaWcbaGaem4zaCgabeaakiabg2da9maaqahabaWaaSaaaeaacqGGJaWidaGadeqaamaadmaabaGaemOray0aa0baaSqaaiabdEgaNbqaaiabigdaXaaakiabcYcaSiabdAeagnaaDaaaleaacqWGNbWzaeaacqaIWaamcqWGIbGyaaaakiaawUfacaGLDbaacqGHLjYScqWGgbGrdaqhaaWcbaGaem4zaCgabaGaeGymaedaaaGccaGL7bGaayzFaaaabaGaemOqaiKaey4kaSIaeGymaedaaaWcbaGaemOyaiMaeyypa0JaeGymaedabaGaemOqaieaniabggHiLdGccqGGUaGlcaWLjaGaaCzcamaabmaabaGaeG4mamdacaGLOaGaayzkaaaaaa@5274@

Given the high-dimensionality of the data set the p-values were adjusted based on the concept of false discovery rate (FDR; [34]). If no probe g is differentially expressed the p-values will follow a U (0,1), where U stands for 'uniform distribution'. The so-called Mixture Distribution Partitioning (MDP) methodology assumes that the distribution of p-values consists of a set of null p0 and alternative p1 components. This partition forms the basis for estimating various quantities as for example the q-values, which were obtained here with the QVALUE software [35]. The problem now is to select a threshold of significance to identify a set of genes likely to be differentially expressed. As an unsupervised criterion we used a q-value cut-off ≤ 0.05 for the P13 vs. P22 planned comparisons, meaning that the maximum expected proportion of false positives incurred when calling a particular gene 'differentially expressed' is 5%.

d) Computer software for statistical analysis

The computer programs used for statistical data analyses were MATLAB algebra program environment (ver. 7.0.4 [77]) together with the collection of tools supplied by the Statistics Toolbox (ver. 5.0.2 [78]). The statistical software packages STATISTICA version 6 [79] and SPSS version 13 [80] were also used.

Mapping of candidate genes

The flies used for physical mapping of candidate genes were collected from a natural population in Bordils (70 Km North-east of Barcelona, Spain; 42° 3' N, 2° 54' E). About 150 males were individually crossed to three or four virgin females from the ch-cu marker strain to help in the identification of polymorphic inversions (the genetic background of this strain is highly homogeneous and fixed for the standard arrangements in all major acrocentric chromosomes but chromosome O, where it is fixed for arrangement O 3+4 .

DNA isolation, DNA amplification, polytene chromosome preparation and in situ hybridization were carried out using standard techniques [81]. The karyotype of D. subobscura consists of five acrocentric chromosomes and a dot chromosome. Following [82] the large chromosomes in this species are traditionally named as A (= X, the sex chromosome), J (= chromosomal element D of Mueller/Sturtevant/Novitski and homologous to arm 3L in Drosophila melanogaster [59]), U (= chromosomal element B and homologous to arm 2L), E (= chromosomal element C and homologous to arm 2R), and O (= chromosomal element E and homologous to arm 3R). The five major acrocentric chromosomes and the dot chromosome are divided into 100 sections (A: 1 – 16; J: 17 – 35; U: 36 – 53; E: 54 – 74; O: 75 – 99; Dot :100), and each section into 3–5 subsections (A, B, ...) [83].

References

  1. Gockel J, Kennington WJ, Hoffmann AA, Goldstein DB, Partridge L: Nonclinality of molecular variation implicates selection in maintaining a morphological cline of Drosophila melanogaster. Genetics. 2001, 158: 319-323.

    PubMed Central  CAS  PubMed  Google Scholar 

  2. de Jong G, Bochdanovits Z: Latitudinal clines in Drosophila melanogaster: body size, allozyme frequencies, inversion frequencies, and the insulin-signalling pathway. J Genet. 2003, 82: 207-223.

    Article  CAS  PubMed  Google Scholar 

  3. Santos M, Brites D, Laayouni H: Thermal evolution of pre-adult life history traits, geometric size and shape, and developmental stability in Drosophila subobscura. J Evol Biol. 2006, 19: 2006-2021. 10.1111/j.1420-9101.2006.01139.x.

    Article  CAS  PubMed  Google Scholar 

  4. Hoffmann AA, Weeks AR: Climatic selection on genes and traits after a 100 year-old invasion: a critical look at the temperate-tropical clines in Drosophila melanogaster from eastern Australia. Genetica. 2007, 129: 133-147. 10.1007/s10709-006-9010-z.

    Article  PubMed  Google Scholar 

  5. Umina PA, Weeks AR, Kearney MR, McKechnie SW, Hoffmann AA: A rapid shift in a classic clinal pattern in Drosophila reflecting climate change. Science. 2005, 308: 691-693. 10.1126/science.1109523.

    Article  CAS  PubMed  Google Scholar 

  6. Levitan M, Etges WJ: Climate change and recent genetic flux in populations of Drosophila robusta. BMC Evol Biol. 2005, 5: 4-10.1186/1471-2148-5-4.

    Article  PubMed Central  PubMed  Google Scholar 

  7. Balanya J, Oller JM, Huey RB, Gilchrist GW, Serra L: Global genetic change tracks global climate warming in Drosophila subobscura. Science. 2006, 313: 1773-1775. 10.1126/science.1131002.

    Article  CAS  PubMed  Google Scholar 

  8. Van Heerwaarden B, Hoffmann AA: Global warming: fly populations are responding rapidly to climate change. Curr Biol. 2007, 17: R16-R18. 10.1016/j.cub.2006.11.035.

    Article  CAS  PubMed  Google Scholar 

  9. Krimbas CB: Drosophila subobscura: Biology, Genetics and Inversion Polymorphism. 1993, Hamburg, Germany: Verlag Dr. Kovac

    Google Scholar 

  10. Prevosti A, Ribó G, Serra L, Aguadé M, Balañá J, Monclús M, Mestres F: Colonization of America by Drosophila subobscura: Experiment in natural populations that supports the adaptive role of chromosomal-inversion polymorphism. Proc Natl Acad Sci USA. 1988, 85: 5597-5600. 10.1073/pnas.85.15.5597.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  11. Huey RB, Gilchrist GW, Hendry AP: Using invasive species to study evolution. Species Invasions: Insights to Ecology, Evolution and Biogeography. Edited by: Sax DF, Gaines SD, Stachowicz JJ. 2005, Sunderland: Sinauer Associates, 139-164.

    Google Scholar 

  12. Balanyà J, Serra L, Gilchrist GW, Huey RB, Pascual M, Mestres F, Solé E: Evolutionary pace of chromosomal polymorphism in colonizing populations of Drosophila subobscura: an evolutionary time series. Evolution. 2003, 57: 1837-1845. 10.1554/02-577.

    Article  PubMed  Google Scholar 

  13. Gilchrist GW, Huey RB, Balanyà J, Pascual M, Serra L: A time series of evolution in action: a latitudinal cline in wing size in South American Drosophila subobscura. Evolution. 2004, 58: 768-780. 10.1554/03-414.

    Article  PubMed  Google Scholar 

  14. Rose MR, Nusbaum TJ, Chippendale AK: Laboratory evolution: the experimental wonderland and the Cheshire Cat syndrome. Adaptation. Edited by: Rose MR, Lauder GV. 1996, San Diego: Academic Press, 221-241.

    Google Scholar 

  15. Santos M, Fernández Iriarte P, Céspedes W, Balanyà J, Fontdevila A, Serra L: Swift laboratory thermal evolution of wing shape (but not size) in Drosophila subobscura and its relationship with chromosomal inversion polymorphism. J Evol Biol. 2004, 17: 841-855. 10.1111/j.1420-9101.2004.00721.x.

    Article  CAS  PubMed  Google Scholar 

  16. Santos M, Céspedes W, Balanyà J, Trotta V, Calboli FCF, Fontdevila A, Serra L: Temperature-related genetic changes in laboratory populations of Drosophila subobscura: evidence against simple climatic-based explanations for latitudinal clines. Am Nat. 2005, 165: 258-273. 10.1086/427093.

    Article  PubMed  Google Scholar 

  17. Rifkin SA, Kim J, White KP: Evolution of gene expression in the Drosophila melanogaster subgroup. Nat Genet. 2003, 33: 138-144. 10.1038/ng1086.

    Article  CAS  PubMed  Google Scholar 

  18. Wray GA, Hahn MW, Abouheif E, Balhoff JP, Pizer M, Rockman MV, Romano LA: The evolution of transcriptional regulation in eukaryotes. Mol Biol Evol. 2003, 20: 1377-419. 10.1093/molbev/msg140.

    Article  CAS  PubMed  Google Scholar 

  19. Khaitovich P, Weiss G, Lachmann M, Hellmann I, Enard W, Muetzel B, Wirkner U, Ansorge W, Pääbo S: A neutral model of transcriptome evolution. PLoS Biol. 2004, 5: 682-689.

    Google Scholar 

  20. Michalak P, Noor MAF: Association of misexpression with sterility in hybrids of Drosophila simulans and D. mauritiana. J Mol Evol. 2004, 59: 277-282. 10.1007/s00239-004-2622-y.

    Article  CAS  PubMed  Google Scholar 

  21. Gibson G, Weir B: The quantitative genetics of transcription. Trends Genet. 2005, 21: 616-623. 10.1016/j.tig.2005.08.010.

    Article  CAS  PubMed  Google Scholar 

  22. Ranz JM, Machado CA: Uncovering evolutionary patterns of gene expression using microarrays. Trends Ecol Evol. 2006, 21: 29-37. 10.1016/j.tree.2005.09.002.

    Article  PubMed  Google Scholar 

  23. Feder ME, Walter J-C: The biological limitations of transcriptomics in elucidating stress and stress responses. J Evol Biol. 2005, 18: 901-910. 10.1111/j.1420-9101.2005.00921.x.

    Article  CAS  PubMed  Google Scholar 

  24. Partridge L, Barrie B, Fowler K, French V: Thermal evolution of pre-adult life-history traits in Drosophila melanogaster. J Evol Biol. 1994, 7: 645-663. 10.1046/j.1420-9101.1994.7060645.x.

    Article  Google Scholar 

  25. Joshi A, Mueller LD: Density-dependent natural selection in Drosophila: trade-offs between larval food acquisition and utilization. Evol Ecol. 1996, 10: 463-474. 10.1007/BF01237879.

    Article  Google Scholar 

  26. Santos M, Borash DJ, Joshi A, Bounlutay N, Mueller LD: Density-dependent natural selection in Drosophila: evolution of growth rate and body size. Evolution. 1997, 51: 420-432. 10.2307/2411114.

    Article  Google Scholar 

  27. Eisen MB, Brown PO: DNA arrays for analysis of gene expresión. Methods Enzymol. 1999, 303: 179-205.

    Article  CAS  PubMed  Google Scholar 

  28. Renn SCP, Aubin-Horth N, Hofmann HA: Biologically meaningful expression profiling across species using heterologous hybridization to a cDNA microarray. BMC Genomics. 2004, 5: 42-10.1186/1471-2164-5-42.

    Article  PubMed Central  PubMed  Google Scholar 

  29. Russo CA, Takezaki N, Nei M: Molecular phylogeny and divergence times of drosophilid species. Mol Biol Evol. 1995, 12: 391-404.

    CAS  PubMed  Google Scholar 

  30. Krimbas CB, Loukas M: The inversion polymorphism of Drosophila subobscura. Evol Biol. 1980, 12: 163-234.

    Google Scholar 

  31. Betrán E, Santos M, Ruiz A: Antagonistic pleiotropic effect of second-chromosome inversions on body size and early life-history traits in Drosophila buzzatii. Evolution. 1998, 52: 144-154. 10.2307/2410929.

    Article  Google Scholar 

  32. Hoffmann AA, Sgrò CM, Weeks AR: Chromosomal inversion polymorphism and adaptation. Trends Ecol Evol. 2004, 19: 482-488. 10.1016/j.tree.2004.06.013.

    Article  PubMed  Google Scholar 

  33. Frischer LE, Hagen FS, Garber RL: An inversion that disrupts the Antennapedia gene causes abnormal structure and localization of RNAs. Cell. 1986, 47: 1017-1023. 10.1016/0092-8674(86)90816-0.

    Article  CAS  PubMed  Google Scholar 

  34. Benjamini Y, Hochberg Y: Controlling the false discovery rate – A practical and powerful approach to multiple testing. J Roy Stat Soc Series B-Methodol. 1995, 57: 289-300.

    Google Scholar 

  35. Storey JD, Tibshirani R: Statistical significance for genomewide studies. Proc Natl Acad Sci USA. 2003, 100: 9440-9445. 10.1073/pnas.1530509100.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  36. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, Harris MA, Hill DP, Issel-Tarver L, Kasarskis A, Lewis S, Matese JC, Richardson JE, Ringwald M, Rubin GM, Sherlock G: Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000, 25: 25-29. 10.1038/75556.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  37. GOToolBox: Functional Investigation of Gene Datasets. Genome Biol. 2004, 5 (12): R101-10.1186/gb-2004-5-12-r101. [http://crfb.univ-mrs.fr/GOToolBox]

  38. Kristensen TN, Sørensen P, Pedersen KS, Kruhøffer M, Loeschcke V: Inbreeding by environmental interactions affect gene expression in Drosophila melanogaster. Genetics. 2006, 173: 1329-1336. 10.1534/genetics.105.054486.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  39. Bochdanovits Z, de Jong G: Experimental evolution in Drosophila melanogaster: interaction of temperature and food quality selection regimes. Evolution. 2003, 57: 1829-1836. 10.1554/02-740.

    Article  PubMed  Google Scholar 

  40. Feder ME, Hofmann GE: Heat-shock proteins, molecular chaperones and the stress response: evolutionary and ecological physiology. Annu Rev Physiol. 1999, 61: 243-282. 10.1146/annurev.physiol.61.1.243.

    Article  CAS  PubMed  Google Scholar 

  41. Parsell DA, Taulien J, Lindquist S: The role of heat-shock proteins in thermotolerance. Phil Trans Roy Soc London Series B. 1993, 339: 279-286. 10.1098/rstb.1993.0026.

    Article  CAS  Google Scholar 

  42. Krebs RA: A comparison of Hsp70 expression and thermotolerance in adults and larvae of three Drosophila species. Cell Stress Chaperones. 1999, 4: 243-249. 10.1379/1466-1268(1999)004<0243:ACOHEA>2.3.CO;2.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  43. Bettencourt BR, Kim I, Hoffmann AA, Feder ME: Response to natural and laboratory selection at the Drosophila hsp 70 genes. Evolution. 2002, 56: 1796-1801. 10.1554/0014-3820(2002)056[1796:RTNALS]2.0.CO;2.

    Article  CAS  PubMed  Google Scholar 

  44. Frydenberg J, Hoffmann AA, Loeschcke V: DNA sequence variation and latitudinal associations in hsp23, hsp26 and hsp27 from natural populations of Drosophila melanogaster. Mol Ecol. 2003, 12: 2025-2032. 10.1046/j.1365-294X.2002.01882.x.

    Article  CAS  PubMed  Google Scholar 

  45. McColl G, Hoffmann AA, McKechnie SW: Response of two heat shock genes to selection for knockdown heat resistance in Drosophila melanogaster. Genetics. 1996, 143: 1615-1627.

    PubMed Central  CAS  PubMed  Google Scholar 

  46. Hoffmann AA, Sorensen JG, Loeschcke V: Adaptation of Drosophila to temperature extremes: bringing together quantitative and molecular approaches. J Therm Biol. 2003, 28: 175-216. 10.1016/S0306-4565(02)00057-8.

    Article  Google Scholar 

  47. Anderson AR, Collinge JE, Hoffmann AA, Kellett M, McKechnie SW: Thermal tolerance tradeoffs associated with the right arm of chromosome 3 and marked by the Hsromega gene in Drosophila melanogaster. Heredity. 2003, 90: 195-202. 10.1038/sj.hdy.6800220.

    Article  CAS  PubMed  Google Scholar 

  48. Lin YJ, Seroude L, Benzer B: Extended life-span and stress resistance in the Drosophila mutant methuselah. Science. 1998, 282: 943-946. 10.1126/science.282.5390.943.

    Article  CAS  PubMed  Google Scholar 

  49. Duvernell DD, Schmidt PS, Eanes WF: Clines and adaptive evolution in the methuselah gene region in Drosophila melanogaster. Mol Ecol. 2003, 12: 1277-1285. 10.1046/j.1365-294X.2003.01841.x.

    Article  CAS  PubMed  Google Scholar 

  50. Goto SG: A novel gene that is up-regulated during recovery from cold shock in Drosophila melanogaster. Gene. 2001, 270: 259-264. 10.1016/S0378-1119(01)00465-6.

    Article  CAS  PubMed  Google Scholar 

  51. Ferrante AW, Reinke R, Stanley ER: Shark, a Src homology 2, ankyrin repeat, tyrosine kinase, is expressed on the apical surfaces of ectodermal epithelia. Proc Natl Acad Sci USA. 1995, 92: 1911-1915. 10.1073/pnas.92.6.1911.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  52. Greenberg AJ, Moran JR, Coyne JA, Wu CI: Ecological adaptation during incipient speciation revealed by precise gene replacement. Science. 2003, 302: 1754-1757. 10.1126/science.1090432.

    Article  CAS  PubMed  Google Scholar 

  53. Costa R, Peixoto AA, Barbujani G, Kyriacou CP: A latitudinal cline in a Drosophila clock gene. Proc R Soc London B. 1992, 250: 43-49. 10.1098/rspb.1992.0128.

    Article  CAS  Google Scholar 

  54. Morgan TJ, Mackay TFC: Quantitative trait loci for thermotolerance phenotypes in Drosophila melanogaster . Heredity. 2006, 96: 232-242. 10.1038/sj.hdy.6800786.

    Article  CAS  PubMed  Google Scholar 

  55. Oakeshott JG, Gibson JB, Anderson PR, Knibb WR, Chambers GK: Alcohol dehydrogenase and glycerol-3-phosphate dehydrogenase clines in Drosophila melanogaster on different continents. Evolution. 1982, 36: 86-96. 10.2307/2407970.

    Article  Google Scholar 

  56. Oakeshott JG, McKechnie SW, Chambers GK: Population genetics of the metabolically related Adh, Gpdh, and Tpi polymorphisms in Drosophila melanogaster. I. Geographic variation in Gpdh and Tpi allele frequencies in different continents. Genetica. 1984, 63: 21-29. 10.1007/BF00137461.

    Article  Google Scholar 

  57. Sezgin E, Duvernell DD, Matzkin LM, Duan YH, Zhu CT, Verrelli BC, Eanes WF: Single-locus latitudinal clines and their relationship to temperate adaptation in metabolic genes and derived alleles in Drosophila melanogaster. Genetics. 2004, 168: 923-931. 10.1534/genetics.104.027649.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  58. Arbona M, De Frutos R, Tanguay RM: Transcriptional and translational study of the Drosophila subobscura hsp83 gene in normal and heat-shock conditions. Genome. 1993, 36: 694-700.

    Article  CAS  PubMed  Google Scholar 

  59. Powell JR: Progress and Prospects in Evolutionary Biology. The Drosophila Model. 1997, New York: Oxford University Press

    Google Scholar 

  60. Hartl DL, Lozovskaya ER: The Drosophila Genome Map: A Practical Guide. 1995, New York: Springer-Verlag

    Google Scholar 

  61. Segarra C, Ribó G, Aguadé M: Differentiation of Muller's chromosomal elements D and E in the obscura group of Drosophila. Genetics. 1996, 144: 139-146.

    PubMed Central  CAS  PubMed  Google Scholar 

  62. Ranz JM, Casals F, Ruiz A: How malleable is the eukaryotic genome?. Genome Res. 2001, 11: 230-239. 10.1101/gr.162901.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  63. Adams MD, Celniker SE, Holt RA, Evans CA, Gocayne JD, Amanatides PG, Scherer SE, Li PW, Hoskins RA, Galle RF, George RA, Lewis SE, Richards S, Ashburner M, Henderson SN, Sutton GG, Wortman JR, Yandell MD, Zhang Q, Chen LX, Brandon RC, Rogers YHC, Blazej RG, Champe M, Pfeiffer BD, Wan KH, Doyle C, Baxter EG, Helt G, Nelson CR, Miklos GLG, Abril JF, Agbayani A, An HJ, Andrews-Pfannkoch C, Baldwin D, Ballew RM, Basu A, Baxendale J, Bayraktaroglu L, Beasley EM, Beeson KY, Benos PV, Berman BP, Bhandari D, Bolshakov S, Borkova D, Botchan MR, Bouck J, Brokstein P, Brottier P, Burtis KC, Busam DA, Butler H, Cadieu E, Center A, Chandra I, Cherry JM, Cawley S, Dahlke C, Davenport LB, Davies A, de Pablos B, Delcher A, Deng ZM, Mays AD, Dew I, Dietz SM, Dodson K, Doup LE, Downes M, Dugan-Rocha S, Dunkov BC, Dunn P, Durbin KJ, Evangelista CC, Ferraz C, Ferriera S, Fleischmann W, Fosler C, Gabrielian AE, Garg NS, Gelbart WM, Glasser K, Glodek A, Gong FC, Gorrell JH, Gu ZP, Guan P, Harris M, Harris NL, Harvey D, Heiman TJ, Hernandez JR, Houck J, Hostin D, Houston DA, Howland TJ, Wei MH, Ibegwam C, Jalali M, Kalush F, Karpen GH, Ke ZX, Kennison JA, Ketchum KA, Kimmel BE, Kodira CD, Kraft C, Kravitz S, Kulp D, Lai ZW, Lasko P, Lei YD, Levitsky AA, Li JY, Li ZY, Liang Y, Lin XY, Liu XJ, Mattei B, McIntosh TC, McLeod MP, McPherson D, Merkulov G, Milshina NV, Mobarry C, Morris J, Moshrefi A, Mount SM, Moy M, Murphy B, Murphy L, Muzny DM, Nelson DL, Nelson DR, Nelson KA, Nixon K, Nusskern DR, Pacleb JM, Palazzolo M, Pittman GS, Pan S, Pollard J, Puri V, Reese MG, Reinert K, Remington K, Saunders RDC, Scheeler F, Shen H, Shue BC, Siden-Kiamos I, Simpson M, Skupski MP, Smith T, Spier E, Spradling AC, Stapleton M, Strong R, Sun E, Svirskas R, Tector C, Turner R, Venter E, Wang AHH, Wang X, Wang ZY, Wassarman DA, Weinstock GM, Weissenbach J, Williams SM, Woodage T, Worley KC, Wu D, Yang S, Yao QA, Ye J, Yeh RF, Zaveri JS, Zhan M, Zhang GG, Zhao Q, Zheng LS, Zheng XQH, Zhong FN, Zhong WY, Zhou XJ, Zhu SP, Zhu XH, Smith HO, Gibbs RA, Myers EW, Rubin GM, Venter JC: The genome sequence of Drosophila melanogaster. Science. 2000, 287: 2185-2195. 10.1126/science.287.5461.2185.

    Article  PubMed  Google Scholar 

  64. Sokal RR, Rohlf FJ: Biometry. 1995, New York: Freeman, 3

    Google Scholar 

  65. Krimbas CB, Loukas M: Drosophila subobscura: lengths of chromosome segments heterozygotes for inversions, to be used in IFR calculations. Eur Dros Pop Biol Group Bull. 1979, 3: 4-11.

    Google Scholar 

  66. Fontdevila A, Zapata C, Alvarez G, Sánchez L, Méndez J, Enríquez I: Genetic coadaptation in the chromosomal polymorphism of Drosophila subobscura. I. Seasonal changes of gametic disequilibrium in a natural population. Genetics. 1983, 105: 935-955.

    PubMed Central  CAS  PubMed  Google Scholar 

  67. Rodríguez-Trelles F: Seasonal cycles of allozyme-by-chromosomal-inversion gametic disequilibrium in Drosophila subobscura. Evolution. 2003, 57: 839-848. 10.1554/0014-3820(2003)057[0839:SCOAGD]2.0.CO;2.

    Article  PubMed  Google Scholar 

  68. Brncic D, Budnik M: Colonization of Drosophila subobscura Collin in Chile. Dros Inf Serv. 1980, 55: 20-

    Google Scholar 

  69. Kerr MK, Churchill GA: Statistical design and the analysis of gene expresión microarray data. Genet Res. 2001, 77: 123-128. 10.1017/S0016672301005055.

    CAS  PubMed  Google Scholar 

  70. David J: A new medium for rearing Drosophila in axenic conditions. Dros Inf Serv. 1962, 36: 128-

    Google Scholar 

  71. Andres AJ, Thummel CS: Methods for quantitative analysis of transcription in larvae and prepupae. Drosophila melanogaster: Practical Uses in Cell and Molecular Biology. Methods in Cell Biology. Edited by: Goldstein L, Fyrberg E. 1994, New York: Academic Press, 44: 565-573.

    Google Scholar 

  72. Canadian Drosophila Microarray Centre. [http://www.flyarrays.com]

  73. Parc Científic de Barcelona. [http://www.pcb.ub.es/homePCB/live/en/p126.asp]

  74. Saeed AI, Sharov V, White J, Li J, Liang W, Bhagabati N, Braisted J, Klapa M, Currier T, Thiagarajan M, Sturn A, Snuffin M, Rezantsev A, Popov D, Ryltsov A, Kostukovich E, Borisovsky I, Liu Z, Vinsavich A, Trush V, Quackenbush J: TM4: A free, open-source system for microarray data management and analysis. Biotechniques. 2003, 34: 374-378.

    CAS  PubMed  Google Scholar 

  75. Edgington ES: Randomization Tests. 1995, New York: Marcel Dekker, 3

    Google Scholar 

  76. Manly BFJ: Randomization, Bootstrap and Monte Carlo Methods in Biology. 1997, London: Chapman & Hall, 2

    Google Scholar 

  77. The MathWorks Inc: : MATLAB, version 7.0.4. The Language of Technical Computing. 2005, [http://www.mathworks.com]

    Google Scholar 

  78. The MathWorks Inc: Statistics Toolbox for use with MATLAB, version 5.0.2. 2005, [http://www.mathworks.com]

    Google Scholar 

  79. StatSoft Inc: STATISTICA (data analysis software system), version 6. 2003, [http://www.statsoft.com]

    Google Scholar 

  80. SPSS Inc: SPSS for Windows. 2004, [http://www.spss.com]

    Google Scholar 

  81. Laayouni H, Santos M, Fontdevila A: Toward a physical map of Drosophila buzzatii: Use of randomly amplified polymorphic DNA polymorphisms and sequence-tagged site landmarks. Genetics. 2000, 156: 1797-1816.

    PubMed Central  CAS  PubMed  Google Scholar 

  82. Mainx F, Koske T, Smital E: Untersuchungen über die chromosomale Struktur europäischer Vertreter der Drosophila obscura Gruppe. Zeitschrift für inducktive Abstammungs- und Vererbungslehre. 1953, 85: 354-372. 10.1007/BF00309673.

    CAS  Google Scholar 

  83. Kunze-Mühl E, Müller E: Weitere Untersuchungen über die chromosomale Struktur und natürlichen Strukturtypen von D. subobscura. Chromosoma. 1958, 9: 559-570. 10.1007/BF02568093.

    Article  PubMed  Google Scholar 

  84. Krimbas CB: The inversion polymorphism of Drosophila subobscura. Drosophila Inversion Polymorphism. Edited by: Krimbas CB, Powell JR. 1992, Boca Raton: CRC Press, 127-220.

    Google Scholar 

Download references

Acknowledgements

We thank Lidia Sevilla from Parc Científic de Barcelona for her invaluable technical assistance with the microarray experiment, Victor Moreno and Miguel Pignatelli for helpful discussions with the experimental design and Gene Ontology classification, and Elias Zintzaras for useful discussions on statistical analysis. We would also like to acknowledge the critical and constructive comments of two anonymous reviewers which have helped to improve the article significantly. This work was supported by grants BOS2003-05904-C02 and CGL2006-13423-01/BOS from the Ministerio de Ciencia y Tecnología (Spain), 2001SGR-00207 from the Direcció General de Recerca (Generalitat de Catalunya) to the GBE, grant GEN2001-4846-C05-02 from the Ministerio de Educación y Ciencia (Spain), and by Fundación Ramón Areces (Spain).

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Mauro Santos.

Additional information

Authors' contributions

HL sampled the thermal populations, made the RNA extractions, participated in the design of the experiment, carried out statistical analysis, Gene Ontology searches, and drafted the manuscript. FG-F, BEC-S and VT designed primers, carried out in situ hybridizations on the polytene chromosomes of D. subobscura, and read all salivary gland squashes from the in situ hybridizations. SB and MC coordinated the microarray experiments for generating the original data. MS set up and maintained the thermal populations, designed the study, carried out statistical analyses and drafted the final version of the manuscript. All authors read and approved the final manuscript.

Electronic supplementary material

12862_2006_336_MOESM1_ESM.xls

Additional file 1: Summary of microarray results from the 4,651 cDNA clones with valid expression observations. The details for the different columns are explained within the file itself. (XLS 2 MB)

12862_2006_336_MOESM2_ESM.xls

Additional file 2: Molecular function GO categories for the differentially expressed genes. Analysis of functional categories defined by the Gene Ontology project [36] using the GOToolBox [37]. (XLS 72 KB)

12862_2006_336_MOESM3_ESM.pdf

Additional file 3: Physical map of differentially expressed genes. Localization by in situ hybridization on the salivary gland chromosomes of Drosophila subobscura of 88 differentially expressed genes. (PDF 77 KB)

Authors’ original submitted files for images

Rights and permissions

Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( https://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Laayouni, H., García-Franco, F., Chávez-Sandoval, B.E. et al. Thermal evolution of gene expression profiles in Drosophila subobscura. BMC Evol Biol 7, 42 (2007). https://doi.org/10.1186/1471-2148-7-42

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/1471-2148-7-42

Keywords