Changes in gene expression linked with adult reproductive diapause in a northern malt fly species: a candidate gene microarray study

Background Insect diapause is an important biological process which involves many life-history parameters important for survival and reproductive fitness at both individual and population level. Drosophila montana, a species of D. virilis group, has a profound photoperiodic reproductive diapause that enables the adult flies to survive through the harsh winter conditions of high latitudes and altitudes. We created a custom-made microarray for D. montana with 101 genes known to affect traits important in diapause, photoperiodism, reproductive behaviour, circadian clock and stress tolerance in model Drosophila species. This array gave us a chance to filter out genes showing expression changes during photoperiodic reproductive diapause in a species adapted to live in northern latitudes with high seasonal changes in environmental conditions. Results Comparisons among diapausing, reproducing and young D. montana females revealed expression changes in 24 genes on microarray; for example in comparison between diapausing and reproducing females one gene (Drosophila cold acclimation gene, Dca) showed up-regulation and 15 genes showed down-regulation in diapausing females. Down-regulation of seven of these genes was specific to diapause state while in five genes the expression changes were linked with the age of the females rather than with their reproductive status. Also, qRT-PCR experiments confirmed couch potato (cpo) gene to be involved in diapause of D. montana. Conclusions A candidate gene microarray proved to offer a practical and cost-effective way to trace genes that are likely to play an important role in photoperiodic reproductive diapause and further in adaptation to seasonally varying environmental conditions. The present study revealed two genes, Dca and cpo, whose role in photoperiodic diapause in D. montana is worth of studying in more details. Also, further studies using the candidate gene microarray with more specific experimental designs and target tissues may reveal additional genes with more restricted expression patterns.


Background
Several insect species in temperate regions develop and reproduce during the summer, but become dormant as the winter approaches to avoid unfavourable conditions for development and reproduction. This dormancy is usually linked with a period of developmental arrest, i.e. diapause, at the embryonic, larval, pupal or adult stage. Adult reproductive diapause, where oogenesis and vitellogenesis are arrested, has been studied in monarch butterflies [1] and in several grasshopper [2] and Drosophila species [e.g. [3]; see also [4]]. In Drosophila species, this kind of diapause has been found to offer resource alloca-tion trade-offs: the females entering reproductive diapause early in the summer will not produce progeny during the ongoing season, but they have a higher stress resistance, age more slowly [5] and have better chances to survive over the winter period and produce progeny during the next summer than the non-diapausing females. In most organisms the onset of reproductive diapause is largely determined by the length of the light period (photoperiodic diapause), which is not surprising, as seasonal changes in the day length in different latitudes are more predictable and consistent than changes in temperature or humidity [4]. Photoperiodic diapause also allows the females to anticipate the coming cold season in advance, which is very important for organisms with a long generation time. From an evolutionary point of view, reproduc-tive diapause is an important biological process which involves many life-history parameters important for survival and reproductive fitness at both individual and population levels.
The discovery of an adult photoperiodic reproductive diapause in Drosophila melanogaster [6,7] opened a possibility to study the role of the circadian clock genes in photoperiodic time measurement and diapause, as well as interactions between genes affecting photoperiodism and hormonal and physiological changes involved in diapause. The finding also gave a chance to search for genes affecting fitness-related traits linked with diapause, e.g. life span, age-specific mortality, fecundity, resistance to cold and starvation stress, lipid content, development time, egg-to adult viability, accessory gland synthesis and mating behaviour [8][9][10]. D. melanogaster females are known to enter adult reproductive diapause in response to short days and low temperatures in Europe [11] and in North America [10,12]. However, Emerson et al., [13] showed lately that in two natural populations representing latitudinal extremes in eastern North America temperature is the main determinant of the diapause. Even though the reproductive diapause in D. melanogaster seems to be a relatively weak response to photoperiod and is only observed at temperatures below 14°C [7,14], the genes found to be linked with sexual maturation and diapause in this species are good candidate genes for diapause studies in more northern Drosophila species with stronger, longer-lasting photoperiodic diapause responses.
Our study species, D. montana, belongs to the D. virilis group, which originated in continental Asia. It has expanded its distribution around the northern hemisphere from 30°N to 70°N in latitude and from 0 to 3000 meters in altitude, adapting to various kinds of environmental conditions [15,16]. The females of this species have the ability to survive through the harsh winter conditions in photoperiodic adult reproductive diapause [17]. They have adapted to environmental conditions in northern Scandinavia with a mean temperature of less than 0°C for up to six winter months, a minimum temperature reaching occasionally -30°C or even below. D. montana flies predict the forthcoming cold season mainly by changes in the day length even though also the temperature seems to have some effect on the percentage of diapausing females [18].
We had two aims in this study. The first aim was to create a custom-made microarray for D. montana with candidate genes known to play a role in traits linked with reproductive diapause in D. melanogaster. The second aim was to determine which of the candidate genes show expression changes, and in which direction, in comparisons among diapausing (i.e. those in reproductive diapause), reproducing and young D. montana females. Our candidate gene microarray proved to be a practical and cost-effective way to use information on gene function obtained largely through mutagenesis in a genetic model species and to study whether the same genes play a role in adaptation to different environmental conditions in an ecologically interesting non-model organism. Out of the 101 genes assayed, 24 showed differences in their expression patterns in comparisons among D. montana females at different reproductive modes.

Study flies
We used D. montana flies from an isofemale line O3F66 originated from Oulanka (Northern Finland; 66°N) and maintained in the laboratory under diapause-preventing conditions (continuous light, 19°C) since summer 2003. In this isofemale line 98-100% of the female flies will reach sexual maturity (i.e., will not enter diapause) in continuous light (L:L) while in a light dark (L:D) cycle of 16:8 97% of them will enter reproductive diapause at 16°C (Salminen et al., unpublished). Female flies for the present study were collected under a light CO 2 anaesthesia within one day after their emergence, i.e. before the photoperiodic induction of diapause occurred (Salminen et al., unpublished). One day old flies (three replicates, each including 13 flies) were transferred straight from the culturing conditions into liquid nitrogen for RNA extractions. The rest of the females were placed in vials (13-15 flies per vial) containing 7 ml of yeast-sucrose-agar medium [19] with about 5 mg of yeast sprinkled on top of the medium. The samples of reproducing females were obtained by maintaining the vials in a climate chamber in continuous light (24:0; diapause preventing conditions) for 14 days. Samples of diapausing females were obtained by keeping the vials in a climate chamber in a L:D cycle of 16:8 (diapause inducing conditions) for the same time period. The temperature in both chambers was 16 ± 1°C; lower temperatures were not used to avoid cold acclimation. Three sets of 13 females representing the two groups (reproducing/diapausing) were snap-frozen in liquid nitrogen immediately after collecting them from the climate chambers during a light period (two to three hours after lights went on in 16:8 light dark cycle). Flies were then stored at -84°C for RNA extractions.

Designing a candidate gene DNA microarray for D. montana
The species-specific candidate gene microarray constructed for D. montana included 101 genes (Additional file 1) known to affect diapause (7 genes), phototransduction (25 genes), courtship behaviour (14 genes) and circadian clock (26 genes), as well as responses to cold (2 genes) and heat (27 genes) in D. melanogaster. The genes were selected from the D. melanogaster genome [20] using FlyBase gene ontology searches on biological response [21]. The most conserved regions of the candidate genes were indentified by aligning the exon regions of the D. melanogaster gene to their orthologs in D. virilis with Blat search in UCSC Genome Bioinformatics database [22]. D. virilis was chosen as a reference species, as it belongs to the same species group as D. montana and has its whole genome sequenced (Agencourt Bioscience Corporation).
The primers for the candidate genes of D. montana were designed for one to three unique 200-500 bp regions with Primer3 program [23] or manually, using information on the gene alignments mentioned above (primer sequences are available from the authors upon request). The gene sequences of D. montana were amplified with PCR and sequenced using both DNA and cDNA as a template to observe possible intron regions. The sequenced gene regions of D. montana were checked against D. melanogaster annotated genes and the genome sequences of D. melanogaster and D. virilis in FlyBase to ensure that they represented right genes and did not have homology to any other regions of the genome. Next, an online program e-Array 5.3 (Agilent) was used to design one to four 60 bp oligonucleotide probes per gene using exon regions. Finally, the probes were synthesized in situ with liquid chemistry and arrayed on an "Agilent 60-mer Multi-Pack Gene Expression Microarray" platform with one-colour system.

Preparing the RNA samples for the microarray
The three samples for the microarray experiments included young females ("young"), 14 days old females cultured under diapause-preventing ("reproducing") and 14 days old females cultured under diapause-inducing conditions ("diapausing"). Total RNA was extracted from the whole body of 13 individuals per sample with three replicates per sample. Qiagen RNA extraction kit with RNase-Free DNase treatment was used according to the manufacturer's protocols (Qiagen). The RNA was purified using Qiagen RNAeasy purification kit and the concentration and the purity of each sample was measured with NanoDrop ND-1000 spectrophotometer (Nano-Drop Technologies, Wilmington, DE, USA). 350 ng of total RNA from each sample was amplified and Cy3labeled with Agilent's Low RNA Input Linear Amplification Kit PLUS (One Color) and processed together with Agilent's One-Colour RNA Spike-in Kit. The RNA/cRNA concentration (and also specific activity for cRNA) was checked with Nanodrop ND-1000 and the RNA/cRNA quality with BioRad's Experion electrophoresis station, both before and after the amplifications. 600 ng of each cRNA sample was hybridized to Agilent's 8 × 15K D. montana custom arrays with Gasket slide at 65°C overnight (17 h) using Agilent's Gene Expression Hybridization Kit and washed with Agilent's Gene Expression Wash Pack and Stabilization and Drying solution. Arrays were scanned with Agilent Technologies Scanner (model G2505B) and numerical results were extracted with Feature Extraction version 9.5.3.1 using 018413_D_F_20071204 grid, GE1-v5_95_Feb07 protocol and GE1_QCM_Feb07 QC metric set.
In addition to the 101 candidate genes we placed on our chip, we used seven genes, Ef1α48D, eIF-4a, Gapdh1, Gapdh2, RpL11, RpL19 and RpL27A, which have been used as control genes in Drosophila whole genome microarrays (GeneChip Drosophila genome 2.0 Array, Affymetrix, 2006). These genes were added on the chip to select the best control genes for microarray validation with qRT-PCR analysis.

Statistical analysis
The microarray contained 197 unique probe sequences representing 101 candidate genes and seven housekeeping genes. Some of the probes had two copies per array, in which case the averages of their signals were calculated. The Agilent control probes included 77 negative control probes and 200 spike-in probes (10 different probe types, each represented by 20 copies per array). The raw images of the microarray were processed with Agilent's Feature Extraction software and further analyzed using R statistical analysis software (R Foundation for Statistical Computing) and related Bioconductor analysis package [24].
Data from small custom arrays cannot be normalized using traditional normalization methods developed for the whole genome arrays, because the common assumption of an unbiased trend of up and down-regulated genes and on the modulation of a small proportion of genes may not hold. There are, however, several alternative methods for small custom arrays [see [25]]. Variance stabilizing normalization method (VSN, [26]) using spike-in controls performs best with the small data sets, as it can be parameterized to allow up to 50% of the genes to be modulated [27]. In our data set this allowed all the experimental probes (representing about half of the probes) to be modulated between sample groups (the spike-in probes not expected to modulate formed the second half). Accordingly, microarray data were normalized within the groups and across all samples using the VSN method. The statistical analyses were carried out using Bioconductor related Limma package [24] and the differentially expressed genes in each comparison were selected requiring a strict selection criteria of fold-change > 2.0 and false discovery rate (fdr) < 0.05. The false discovery rates were calculated using the Benjamini and Hochberg's method as implemented in the Limma package.

Quantitative real time PCR (qRT-PCR)
Expression changes in six of the 24 genes that showed upor downregulation on the microarray were validated with qRT-PCR. The genes represent all the main categories based on their function in D. melanogaster (diapause, phototransduction, courtship behaviour, circadian clock, and cold and heat tolerance, Additional file 1). Primers for these genes (Table 1) were the same as those used in the microarray experiment, except in the case of Dca and cpo genes for which the original PCR product was too long to be used in qRT-PCR analysis. Samples of diapausing, reproducing and young females were collected as described earlier and cDNA was generated by reverse transcription with iScript cDNA kit (Bio-Rad) from newly made RNA extractions. RNA was exracted as described earlier and the purity of each sample was measured with NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). Ef1α48D(Elongation factor 1 aplha 48D) was used as reference gene for comparisons of diapausing and reproducing flies and eIF-4a (Eukaryotic initiation factor4a) gene was used in comparison of diapausing or reproducing females to young females. These genes were selected on the basis of their expression stability among all the control genes (i.e. housekeeping genes) tested in the microarray analysis. Amplification efficiency of each primer pair was checked with serial dilutions of cDNA (data not shown) and the real efficiency values (all with E-values > 90% and R 2 -values > 0.980) were used in the analysis. Real time PCR analysis was performed in 25 μl reactions containing SYBR Green SuperMix (BioRad), 0.3 μM of each gene-specific primer and 1 μl of cDNA solution. Cycling conditions in BioRad CFX96 instrument were 95°C for 3 min, followed by 40 cycles of denaturing at 95°C for 10 sec and annealing at 54°C for 30 sec and extension at 72°C for 30 sec. Fold change relative to the expression of the sample genes and reference genes was calculated using normalized expression (ΔΔ(Ct)) method with default threshold values using CFX Manager Software which will give the same result as the Pfaffl method [28] when using one target and one reference gene [29]. Two to three biological replicates (each mixture of 13 individual flies) were analyzed in each of the comparisons and the purity of the PCR reactions was checked with the melting curve analysis after amplification. Independent-samples t-test (SPSS 13 for windows, SPSS Inc.) was used to test the statistical significance of normalized gene expression in different comparisons.

Results
The microarray protocol used in the present study proved to be both reliable and repeatable. Experimental probe intensities of the raw data on the microarray were highly correlated between the three replicates of each of the three female samples (median Pearson's correlation = 0.885). Correlations between the 'technical replicates' (RNA from three samples of diapausing females hybridized on arrays at two different hybridization times) were even higher (Pearson's correlations > 0.990), indicating high repeatability and low technical variation. The microarray data has been deposited in NCBI's Gene Expression Omnibus [30] and are accessible through GEO Series accession number GSE18475 and GPL accession number GPL9293.

Up-and down-regulation of candidate genes in microarray
Comparison between diapausing, reproducing and young D. montana females revealed significant expression changes in 24 of the 101 candidate genes on the microarray (Fig. 1). We were most interested in genes showing differential expression in diapausing females versus reproducing females. In this comparison the diapausing females showed up-regulation in only one gene, Drosophila cold acclimation gene (Dca), and down-regulation Table 1

Hsp26
Heat tolerance 5' -CATCGGAGGACAGAGAGGAG -3' 5' -ATTGTCCACCACCTTGACGT -3' in 15 genes. Down-regulation of seven of the genes (cac, nonA, CkIIα, dco, Hop, Hsc70-3, Hsp83) was specific to diapause state (they only appeared in the comparison between diapausing and reproducing females) while in five genes (shakB, Ddc, e, slo, so) the expression changes were linked with the age of the females rather than with their reproductive status (Fig. 1). One additional gene, couch potato (cpo), which is known to affect diapause in D. melanogaster (see Discussion), showed a significant down-regulation in diapausing and reproducing females compared to young females and a marginally significant (FC = 1.8 and p = 0.055) up-regulation in diapausing females compared to reproducing females. Eigth of the 15 genes (CG7650, FKBP59, tilB, tws, Fmr1, slmb, Prp5, Cdk9) that were down-regulated in diapausing females when compared to reproducing females showed downregulation also in young females when compared to reproducing females. The last two genes, Hsp20 and Hsp26, were expressed only when reproducing females were compared to young females (Fig. 1).

Confirmation of microarray results with qRT-PCR
Changes in the function of six genes, cpo (diapause), FKBP59 (phototransduction), tilB (courtship behaviour), Fmr1 (circadian clock), Dca (cold tolerance) and Hsp26 Normalized fold expression of the candidate genes in microarray analysis (heat tolerance), showing differential expression between at least two of the three female developmental stages in the microarray analysis were confirmed with qRT-PCR (we made altogether 12 comparisons). All the tested genes showed expression changes in the same direction in qRT-PCR as in the microarray analysis, even though the magnitude and statistical significance of changes differed between the two methods in half of the cases (Additional file 2). Comparisons between the diapausing and reproducing females showed significant up-regulation in diapausing females in cpo and Dca and down-regulation in tilB, while in Fmr1 and FKBP59 the expression difference remained non-significant (Additional file 2, Fig. 2).
In the case of cpo it is worth to note that qRT-PCR revealed significant differences both in comparisons of diapausing and reproducing and reproducing and young females while in microarray analysis the difference between diapausing and reproducing females was marginally non-significant (p = 0.055).

Feasibility of a candidate gene array for detecting evolutionary interesting genes in non-model organisms
Gene expression studies have previously been performed for genetically less well-known Drosophila species using microarrays constructed for D. melanogaster, but due to high levels of genetic divergence between D. melanogaster and the species under investigation (usually several My), these kinds of data may be difficult to interpret and may lack reliability [e.g. [31]]. Another solution with nonmodel organisms is to prepare custom microarrays or to trace changes in gene expression via transcriptome analy-ses [e.g. [32]]. Both of these methods are feasible, but they are still quite expensive and difficult to analyze when genome information is incomplete. We decided to create a candidate gene microarray for D. montana to survey the genes underlying the life-history traits important in reproductive diapause in a model Drosophila species.
The used multi-pack gene expression microarray with one-colour system provides a possibility to compare several (up to eight) samples with each other on one chip, and thus is cost-effective and appropriate for ecological study designs. The most substantial restriction in a candidate gene microarray is that it gives information only on a limited number of genes. However, new genes can be added to new versions of the array once detected in model organisms or found in the study species e.g. through transcriptome or QTL analyses. Recently published work on gene expression during D. melanogaster egg development before and after reproductive diapause indicated several genes connected to diapause behaviour [33] worth of adding to our array.

Changes in gene expression pattern linked with reproductive diapause
Emerson et al. have argued that the incidence of diapause is determined by genes in three successive modules, photoperiodism, hormonal events and the diapause itself and an extra module, circadian clock, which can potentially affect diapause through any of the other three [34]. The present study design has allowed us to detect gene expression changes mainly in the hormonal and diapause modules.
The most interesting findings in our study were the upregulation of Dca gene in diapausing females when com-Normalized fold expressions of candidate genes in qRT-PCR analysis pared to reproducing and young females and the expression changes in cpo gene between the three female types. Dca (also called smp-30, according to its mammalian homologue) has been suggested to be associated with Ca 2+ regulation [35], and it may play a role in maintaining Ca 2+ content during cold exposure [36]. In D. melanogaster the expression level of this gene has been found to increase during 1-2 day cold acclimation in 15°C [36] and to decrease during recovery after cold treatment in 0°C [37,38]. Expression level of Dca has also been found to increase with age in D. melanogaster [36], contrary to the situation in some mammal species [35]. In our study, the high expression level of Dca in diapausing females may be connected to female cold tolerance, which is known to increase during diapause in D. montana (L. Vesala, unpublished). The fact that Dca showed no expression changes between the reproducing and young females suggests that the detected transcriptional changes may not be connected to fly aging. The cpo gene has been suggested to play an important role in diapause in D. melanogaster [e.g. [39]]. This gene could provide a link between the insulin signaling pathway and the downstream hormones involved directly in the regulation of vitellogenesis [34]. cpo is expressed in the peripheral and central nervous system of Drosophila embryos, larvae, and adults [40] and was originally named "couch potato" because several partial loss-of-function alleles of this gene cause hyperactive behavior in D. melanogaster [40]. cpo is also known to be involved in the latitudinal cline of diapause behavior in D. melanogaster [39]. In our study, cpo was highly down-regulated in diapausing and reproducing females, when compared to young females and it also showed significant up-regulation in diapausing females when compared to reproducing females in qRT-PCR analysis (in microarray analysis the last comparison gave a marginally significant p-value). These findings suggest that cpo could play a role in diapause in D. montana.
Another interesting group of genes in the scope of the present study are the seven genes down-regulated in diapausing females, but not in young females, when compared to reproducing females. Two of these genes, Casein kinase alpha II subunit (CkIIα) and discs overgrown (dco) have been listed as circadian rhythm genes in FlyBase (see Additional file 1) and they may be involved in several biological processes such as locomotor activity or mating behavior [41], which might explain their low expression in diapausing females. Two genes, no on or off transient-A (nonA) and cacophony (cac), have been found to affect courtship behavior and male courtship song in D. melanogaster [42]; nonA gene is also known to regulate other genes and have pleiotropic effects on fly vision [43], whereas cac is known to be involved in adult locomotor behavior [44], phototransduction and visual behavior [45]. Finally, the three last genes down-regulated only in diapausing females were the heat shock protein genes Hsp70/Hsp90 organizing protein homolog (Hop), Heat shock protein cognate 3 (Hsc70-3)and Heat shock protein 83 (Hsp83). Some studies have indicated that changes in the expression level of heat shock genes could be connected with insect diapause, but the data are not congruent [46,47]. It is worth noting that while other heat shock genes in our study showed significant down-regulation in diapausing females, Heat shock protein 20(Hsp20)and Heat shock protein 26 (Hsp26) showed expression changes only between reproducing and young females (Hsp20 was down-regulated and Hsp26 up-regulated in reproducing females when compared to young females).
Some genes that showed expression variation on our array may not play an important role in the regulation of reproductive diapause itself, but instead may exhibit pleiotropic effects on life-history traits associated with diapause. On the other hand, some genes like the insulinpathway genes (Dp110)and Insulin-like receptor(InR) did not show expression changes on our array, even though insulin signaling has been shown to play a role in the diapause in nematodes [48], mosquitoes [49] and Drosophila [50]. In particular, Dp110 gene has been shown to be necessary for ovarian development of D. melanogaster and allelic variation in this gene has been found to be associated with latitudinal increase in the incidence of diapause [50]. Also some genes important in photoperiodic time measurement, like period (per) and timeless (tim), were lacking from the list of genes showing expression variation in our study. Expression changes in these genes need to be studied with specific experimental design, taking into account e.g. the daily cycles in the expression level of the genes and targeting the studies on specific tissues like the heads of the flies.

Conclusions
The present study shows that a candidate gene microarray provides a practical and cost-effective platform for detecting genes that play a role e.g. in adaptation into different kinds of environmental conditions in ecologically interesting but genetically less well-known organisms. The value of this kind of arrays depends on how well the chosen candidate genes represent the genes and genetic pathways affecting the traits of interest, how repeatable the results are and how well the experimental design supports the asked questions. New genes can be added on array once detected in other species (see [33]) or found in the study species e.g. through transcriptome analyses. Candidate gene arrays will be especially useful when combined with other methods of the candidate gene approach; using SNPs of the array genes as markers in association studies will, for example, help to trace the link between phenotype, genotype and gene expression. The genes that show expression changes on array, like Dca and cpo in the present study, will be good targets of future studies on molecular evolution. Since traits like diapause incidence vary between populations, especially along latitudinal clines, the frequency of any putative mutation changing the function of genes affecting these traits would be predicted to show variation linked with phenotypic variation. Also, it will be interesting to trace signs of selection on the coding region of the genes.

Additional material
Authors' contributions MK constructed the microarray chip, analyzed the results, performed the qRT-PCR analysis, and wrote the manuscript. TS participated in the microarray construction, study design and manuscript writing and cultured the flies and extracted the RNA for the microarray study. AL performed the statistical analysis for microarray results and participated in the manuscript writing. LV made the RNA extractions for the qRT-PCR analysis and participated in the manuscript writing. AH participated in the design of the study and coordination and helped to write the manuscript. All authors read and approved the final manuscript.