Multiple paths to cold tolerance: the role of environmental cues, morphological traits and the circadian clock gene vrille

Abstract

Background

Tracing the association between insect cold tolerance and latitudinally and locally varying environmental conditions, as well as key morphological traits and molecular mechanisms, is essential for understanding the processes involved in adaptation. We explored these issues in two closely-related species, Drosophila montana and Drosophila flavomontana, originating from diverse climatic locations across several latitudes on the coastal and mountainous regions of North America. We also investigated the association between sequence variation in one of the key circadian clock genes, vrille, and cold tolerance in both species. Finally, we studied the impact of vrille on fly cold tolerance and cold acclimation ability by silencing it with RNA interference in D. montana.

Results

We performed a principal component analysis (PCA) on variables representing bioclimatic conditions on the study sites and used latitude as a proxy of photoperiod. PC1 separated the mountainous continental sites from the coastal ones based on temperature variability and precipitation, while PC2 arranged the sites based on summer and annual mean temperatures. Cold tolerance tests showed D. montana to be more cold-tolerant than D. flavomontana and chill coma resistance (CT_min) of this species showed an association with PC2. Chill coma recovery time (CCRT) of both species improved towards northern latitudes, and in D. flavomontana this trait was also associated with PC1. D. flavomontana flies were darkest in the coast and in the northern mountainous populations, but coloration showed no linkage with cold tolerance. Body size decreased towards cold environments in both species, but only within D. montana populations largest flies showed fastest recovery from cold. Finally, both the sequence analysis and RNAi study on vrille suggested this gene to play an essential role in D. montana cold resistance and acclimation, but not in recovery time.

Conclusions

Our study demonstrates the complexity of insect cold tolerance and emphasizes the need to trace its association with multiple environmental variables and morphological traits to identify potential agents of natural selection. It also shows that a circadian clock gene vrille is essential both for short- and long-term cold acclimation, potentially elucidating the connection between circadian clock system and cold tolerance.

Background

Species geographical distribution is largely defined by their ability to tolerate stressful conditions and to respond to daily and seasonal temperature changes (e.g. [1,2,3,4,5,6]). Accordingly, ectothermic species, especially the ones living at high latitudes, may use a variety of physiological and behavioural strategies to increase their cold tolerance [7, 8]. To understand species’ cold adaptation in these regions, it is essential to identify latitudinal and bioclimatic selection pressures driving these adaptations [2, 9,10,11] and to trace associations between cold tolerance and adaptationally important morphological traits. Also, it is also important to discover molecular mechanisms underlying cold tolerance, but a functional link between candidate genes and cold tolerance has relatively rarely been established (but see e.g. [12, 13]).

Insect cold tolerance is often measured with two ecologically relevant methods. Insects’ chill coma resistance can be assessed by measuring their critical thermal minimum (or chill coma temperature, CT_min), where they lose neuromuscular coordination and fall into chill coma during a gradual decrease in temperature [14]. Chill coma recovery time (CCRT), on the other hand, is based on time required for an individual to recover from the coma after removal of the chill coma inducing temperature [15, 16]. Even though the ability to resist chill coma and recover from it are thought to be somewhat of a continuum, they are affected by different physiological and molecular mechanisms (e.g. [17, 18]). The onset of chill coma is associated with depolarisation of muscle resting membrane potential due to low temperature (e.g. [19, 20]), while the recovery process involves energetically costly restoration of ion gradients and upregulation of genes involved in repairing cold injuries [21,22,23,24,25]. Both the high resistance to chill coma and the fast recovery from it have several implications on insect fitness in cold environments [15, 18, 26]. For example, fast recovery from chill coma after winter, night or otherwise sudden temperature drops ensures finding high-quality mates on feeding and breeding sites, as well as fast escape from predators. High chill coma resistance covers largely the same advantages, but it may also save insects from falling into energetically costly chill coma. Thus, these traits are likely to be under strong selection.

Latitudinal clines in cold tolerance traits, detected in a variety of species, including Drosophila flies [16, 27, 28], Myrmica ants [29], Porcellio laevis woodlouse [30], highlight adaptive genetic variation in these traits. However, it is challenging to distinguish whether such variation has evolved in response to changes in photoperiod (day length), temperature or their combination [7]. One possible approach is to use GIS (geographic information system)-based environmental data to trace correlations between evolutionary changes in the studied traits and the environments the populations or species experience [31]. Studies on clinal variation in insect cold tolerance can also be complicated by the fact that many species spend the coldest time of the year in reproductive diapause, which induces various kinds of changes in their physiology, metabolism and behaviour in addition to increasing their cold tolerance [32]. Thus, it is important to investigate genetic variation in cold tolerance in temperature and light conditions, where insects’ reproductive status is controlled.

In addition to cold tolerance, morphological traits, like body colour (the degree of melanism) and body size, often show latitudinal variation. Thus, to understand insect adaptation to cold environments, it would be important to investigate cold tolerance traits together with other adaptationally important traits and consider correlations between them. In insects, an increase in melanism towards higher latitudes has been detected both between and within species [33,34,35]. According to thermal melanin hypothesis, this can be explained by an increased ability of dark individuals to absorb solar radiation and warm up fast in cold environments with low solar radiation (reviewed in [36]). However, body colour may also be affected by other selection pressures, including protection against UV-radiation [37, 38], desiccation [39, 40] and pathogens [41]. Body size is one of the most important quantitative traits associated with insect metabolism, fecundity, mating success and stress tolerance, and thus it is likely under several selection pressures in nature (e.g. [42, 43]). Body size of ectothermic species has been either shown to increase (Bergmann’s rule) or decrease (converse Bergmann’s rule or U-shaped cline) towards higher latitudes and cooler climates (e.g. [42, 44, 45]). Converse Bergmann’s rule applies especially for insects with a long generation time, for which the short growing season on high latitudes limits the time available for development, growth and foraging [45,46,47]. Overall, parallel clines in cold tolerance and morphological traits give only indirect evidence on the functional linkages between these traits. One way to obtain more direct evidence on their association would be to measure cold tolerance and morphological traits for the same individuals. For example, it has been shown that melanistic wood tiger moths (Parasemia plantaginis) warm up more quickly than the less melanistic ones [48] and Drosophila montana males’ overwinter survival increases along with an increase in body size in nature [49].

At high latitudes with clear seasonal and diurnal temperature variation, insects can enhance their cold tolerance through long- and short-term cold acclimation [50, 51]. In nature, insects can anticipate the forthcoming cold season from a decreasing temperature and/or day length and adapt to these changes through a gradual increase in cold tolerance [52, 53]. Insects can also be experimentally cold-acclimated by maintaining them in relatively low temperature or short day conditions for a few days to weeks prior to cold shock [54, 55]. Short-term cold acclimation of minutes to hours (rapid cold hardening, RCH) has been suggested to share mechanisms with longer-term acclimation [56] and allow organisms to cope with sudden cold snaps and to optimize their performance during diurnal cooling cycles [24, 53]. Here the circadian clock system, which monitors daily and seasonal light and temperature cycles and entrains behavioural and physiological rhythms to match with them [57], could play a vital role. In the central circadian clock, described most thoroughly in D. melanogaster, daily rhythms are driven via one or several transcriptional feedback loops involving changes in the expression of core circadian clock genes (reviewed in [58,59,60]). Intriguingly, in Drosophila montana flies cold tolerance traits (CT_min and CCRT) are under photoperiodic regulation [55, 61], and hence could be associated with circadian clock system. Also, the expression levels of circadian clock genes have been found to change during long-term cold acclimation e.g. in several Drosophila species [62,63,64,65] and Gryllus pennsylvanicus cricket [66]. In plants, the linkage between circadian clock and cold acclimation has been established [67, 68], but in insects direct evidence on this link and its extent is still missing.

Northern Drosophila virilis group species possess a very high cold tolerance compared to most other species of the genus [1, 2]. Our study species, D. montana and D. flavomontana, belong to this group, and they live in diverse climatic conditions in the low-altitude western coast and in the high-altitude Rocky Mountains of North America across several latitudes (Fig. 1). D. montana lives generally at higher latitudes and altitudes than D. flavomontana, but in some sites the species occur sympatrically [69,70,71]. Moreover, the body colour of D. montana is almost black, while that of D. flavomontana varies from light to dark brown [69]. These features make D. montana and D. flavomontana an ideal species pair to investigate genetic adaptation to cold environments, complementing latitudinal studies performed on less cold-tolerant southern species. They also enable us to study connections between fly cold tolerance and morphological traits and to trace selection pressures underlying cold adaptation. Another interesting point is that the circadian clock system of D. virilis group species differs from that of D. melanogaster, and shows features that have helped them to adapt to high latitudes [72,73,74,75,76]. Here vrille, which is one of the core genes in the central circadian clock system and a key regulator of circadian behavioural rhythms in D. melanogaster [58,59,60], is of special interest, as it has been found to be highly upregulated during cold acclimation in D. montana [64, 65, 77].

Fig. 1

Fly collection sites. Table shows fly collecting sites and years, and the coordinates for each site. Map contains information on whether we have samples from one or both species in each site on the western coast and in the Rocky Mountains (brown area on the map) in North America (detailed information given in Additional file 1: Table S1). The map template obtained from https://d-maps.com/carte.php?num_car=5082&lang=en

In this study, we addressed two main objectives. (1) We traced genetic variation in cold tolerance traits (CT_min and CCRT) across distribution ranges of D. montana and D. flavomontana, and identified environmental selection pressures and morphological (body colour and size) correlations driving these adaptations. (2) We compared vrille sequence variation to cold tolerance traits among clinal populations of D. montana and D. flavomontana, and tested the role of this gene in cold tolerance and cold acclimation ability in D. montana using RNA interference (RNAi). In the first part of the study, we quantified environmental variation across fly collecting sites by performing a principal component analysis (PCA) on several bioclimatic variables, and investigated whether latitude, as a proxy of photoperiod, and latitudinally or locally varying climatic conditions have shaped fly cold tolerance (CT_min and CCRT). Here we predicted that genetic variation in cold tolerance traits shows association with latitudinally varying photoperiods, which serves as the most reliable cue for seasonal temperature changes at given localities [7]. We also quantified genetic variation in fly body colour and body size across species’ distribution range, identified the likely selection pressures affecting these traits and measured direct correlations between the morphological and cold tolerance traits from the same individuals. Body colour could be associated with cold adaptation if it shows latitudinal cline, becoming darker towards northern latitudes, and if it correlates directly with cold tolerance measures (thermal melanin hypothesis, [36]). If large body size and better cold tolerance are associated with each other (e.g. [49]), the size could be expected to increase towards cold environments (Bergmann’s rule, [47, 78]) and the traits should be directly correlated. We studied these adaptations in non-diapausing individuals in a common environment to eliminate plastic responses. In the second part of the study, we evaluated the adaptive role of vrille in fly cold tolerance by comparing its nucleotide and amino acid variation with variation in two cold tolerance traits (CT_min and CCRT) among D. montana and D. flavomontana cline populations. We also investigated the role of this gene in regulating females’ cold tolerance traits and cold acclimation ability by silencing it with RNAi in D. montana. We expected vrille to be involved at least in females’ ability to become cold-acclimated, as it shows expression changes during cold acclimation in D. montana [64, 77].

Results

Variation in the climatic conditions at fly collecting sites

We investigated macroclimatic variability among the sites, where D. montana and D. flavomontana strains originated from, by performing a principal component analysis (PCA) on 19 bioclimatic variables, growing season length (days) and altitude (Additional file 1: Tables S2, S3). PCA revealed four principal components (PCs) with eigenvalues > 1 (Additional file 1: Table S4). The first two PCs explained more than 83% of the total variation (Fig. 2; Additional file 1: Table S4) and were included in further analyses.

Fig. 2

Climatic conditions in the fly collecting sites. Principal component analysis (PCA) was performed on 21 variables describing environmental conditions in fly collecting sites (Additional file 1: Tables S2, S3). Big black arrows show the change in given conditions

PC1 clearly separated the Rocky Mountains sites from the ones on the western coast. Variables with the highest contribution on this separation included altitude and the ones describing daily and seasonal temperature variation (bio2, bio4, bio7), the minimum temperature of the coldest month (bio6), and the mean temperature of the coldest quarter (bio11) and precipitation (Fig. 2; Additional file 1: Table S5). Together they showed the high-altitude Rocky Mountains sites to have higher temperature variation and colder winters than the ones on the western coast sites. On the other hand, the western coast sites had higher precipitation throughout the year than the Rocky Mountains sites.

PC2 arranged the sites on the basis of the growing season length, annual mean temperature (bio1), the mean temperature of the warmest quarter (bio10), the maximum temperature of the warmest month (bio5) and isothermality (bio3, i.e. how large day-to-night temperatures oscillate relative to the summer-to-winter oscillations; Fig. 2; Additional file 1: Table S5). They showed that while some fly collecting sites differ in latitude (photoperiod), they resemble each other in growing season length and summer temperatures due to their variability in altitude and closeness to sea.

The effects of latitude, climatic conditions and morphological traits on fly cold tolerance

Studying the effects of photoperiod, temperature-related factors (PC1, PC2), body colour and size (measured as weight) on fly cold tolerance enabled us to identify selection pressures affecting this trait. In these experiments, we used three isofemale strains for majority of the populations, but for four populations we succeeded to establish only one or two strains due to the rarity of the species at the collection sites (Additional file 1: Table S1). All traits were measured for non-diapausing D. flavomontana and D. montana females, reared in constant light at 19 °C, to minimize plastic responses. Contrary to other traits, weight was analysed only for the flies used in CCRT tests (see ‘Methods’). The simplest model, which enabled us to distinguish between latitudinally varying photoperiods and temperatures, included latitude and PC2 as explanatory factors (see Fig. 2; Additional file 1: Table S6). The more complicated models included macroclimatic conditions varying between the western coast and the Rocky Mountains (PC1), different interaction terms and fly body colour and size (Additional file 1: Table S6).

The best-fit models explaining the cold tolerance (CT_min = chill coma temperature, CCRT = chill coma recovery time), body colour and body size of D. flavomontana and D. montana females are presented in Fig. 3A–D (model comparisons are given in Additional file 1: Table S6). The models show that the selection pressures driving the evolution of these traits vary between the species. Most of the pairwise correlations between fly cold tolerance traits, colour and weight were non-significant, and were not included in the best-fit models; Figs. 3, 4; Additional file 1: Table S6).

Fig. 3

Cold tolerance, morphological traits (body colour and body size) and environmental conditions. The best-fit models are presented in tables. Relationship between latitude (as a proxy of photoperiod) and PC2 (as a proxy of latitudinally varying temperature) and A chill coma temperature (CT_min), B chill coma recovery time (CCRT), C fly body colour (measured as colour intensity) and D fly body size (measured as weight) in D. flavomontana (D. fla) and D. montana (D. mon) populations. The effects of local climatic conditions (PC1; see Fig. 2) on the western coast and in the Rocky Mountains are illustrated in grey scale (lighter colours represent the western coast populations and darker ones the Rocky Mountains populations). Error bars represent bootstrapped 95% confidence intervals (Mean ± CI). Significant regression lines for latitude or PC2 with standard errors are shown. Significance levels were obtained from GLMMs or LMM: ^NS non-significant, *P < 0.05, **P < 0.01 and ***P < 0.001

Fig. 4

Correlations between fly cold tolerance and morphological traits (body colour or body size). Correlation between A D. flavomontana CT_min and body colour, B D. flavomontana CCRT and body colour, C D. flavomontana CCRT and body size (measured as weight), D D. montana CT_min and body colour, E D. montana CCRT and body colour, F D. montana CCRT and body size (measured as weight). Error bars represent bootstrapped 95% confidence intervals. Significant regression lines are shown with standard errors (best-fit models are presented in Fig. 3 and Additional file 1: Table S6)

CT_min of D. flavomontana showed only low variation and was not significantly explained by any of the variables (Fig. 3A; Additional file 1: Table S7). However, in D. montana this trait showed significant association with PC2 (Fig. 3A; Additional file 1: Table S7), which suggests that the chill coma resistance of D. montana flies is highest in northern populations with a short growing season and cold summer and low annual mean temperatures (Fig. 2). CCRT tests showed D. flavomontana flies’ cold tolerance to be significantly associated with latitude and to improve towards North (Figs. 2, 3B; Additional file 1: Table S7). Moreover, this trait was affected by PC1, especially on latitudes around 50–55° N, suggesting that fly cold tolerance is higher in the humid, low-altitude western coast populations than in the high-altitude Rocky Mountains populations with colder temperatures and higher temperature variation (Figs. 2, 3B; Additional file 1: Table S7). CCRT of D. montana was significantly associated with latitude, improving towards North (Figs. 3B; Additional file 1: Table S7). Moreover, in this species large flies recovered faster from chill coma than the small ones (Figs. 3B, 4G; Additional file 1: Table S7).

Drosophila flavomontana flies’ body colour was significantly affected by latitude, PC1 and an interaction between PC2 and PC1 (Fig. 3C; Additional file 1: Table S7). In the Rocky Mountains, the flies became darker (their colour intensity decreased) towards North, while in the western coast populations flies were equally dark and darker than the ones from the Rocky Mountains populations on similar latitudes. D. montana body colour showed only minor population-level variation and no significant association with latitude, PC2 or PC1. The body size of D. flavomontana was significantly associated with PC1 and an interaction between PC1 and PC2 (Fig. 3D; Additional file 1: Table S7), increasing towards warmer winters and summers and longer growing season. The body size of D. montana was significantly associated with PC2, being highest in sites with high summer temperatures and a long growing season (Fig. 3D; Additional file 1: Table S7).

Association between the nucleotide and amino acid variation in vrille and variation in fly cold tolerance traits among clinal D. montana and D. flavomontana populations

To elucidate the adaptive potential of vrille in cold tolerance among D. montana and D. flavomontana populations, we traced the association between pairwise differences in vrille sequence variation (the number of nucleotide and amino acid differences) and differences in fly cold tolerance traits among the studied populations. The number of amino acid differences and the differences in mean CT_min in D. montana populations were significantly correlated (Mantel test: r = 0.60, P = 0.007; Fig. 5B; Additional file 1: Table S8), i.e. the greater the amino acid differences were between D. montana populations, the greater were the differences in their chill coma resistance. Other measured correlations were non-significant (Fig. 5; Additional file 1: Table S8).

Fig. 5

The association between vrille sequence variation and cold tolerance among D. montana and D. flavomontana populations. A Pairwise mean CT_min differences and the number of pairwise nucleotide differences among D. montana and D. flavomontana populations. B Pairwise mean CT_min differences and the number of pairwise amino acid differences among D. montana and D. flavomontana populations. C Pairwise mean CCRT differences and the number of pairwise nucleotide differences among D. montana and D. flavomontana populations. D Pairwise mean CCRT differences and the number of pairwise amino acid differences among D. montana and D. flavomontana populations

Daily rhythm of vrille and the expression levels of vrille after RNAi in D. montana

Since the expression levels of circadian clock genes are known to fluctuate during the day, we first investigated daily rhythm in vrille expression in LD (light:dark cycle) 18:6, where the females of the study population of D. montana (Seward, Alaska, USA) can be expected to be non-diapausing (and which was also verified, see methods). Our results showed that vrille has a clear daily rhythm and that its highest expression corresponds to ZT16 and ZT20 (ZT = Zeitgeber Time which refers to the number of hours after the lights are switched on; see Fig. 6A). From these two time points we chose to perform RNAi-injections at ZT16, when the lights were still on.

Fig. 6

vrille expression during the day and after RNA interference (RNAi). A Relative normalised expression of naturally cycling vrille at six time points, starting at ZT0 (Zeitgeber Time = sampling time every 4 h over a 24 h period), in LD 18:6 and 19 °C. Yellow represents the time of the day when the lights were on and dark grey when the lights were off. B Relative normalised expression levels of vrille in LacZ-injected and no-injection control females, and in females injected with dsRNA targeting on vrille 48 h after the injections at ZT16. Error bars represent bootstrapped 95% confidence intervals. Significance levels were obtained from linear model and only significant observations are shown: *P < 0.05, **P < 0.01 and ***P < 0.001

We next compared vrille expression levels in LacZ-injected control females to those of vrille-injected and no-injection females 12, 24 and 48 h after the injections (Additional file 1: Fig. S1, Fig. 6B). This enabled us to measure the effects of vrille-RNAi on the expression level of this gene controlling possible effects of immune responses and physical damage. Differences between LacZ- and vrille-injected females were most pronounced 48 h after the RNAi-injections (Fig. 6B; Additional file 1: Fig. S1), where vrille-injected females had approximately 56% lower vrille expression compared to LacZ controls (Fig. 6B; Additional file 1: Table S9). Accordingly, all RNAi experiments were performed 48 h after the injections at ZT16 in LD 18:6.

The effects of vrille-RNAi on D. montana females’ cold tolerance and cold acclimation ability

The effects of vrille-RNAi on female cold tolerance and cold acclimation ability were studied by quantifying these traits in vrille-injected, LacZ-injected and no-injection females. Comparisons between LacZ- and vrille-injected females enabled us to determine whether a expression of vrille increases D. montana females’ ability to resist chill coma (CT_min), to recover from it (CCRT) and/or to achieve higher cold-tolerance after cold acclimation. Comparisons between LacZ-injected and no-injection females, on the other hand, revealed possible immune responses to dsRNA and/or physical damage caused by the injection itself. All cold tolerance tests were started 48 h after the injections at ZT16, when the effects of vrille-RNAi were at a highest level (Fig. 6, Additional file 1: Fig. S1). Cold tolerance tests (CT_min and CCRT) were conducted the same way as for the females of cline populations, except that now the females were maintained in LD 18:6 instead of LL throughout the experiments (also during the acclimation treatment). Moreover, to measure females’ cold acclimation ability, half of the females (cold acclimation group) were maintained in 6 °C and the other half (non-acclimation group) in 19 °C for 5 days prior to cold tolerance tests. RNAi-injections were performed 2 days (48 h) before finishing acclimation treatment and performing the tests.

We first investigated whether cold acclimation in 6 °C (cold-acclimated females) had improved female cold tolerance compared to the flies kept in 19 °C (non-acclimated females) within the three experimental groups (LacZ-injected females, no-injection females and vrille-injected females). In CT_min tests, cold acclimation improved the chill coma resistance in LacZ-injected and no-injection females by decreasing their CT_min on average by 0.7 °C and 0.6 °C, respectively, while acclimation had no significant effect on the chill coma resistance of vrille-injected females (Fig. 7A; Additional file 1: Table S10). In CCRT tests, cold-acclimated no-injection females recovered from chill coma on average 1.5 min faster than the non-acclimated ones (Fig. 7B; Additional file 1: Table S10), i.e. their ability to recover from coma was faster, as could be expected. However, cold acclimation had no significant effect on the recovery time of LacZ-injected females, which suggests that either immune responses or physical damage had overridden the positive effects of the acclimation (Fig. 7B; Additional file 1: Table S10). Finally, CCRT tests for vrille-injected females showed that cold acclimation had slowed down their recovery time by ~ 3.5 min instead of fastening it (Fig. 7B; Additional file 1: Table S10). Such a significant effect in females’ ability to recover from chill coma cannot be explained solely by immune responses or physical damage, which means that expression of vrille is essential for females’ cold acclimation ability.

Fig. 7

D. montana cold tolerance after RNA interference (RNAi). A Chill coma resistance (CT_min) and B chill coma recovery time (CCRT) for LacZ-RNAi, no-injection, and vrille-RNAi females that were kept at + 19 °C for the whole time, or at first at + 19 °C and the last 5 days at + 6 °C (cold acclimation period). Dashed lines indicate significant effects of the acclimation in each treatment, and solid lines significant differences between the LacZ control and the other treatments in females that were or were not acclimated. Significance levels were obtained from GLMMs and only significant observations are shown: *P < 0.05, **P < 0.01 and ***P < 0.001. Numbers below boxplots refer to sample sizes. Whiskers represent ± 1.5 × IQR

Next, we compared CT_min and CCRT of LacZ-injected females to those of the no-injection and vrille-injected females separately among the non-acclimated and cold-acclimated females. In CT_min tests, chill coma resistance did not differ significantly between LacZ-injected and no-injection females in either acclimation groups (Fig. 7A; Additional file 1: Table S11). On the other hand, vrille-injected females entered chill coma in 1 °C (non-acclimated females) and 2.5 °C (cold-acclimated females) higher temperature than LacZ-injected females of the same groups (Fig. 7A; Additional file 1: Table S11). These results show that low expression levels of vrille significantly decreased females’ ability to resist low temperatures. In CCRT tests, both the non-acclimated and cold-acclimated LacZ-injected females recovered from chill coma ~ 3 min more slowly than the respective no-injection females (Fig. 7B; Additional file 1: Table S11), which suggests that immune and/or injection effects might have played a role in chill coma recovery in both groups. The recovery time of non-acclimated vrille-injected females was on the same level as that of the LacZ-injected females (Fig. 7B; Additional file 1: Table S11). However, cold-acclimated vrille-injected females recovered from chill coma ~ 3.5 min more slowly than LacZ-injected females (Fig. 7B; Additional file 1: Table S11), which again brings up the poor acclimation ability of vrille-injected females.

Discussion

Numerous studies have found inter- and intraspecific latitudinal variation in insect cold tolerance (e.g. [1,2,3,4, 6, 16, 27,28,29,30]) and identified candidate genes for it (e.g. [62,63,64, 66]). Spatially varying selection on insect fitness along latitudinal clines can be based on photoperiod and/or climatic factors, and sometimes it is difficult to detect the actual selection pressures driving the evolution of cold tolerance. To deepen our understanding on clinal adaptation, it is important to gather several independent sources of evidence, including sibling species, multiple populations or geographic regions and environmental correlations that account for population structure [79], as well as to investigate associations between cold tolerance and morphological traits (e.g. body colour and size) that often show parallel clines. Here, we explored factors affecting cold tolerance of non-diapausing females of D. montana and D. flavomontana originating from diverse climatic environments across different latitudes on the western coast and the Rocky Mountains of North America. Moreover, tracing sequence variation in both species and performing RNAi experiments in D. montana on one of the key regulators of the circadian behavioural rhythms, vrille, enabled us to investigate the role of this gene on female cold tolerance and cold acclimation ability.

Effects of latitude and bioclimatic variables on fly cold tolerance vary even between closely-related species

The principal component analysis (PCA), which we performed on macroclimatic variables on fly collection sites, showed PC1 to separate the low-altitude coastal sites from the high-altitude mountainous ones based on differences in winter temperatures, temperature variability and precipitation. PC2 further arranged the clinal populations according to their summer and annual mean temperatures. Thus, in our dataset latitude represents photoperiodic differences (day length) between the sites, while PC2 corresponds to latitudinally varying mean temperatures and PC1 to macroclimatic variation between the coast and mountains.

Overall, genetic variation was lower in chill coma resistance (CT_min) than in chill coma recovery time (CCRT), and CT_min showed significant variation only in D. montana. There are some plausible explanations for this pattern. CT_min may not be as suitable as CCRT for studying insects’ inherent cold tolerance, as rapid cold hardening (RCH) during the gradual cooling period in CT_min tests can change the composition of membrane phospholipids, improving the chill coma resistance [19, 20, 80, 81]. A trade-off between inherent cold tolerance and RCH, detected in several studies [26, 56, 82, 83], could lead to skewed or unapparent population-level variation in CT_min. For example, sub-Antarctic Ectemnorhinus weevils and Glossina pallidipes tsetse flies show low inherent genetic variation and greater plastic responses in CT_min between populations [84, 85]. In contrast, CCRT has high heritability [27, 56, 86,87,88], and fast recovery from chill coma may be under strong selection due to energetically costly repair of cold injury and high fitness advantages involved in it (e.g. earlier recovery ensures high-quality territories, mates, feeding and breeding sites and escape from predators) [15, 18, 21,22,23,24,25]. These characteristics make CCRT a good indicator of climatic adaptation. Latitudinal variation in CCRT, but not in CT_min, has also been observed in two widespread ant species, Myrmica rubra and Myrmica ruginodis [29]. Although it is possible that the low variation in cold tolerance in our study is a consequence of relatively low number of studied strains or populations, a recent study found similar trends when using mass-bred D. montana populations for a limited number of places from North America [89].

We predicted that cold tolerance traits have evolved in response to photoperiod, because it serves as a more reliable cue for seasonal temperature changes than environmental temperature itself [7]. In contrast to our expectations, D. montana CT_min was associated with latitudinally varying summer and annual mean temperatures (PC2) instead of photoperiod (latitude). This finding could indicate that the non-diapausing D. montana flies have adapted to cope with sudden temperature drops and diurnal temperature cycles at their home sites during the late spring, summer and early autumn [53]. CCRT of both species, on the other hand, was associated with photoperiod (latitude), as hypothesized. In D. montana, northern females recovered faster from chill coma than the southern ones independent of local climatic conditions (PC1 or PC2), while in D. flavomontana CCRT was also associated with macroclimatic conditions varying between the coastal and mountain sites (PC1). Surprisingly, D. flavomontana females from the colder high-altitude mountain populations recovered more slowly from chill coma than the ones from the humid low-altitude coastal populations. The most logical explanation for this is that high daily and seasonal temperature variation in the mountains create opposing selection pressures, which select for plastic responses. Accordingly, the flies could compensate their low inherent cold tolerance with high cold acclimation ability, as observed in several studies [26, 56, 82, 83, 90, 91]. For example, non-diapausing D. montana females from the coast have been shown to recover faster from chill coma than those from the mountains, but diapausing females from both regions recover at similar rates [91]. Mountain D. flavomontana females could also enhance their cold tolerance by entering reproductive diapause at earlier time of the year than females from the coastal populations on the same latitude, like D. montana females do [92]. Finally, the mountain flies could occupy lower mountain slopes, or thick snow cover could protect them during the cold season.

Fly body colour and size and the effects of these traits on fly cold tolerance

Dark cuticle pigmentation (high melanism) can be expected to offer an advantage in cold environments, as it increases insects’ ability to absorb solar radiation and enables them to warm up faster in cold environments (thermal melanin hypothesis, [36]). This assumption has received support e.g. from latitudinal variation in the degree of melanism in 473 European butterfly and dragonfly species [33]. However, clinal variation in body colour can be induced also by selection favouring light individuals in the south due to their higher avoidance of UV radiation (UV protection hypothesis, [37, 38]). Overall, body colour genes are highly pleiotropic and can play a role also in a variety of other processes, including immunity, camouflage and mate choice [93, 94]. In our study, D. montana flies were darker and more cold-tolerant than D. flavomontana flies, which could be expected as D. montana is found on higher latitudes and altitudes than D. flavomontana [69, 70]. Body colour of D. montana showed no significant variation between populations, while that of D. flavomontana showed two trends. In the Rocky Mountains, D. flavomontana flies became darker towards North, as predicted by thermal melanin hypothesis [36], but the lack of association between fly colour and cold tolerance reduced the support to this theory (see [36]). Rocky Mountains cline could also be explained by UV protection hypothesis [37, 38], if the light flies collected from the southern populations proved to have higher UV resistance. The second trend, detected in D. flavomontana body colour, was increased melanism in the coastal populations. Also this trend could be explained by UV protection hypothesis, as the flies of the misty coastal populations receive less UV radiation than the ones living on high mountains. Desiccation resistance, linked with high melanism in many systems [39, 40, 95], is not likely to play an important role in the formation of either trend, as dark D. flavomontana individuals inhabit humid rather than arid regions. Sexual selection could play part in D. flavomontana colour variation, as cuticular hydrocarbons (CHCs), which can be regulated by pigmentation genes [96], are under sexual selection in this species [71]. Moreover, the darkness of D. flavomontana flies throughout the western coast could be explained by a founder effect, since this species has only recently distributed to the coastal area through British Columbia in Canada, where its body colour is dark [69, 71]. Finally, D. flavomontana may hybridise with D. montana to some degree [69,70,71], which could potentially have led to introgression of dark body colour from D. montana to D. flavomontana.

We hypothesized that body size (measured as weight) of our study species increases towards cold environments like e.g. in Drosophila serrata [27], consistent with Bergmann’s rule [45,46,47], and that the body size is correlated with cold tolerance traits. Contrary to our expectations, the body size of both species was largest in sites with warm summers and winters, which gives support for the converse Bergmann's rule [45,46,47]. This likely results from the short growing season in cold regions and long generation time of D. montana and D. flavomontana (~ 7 weeks) when the time available for development, growth and foraging is limited, resulting in smaller individuals [45,46,47]. However, D. montana body size is clearly affected by multiple selection pressures. While the body size of these flies decreased towards cold environments at population level, the large individuals recovered faster from chill coma than the smaller ones within populations. The latter finding is consistent with a previous study, where the overwinter survival of D. montana males was found to increase along with an increase in body size in nature [49]. In D. flavomontana, the largest individuals came from the southern Rocky Mountains and from the coastal region with warm summers and mild winters. Lack of correlation between body size and cold tolerance, detected in this species, resembles the situation in several other insect species [97,98,99,100].

Circadian clock gene vrille plays an essential role in D. montana females’ short- and long-term cold acclimation

Insects’ circadian clock system monitors changes in daily light and temperature cycles and entrains behavioural and physiological rhythms to match with them [57]. This clock system shows adaptive divergence between northern and southern species [72,73,74,75], and several clock genes have also been found to be differentially expressed during cold acclimation in studied species [62, 63, 66], suggesting an association between the circadian clock system and cold adaptation. Here, we traced correlation between the number of nucleotide or amino acid differences in vrille and variation in the mean chill coma resistance (CT_min) or chill coma recovery time (CCRT) in D. montana and D. flavomontana populations, and found a significant correlation only between amino acid variation and CT_min in D. montana. We then used RNAi to silence vrille and to investigate its effects on D. montana females’ cold tolerance traits and cold acclimation ability, and found this gene to play an important role in CT_min and in females’ cold acclimation ability, but not in CCRT. These findings together highlight the importance of vrille in enhancing females’ cold tolerance both during the rapid cold hardening occurring in CT_min test [81] and during the longer-term cold acclimation. Both traits, short- and long-term cold acclimation, have important ecological implications in cold environments with high daily and seasonal temperature changes.

Although we found a link between vrille and short- and long-term cold acclimation, it is not clear what are the genetic mechanisms underlying this link and to what extent the link exists. First, vrille might contribute to cold acclimation ability via circadian clock system, which may either directly contribute to cold acclimation, or regulate other pathways contributing to it. Here, any disruptions in the clock genes could cause similar effects. Alternatively, vrille might have more direct role on cold acclimation ability. For example, circadian clock system has been linked to Na⁺/K⁺ ion channel activity in D. melanogaster [101], which is also involved in avoiding cold injury, as observed in D. montana [102]. The second open question is whether the direct or indirect effects of vrille are limited to D. montana, D. virilis group species, or insects or organisms in general. D. virilis group species are among the most cold-tolerant Drosophila species [1, 2], their circadian clock system is adapted to northern habitats [72,73,74,75] and vrille expression is highly upregulated during cold-acclimation [64, 65], which gives indirect support for the role of vrille or the circadian clock system in cold adaptation in this group.

Although becoming more feasible, pinpointing the functions of genes affecting ecologically relevant traits is still rare. Previous RNAi experiments have verified e.g. Hsp22 and Hsp23 genes to contribute to CCRT in D. melanogaster [12], myo-inositol-1-sphosphate synthase (Inos) to contribute to D. montana flies’ survival during cold stress (5 °C), but not to affect their CCRT [13], and Gs1-like (Gs1l) gene to be involved in D. montana cold acclimation [103]. vrille is an interesting addition to this list, and while the details on the specific molecular mechanisms remain unclear, it gives new insights on the genetics underlying insect cold tolerance and cold acclimation.

Conclusions

Studying the mechanisms that generate and maintain variation in species stress tolerances is essential for understanding adaptation processes. We show that insect cold tolerance may rely on different environmental cues and morphological traits even in closely-related species, and that insects’ chill coma recovery time may be affected by stronger selection pressures in nature than chill coma resistance. These kinds of studies may have applications e.g. in predicting the likely outcomes of climate change or invasion biology. Species, whose cold tolerance is tightly linked with photoperiod, may encounter more difficulties in adapting to changing temperature conditions than the ones whose tolerances are driven by local climatic conditions. We also propose that vrille, and possibly the whole circadian clock system, play an essential role in molecular mechanisms underlying short- and long-term cold acclimation, both of which are ecologically important traits on high latitudes with high daily and seasonal temperature variation. In the future, studies on the relationship of reproductive stage and cold tolerance traits in individuals originating from diverse climatic conditions could deepen our understanding on adaptation to seasonally varying environments. Also, investigating the effects of silencing the other core circadian clock genes under cold environment and in other organisms would give valuable insights on the role of the circadian clock system on cold acclimation in general.

Methods

Genetic variation in cold tolerance, body colour and body size in D. montana and in D. flavomontana populations

Study species and populations

Drosophila montana and D. flavomontana belong to the montana phylad of the Drosophila virilis group, and our recent whole genome analyses have shown their divergence time to be ~ 1 mya (Poikela et al., in preparation). D. montana is distributed around the northern hemisphere across North America, Asia and Europe [70], while the distribution of D. flavomontana is restricted to North America [69, 70]. In the central Rocky Mountains, D. montana is found at altitudes from 1400 m to well over 3000 m, while D. flavomontana is found mainly below 2000 m. In the western coast, where D. flavomontana has probably invaded only during the last decades, both species live at much lower altitudes (see [69, 71]).

We investigated cold tolerance traits and measures of body colour and body weight using females from 23 D. montana and 20 D. flavomontana isofemale strains, which were established from the progenies of fertilized females collected from several sites in North America between 2013 and 2015 (Fig. 1). Each site was represented by three isofemale strains per population per species, when possible (Fig. 1; Additional file 1: Table S1). D. flavomontana is particularly rare in the western coast of North America [71] and we succeeded to collect only 1–2 strains from the study populations. All the strains were maintained in continuous light at 19 ± 1 °C since their establishment (15–30 generations) to eliminate plastic responses and prevent females from entering reproductive diapause. For the experiments, we sexed emerging flies under light CO₂ anaesthesia within 3 days after emergence. Females were changed into fresh malt-vials once a week and used in experiments at the age of 20 ± 2 days, when they all had fully developed ovaries [104].

Cold tolerance traits

We investigated variation in female cold tolerance using two well-defined and ecologically relevant methods: chill coma temperature (CT_min; also called critical thermal minimum) and chill coma recovery time (CCRT). CT_min corresponds to the temperature, at which the fly resists cold until it loses all neurophysiological activity and coordination and falls into a chill coma during a gradual decrease in temperature (reviewed in [14]). In this test, we placed the females individually in glass vials, which were submerged into a 30% glycol bath. We then decreased the bath temperature from the starting temperature of 19 °C at the rate of 0.5 °C/min and scored the CT_min for each fly. The second method, CCRT, measures the time taken from a fly to recover from a standardized exposure time at a chill-coma-inducing temperature (reviewed in [14]). In this test, we submerged the females individually in glass vials into a 30% glycol bath for 16 h at − 6 °C [65]. After returning the vials into 19 ± 1 °C in fly laboratory, we measured the time required for each female to recover from the chill coma and stand on its legs. CT_min tests were replicated 21 times and CCRT tests 20 times with Julabo F32-HL Refrigerated/Heating Circulator and Haake k35 DC50 Refrig Circulating Bath Chiller, respectively. To account for possible variation between replicates, each test included 1–3 females from each strain.

Fly body colour and body weight

We analysed variation in the body colour and body weight (as a proxy of body size) of the same females that had been phenotyped in CT_min or CCRT tests. Immediately after the cold tolerance tests, the females were put individually into tightly sealed PCR plates and kept in -20 °C freezer until measuring their body colour and weight. For body colour measurements, the females were photographed under Leica L2 microscope with × 5 magnification, using Leica DFC320 R2 camera together with the Leica Application Suite software v4.3.0. Exposure time, zoom factors and illumination level were kept constant, and the photographing arena was surrounded by a plastic cylinder to eliminate glares on the chitinous surface of the fly. All photographs were taken within 3 months after the cold tolerance tests. Images were saved in 24-bit (RGB) format and the colour intensity was measured using grayscale image analysis (ImageJ, [105]); linearly scaling from 0 to 255 (0 = black, 255 = white). We took colour measurements from part of thorax (scutum), as our preliminary tests showed that it best incorporates the colour variation among flies (Additional file 1: Fig. S2). Body weight was measured by weighting the females with Mettler Toledo™ NewClassic Balance (model ME). The weight of the females used in CT_min tests was measured after females had been frozen for 2–4 months, which appeared to be problematic as the females had started to dry and their weight correlated with the freezing time (Additional file 1: Fig. S3). The weight of the females used in CCRT tests was measured after 6–17 months in freezer. These females had lost most of their body liquids so that their weight was close to dry weight and freezing time had no effect on it (Additional file 1: Fig. S3). Accordingly, only the latter dataset was used in the analyses.

Statistical analyses

We used different statistical models to investigate whether variation in fly cold tolerance, body colour and body weight was associated with latitude, as a proxy of photoperiod, and/or local climatic variables, and whether cold tolerance traits and morphological traits showed a correlation with each other. In these models, we used either CT_min data (chill coma temperatures in Celsius degrees + 10 °C to prevent negative values from affecting the analysis), CCRT data (in minutes), body colour (measured from CCRT flies) or body weight (mg; measured from CCRT flies) as response variables. D. flavomontana CT_min data were normally distributed (Additional file 1: Fig. S4) and they were analysed with generalized linear mixed model (GLMM) with Gaussian distribution (equivalent of linear mixed model). Other datasets showed deviation from the normality (Additional file 1: Fig. S4), and they were analysed using generalized linear mixed model (GLMM) with gamma distribution, using glmmTMB function from glmmTMB package [106]. Technical replicates and isofemale strains were handled as crossed random effects.

In our dataset latitude and altitude were negatively correlated (Pearson correlation coefficient = − 0.82) and to prevent this multicollinearity from affecting the analysis, we used only latitude as an explanatory factor. However, we considered the effect of altitude and climatic variables through a principal component analysis (PCA). We downloaded climatic information from WorldClim database v2.1 (2.5 min spatial resolution; current data 1970–2000; [107]; http://www.worldclim.org) and extracted 19 bioclimatic variables for each site using their latitudinal and longitudinal coordinates (Additional file 1: Table S3; Fig. 1) with raster package v. 2.8-19 [108]. Moreover, we obtained growing season length (days) for each site from article by [92] and www.weatherbase.com. Growing season is defined as the average number of days per year when the average daily temperature is at least 5 °C. We performed the PCA on the 19 bioclimatic variables, altitude and growing season length (Additional file 1: Table S3) to summarize climatic differences on temperature and precipitation in each site using “FactoMineR” package [109].

We included the first two PCs in the model comparison of cold tolerance, body colour and body weight (see ‘Results’). The simplest model included latitude and PC2, which enabled us to distinguish between the latitudinally varying photoperiods and temperatures. Moreover, the effect of macroclimatic conditions varying between the western coast and the Rocky Mountains (PC1), interaction terms between latitude/PC2 and PC1, and body colour (divided by 100 to scale the variables) and body weight were included in the model selection (Additional file 1: Table S6). The best-fit model for CT_min, CCRT, body colour and body weight of both species was chosen for further analysis based on Akaike information criterion (AIC) value and Akaike weight (Additional file 1: Table S6) using aictab function from AICcmodavg package [110]. All the analyses were conducted in R (v1.2.1335-1) and R studio (v3.6.1).

Investigating circadian clock gene vrille and its effects on fly cold tolerance using sequence variation and RNA interference (RNAi)

The candidate gene

Drosophila virilis species group shows several distinct features in their circadian clock system compared to D. melanogaster, which has enabled their adaptation in high-latitudes [72,73,74]. Interestingly, vrille was the only core circadian clock gene that was found to show both population-level differentiation in clinal cold tolerance study [89] as well as expression level differences when cold-acclimated and non-acclimated D. montana flies were compared [77]. This gene was also one of the differentially expressed genes in our earlier microarray study investigating cold tolerance of D. montana and D. virilis [65]. Accordingly, we investigated the sequence variation in vrille and traced its association with cold tolerance traits (CT_min and CCRT) in D. montana and D. flavomontana populations. We also used RNAi to silence the expression of vrille and investigated its effects on cold tolerance traits (CT_min and CCRT) and cold acclimation ability in D. montana, where this gene showed correlation with CT_min at the amino acid level (see Fig. 5B).

Association between the number of nucleotide or amino acid differences and differences in cold tolerance traits in D. montana and D. flavomontana populations

To evaluate the role of vrille in cold adaptation, we investigated the association between vrille sequence variation and variation in cold tolerance traits (CT_min and CCRT) in D. montana and D. flavomontana populations. We obtained published D. montana genomic sequences from NCBI under accession number LUVX00000000 and the respective genome annotation from Dryad: https://0-doi-org.brum.beds.ac.uk/10.5061/dryad.s813p55 [111]. Based on D. virilis and D. melanogaster orthologs in the annotation file, vrille sequence is found in scaffold2327-size32305 at position 13,933–17,558. We used Illumina 150 bp paired-end data from single wild-caught females or their F₁ female progeny, i.e. the founder females of the studied isofemale strains, from all study populations of both species (one individual per population; Additional file 1: Table S1; Poikela et al., in prep.).

Before mapping, Illumina paired-end reads were quality-checked with FastQC v0.11.8 (Andrews 2010: FastQC: a quality control tool for high throughput sequence data. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc) and trimmed for adapter contamination and low-quality bases using fastp v0.20.0 [112]. Each trimmed Illumina sample was mapped against vrille D. montana reference vrille gene using BWA mem v0.7.17 [113]. The alignments were sorted with SAMtools v1.10 [114] and PCR duplicates marked with sambamba v0.7.0 [115]. The resulting BAM-files were used for variant calling with freebayes v1.3.1-dirty [116], and the variant calling (VCF) file was filtered for minimum depth and quality of 20 using bcftools view v1.8 (SAMtools. 2018. https://github.com/samtools/bcftools/releases/tag/1.8). Finally, vrille FASTA sequences for each sample were obtained using bcftools consensus.

The vrille sequences together with information on exons from the annotation file (see Additional file 1: Figs. S5, S6) were imported to Geneious v9.1.8. (Biomatters ltd). For exons, D. montana and D. flavomontana nucleotide sequences and translated amino acid sequences were aligned using muscle algorithm. The number of pairwise nucleotide and amino acid differences were then calculated for D. montana and D. flavomontana populations (Additional file 1: Tables S12, S13). Similarly, pairwise differences in mean CT_min (in Celsius degrees) and CCRT (in minutes) were calculated for D. montana and D. flavomontana populations (Additional file 1: Tables S14, S15). The correlations between distance matrices were statistically tested with a Mantel test with 1000 permutations using mantel.rtest function from ade4 package [117]. The correlations were tested between (i) the number of nucleotide differences and differences in mean CT_min (ii) the number of amino acid differences and differences in mean CT_min, (iii) the number of nucleotide differences and differences in mean CCRT, and (iv) the number of amino acid differences and differences in mean CCRT among D. montana and D. flavomontana populations (only within-species correlations).

Study material for RNAi

We performed RNAi study using only the more cold-tolerant D. montana, where vrille showed correlation with CT_min at the amino acid level (see Fig. 5B) and which has become a model species for studying cold adaption at phenotypic and genetic level (e.g. [61, 64, 92]). Here, we used females of a mass-bred cage population originating from Seward (Alaska, USA; see Fig. 1) to increase genetic diversity. The population cage has been established from the 4th generation progenies of 20 isofemale strains (total of 400 flies) in 2013 and maintained in continuous light at 19 ± 1 °C since its establishment (30 generations). Malt bottles with freshly laid eggs were transferred from the population cage into a climate chamber in LD 18:6 (light:dark cycle) at 19 ± 1 °C to reinforce flies’ circadian rhythmicity prior to the experiments and ensure they are non-diapausing. Critical day length for the induction of reproductive diapause (CDL; 50% of the females of given population enter diapause) in 19 °C is LD 17:7 in Seward population [92], and thus the females emerging in LD 18:6 can be expected to develop ovaries. After ~ 4 weeks, we collected newly emerged females (≤ 1 day old) using light CO₂ anaesthesia and placed them back into the above-mentioned conditions in malt-vials. Females were changed into fresh vials once a week until they were used in the experiments at the age of 21 days.

Defining daily expression rhythm of vrille for RNAi study

Expression levels of circadian clock genes are known to show daily fluctuations, and in order to perform RNAi experiments at a right time of the day, we defined the time when the expression of vrille is highest at LD 18:6. To do this, we collected females every four hours (i.e. at six time points/day), starting at 6 am (ZT0) when lights were switched on. At each time point, we stored samples of females into – 80 °C through liquid N₂ and transferred them later on into RNAlater®-ICE frozen tissue transition solution (Thermo Fisher Scientific). We then checked the size of female ovaries (see [104]) and used only the females with fully developed ovaries (> 95% of the females). For each time point, RNA was extracted from three pools, each of which consisted of three whole females, using ZR Tissue & Insect RNA MicroPrep kit with DNAse treatment (Zymo Research®). RNA purity was analysed with NanoDrop® ND-1000 spectrophotometer (Thermo Fisher Scientific) and concentration with Qubit® 2.0 Fluorometer and RNA BR kit (Thermo Fisher Scientific). cDNA was synthesized using equal quantities of RNA (200 ng) with iScript Reverse Transcription kit (Bio-Rad Laboratories®).

We measured the expression levels of vrille with quantitative real time PCR (qPCR). qPCR primers for vrille and reference genes were designed based on D. montana genomic sequences under NCBI accession number LUVX00000000 [111], together with information from D. virilis exons using Primer3 (primer3.ut.ee) and NetPrimer (www.premierbiosoft.com/netprimer) programs (gene accession numbers and primer sequences in Additional file 1: Table S16). qPCR mix contained 10 µl of 2 × Power SYBR Green PCR Master Mix (Bio-Rad Laboratories), 0.3 µM of each gene-specific primer and 1 µl of cDNA solution. Cycling conditions in Bio-Rad CFX96 instrument were: 3 min at 95 °C, 10 s at 95 °C, 10 s at annealing temperature of 53 °C (for reference genes 56 °C) and 30 s at 72 °C (repeated 40 ×), followed by melting curve analysis (65–95 °C) for amplification specificity checking. Each run included three technical replicates for each sample and the final threshold value (Cq) was defined as a mean of the technical replicates that produced good quality data. The relative qPCR data was normalised with ∆∆(Ct) normalisation method [118] using two reference genes, Tubulin beta chain (Tub2) and Ribosomal protein L32 (RpL32), that showed equal expression levels in all samples (data not shown). Real efficiency values of the genes used in the qPCR are given in Additional file 1: Table S16.

Synthesis of double-stranded RNA for RNAi

LacZ, which codes a part of a bacterial gene, was used as a control for dsRNA injections. We generated fragments of vrille and LacZ genes, with the length of 347 and 529 bp, respectively, with PCR (primer information given in Additional file 1: Table S16). PCR products were purified with GeneJET Gel Extraction kit (Thermo Fisher Scientific) and cloned using CloneJET PCR Cloning kit (Thermo Fisher Scientific). PCR products were ligated into the vector (pJET1.2/blunt Cloning vector), transformed into E. coli Zymo JM109 (Zymo Research) cells, which were grown on Luria Broth (LB) ampicillin plates. Individual colonies were picked up after 16 h and cultivated overnight in LB solution with ampicillin using a Unimax 1010 shaker with incubator (Heidolph Instruments). The samples from the extracted bacterial solutions were analysed for the size of the products in the second PCR, which was carried out with pJET primers from the cloning kit, followed by agarose gel runs. We then selected the colonies with the right size products for the third PCR using pJET primers, where the R primer contained T7 promoter sequence at the 5′ end of the primer (primer sequences in Additional file 1: Table S16). PCR products were first purified with GeneJet Gel Extraction kit and then used in transcription synthesis of the double-stranded RNA (dsRNA), using the TranscriptAid T7 High Yield Transcription kit (Thermo Fisher Scientific). Finally, we purified and precipitated the synthesized products with ethanol, suspended them in salt buffer and quantified them using NanoDrop and agarose gel.

RNAi microinjection procedure, response time screening and vrille expression at the chosen response time

Injecting dsRNA targeting vrille gene is expected to cause gene-specific effects, but it may also cause immune responses and physical damage (injection) in the flies. Accordingly, we used LacZ (encoding for bacterial gene) injections as a control for both immune response to dsRNA and to physical damage of injections, and no-injection as a baseline control. We checked the effectiveness of RNAi treatment on vrille expression 12 h, 24 h and 48 h after injections, at a time of the day when it shows highest expression (ZT16, see ‘Results’). The females were injected into thorax with 138 nl of ~ 20 µM dsRNA targeting vrille or LacZ using Nanoject II Auto-nanoliter Injector with 3.5″ Drummond glass capillaries (Drummond Scientific) pulled with P-97 Flaming/Brown Micropipette Puller (Sutter Instrument). No-injection control females were not injected, but were otherwise handled in the same way as the injected females. To prevent CO₂ anaesthesia from inducing variation between these groups, we injected six females at a time. For each response time (12 h, 24 h, 48 h) all three treatments (vrille, LacZ, no-injection) were performed. Each treatment consisted of three pools, each pool containing 10 whole females. After reaching the response time, the females were transferred into − 80 °C through liquid N₂. Then, the females were transferred into RNAlater®-ICE solution and their ovaries were checked as explained above. RNA was extracted from the pools using TRIzol® Reagent with the PureLink® RNA Mini kit (Thermo Fisher Scientific). RNA purity was analysed with NanoDrop and concentration with Qubit and RNA BR kit. cDNA was synthesized using equal quantities of RNA (143 ng) using SuperScript IV First-Strand Synthesis System (Thermo Fisher Scientific). Expression levels of vrille 12, 24 and 48 h after the injections were quantified with qPCR, as described above. The response time of 48 h was the most effective (Additional file 1: Fig. S1), and its effect was double-checked with another round of cDNA synthesis (200 ng) and qPCR (Fig. 6B).

Experimental design for studying the functional role of vrille in cold tolerance using RNAi

We investigated the role of vrille gene in resisting chill coma (CT_min) and recovering from it (CCRT) in non-diapausing D. montana females originating from Seward, Alaska, USA. Moreover, we considered the plastic effects of vrille during cold acclimation. Prior to performing cold tolerance tests, females were maintained for 16 days in LD 18:6 at 19 °C, corresponding to summer conditions at their home site (LD 18:6 was used throughout the experiment). Then, half of these females were subjected to cold acclimation treatment at 6 °C for 3 days (cold-acclimated females; [65]), while the other half was maintained at 19 °C. At the age of 19 days, both cold-acclimated and non-acclimated females were collected from the chambers, anesthetized with CO₂ and injected as described above. They were then placed back to either cold acclimation (6 °C) and non-acclimation (19 °C) conditions for two more days, as the expression levels of target genes had been found to be lowest 48 h after RNAi treatment (Additional file 1: Fig. S1). At the age of 21 days, females’ cold tolerance was quantified by measuring their chill coma temperature (CT_min) or chill coma recovery time (CCRT) using Julabo F32-HL Refrigerated/Heating Circulator. Sample sizes for CT_min and CCRT tests were 26–32 and 25–30 females per treatment, respectively.

Statistical analyses of RNAi experiment

For investigating the effects of RNAi, expression levels of vrille in LacZ-injected females were compared to no-injection control females and females injected with dsRNA targeting on vrille. The relative normalised expression values were analysed using a linear model (ANOVA) in base R.

To test the cold acclimation effect within each treatment group (vrille-RNAi, LacZ-RNAi, no-injection control), we compared the cold tolerance of the females that had or had not been cold-acclimated. We then investigated the gene specific effects on cold tolerance by comparing LacZ-injected females to vrille-injected females and traced possible immune and physical effects of the injections by comparing LacZ-injected females to no-injection females. These analyses were performed separately for the cold-acclimated and non-acclimated females. All data showed deviation from normality (Additional file 1: Fig. S7) and were analysed with generalized linear mixed model using gamma distribution (GLMM; glmmTMB function from glmmTMB package [106]. In the models, response variables were either CT_min (Celsius degrees + 10 to prevent negative values from affecting the results) or CCRT (minutes) data, and the test replicates were used as a random effect. All the analyses were conducted in R (v1.2.1335-1) and R studio (v3.6.1).