Selection for reproduction under short photoperiods changes diapause-associated traits and induces widespread genomic divergence

Changes in day length often act as token stimuli, which help organisms to anticipate the forthcoming cold season (Tauber et al., 1986). In northern arthropod species with a facultative reproductive diapause, emerging females enter reproductive diapause when the day length decreases below a critical point. The timing of diapause, described as a critical day length (CDL), critical night length (CNL) or critical photoperiod, is under strong local selection pressure largely because of restriction in the length of the growing period and stressful overwintering conditions. Several insect species have been found to show latitudinal variation in this trait (e.g. Bradshaw and Lounibos, 1977; Lankinen, 1986a; Schmidt et al., 2005; Tyukmaeva et al., 2011), as well as in the frequency of diapausing individuals (Schmidt et al., 2005), cold tolerance (Božičević et al., 2016; Sörensen et al., 2016) and the period and damping rate of their free-running rhythms (Allemand and David, 1976; Lankinen, 1986a; Pittendrigh and Takamura, 1989). However, the extent to which these multiple responses are correlated or independent is still poorly understood.

The circadian clock systems of most multicellular organisms consist of at least one specialized pacemaker oscillator (the central clock) that responds to environmental signals and coordinates rhythmic output in peripheral oscillators. In Drosophila melanogaster, the peripheral oscillators in different tissues can also be directly entrained by light (Bell-Pedersen et al., 2005), and daily (circadian) rhythms such as locomotor activity and eclosion are controlled by different oscillators, which may differ in period length and light sensitivity (Engelmann and Mack, 1978). Furthermore, circadian rhythms may differ in the peptidergic circuit that links the clock to motor outputs modulating these rhythms (King et al., 2017; King and Sehgal, 2018). Photoperiodic control of seasonal (circannual) rhythms is even more complicated, as these rhythms can be regulated either by a circadian clock network or a photoperiod timer and/or by their interaction (Hut et al., 2013; Denlinger et al., 2017). The photoperiodic timer helps organisms to distinguish short photoperiods from long ones and the traits it regulates are typically on/off traits such as diapause response or migration (Saunders, 1981; Bradshaw and Holzapfel, 2010). Genetic and physiological mechanisms behind photoperiodic time measurement are still largely unknown, some models stressing their independence from, and others their integration with, the circadian clock. In hourglass-type models, photoperiodic time measurement is suggested to function independently and be based on the accumulation of a hypothetical chemical substance during the dark period, so that a photoperiodic response is induced after a certain number of short-night cycles have accumulated in the counter system (Lees, 1973, 1986). A circadian-based model for photoperiodic time measurement, originally proposed by Bünning (1936), has been developed into several robust models involving one or more circadian oscillators that may be coupled to each other and which may show a certain degree of damping (see Vaz Nunes and Saunders, 1999). An hourglass-like clock and the damped oscillator model could, in fact, be based on the same clock mechanism and differ only in the damping coefficient of the oscillator(s) concerned (Saunders, 2005).

Insect photoperiodic responses involve a sequence of events from photoreception and the measurement of night or day length to the accumulation of a hypothetical diapause titer in a counter mechanism and the downstream regulation of events leading to release or retention of neurohormones regulating diapause (Saunders, 2014). Both the measurement of successive photoperiods and their accumulation by the counter mechanism are likely to function under the circadian system (Saunders, 2012), while the downstream cascading system governing circadian behaviors and photoperiodic responses may be completely different (Goto, 2013). Denlinger et al. (2017) suggested that rapid adaptive response in traits regulated by photoperiodic timer could occur without disrupting daily circadian organization even if one or a few circadian clock genes function pleiotropically as a phase reference point for the photoperiodic timer. Circadian clock output genes, whose dysregulation disrupts behavioral rhythms without affecting oscillations in the molecular clock itself, have been identified in several species (Omura et al., 2016; King et al., 2017; King and Sehgal, 2018). It is also becoming obvious that neurotransmitters and neuropeptides play a central role in regulating the production and/or release of hormones involved in diapause induction (Sim and Denlinger, 2013; Nässel et al., 2013; Andreatta et al., 2018). Reproductive diapause can further trigger a hormonal cascade, which induces changes in metabolism (Kubrak et al., 2014), increases cold tolerance (Vesala et al., 2012; Wallingford et al., 2016) and reduces insects' mating drive and pheromone production (Ala-Honkola et al., 2018). However, there is still debate on whether changes in CDL and other photoperiodically regulated traits occur through local selection directed on the circadian clock system (Pittendrigh and Takamura, 1987) or directly on the traits involved (Bradshaw et al., 2003).

Quasinatural selection, where the study organisms are transferred into altered environmental conditions over successive generations, offers an efficient way to examine the role of environmental factors in trait evolution and to study trait correlations and trade-offs (Fry, 2003). For example, Agrawal (2000) performed such selection inducing a host-shift in a polyphagous spider mite, Tetranychus urticae, and found evidence of a correlation between mites' performance on a novel host and their host-plant preference. The ‘evolve and resequence’ approach goes one step further by allowing evolutionary responses to selection to be examined at the genome level to identify diverged genomic regions and functional pathways, as well as localizing genes diverging under the altered selection regime (Turner et al., 2011; Burke, 2012; Wiberg et al., 2017). Martins et al. (2014) used this technique and exposed Drosophila melanogaster flies to Drosophila C virus for 20 generations. After selection, they pool-sequenced treated and control lines and found two genome regions containing several candidate genes connected with an increased survival of the flies under different viral treatments.

We have used quasinatural selection, combined with the ‘evolve and resequence’ approach, to find out whether selection for reproduction under short photoperiods changes CDL and other traits associated with reproductive diapause, and whether it affects these traits directly or by modifying the action of the circadian clock network. The selection experiment was performed on a northern fly species, Drosophila montana, where the photoperiod is a key cue for the induction of adult reproductive diapause. This species shows robust latitudinal clines in CDL (Tyukmaeva et al., 2011; Lankinen et al., 2013; Venera Tyukmaeva, P.L., J.K, H.K. and A.H., unpublished data) and in the frequency of diapausing females under light:dark (LD) cycles that deviate markedly from 24 h photoperiodic cycles in Nanda–Hamner (N–H) experiment (P.L., unpublished results). D. montana is one of the most cold-tolerant Drosophila species (Kellermann et al., 2012), and its cold-tolerance has been found to show variation between latitudinally and altitudinally distinct populations (Vesala and Hoikkala, 2011; R.A.W.W., unpublished results), as well as seasonal photoperiodic cold acclimation (Vesala et al., 2012). We hypothesize that if our selection affects CDL directly, and not through changes in the central circadian clock network, it should affect CDL and other diapause-associated traits, but it should not disrupt fly locomotor activity rhythm, which is regulated by the central circadian clock (King et al., 2017). We also expect the genetic divergence of selection and control lines to resemble divergence detected among latitudinally distinct wild D. montana populations (R.A.W.W., unpublished data), and to include SNP variation associated with genes important in the photoperiodic control of CDL and other traits which diverge during selection.

We established a genetically variable base population of Drosophila montana Stone, Griffen and Patterson 1941 by collecting 102 fertilized females from a northern Oulanka (Finland; 66.4°N) population in summer 2013. Species identification was performed by sequencing part of the mtDNA COI region of one progeny fly per female as described in Simon et al. (1994), using forward and reverse primers: COI_1F 5′-ATCTATCGCCTAAACTTCAGCC-3′ and COI_1R 5′-ACTTCAGGGTGACCAAAAAATC-3′. F₃ progenies of wild-caught D. montana females (about 10 virgin sexually mature females and males per progeny) were transferred into a population cage containing eight bottles with malt medium (Lakovaara, 1969). The flies were allowed to mate and produce progeny in the cage for another 3 generations before starting selection (Fig. 1A). During selection, the flies were maintained in malt vials in wooden light-insulated cabinets, illuminated with compact fluorescent lamps (LedEnergie) with light intensity of 600–1100 lx during the photophase, in 17±1°C. Length of the light period was decreased by 0.5 h in each generation in the chambers containing selection line replicates until the incidence of female diapause exceeded 50%. This helped us to keep the strength of selection close to 50% (keeping the flies continuously in the same day length would have lowered the diapause percentage and reduced selection strength in each generation; Roff, 1994). Control line replicates were maintained through the experiment in constant 24 h light (LL) mimicking the light condition at Oulanka during flies' breeding season in early summer.

In generation F₆, 200 females and males were separated within 3 days of emergence. When sexually mature (∼3 weeks old), the flies were allowed to mate in malt vials for 48 h in groups of 10 flies of each sex. They were then transferred into fresh malt vials for further mating and egg laying and into new vials 3 times at 3 day intervals to prevent larval crowding. Half of the egg vials were transferred into 19 h:5 h light:dark (LD) (selection line; CDL of the base population) and the other half kept in LL (control line) before the flies of next generation flies started to emerge (Fig. 1B). The same procedure was repeated in generation 7 in respective photoperiods, using 100 sexually mature females and males per line, to establish 3 replicates for both lines (Fig. 1C). In generations F₈–F₁₅, selection was continued with 100 females and 50 males per replicate for selection line replicates (half of the females were expected to enter diapause), and 50 females and 50 males for control line replicates (Fig. 1D). In each generation, the flies were allowed to mate for 48 h in the same LD conditions in which they had emerged, after which they were transferred into malt vials for egg laying as described above. The egg vials of selection line replicates were transferred in each generation into 0.5 h shorter photoperiods (in generation F₈ progeny flies were kept in the same LD conditions as their parents because there were few flies, but this was compensated for by transferring the progenies into 1 h shorter photoperiods in the next generation).

The proportion of diapausing females was calculated in each generation by dissecting 77–100 females for each selection line replicate and about 50 females for the control line replicates after progeny production. Selection was completed in generation F₁₅, where the frequency of diapausing females increased above 50% in all selection line replicates. After finishing active selection, the selection line replicates were maintained in LD 15:9 (the final LD conditions used in the selection experiment) and the control line replicates in LL at 16°C in mixed generations.

To determine whether selection for reproduction under short photoperiods induced changes in CDL and other photoperiodically controlled traits, and whether it also affected traits thought to be regulated by circadian clocks, we assayed several phenotypic traits at one or more stages of selection experiment. The first set of phenotypic assays, including CDL and female locomotor activity rhythm, was performed in the base population (Fig. 1B). CDL, cold tolerance and locomotor activity rhythm were analyzed for all selection and control line replicates 2–4 generations after selection had finished (Fig. 1E). Prior to performing the latter assays, all line replicates were maintained in constant light (LL) for two to three generations to exclude possible nongenetic components of line differences according to the protocols of quasinatural selection (Fry, 2003). Generations where the assays were performed are given in Fig. 1.

The frequency of diapausing females under photoperiodic cycles that deviate from daily 24 h cycles in N–H experiments was studied using the flies that had been maintained in LD 15:9 (selection line replicates) and LL (control line replicates) in mixed generations in 16°C for 8 to 10 months after finishing active selection. These experiments, as well as the last CDL measurements, were performed in the University of Oulu (all other experiments were done in the University of Jyväskylä).

Critical day length for diapause induction (CDL) was estimated from the photoperiodic response curves as the point where 50% of females enter diapause. Emerged flies were maintained under different LD cycles for 21 days in malt vials in similar cabinets as the ones used for fly maintenance during selection (compact LedEnergie fluorescent lamps with light intensity of 600-1100 lux during the photophase). After this, the reproductive stage of the females was determined by dissecting their abdomen and recording whether their ovaries contained at least one fully developed egg or whether they had entered diapause (characterized by small ovaries with no egg yolk; Tyukmaeva et al., 2011). The photoperiodic response curve of the base population females was plotted for the LDs 20:4, 19:5, 18:6, 17:7 and 16:8 in 17±1°C (N=59–95 females/LD; exact sample sizes in Table 1). Respective curves for the control and selection line replicates after the experiment were drawn in LDs 21:3, 20:4, 19:5, 18:6, 17:7, 16:8, 15:9, 14:10 and 13:11 in 17±1°C (N=69–97 females/LD; more detailed information about the sample size can be found in Table 1) in generation F17.

CDL of selection and control line replicates was measured a second time 8 to 10 months after finishing active selection using the same protocol and similar cabinets. Here, the cabinets were illuminated with one white fluorescent lamp/chamber (9 W, Megaman, Germany) with light intensity of 300-1000 lux, and the photoperiodic response curves were drawn in LDs 19:5, 18:6, 17:7 and 16:8 for the selection line replicates and in LDs 22:2, 21:3, 20:4, 19:5 and 18:6 for the control line replicates in 16±0.3°C (N=131–507 females/LD; exact sample sizes in Table 1).

The N–H experiment involves sets of experiments where a short light period (10–12 h) is usually coupled with varying lengths of dark periods so that the total period of the LD cycle varies between experiments (Nanda and Hamner, 1958). If the photoperiodic time measurement relies on the function of circadian clock(s), the peak amplitude of insect photoperiodic responses should fluctuate in a ∼24 h rhythm reflecting the period of the circadian oscillator. Based on our earlier N–H experiment (Kauranen et al., 2013), we know that D. montana females do not show cycling amplitude peaks in their diapause response under changing night lengths, and in the present study, our main interest was focused on the effects of selection on females' overall tendency to enter diapause under these conditions. N–H experiments were performed on selection and control line replicates, which had been maintained in LD 15:9 and LL, respectively, in mixed generations in 17°C for 4 to 5 generations after finishing active selection. The study consisted of five separate experiments, performed in LD 12:6, 12:12, 12:24 12:36 or 12:72 (choice of these photoperiods was based on a previous N–H experiment; Kauranen et al., 2013). In addition, a sample of flies was maintained in constant darkness (DD) to see if the females require light input to enter diapause. Experimental flies (87–388 females/each LD cycle or DD; specific information about the sample sizes can be found in Table 2) were transferred into the climate chambers at the pupal stage in bottles containing malt medium, so that the bottles containing flies of one selection and one control line replicate were always in the same chamber. Cabinets were illuminated with compact fluorescent lamps (9 W, Megaman, Germany) with light intensity of 300–1000 lx during the photophase at 16±0.3°C. Flies were kept in the chambers for 21±2 days until they reached maturity and then transferred into a freezer until the diapause state of all females was determined as explained above.

Cold tolerance of D. montana flies increases in response to a decreasing photoperiod in late summer and autumn (photoperiodic cold acclimation; Vesala et al., 2012), and to take this phenomenon into account, we measured the cold tolerance of selection and control line for the flies reared in LL and in LD 19:5, 17:7, 15:9 and 13:11. Parental females were allowed to lay eggs in the given LD conditions in successive 3 h periods and the emerging flies were collected at 12 h intervals, which allowed us to collect material evenly along the modal distributions of flies' pre-adult development time of each replicate. Each of the five samples per line replicate consisted of 16 males and 16 females, which were used in the experiments at the age of 20–27 days. Female reproductive stage was determined after the experiment as before. Flies' cold tolerance was measured using the critical thermal minimum (CT_min) method (Overgaard et al., 2011); the flies were transferred into a water-glycol bath in sealed vials in 16°C, after which temperature was decreased by 0.5°C min⁻¹. The temperature when the flies entered chill coma (i.e. were no longer able to stand on their feet) was recorded as their CT_min. The total number of experiments was 64, and the males and females from different selection regimes and LDs were randomly assigned.

Locomotor activity rhythm of the flies was traced in Trikinetics Drosophila High Resolution Activity Monitors (Waltham, MA, USA) (N=27–32 females and males per replicate; see Table 3) with 9 consolidated infrared beams per tube. Flies were collected 1–2 days after emergence and entrained under LD 19:5 for 6 days, after which their free-running rhythmicity was measured in LL for 8 days at 17°C (see Kauranen et al., 2012).

50 females of each selection and control replicate were frozen in generation F₁₆ and stored at −20°C. DNA was extracted from 10 females at a time and the extractions pooled into 6 samples of 50 flies (one for each replicate of the selection and control line). DNA extraction was made using CTAB solution (Applichem, BioChemica, Germany) with Proteinase K (Qiagen, Germany) treatment, followed by phenol–chloroform–isoamyl alcohol purification with RNAaseA (Qiagen, Germany) treatment. DNA concentration was checked using a Qubit spectrophotometer (Thermo Fisher Scientific, MA, USA) and the purity with TapeStation (Agilent Technologies, Santa Clara, CA, USA). Extracted DNA samples were sequenced at the Edinburgh genomics (Scotland, UK) facility using Illumina HiSeq X platform with 350 bp insert size and by 150 bp pair-end sequencing, achieving a coverage of about 200×.

Reads were trimmed to remove remaining adapter sequences and low-quality bases, using Trimmomatic v.0.32 (Bolger et al., 2014). Because of a general drop-off in DNA quality at the end of reads, all reads were first cropped to 145 bp and then clipped from the leading and trailing edges if base quality fell below 20. Reads were also clipped if the base quality fell below 20 in a sliding window of 5 bp, and finally any reads shorter than 100 bp were removed. Reads were mapped to the D. montana reference genome (Parker et al., 2018) using bwa mem (v.0.7.7) (Li, 2013). Samtools (v.1.3) (Li et al., 2009) was used to keep only properly paired reads and reads with a mapping quality greater than 30 (Schlötterer et al., 2014). Additionally, regions around indels were locally re-aligned with Picard v 2.14.1 (Broad Institute, http://broadinstitute.github.io/picard) and GATK (McKenna et al., 2010; DePristo et al., 2011). Finally, PCR duplicates were removed with samtools.

SNPs were called first with samtools mpileup and a newly developed heuristic SNP caller, PoolSNP (Kapun et al., 2018 preprint). The minimum base quality to consider a read to be used for base calling was 15. An allele was considered if the number of reads supporting it was >4 across all populations. SNPs were called if the coverage was greater than 33 and less than the 95th percentile of coverage (as calculated for each contig separately). In total, 4,297,404 SNPs were called for further analysis. Furthermore, 1,312,706 SNPs could be placed on 4 autosomal linkage groups and the X-chromosome using an available linkage map (Parker et al., 2018).

All statistical analyses were performed with R 3.5.1 (https://cran.r-project.org). For all regression type models mentioned below model selection was done by assessing the best-fitting models with Akaike information criterion using the AIC function and log-likelihood tests using the anova function in R.

For the analysis of the CDL responses the data were made binary, with values of 0 for diapausing female and 1 for non-diapausing female. To determine the CDL for the base population and the selection and control line replicates, the photoperiodic response curves (females' diapause percentages under different LDs) were analyzed by a dose response model using the drm function from the drc R package (Ritz et al., 2015). This package fits a range of models with sigmoid or biphasic distributions, and here we are predicting the proportion of diapausing females as the response that occurs per change (i.e. ‘dose’) in hours of light in the LD cycle. The best-fitting model was chosen using a lack-of-fit test and AIC scores using the mselect and modelFit functions (Ritz et al., 2015). This lack-of-fit test uses an approximate F-test to compare dose-response models with different parameters to a more general ANOVA model to see which is the best fitting (Bates and Watts, 1988). This model was fitted with a three-parameter Weibull function (coded in the function as ‘fct=W2.3’ in R), where upper limit of the model is set equal to 1.

Differences in the percentage of diapausing females in the selection and control line replicates under different photoperiods (12 h light+6, 12, 24, 36 or 72 h dark) or constant darkness (DD) was tested for each LD or DD period with a generalised linear mixed model (GLMM) with a bimodal distribution and a logit link function. Here, selection regime (selection or control) was used as a fixed factor and replicate as a random factor. As the six LD/DD cycles were tested separately, P-values for comparisons between selection and control lines were corrected using Bonferroni correction for multiple testing.

Differences in the CT_min values of the flies were tested using a linear mixed model (LMM) and fitted using maximum likelihood (ML) and lmer function in R. Selection regime (selection or control line), LD (expressed as hours of light) and fly reproductive type (male, diapausing female and non-diapausing female), and their interactions, were used as fixed factors and replicate as a random factor. For interpretation, this model was fitted to a type II ANOVA with P values computed with a Wald test and interpretation aided by interaction plots produced from the LMM model.

The locomotor activity data were analyzed with the ActogramJ program (Schmid et al., 2011; available at https://imagej.net/ActogramJ) to trace the rhythmicity of the flies and to determine the length of the period (τ), i.e. the length of the fly's intrinsic day in LL. Rhythmicity of the flies was determined from females' daily activity rhythm profiles by displaying the data as double-plotted actograms (48 h plots) both under entraining (LD) and free-running conditions (LL) and using the Lomb–Scargle periodogram method with a significance level of 0.05. The power of the Lomb–Scargle periodogram analysis was defined as the amplitude of the peak only for the rhythmic flies from Lomb–Scargle periodogram with significance level P<0.05. Flies that did not survive throughout the experiment (both entrained and free-running conditions) were excluded from the analysis. For more detailed information on the locomotor activity analysis, see Kauranen et al. (2012).

The data were made binary with values of 0 for arrhythmic flies and 1 for rhythmic flies to determine whether selection had affected flies' ability to retain their activity rhythm in constant light (LL). A fly was determined to be rhythmic if the periodogram analysis detected significant periodicity in its activity rhythm measured across several consecutive days (the significance level in the Lomb–Scargle periodogram analysis was adjusted to 0.05). Only the flies possessing a clear entraining rhythm in 19:5 LD were used in the analysis. We analyzed the instance of rhythmicity with generalized linear models (GLM) with bimodal distribution and a logit link function. Here, selection treatment (baseline, selection or control) and sex were used as fixed factors. As there were no significant differences between the models of control and selection line data with or without replicates (likelihood ratio test, d.f.=1, χ²=1.22, P=0.26) or sex (likelihood ratio test, d.f.=1, χ²=2.5, P=0.26), we dropped replicate and sex from the model and compared the selection and control lines to the base line directly. To examine whether selection had affected the length of the period of flies' activity rhythm, the difference in the length of the period between LD and LL treatment was first calculated for each fly. The data were then analyzed with a linear mixed model (LMM) using restricted maximum likelihood estimation (REML), where selection regime (selection or control) and sex were used as fixed factors and replicate as a random factor. Only flies possessing a clear entraining rhythm in 19:5 LD were used in analysis.

We also performed a transcription factor (TF) binding motif enrichment analysis using the AME tool (with default options) in the MEME suite and the Fly Factor Survey (http://mccb.umassmed.edu/ffs/) database of TF binding motifs (Bailey et al., 2009; Buske et al., 2010). We used the regions 30 bp up- and downstream of top SNPs as input.

Finally, we ran functional annotation clustering with DAVID (v.6.8) (Huang et al., 2009a,b) using the Drosophila virilis annotation as a background set. Additionally, phenotypic enrichment analysis was performed with DroPhEA (Weng and Liao, 2011) to test for enrichment of known phenotypes within the phenotype levels 2 through 7.

The percentage of diapausing females varied over the course of selection between 22 and 53% in selection line replicates (it only remained below this in generation F₈), while in control line replicates, respective percentages varied between 0 and 8%. Selection was stopped in generation F₁₅, when more than 50% of the females of each replicate entered diapause under 15:9 LD. CDLs of the base population and the selection and control line replicates were determined from the photoperiodic response curves as the point where 50% of the females enter diapause (Fig. 2). In the base population, CDL was 19.4:4.7 LD, and two generations after the end of selection (generation F₁₇) the mean CDL over the three selection line replicates was 17.5:6.5 LD and the three control line replicates 18.7:5.3 LD. CDL of the selection line was significantly shorter than those of the base population (1.82 h; Est.=1.22, s.e.=0.06, t=20.45, P<0.001) or the control line (1.27 h; Est.=1.19, s.e.=0.03, t=36.09, P<0.001), while the difference between the base population and the control line was not significant (0.55 h; Est.=0.03, s.e.=0.05, t=0.66, P=0.509). Selection had no effect on female diapause propensity under short-day conditions (day lengths of 16 h or shorter), where the diapause percentages were close to 100% in both lines.

Photoperiodic response curves and CDLs of the control and selection line replicates were measured again after 8–10 months of relaxed selection, prior to performing N–H experiments. Here, the CDLs of the three selection line replicates were 17.5:6.5, 17.4:6.6 and 17:7 LD and for the three control line replicates 20.5:3.5, 20.2:3.8 and 20.2:3.8 LD, with a strongly significant difference in the mean dose response of two lines (Est.=2.62, s.e.=0.07, t=36.7, P<0.001).

N–H experiments revealed a clear connection between CDL and females' diapause frequency under the photoperiods that deviate from 24 h cycles: females' tendency to enter diapause was significantly lower in selection line replicates (short CDL) than in control line replicates (long CDL) in all unnatural LDs and DD (P<0.001), but not in 12:12 LD conditions (Table 4, Fig. 3). The difference was particularly large in the shortest scotophase (12:6 LD), where CNL (6 h) was shorter than the CNL of the selection line replicates in critical photoperiod (CDL 17.7; CNL 6.3 h). D. montana females with a long northern-type CDL have also been found to enter diapause under photoperiods that deviate from 24 h cycles in N–H experiments more frequently than those having a short southern-type CDL in wild clinal populations (P.L., unpublished data).

Cold tolerance of the selection and control line males and the non-diapausing females of selection line increased (CT_min value decreased) towards shorter day lengths, while the diapausing females of both lines showed no changes in their cold tolerance and the non-diapausing females of control line even showed a slight decrease in this trait (Fig. 4). Control line males had about 0.5°C lower CT_min values than the selection line males in all LDs, and the cold tolerance of diapausing females in all LDs was as low as that of the males of the same line under the shortest photoperiod (Fig. 4). The replicates explained only 3.77% of variance in flies' cold-tolerance. Changes in CT_min values were significant across both of the LD conditions and between the fly reproductive types (males and non-diapausing and diapausing females; Table 5), while the difference between the selection conditions remained on the borderline of significance (Table 5). Furthermore, there were strongly significant interactions between fly reproductive type and LD and selection regime, but the three-way interaction was not significant (Table 5). Photoperiodic cold acclimation was most clear in males, where the cold tolerance increased by an estimated 0.15°C per 1 h decrease in photoperiod.

The locomotor activity rhythm of the flies was studied under both entraining (19:5 LD) and constant light (LL) conditions, the latter representing flies' free-running rhythmicity. In 19:5 LD, all base population females were rhythmic (males were not studied), while in selection and control line replicates, the percentage of rhythmic flies varied 86.4–100% and 63.9–86.7%, respectively (Table 3). In LL, the selection line flies showed higher free-running rhythmicity than the control line flies (GLM: Est.= 0.62, s.e.=0.28, z= 2.22, P=0.027), while the differences between the base population and the selection (GLM: Est.=0.30, s.e.=0.32, z=0.9, P=0.367) or the control line (GLM: Est.=−0.25, s.e.=0.32, z=−0.75, P=0.452) were not significant (see Table 3).

Rhythmic flies possessed a clear evening-activity peak in 19:5 LD, and most flies also maintained this peak in LL. However, the length of the free-running period (τ) in LL did not show significant differences between the selection and control lines (LMM: Est.=−7.53, s.e.=24.39, t=−0.309, P=0.77) or between the sexes (LMM: Est.=11.41, s.e.=23.54, t=0.48, P=0.63).

Selection and control line replicates showed significant differences (q value <0.05) in allele frequencies of 4496 SNPs (hereafter called ‘divergent SNPs’), 1445 of which could be localized to the main chromosome arms of D. montana using the existing anchored genome. A Manhattan plot of these SNPs (Fig. 5) indicates that instead of being randomly distributed across the main chromosomes, they are enriched on chromosomes 3 and 4 (Table S1; χ²=63,552, d.f.=4, P<0.01). This pattern holds when using the linkage group length as a proportion of the total genome length or the proportion of SNPs out of all SNPs as the expected proportions of divergent SNPs. Thus, there seems to be an excess of divergent SNPs on the 3rd and 4th chromosomes (Table S1).

Allele frequencies across the control and selection lines in the two treatments showed a combination of large allele frequency differences, as well as smaller differences with a very consistent effect between the lines. The regions around top SNPs were enriched for a single TF binding motif, Adult enhancer factor 1 (Aef1). Additionally, top SNPs lay within 10 kb of 1756 (Table S2) D. montana gene models, which have an identifiable ortholog in D. virilis. These were tested for functional annotation clustering with DAVID (v.6.8) (Huang et al., 2009a,b), using the D. virilis annotation as a background set. This found 100 functional clusters, 16 of which were significantly enriched (Table S3). Most of these clusters consisted of factors linked with signaling, cluster 1 involved factors linked with membrane/transmembrane changes, clusters 4, 11 and 14 with ion transport, clusters 9, 10 and 16 with signal transduction and clusters 3, 6 and 8 with neuropeptide/neurotransmitter signaling. The remaining clusters involved factors affecting mainly embryogenesis, differentiation and growth: cluster 2 involved immunoglobulins, cluster 5 epidermal growth factor, cluster 7 homeobox, clusters 12 and 15 protein kinase and tyrosine activity and cluster 13 transcription regulation.

Genes located within 10 kb of significant SNPs are listed in Table S4, and the list of most-divergent 100 genes in Table S5. Comparing the gene list with genes known to be involved in photoreception/signal transduction, circadian rhythms and neurohormonal control of reproductive diapause revealed several SNP-associated genes playing a role at different steps of diapause induction and CDL determination. A large number of genes associated with diverged SNPs were involved in photoreception (pinta, Cpn, w, st, Arr1, crb, cry and Pld) and/or ion channel function (TrapA1, Trpm, Trpγ, norpA, stops, Shab, brv3, dysc and slob). Several genes were also involved in time measurement and activity rhythms through the circadian clock system (e.g. Shab, dysc, cry, sgg, nf1, Fmr1, dh31, Lkr, SIFaR, Oamb, DAT, DopR1, DopR2, Hmgr, AstA-R2 and Inr). Furthermore, at least 12 of the genes were G protein-coupled receptors that are active in neuropeptide signaling and/or neurotransmitter (mainly octopamine and dopamine) synthesis or reception (Oamb, Octbeta3R, Ple, Ddc, DAT, DopR1, DopR2, Hmgrc, AstA-R2, Inr, Pdk1, Ide, sNPF-R, slob, GABA-B-R1, GABA-B-R2, SIFaR and Tor). Several of the latter genes, as well as Inr, Pi3K59F, jhamt, dpp and Jheh1, are known to affect fly development and diapause through insulin signaling and 20-hydroxyecdysone, and on juvenile hormone metabolism. More information on the function of these, as well as some of the key references, is given in Table S6.

The correct timing of reproductive diapause increases female survival rate during the unfavorable season and their progeny production in the following warm period, and thus it is highly adaptive, especially in temperate regions (Hahn and Denlinger, 2007). Evolution of longer CDLs at higher latitudes has been suggested to result from selection from longer summer days (photic environment) on the circadian clock system (Pittendrigh and Takamura, 1987) or directly on the critical photoperiod/CDL (Bradshaw et al., 2003). In our study, divergence of CDL and other diapause-associated traits between selection and control lines occurred without disrupting circadian rhythm in fly locomotor activity, which gives support to the latter view.

Diapause propensity and CDL have been shown to exhibit high genetic variation within populations, and thus these traits can be expected to be rapidly altered by selection (Tauber et al., 1986; Lankinen et al., 2013). Our earlier studies have revealed about 2 h difference between the CDLs of the most northern (67°N) and southern (62°N) D. montana populations in Finland, as well as high variation (up to 3 h) within the Oulanka (66.4°N) population used in the present study (Tyukmaeva et al., 2011; Lankinen et al., 2013). CDL of selection line replicates decreased by 1.18 h, on average, compared with control line replicates over 9 generations of selection, which corresponds to 2.5–5 deg change on a latitudinal scale (Danilevsky et al., 1970; Tyukmaeva et al., 2011). It also suggests that at the flies' home site (66.4°N), the base population and control line females would enter diapause about 2 weeks earlier than the selection line females (July 31st and August 15th, respectively). Similar shifts towards shorter (southern type) critical photoperiod have been detected in the wild as a response to climate warming, for example, in Wyeomyia smithii mosquitoes (Bradshaw and Holzapfel, 2001) and Phratora vulgatissima leaf beetles (Dalin, 2011).

CDLs of the selection and control line replicates were measured for a second time after maintaining the flies in mixed generations in LD 15:9 (selection line) and LL (control line) at 16°C for 4 to 5 generations. Here, the mean CDL of the selection line replicates (17.6) was about the same as that measured in generation F₁₇ (17.7), while that of the control line replicates increased from 18.9 to 20.5 h. The second set of CDLs was measured at a temperature ∼1°C lower than the first set, and thus the two measurements are not strictly comparable (CDL of D. montana becomes longer in lower temperatures; V. Tyukmaeva, personal communication). A plausible explanation for why the effect of decreased temperature is evident in the CDL of the control line, but not in that of the selection line, is that in the latter line the temperature effect has been compensated by continued selection for shorter CDL in LD 15:9. Oikarinen and Lumme (1979) have previously succeeded in shortening CDL of Drosophila littoralis females by 1.1–1.8 h by maintaining the flies for seven generations in LD 18:6, which shows that this kind of selection can be quite effective.

In the N–H experiment, circadian clock oscillations are suggested to be involved in the photoperiodic time measurement, if the insects' short- or long-day responses show amplitude peaks under LD cycles whose period (T) is close to 24 h and its multiples (Nanda and Hamner, 1958). We have previously shown that D. montana females measure night rather than day length for diapause timing and that their diapause frequency shows a peak only when T was close to 24 h (12 h light+12 h dark; Kauranen et al., 2013). In the present study, selection line females showed a clear diapause peak in LD 12:12 (T=24 h), while the control line females with a long ‘northern’ CDL entered diapause at a high frequency under all photoperiods studied. All these findings can be explained by the damping version of the circadian external coincidence model (Lewis and Saunders, 1987), if the damping of circadian oscillator(s) increases towards North (note that in the present study free-running locomotor activity rhythm dampened more quickly in the ‘northern type’ control than in the ‘southern type’ selection line). In this model, females are expected to enter diapause when the night length exceeds their CNL (24 h CDL), which could also explain why the diapause frequency of the selection line females (CNL 6.5), but not that of the control line females (CNL 5.3), dropped drastically in LD 12:6 (CNLs are calculated from the CDLs measured in Oulu University during the same time period and temperature as the N–H experiments). Line differences in diapause frequency in DD resembles a latitudinal cline detected in the diapause frequency of D. littoralis females in this condition (Lumme and Oikarinen, 1977).

D. montana flies have been found to show some variation in their cold tolerance between populations from different latitudes and altitudes (Vesala and Hoikkala, 2011; R.A.W.W., unpublished data), as well as notable photoperiodic cold acclimation towards the cold season (Vesala et al., 2012). In the present study, selection decreased the cold tolerance of males and diapausing females in all photoperiods examined, while the results for non-diapausing females were less clear. Male cold tolerance increased linearly towards shorter photoperiods, while that of the diapausing females showed practically no changes under different LDs, which suggests that female cold tolerance is phenotypically associated with diapause, while male cold tolerance shows clearer photoperiodic cold acclimation. D. montana males are known to enter reproductive diapause (become reproductively inactive) under the same day lengths that induce diapause in females (Ala-Honkola et al., 2018), but the present study shows that this does not increase their cold tolerance. Vesala and Hoikkala (2011) have previously shown that the strength of association between female diapause and cold tolerance varies between D. montana populations, and Vesala et al. (2012) have demonstrated that the females of a D. montana strain, which lacks diapause, are less cold tolerant than those of the other strains of this species. Furthermore, Teets and Denlinger (2013) have argued that seasonal cold adaptation can be a component of overwintering diapause, but the two processes can also occur independently, and Pegoraro et al. (2014) have suggested that differences between diapause and cold tolerance phenotypes may represent two separate photoperiodic circuits that use different genetic networks.

Decrease in CDL and other diapause-associated traits in the selection lines was not accompanied by changes in the period of fly free-running locomotor activity rhythm, which shows that the traits associated with reproductive diapause can occur without disturbing mechanisms of the circadian clock. Another species of the D. virilis group, D. littoralis, has been found to show correlated latitudinal variation in CDL and flies' circadian eclosion rhythm (Lankinen, 1986b), but in this species, the traits also diverged in a selection experiment (Lankinen and Forsman, 2006). Lack of correlation between the circadian clock-regulated traits (pupal eclosion, egg hatch and oviposition rhythms) and photoperiodic responses has been verified also in the pink boll worm moth Pectinophora gossypiella (Pittendrigh and Minis, 1971) and between circadian eclosion rhythm and CDL in Drosophila auraria (Pittendrigh et al., 1984; Pittendrigh and Takamura, 1987). Furthermore, Bradshaw et al. (2003) succeeded in changing both the critical photoperiod and the amplitude, but not the period, of diapause response rhythmicity of W. smithii mosquitoes under extended days in an N–H selection experiment.

Experimental evolution combined with genome sequencing has proved to be an effective approach to study the connection between phenotypic responses and genetic changes during selection (Tobler et al., 2013; Schlötterer et al., 2014; Graves et al., 2017). In total, 4496 SNPs show a significant difference between the control and diapause selection lines, 1445 of which could be placed on linkage groups. There was a significant deviation from a random distribution of SNPs. It seems likely that the number of SNPs involved could be overestimated because of linkage and hitchhiking effects. Unfortunately, little is known of the recombination map in D. montana. In a functional enrichment analysis, we detected 16 functional gene clusters that diverged between selection and control line replicates; 12 of these clusters contained factors involved in signal transduction, including membranes, ion channels and neuropeptide/neurotransmitter signaling. Changes in the function of genes within these clusters could affect signal transduction at any level, from photoreception and the passage of visual and other signals in fly brain to the regulation of various hormones through neurotransmitters. The remaining 4 clusters involved factors affecting embryogenesis, differentiation and growth. Notably, many of above-mentioned gene clusters, especially those in the first group of clusters, have been shown to also differ between wild populations of this species (Parker et al., 2018; R.A.W.W. et al., unpublished data).

A transcription factor motif enrichment analysis found that a motif of the transcription factor Adult enhancer factor 1 (Aef1) was over-represented in the regions surrounding top SNPs. In D. melanogaster, Aef1 is a transcription factor that is a repressor of Adh in the adult fat body (Falb and Maniatis, 1992). Aef1 is also an inhibitor of yolk protein 1 and 2 in males (Tarone et al., 2012) and seems to influence olfactory behavior (Tunstall et al., 2012).

Recent studies in insects have discovered genes and gene networks involved in different steps of diapause induction, including photoreception and transduction (Fowler and Montell, 2013; Kistenpfenning et al., 2018), circadian regulation (Zhang et al., 2010; Dubowy and Sehgal, 2017; Liang et al., 2017) and neurohormonal control of reproductive diapause (Nässel and Winther, 2010; Diniz et al., 2017). Furthermore, Bryon et al. (2013) and Zhao et al. (2017) detected upregulation of several genes encoding G protein-coupled receptors, especially those for octopamine, neuropeptide F, proctolin and tachykinin, during early diapause in the two-spotted spider mite, T. urticae. In our study, several genes associated with SNPs diverging between selection and control line replicates were involved in photoreception, ion channel function and/or in the function and output of the circadian clock system. The list of genes also included G protein-coupled receptors that function in neuropeptide signaling and/or neurotransmitter synthesis or reception. Finally, several genes functioned in insulin signaling and on 20-hydroxyecdysone and juvenile hormone pathways that play a central role in diapause induction. More information on the function of implicated genes, as well as key references, are provided in Table S5. Of course, further studies are required to examine the functional involvement of any of these genes in photoperiodic adaptation, but our results confirm that changes in CDL involve extensive genome-wide changes in genes involved in signaling.

Photoperiodic reproductive diapause has evolved multiple times in arthropod species, including within the genus Drosophila (Hand et al., 2016), and the interaction between PTM and circadian clock also seems to vary between species (Meuti and Denlinger, 2013). D. montana is a northern Drosophila species with a robust reproductive diapause, unimodal daily activity rhythm and an ability to maintain free-running locomotor activity rhythm in constant light, but not in constant darkness (Kauranen et al., 2012). This and other D. virilis group species also differ from a more southern species, D. melanogaster, at the neuronal level, for example, in the expression pattern of the blue light photopigment cryptochrome (CRY) and the neuropeptide pigment-dispersing factor (PDF) (Kauranen et al., 2012; Hermann et al., 2013). D. virilis group species (including D. montana) probably lost CRY and PDF expression in their long (l-LNvs) and small (s-LNvs) ventrolateral clock neurons, respectively, but gained PDF expression in their central brain, under selection from environments subject to major photoperiodic changes throughout the year with extremely long photoperiods in summer (Menegazzi et al., 2017; Beauchamp et al., 2018). Beer et al. (2017) give Acyrthosiphon pisum aphids and northern D. virilis group species, D. montana, Drosophila ezoana and D. littoralis, as examples of the species exhibiting weak circadian rhythms and robust seasonal responses to shortening photoperiod. Weak clocks have been suggested to be generally more plastic than strong clocks, which may partly explain why they can easily adapt to extreme changes in day length that prevail in the North (Abraham et al., 2010). Even though anticipating daily reoccurring events is an advantage, an oscillator that drives strong circadian rhythms over a long time may be potentially less flexible in resetting to seasonal changes (Beer et al., 2017).

Short selection experiments in the laboratory differ in many fundamental ways from the long-term selection occurring in the wild, but nevertheless, our study revealed important aspects of life history evolution. First, it showed that adaptive divergence of CDL and other diapause-related traits can occur without disturbing circadian clock-regulated locomotor activity rhythm, and that CDL also shows plasticity for reproduction under short-day conditions. Second, resequencing the selection and control line replicates revealed widespread genomic divergence, including significant divergence in SNPs located in or near several genes important for photoreception, circadian rhythms, signal transduction, aminergic signaling and hormone secretion. It would be interesting to find out how these genes function under different LDs and temperatures, for example, through alternative splicing, and whether they show latitudinal variation in D. montana populations. Such information would help us to trace the steps leading to changes in CDL and other photoperiodically regulated traits in the wild and to understand how adaptation to new or changing environmental conditions occurs at the genetic level. As Bergland et al. (2014) have stated, environmental heterogeneity can promote the long-term persistence of functional polymorphisms within populations that can fuel a fast directional adaptive response in the future.

Genome sequencing was performed at Edinburgh genomics (Scotland, UK). We are grateful to Professor David S. Saunders for his valuable comments on an earlier version of this manuscript.