Puig Giribets et al. BMC Evolutionary Biology (2020) 20:17 https://doi.org/10.1186/s12862-020-1584-z RESEARCH ARTICLE Open Access Basal hsp70 expression levels do not explain adaptive variation of the warm- and cold-climate O3 + 4 + 7 and OST gene arrangements of Drosophila subobscura Marta Puig Giribets, Mauro Santos and María Pilar García Guerreiro* Abstract Background: Drosophila subobscura exhibits a rich inversion polymorphism, with some adaptive inversions showing repeatable spatiotemporal patterns in frequencies related to temperature. Previous studies reported increased basal HSP70 protein levels in homokaryotypic strains for a warm-climate arrangement compared to a cold-climate one. These findings do not match the similar hsp70 genomic organization between arrangements, where gene expression levels are expected to be similar. In order to test this hypothesis and understand the molecular basis for hsp70 expression, we compared basal hsp70 mRNA levels in males and females, and analysed the 5′ and 3′ regulatory regions of hsp70 genes in warm- and cold-climate isochromosomal O3 + 4 + 7 and OST lines of D. subobscura. Results: We observed comparable mRNA levels between the two arrangements and a sex-biased hsp70 gene expression. The number of heat-shock elements (HSEs) and GAGA sites on the promoters were identical amongst the OST and O3 + 4 + 7 lines analysed. This is also true for 3′ AU-rich elements where most A and B copies of hsp70 have, respectively, two and one element in both arrangements. Beyond the regulatory elements, the only notable difference between both arrangements is the presence in 3′ UTR of a 14 bp additional fragment after the stop codon in the hsp70A copy in five O3 + 4 + 7 lines, which was not found in any of the six OST lines. Conclusions: The equivalent hsp70 mRNA amounts in OST and O3 + 4 + 7 arrangements provide the first evidence of a parallelism between gene expression and genetic organization in D. subobscura lines having these arrangements. This is reinforced by the lack of important differential features in the number and structure of regulatory elements between both arrangements, despite the genetic differentiation observed when the complete 5′ and 3′ regulatory regions were considered. Therefore, the basal levels of hsp70 mRNA cannot account, in principle, for the adaptive variation of the two arrangements studied. Consequently, further studies are necessary to understand the intricate molecular mechanisms of hsp70 gene regulation in D. subobscura. Keywords: Drosophila subobscura, hsp70 expression, Chromosomal arrangements, Thermal adaptation * Correspondence: mariapilar.garcia.guerreiro@uab.es Grup de Genòmica, Bioinformàtica i Biologia Evolutiva (GGBE), Departament de Genètica i de Microbiologia, Universitat Autònoma de Barcelona, Bellaterra, 08193 Barcelona, Spain © The Author(s). 2020 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Puig Giribets et al. BMC Evolutionary Biology (2020) 20:17 Page 2 of 9 Background +5.4σ from the expected average frequency (in close Inversions can foster local adaptation as they protect agreement with the temperature anomaly), and high selec- chromosomal inverted sequences from recombination in tion coefficients are needed to explain this shift [18]. heterokaryotypes. Drosophila subobscura exhibits a rich This interpretation is not without caveats. A trend to a inversion polymorphism [1], with the O chromosome decrease in basal HSP70 protein expression was ob- (Muller element E that also corresponds to chromosome arm 3R in Drosophila melanogaster ([2], p., 307) having served in strains of D. melanogaster [19, 20] and Dros- the highest number of inversions that form complex ophila buzzatii [21] continually exposed to high gene arrangements of overlapping and non-overlapping temperatures, probably to avoid the deleterious effects of inversions [3]. The overlapping inversions cluster in the increased HSP70 levels. On the other hand, increased distal segment of the chromosome and form the O3þ4 basal HSP70 protein and mRNA levels have been ob- phylad, whose different gene arrangements have higher served in certain strains of Drosophila virilis with in- frequencies at lower latitudes [1, 4]. Thus, the gene ar- creased thermotolerance [22]. Therefore, basal levels of rangements on chromosome O can be divided into two groups: group 1 (cold-climate) includes OST that shows HSP70 protein expression might be contingent on the a negative correlation coefficient with maximum temper- thermal history and the species under study. A positive atures along the latitudinal clines, or a positive correl- link between hsp70 copy number, expression levels, and ation coefficient with latitude; and group 2 (warm- thermal resistance in artificial [23–25] and natural [26] climate) that includes gene arrangements showing the Drosophila strains has been reported, but this does not reverse pattern (e.g., O3þ4 , O3þ4 þ 7 and O3þ4þ8 ). In apply to D. subobscura. In this species, we have recently agreement with the clinal patterns, warm-climate inver- shown that the genomic sequence organization of the sions increase in frequency during summer and decrease during winter [5, 6]. Furthermore, they show long-term hsp70 gene family (located in the 94A subsection of the trends of increasing frequency related to global warming O chromosome, which is inside the region covered by [7, 8]. Overall, these spatiotemporal patterns strongly O3 and O4 inversions [27];), is essentially similar in cold- suggest that temperature is the main selective factor (OST) and warm-climate ( O3þ4 , O3þ4þ8 and O3þ4þ16 ) driving inversion frequencies, but the underlying genetic gene arrangements, and consists of only two genes, each causes are still poorly understood. The heat-shock HSP70 protein is involved in protein with four heat-shock elements and three GAGA sites on folding, assembly, and can be induced by several its promoter [27]. To date, no clear evidence has been stressors such as a sudden increase in temperature [9, found to link the genomic sequence organization of the 10]. It appears to be the primary protein involved in hsp70 gene family with the different HSP70 protein ex- thermotolerance in D. melanogaster [11, 12], and allele pression patters previously detected [14]. Knowing that proteins are produced from mRNA frequencies show latitudinal clines and change in re- which, in turn, depends on gene regulatory regions, we sponse to thermal evolution in the laboratory [13]. In a extended the former approach to gauge if hsp70 mRNA previous work with D. subobscura, expression levels of HSP70 protein in flies carrying cold- and warm-climate levels are similar in cold- and warm-climate arrange- gene arrangements were analysed [14]. It was found that ments as expected by their similar genomic organization. homokaryotypic O3þ4 flies had increased HSP70 basal We quantified basal hsp70 mRNA expression levels in males and females, and included a higher number of iso- protein levels that did not boost after a heat-shock, in genic OST and O3þ4 þ 7 strains than in our previous contrast to what happened to their OST homokaryotypic counterparts. This result was somehow surprising because study [27]. We expanded the study to the 5′ and 3′ HSP70 protein expression entails fitness costs [15–17]. hsp70 regulatory regions, which are known to be in- The interpretation was that low-latitude populations with volved in transcriptome variation in D. melanogaster relatively high O3þ4 frequencies are more likely to be ex- [28]. A clarification is in order. The arrangement O3þ4 posed to occasionally higher temperatures than popula- þ 7 has the overlapping inversions O3þ4 in the distal tions located at higher latitudes, and the benefits of segment and the single O7 inversion in the proximal seg- mounting the stress response might outweigh the costs. ment of chromosome O. Because our lines derived from Actually, an extreme heat wave in Europe in spring 2011 a natural population where both gene arrangements caused frequencies of warm-climate O3þ4 þ 7 to shift (OST and O3þ4 þ 7 ) have a relatively high frequency Puig Giribets et al. BMC Evolutionary Biology (2020) 20:17 Page 3 of 9 (Berbikiz, north of Spain [18]) the heterokaryotype O3þ4 and proximal (HSE1, HSE4) and four in HSE2 and HSE3 =O3þ4 þ 7 freely recombines in the distal segment and, (Fig. 1). Within each HSE, polymorphisms appear to be overrepresented in the first and fifth position of therefore, the gene content in the O3þ4 region is ex- nGAAn/nTTCn pentanucleotides, which are the least pected to be the same in both arrangements. The reason constrained [29]. Some insertion/deletion (INDEL) poly- for using O3þ4 þ 7 was because the integrity of the wild morphisms were also detected in the presumably uncon- chromosome is better preserved during the isogenization strained regions separating the more distal HSEs. These polymorphisms tend to be shared between the two process with the available balancer stock in D. subobs- hsp70 paralogs of the same line as a result of ectopic cura (see below). gene conversion [27]. We did not observe differences be- We found equivalent mRNA amounts in OST and tween isogenic lines either in the two core proximal pro- O3þ4 þ 7 arrangements, a sex-biased hsp70 mRNA gene moter elements (TATA box and TSS), or in the total expression compatible with the high degree of conserva- number of GAGA sites in the proximal promoters. The tion of the cis-regulatory elements observed, and an latter is known to be recognized by the transcription fac- identical genomic organization of the gene family previ- tor GAF and to participate, like other promoters, in the ously reported in our study. regulation of global gene expression [30]. The distal regulatory region corresponds to the 3′ un- Results translated mRNA region (3’UTR), which includes a re- High degree of conservation in the hsp70 proximal and gion spanning 198 ± 4 bp in all hsp70A and hsp70B distal CRE elements between OST and O3 + 4 + 7 copies, except in one hsp70B copy from OST (4th line) arrangements that is 381 bp long. The in silico screening let us identify We have sequenced the proximal and distal regulatory conserved regulators consisting in AU-rich elements regions of the hsp70 locus in 12 isogenic lines (6 lines (AREs), known to be involved in the post-transcriptional for OST and 6 for O3þ4 þ 7 gene arrangements). The regulation of gene expression in Drosophila [31] by proximal 5′ region consists of a ~ 550 bp, of which ~ recruiting ARE-binding proteins that trigger mRNA 290 bp correspond to the region upstream of the tran- deadenylation and turnover. In most lines we identified scription start site (TSS) in the two hsp70 gene copies two and one ARE motifs, consisting in a single AUUUA (hsp70A and hsp70B). Multiple sequence alignments in- pentamer in hsp70A and hsp70B copies, respectively. dicated high CRE (Cis-Regulatory Element) conservation There are a few exceptions; namely, the 1st and 6th OST between OST and O3þ4 þ 7 . As previously found [27], lines that have an unique motif in hsp70A (Fig. 2 and four highly conserved heat shock elements (HSEs) were Additional file 3) and the 4th OST line, where the 3’UTR identified (Fig. 1), and no fixed nucleotide differences region of the hsp70B copy is longer and has a second between arrangements were observed (Additional files 1 ARE motif (Fig. 2 and Additional file 4). We found that and 2). The genomic organization of HSEs resembles ARE motifs are, in general, highly conserved in both that of other Drosophila species [29] and consists of at copies of OST and O3þ4 þ 7, but the hsp70A copy in five least three continuous alternate nGAAn/nTTCn penta- O3þ4 þ 7 lines has a 14 bp additional fragment after the nucleotide subunits representing binding sites for the stop codon, which was not found in OST lines (Add- heat shock factor (HSF) monomer. Although the num- itional file 3). ber of repetitions is identical in the 12 isogenic lines, it Table 1 summarizes the level of genetic differentiation varies between HSE sequences: three in the most distal between chromosomal rearrangements for the complete Fig. 1 Representation of hsp70 5′ cis-regulatory elements of OST and O3 + 4 + 7 lines. The two copies of the gene (A and B) are included. Green rectangles: HSEs 1–4, where numbers in brackets represent the number of nGAAn or nTTCn pentanucleotide units. Red ellipse: GAGAG (G+) repeat motifs. Blue ellipses: CTCTC (G-) repeat motifs. Orange triangle: TATA box. The scale represents nucleotides relative to the transcription start site (+ 1), indicated by a forward arrow Puig Giribets et al. BMC Evolutionary Biology (2020) 20:17 Page 4 of 9 Fig. 2 Schematic representation hsp70 3′ UTR region of OST and O3 + 4 + 7 lines. The two copies of the gene (A and B) are included. The region represented spans from the stop codon to the polyadenylation signal. Green rectangles: ARE motifs. In red: polyadenylation signal. Square brackets: number of nucleotides from the stop codon and between ARE motifs nucleotide sequence of 5′ and 3′ regulatory regions of A in Table 2, and results are plotted in Fig. 3. Global hsp70 and B copies. Regulatory regions in both copies, except mRNA levels were about the same between arrange- 3’B, exhibited significant genetic differentiation between ments (P = 0.473) but were higher in females than in the two arrangements studied. However, whereas differ- males (P = 0.005), particularly for the O3 + 4 + 7 gene ar- entiation levels were similar in 5′ regulatory regions for rangement; but the arrangement × sex interaction was A and B copies, they were significantly higher in 3′ regu- not statistically significant (P = 0.296). latory regions from A copy, where the average number of differences per site (Dxy) was much higher than in the B copy. Discussion Here we have sequenced the proximal and distal regula- tory regions of the hsp70 loci and have analysed mRNA Basal hsp70 mRNA expression levels expression levels in males and females from 6 isochro- hsp70 mRNA levels have been quantified in 67 (out of mosomal lines of warm- (O3þ4 þ 7) and 6 isochromoso- 72) D. subobscura samples consisting of 3 biological rep- mal lines of cold-climate (OST) gene arrangements in D. licates per isogenic line (6 OST and 6 O3þ4 þ 7) and sex, subobscura. Our aim was to shed some light on the basal except five lines where only two replicates for females amounts of hsp70 mRNA in these arrangements that and/or males could be analysed. The ANOVA is shown show an almost identical genetic organization in the Table 1 Genetic differentiation between O3 + 4 + 7 and OST arrangements for 5′ and 3′ regulatory regions Region Comparisons Gen copy SPSTSF347 SP347SFST Shared Dxy Fst P 5’ OST vs O3 + 4 + 7 Conc 17 35 2 0.031 0.344 0.000*** A (6) 13 33 2 0.034 0.306 0.008** B (6) 10 28 2 0.030 0.330 0.045* 3’ OST vs O3 + 4 + 7 Conc 17 38 7 0.127 0.045 0.000*** A (6) 16 5 0 0.075 0.713 0.007** B (6) 1 9 0 0.014 0.081 0.068 Sample sizes are in parenthesis. Fst: Proportion of nucleotide diversity due to differences between populations; Dxy: Average number of nucleotide substitutions per site between populations; SPSTSF347: Sites polymorphic in OST and fixed in O3 + 4 + 7; SP347SFST: Sites polymorphic in O3 + 4 + 7 and fixed in OST; Conc: concatened data set; Statistical significance was assessed using the P-value of Snn statistic; * P < 0.05; ** P < 0.01; *** P < 0.001. Puig Giribets et al. BMC Evolutionary Biology (2020) 20:17 Page 5 of 9 Table 2 ANOVA for the effects of chromosome gene arrangements (O3 + 4 + 7 and OST), isochromosomal line, and sex on mRNA relative expression (2−ΔC T ). Data plotted in Fig. 3 Source d.f. MS F P (parametric) P (permutation) Arrangement 1 0.000084 0.56 0.473 0.471 Line (Arrangement) 10 0.000152 1.36 0.223 0.222 Sex 1 0.000955 8.55 0.005** 0.004** Arrangement × Sex 1 0.000125 1.11 0.296 0.294 Error 53 0.000112 Arrangement was nested within lines. MS: Mean squares, d.f.: degrees of freedom; ** P < 0.01. hsp70 locus [27]. We have now performed for the first Our main aim here was to correlate the genome time a thorough analysis of 5′ and 3′ hsp70 regulatory organization of hsp70 loci between cold- (OST) and a regions together with the quantification of mRNA levels warm-climate ( O3þ4 þ 7 ) gene arrangements with in OST and O3 + 4 + 7 arrangements. We found similar mRNA basal expression patterns. The use of isochromo- hsp70 mRNA levels in both arrangements, in accordance somal lines precludes any extrapolation to wild, outbred with the identical gene configuration found in our previ- flies. mRNA and protein levels are not always correlated ous work [27], but also found differences due to sex [32], and we did not subject our flies to a thermal stress because inbreeding can interfere with the heat-shock re- (Table 2). Since transcript levels are under the influence sponse [33]. Besides, inbred individuals might express of 5′ and 3′ regulatory elements, we analysed these re- conditional lethality under heat stress [34]. Therefore, gions and found that the number of 5′ HSE elements in the present results cannot be directly comparable to both O3þ4 þ 7 and OST arrangements was also identical, those obtained by Calabria et al. [14] because these au- as previously observed in other chromosomal arrange- thors measured both basal and heat-induced HSP70 pro- ments as O3 + 4 + 8 and O3 + 4 + 16, where these elements tein levels in outbred strains (obtained from cyclically evolve in concert through gene conversion [27]. Analysis permuted reciprocal crosses between isochromosomal of 3′ UTR regions also showed no noticeable differences lines). The differences observed here in the number of between chromosomal arrangements in the number of ARE motifs among the two hsp70 copies could contrib- ute to ensure basal amounts of HSP70 protein, basically ARE elements, considered to be critical regulators of from one of the gene copies. Moreover, the existence in mRNA turnover [31]. However, differences were ob- D. melanogaster of 3 and five ARE motifs, respectively, served between paralogous genes: in most lines, two in the two hsp70 gene clusters compared to D. subobs- ARE motifs were found in hsp70A but only one in cura could contribute to explain why in the previous hsp70B, which could point to post-transcriptional work [14] certain basal protein amounts were detected. regulation events responsible for RNA shortening and In this regard, a recent study suggested that differences degradation. Interestingly, the 3’UTRs were already con- in acclimation capacity between D. melanogaster and D. sidered as more divergent between paralogs than the subobscura correlate with structural differences in the coding and proximal promoter regions [27]. hsp70 regulatory regions [35]. However, in view of the highly conserved regulatory elements between chromo- somal arrangements found here, we cannot easily ex- plain the differences observed in protein amounts between arrangements [14]. These regulatory elements seem to be particularly highly constrained even though significant genetic differentiation between arrangements was detected when regulatory 5′ and 3′ were considered as a whole (Table 1). This is in concordance with a pre- vious study where differences between arrangements at non-coding regions were also observed [36]. Previous findings pointed out that changes in the ex- pression of particular genes adjacent to inversion break- points could contribute to an inversion evolutionary success [37, 38]. However, studies using genome wide transcriptome analyses in D. melanogaster [39] and D. pseudoobscura [40] concluded that, whereas inversions affect the expression of numerous genes, little evidence Fig. 3 Hsp70 expression in females (black) and males (grey) carrying O3 + 4 + 7 of direct position effects was globally detected at the and OST arrangements. Error bars are 95% of confidence intervals chromosomal level. In D. subobscura the hsp70 loci are Puig Giribets et al. BMC Evolutionary Biology (2020) 20:17 Page 6 of 9 located in subsection 94A of the hypothetical O3 ances- equivalent. This provides some evidence of a parallelism be- tral arrangement. This subsection is inside the derived tween hsp70 mRNA expression and its genetic O4 and OST inversions that gave rise to present day organization. Therefore, the basal levels of hsp70 mRNA O3þ4 and OST gene arrangements, respectively [27]. cannot, in principle, explain the spatiotemporal patterns Using the recently assembled genome of D. subobscura observed for these two gene arrangements. A potential [41] we have estimated that the hsp70 loci are around 0.8 drawback of our study is that we cannot extrapolate the re- Mb from the closest O4 distal breakpoint and likely too far sults to wild, outbred flies. Further experiments are neces- away as to be directly affected by this breakpoint, an infer- sary to understand the actual role (if any) played by hsp70 ence that is in agreement with the similar RNA amounts genes in the adaptive value of the rich chromosomal poly- found in the two chromosomal arrangements analysed morphism on the distal segment of chromosome O, where here. Four genes (Pxd, Set 8, CG5225 and Acf) are located the hsp70 genes are located. near O4 inversion breakpoints [42], but at present we do not know whether potential position effects affecting these Methods genes play any role in the in the adaptive value of this inver- Base stocks and fly handling sion. In any case, fitness interactions among loci (epistasis) Wild D. subobscura males were collected by net- seem to be important because a recombination load (i.e., sweeping over banana baits fermented with baker’s yeast the loss of fitness caused by recombination between benefi- in Berbikiz (north of Spain; 43°11′20.31″ N, 3°5′23.74″ cial combinations of interacting alleles) was detected in W) in 2012. They were individually crossed to three or OST =O3þ4 þ 7 heterokaryotypes [43]. four virgin females from the ch-cu marker strain, which A potential limitation of our study is that we cannot ex- carries the recessive markers on the O chromosome trapolate the results to mRNA expression levels during the cherry eyes (ch) and curled wings (cu) and has a highly larval stage. Insects with complex life cycles as D. subobs- homogeneous genetic background, to identify the gene cura might evolve different thermal sensitivities across life arrangements of one set of the five major wild-type stages, particularly between mobile adults (with a higher chromosomes (the other set of homologous chromo- capacity for thermoregulatory behaviour) and more sessile somes coming from the ch-cu strain). Thereafter, off- (eggs, larvae, pupae) life stages. A previous study assessing spring males from each line were crossed to females the thermal evolution of total mRNA expression levels in from the Va/Ba (Varicose/Bare) balancer stock [49], D. subobscura larvae did not detect differential gene expres- which has the same genetic background than the ch-cu sion for hsp70 among thermal environments, though the strain, to obtain isochromosomal lines for the O highest temperature used in the experiments was likely not chromosome in an otherwise highly homogeneous gen- high enough as to induce a heat-shock response [44]. We etic background. cannot unequivocally generalize our present results to other A total of 6 OST and 6 O3þ4 þ 7 isogenic lines were life stages, but we think that there are not compelling rea- used in the experiments. All 12 lines were reared at sons to believe that the basal mRNA expression of hsp70 18 °C (the optimum temperature [50]) with a 12:12 light/ would have been different between O3þ4 and OST gene ar- dark cycle on David’s killed-yeast Drosophila medium rangements had we used larvae instead of adults. [51]. To obtain the experimental flies, eggs from the dif- An important finding is the trend towards higher mRNA ferent lines were placed in 2 × 8 cm vials (45 eggs/vial) amounts in females compared to males. Differences be- with 6 mL of food, which guaranteed low or negligible tween sexes are, in principle, not surprising because in larval crowding that could have triggered the hsp70 Drosophila a substantial fraction of the transcriptome stress response [20]. For each isogenic line tested, two shows sex-dependent regulation [45]. In our case, the dif- independent sets of an average of three biological repli- ferences could be somehow understandable if the HSP70 cates, each containing three 7-day old adult males and protein played sex-specific functions. In D. melanogaster, females, were used for mRNA analyses. this protein is expressed during spermatogenesis under heat-shock conditions in adults [46], and at normal rearing temperatures in the larval stage and onwards [16, 47]. In Characterization and sequence analysis of the hsp70 5′ the case of females, a novel role for the molecular and 3′ regulatory regions chaperone HSP70 in the regulation of border cell migration For each isogenic line, two genomic regions expanding during Drosophila oogenesis has also been reported [48]. ~ 300 bp upstream the transcription starting site and ~ 205 bp downstream the stop codon were sequenced for the two hsp70 gene copies (hsp70B and hsp70A; see Conclusions [27]). The primers used for amplification of 3′ hsp70B The mRNA levels from isogenic lines for OST and O3 + 4 + 7 and hsp70A regulatory regions were, respectively: P154- gene arrangements in D. subobscura are quantitatively proximal (5′- ACC CGC AAA ATT GAA CCC AA-3′), Puig Giribets et al. BMC Evolutionary Biology (2020) 20:17 Page 7 of 9 CD_1_F (5′-ACC ATA CAG AAC GAC AAG GGT- 3′) lightcycler at 59 °C using hsp70 specific primers. All as- and DMT_proximal (5′-AGT CGG AAT TGT GAA says were performed in three technical replicates, and GCC TT-3′). Due to the SNPs (Single nucleotide poly- the rp49 housekeeping gene was used as an endogenous morphisms) found in 5’regulatory regions between lines, control because this gene showed to be expressed different primer sets were used: a) hsp70B: ARN_1R (5′- equally in all samples analysed (P = 0.772, see Add- TGT GAT GCT TTG GCC CAG AT-3′) and D5_NP itional file 5). Relative expression of hsp70 was then cal- (5′-GAT CCT CAA GGA ATG CGA AA-3′)/D7_ culated with the comparative CT method [59] and ΔCT NP(5′-CTG TAG ATC TTC AGG CAT AAG CTG-3′)/ values for the set of samples are provided in NP-1-FW (5′ GTA TTG TGG CTT CTT AAC GAG Additional file 6. GTTC C-3′); b) hsp70A: ARN_1R and and D1_NP (5′- For each isogenic line, we amplified a conserved hsp70 CTT CGG TGG GTT GGA TTG TAG ATT C-3′)/ D3_ coding fragment of 226 bp, common to both hsp70 cop- NP(5′-CTA CTT TCT CTC CTG TGT ATG TCT GG- ies described in D. subobscura [27]. The primers used 3′)/D4_NP(5′-CTC TCC TGT GTA CAT ATG TCT for amplification were: HSP70_1L (5′-CAC AGT CTT TGC G-3′). To identify several CREs (Cis-Regulatory El- TGA CGC CAA GC-3′) and HSP70_1R (5′-TGT GAT ements) that are typically found in the proximal 5′ GCT TTG GCC CAG AT-3′). For the rp49 housekeep- flanking region, the sequences of the two hsp70 genes ing gene, the primers RP49F (5′-ACA TCG GTT ATG were first separately aligned using the MAFFT E-INS-i GCT CCA C-3′) and RP49R (5′-GAT TTC CTT GCG algorithm [52] available at: https://mafft.cbrc.jp/align- CTT CTT TG-3′) were designed from a D. subobscura ment/server/. Next, we used the information available GenBank sequence (Accession: AJ228921), amplifying a on the hsp70 proximal promoters of D. subobscura [27] 164 bp segment from the gene’s second exon. Primer ef- and other Drosophila species [29] to identify the con- ficiencies of hsp70 and rp49, calculated on a standard served heat shock elements in the sequences of the 12 curve slope using serially diluted templates, were 98 and isogenic lines (HSEs 1–4), and the Eukaryotic Neural 99.6%, respectively. Network Promoter Prediction software (http://www. fruitfly.org/seq_tools/promoter.html) to find other CREs such as GAGA binding sites and some core promoter el- Statistical methods ements, including the transcription start site (TSS) and Hsp70 mRNA expression levels were analysed using the the TATA box. 3′ UTR sequences were aligned using linear model. the MAFFT E-INS-i algorithm [52] and manually refined yijkl = μ + Ai + Lj(i) + Sk + ASik + εijkl, (1) using the Genedoc alignment editor (version 2.7) [53]. where μ is the overall grand mean, Ai is the fixed effect The web-based ARE score application (http://arescore. of gene arrangement (OST, O3þ4 þ 7), Lj(i) is the random dkfz.de/arescore.pl) was used to identify post- effect of the jth isochromosomal line (j = 1, 2, ⋯, 6) transcriptional regulators of mRNA levels such as AU- within arrangement i, Sk is the fixed effect of sex, and rich elements [54]. The DNASp version v6.12 program εijkl is the residual error associated with the mRNA ex- [55] was used for genetic differentiation analyses between pression levels of the lth replicate from the kth sex with arrangements, assessed using Dxy [56] and Fst [57], and the ith gene arrangement from the jth isochromosomal the significance for Snn [58] was obtained with 1000 repli- line. Following Schmittgen and Livak [59], statistical cates. Analyses were performed excluding gaps. tests for mRNA relative expression were performed after the 2−ΔC T transformation. Quantification of hsp70 transcripts by qRTPCR To avoid potential problems with the standard as- Total RNA extraction was performed separately in adult sumption that errors are independent, and identically males and females using the Trizol Lysis reagent distributed as normal random variables with common (Thermo Fisher Scientific Inc.) following manufacturer’s variance and zero expectation, we also performed per- instructions, and then treated with 1 μl (2 U/μl) of mutation tests [60, 61] to check the validity of paramet- Ambion DNase I (Thermo Fisher Scientific Inc.) 4 h at ric P-values. In model (1) the permutation test is a four- 37 °C. Concentration of DNA-free samples was checked stage process. First, random permutation among lines, in a Nanodrop-2000 spectrophotometer (Thermo Scien- within arrangements, for the among-line Fline(arrangement); tific Nanodrop) and adjusted to ~ 1 μg/μL and then second, random permutation among lines and arrange- cDNA synthesis was carried out with anchored-Oligo ments for the between-arrangements Farrangement; third, (dT)18 primers using the Transcriptor First Strand random permutation between sexes for the between-sex cDNA Synthesis kit (Roche Diagnostics SL, SPAIN). Fsex; and fourth, random permutation between arrange- Gene expression was quantified by measuring fluores- ments and sexes for the arrangement × sex interaction cence intensity using iQ SYBR Green Supermix (BioRad, Farrangement × sex. In each case, 9999 random allocations Hercules CA, USA) on a CFX96 BioRad Real-Time plus the observed F-values were obtained, which gave us Puig Giribets et al. BMC Evolutionary Biology (2020) 20:17 Page 8 of 9 10,000 values of each F-statistics for each dependent Competing interests variable. We then calculate the percentage of replica- The authors declare no competing interests. tions under the null hypothesis that provide as large a Received: 4 October 2019 Accepted: 16 January 2020 value of each F as the obtained F. The equal expression of rp49 gene was tested on Ct values using one-way ANOVA with four sample groups: References females OST, males OST, females O3þ4 þ 7 and males 1. Menozzi P, Krimbas CB. The inversion polymorphism of D. subobscura revisited: Synthetic maps of gene arrangement frequencies and their O3þ4 þ 7 . The Levene’s test was used to test the homo- interpretation. J Evol Biol. 1992;5:625–41. geneity of variances between groups. 2. . Powell JR: Progress and Prospects in Evolutionary Biology. The Drosophila Model, New York: Oxford Univ Press; 1997. Analyses were performed using the statistical software 3. Krimbas CB. Drosophila subobscura, Biology, Genetics and Inversion STATISTICA 8.0 [62]. Permutation tests were imple- Polymorphism. Hamburg: Verlag Dr. Kovač; 1993. mented in the algebra environment MATLAB using 4. Balanyà J, Serra L, Gilchrist GW, Huey RB, Pascual M, Mestres F, Solé E. evolutionary pace of chromosomal polymorphism in colonizing populations tools supplied by the Statistics Toolbox [63]. of Drosophila subobscura: an evolutionary time series. Evol. 2003;57:1837–45. 5. Fontdevila A, Zapata C, Alvarez C, Sánchez L, Méndez J, Enriquez I. Genetic coadaptation in the chromosomal polymorphism of Drosophila subobscura. Supplementary information I. Seasonal changes of gametic disequilibrium in a natural population. Supplementary information accompanies this paper at https://doi.org/10. Genetics. 1983;105:935–55. 1186/s12862-020-1584-z. 6. Rodríguez-Trelles F, Alvarez G, Zapata C. time-series analysis of seasonal changes of the O inversion polymorphism of Drosophila subobscura. Additional file 1. Multiple sequence alignment of 5’proximal promoters Genetics. 1996;142:179–87. of hsp70A. 7. Rodríguez-Trelles F, Rodríguez MA. rapid micro-evolution and loss of Additional file 2. Multiple sequence alignment of 5’proximal promoters chromosomal diversity in Drosophila in response to climate warming. Evol of hsp70B. Ecol. 1998;12:829–38. 8. Balanyá J, Oller JM, Huey RB, Gilchrist GW, Serra L. global genetic change Additional file 3. Multiple sequence alignment of hsp70A 3’UTR. tracks global climate warming in Drosophila subobscura. Sci. 2006;313:1773– Additional file 4. Multiple sequence alignment of hsp70B 3’UTR. 5. Additional file 5. Statistical results for the homogeneous expression of 9. Lindquist S, Craig EA. the heat-shock proteins. Ann Rev Genet. 1988;22:631– rp49 housekeeping gene. 77. 10. Lindquist S. the heat-shock response. Annu Rev Biochem. 1986;55:1151–91. Additional file 6. Raw gene expression data. 11. Parsell DA, Taulien J, Lindquist S. the role of heat shock proteins in thermotolerance. Phil Trans Royal Soc London Series B. 1993;339:279–86. Abbreviations 12. Feder ME, Blair N, Figueras H. natural thermal stress and heat-shock protein ARE: AU-rich element; CREs: Cis-regulatory elements; HSEs: Heat shock expression in Drosophila larvae and pupae. Funct Ecol. 1997;11:90–100. elements; HSF: Heat shock factor; SNPs: Single nucleotide polymorphisms; 13. Bettencourt BR, Kim I, Hoffmann AA, Feder ME. response to natural and UTR: Untranslated region laboratory selection at the Drosophila hsp70 genes. Evol. 2002;56:1796–801. 14. Calabria G, Dolgova O, Rego C, Castañeda LE, Rezende EL, Balanyà J, Pascual Acknowledgements M, Sørensen JG, Loeschcke V, Santos M. Hsp70 protein levels and We thank two anonymous reviewers for comments on the manuscript. Rosa thermotolerance in Drosophila subobscura: a reassessment of the thermal Tarrío and Francisco Rodríguez-Trelles supplied the isogenic D. subobscura co-adaptation hypothesis. J Evol Biol. 2012;25:691–700. lines used in this work. 15. Krebs RA, Loeschcke V. acclimation and selection for increased resistance to thermal stress in Drosophila buzzatii. Genetics. 1996;142:471–9. 16. Krebs RA, Feder ME. tissue-specific variation in Hsp70 expression and Authors’ contributions thermal damage in Drosophila melanogaster larvae. J Exp Biol. 1997;200: MPGG designed the research; MPG and MPGG performed the research and 2007–15. sequence analysis, MS performed statistical analysis. MPG, MPGG and MS 17. Feder JH, Rossi JM, Solomon J, Solomon N, Lindquist S. the consequences wrote the article. All authors have read and approved the manuscript. of expressing hsp70 in Drosophila cells at normal temperatures. Genes Dev. 1992;6:1402–13. Funding 18. Rodriguez-Trelles F, Tarrio R, Santos M. Genome-wide evolutionary response This work was supported by the research grants CGL2013-42432P and to a heat wave in Drosophila. Biol Lett. 2013;9:20130228. CGL2017-89160P from the Spanish Ministerio de Ciencia e Innovación, and 19. Bettencourt BR, Feder ME, Cavicchi S. experimental evolution of Hsp70 grants 2014SGR 1346 and 2017SGR 1379 from Generalitat de Catalunya to expression and thermotolerance in Drosophila melanogaster. Evol. 1999;53: the Grup de Genòmica, Bioinformàtica i Biologia Evolutiva (GGBE). MPG was 484. supported by a PIF predoctoral fellowship from the Universitat Autònoma de 20. Sørensen JG, Loeschcke V. larval crowding in Drosophila melanogaster Barcelona (Spain). The founding bodies played no role in the design of the induces Hsp70 expression, and leads to increased adult longevity and adult study and collection, analysis, and interpretation of data and in writing the thermal stress resistance. J Insect Physiol. 2001;47:1301–7. manuscript. 21. Sørensen JG, Michalak P, Justesen J, Loeschcke V. expression of the heat- shock protein HSP70 in Drosophila buzzatii lines selected for thermal Availability of data and materials resistance. Hereditas. 1999;131:155–64. The data analysed in this study are included in this article and its additional 22. Garbuz D, Evgenev MB, Feder ME, Zatsepina OG. Evolution of files. All newly reported sequences in this study were deposited in GenBank/ thermotolerance and the heat-shock response: evidence from inter/ EMBL/DDBJ database under accession numbers: MK862365- MK862404. intraspecific comparison and interspecific hybridization in the virilis species group of Drosophila. I. Thermal phenotype. J Exp Biol. 2003;206:2399–408. Ethics approval and consent to participate 23. Gong WJ, Golic KG. loss of Hsp70 in Drosophila is pleiotropic, with effects Not applicable. on thermotolerance, recovery from heat shock and neurodegeneration. Genetics. 2006;172:275–86. Consent for publication 24. Feder ME, Cartaño NV, Milos L, Krebs RA, Lindquist SL. effect of engineering Not applicable. Hsp70 copy number on Hsp70 expression and tolerance of ecologically Puig Giribets et al. BMC Evolutionary Biology (2020) 20:17 Page 9 of 9 relevant heat shock in larvae and pupae of Drosophila melanogaster. J Exp 48. Cobreros L, Fernández-Miñán A, Luque CM, González-Reyes A, Martín- Biol. 1996;199:1837–44. Bermudo MD. a role for the chaperone Hsp70 in the regulation of border 25. Bettencourt BR, Hogan CC, Nimali M, Drohan BW. Inducible and constitutive cell migration in the Drosophila ovary. Mech Dev. 2008;125:1048–58. heat shock gene expression responds to modification of Hsp70 copy 49. Sperlich D, Feuerbach-Mravlag H, Lange P. genetic load and viability number in Drosophila melanogaster but does not compensate for loss of distribution in central and marginal populations of Drosophila subobscura. thermotolerance in Hsp70 null flies. BMC Biol. 2007;6:5. Genetics. 1977;86:835–48. 26. Evgen'ev MB, Zatsepina OG, Garbuz D, Lerman DN, Velikodvorskaya V, 50. Rego C, Balanyà J, Fragata I, Matos M, Rezende EL, Santos M. Clinal patterns Zelentsova E, Feder ME. evolution and arrangement of the hsp70 gene of chromosomal inversion polymorphisms in Drosophila subobscura are cluster in two closely related species of the virilis group of Drosophila. partly associated with thermal preferences and heat stress resistance. Evol. Chromosoma. 2004;113:223–32. 2010;64:385–97. 27. Puig Giribets M, García Guerreiro MP, Santos M, Ayala FJ, Tarrío R, 51. David J. A new medium for rearing Drosophila in axenic conditions. Dros Inf Rodríguez-Trelles F. chromosomal inversions promote genomic islands of Serv. 1962;36:128. concerted evolution of Hsp70 genes in the Drosophila subobscura species 52. Katoh K, Standley DM. MAFFT multiple sequence alignment software subgroup. Mol Ecol. 2019;28:1316–32. version 7: improvements in performance and usability. Mol Biol Evol. 2013; 28. Genissel A, LM MI, Wayne ML, Nuzhdin SV. Cis and trans regulatory effects 30:772–80. contribute to natural variation in transcriptome of Drosophila melanogaster. 53. Nicholas KB, Nicholas HBJ, DWI D. GeneDoc: Analysis and visualization of Mol Biol Evol. 2008;25:101–10. genetic variation. EMBNEW NEWS. 1997;4:14. 29. Tian S, Haney RA, Feder ME. Phylogeny disambiguates the evolution of 54. Spasic M, Friedel CC, Schott J, Kreth J, Leppek K, Hofmann S, Ozgur S, heat-shock cis-regulatory elements in Drosophila. PLoS One. 2010;5:e10669. Stoecklin G. Genome-wide assessment of AU-rich elements by the AREScore 30. Omichinski JG, Pedone PV, Felsenfeld G, Gronenborn AM, Clore GM. the algorithm. PLoS Genet. 2012;8:e1002433. solution structure of a specific GAGA factor-DNA complex reveals a modular 55. Rozas J, Ferrer-Mata A, JC S-DB, Guirao-Rico S, Librado P, Ramos-Onsins SE, binding mode. Nat Struct Biol. 1997;4:122–32. Sánchez-Gracia A. DnaSP 6: DNA sequence polymorphism analysis of large 31. Cairrao F, Halees AS, KSA K, Morello D, Vanzo N. AU-rich elements regulate data sets. Mol Biol Evol. 2017;34:3299–302. Drosophila gene expression. Mol Cell Biol. 2009;29:2636–43. 56. Nei M. Molecular Evolutionary Genetics. New York: Univ Press; 1987. 32. Becker K, Bluhm A, Casas-Vila N, Dinges N, Dejung M, Sayols S, Kreutz C, 57. Hudson RR, Slatkin M, Maddison WP. estimation of levels of gene flow from Roignant JY, Butter F, Legewie S. Quantifying post-transcriptional regulation DNA sequence data. Genetics. 1992;132:583–9. in the development of Drosophila melanogaster. Nat Commun. 2018;9:4970. 58. Hudson RR. a new statistic for detecting genetic differentiation. Genetics. 33. Franke K, Fischer K. inbreeding interferes with the heat-shock response. 2000;155:2011–4. Heredity. 2015;114:80–4. 59. Schmittgen TD, Livak KJ. analyzing real-time PCR data by the comparative 34. Vermeulen CJ, Bijlsma R. characterization of conditionally expressed mutants CT method. Nat Protoc. 2008;3:1101–8. affecting age specific survival in inbred lines of Drosophila melanogaster lethal 60. Edgington ES. Randomization Tests. 3rd ed. New York: Marcel Dekker; 1995. conditions and temperature-sensitive periods. Genetics. 2004;167:1241–8. 61. BFJ M. Randomization, Bootstrap and Monte Carlo Methods in Biology. 3rd 35. Sørensen JG, Giribets MP, Tarrío R, Rodríguez-Trelles F, Schou MF, Loeschcke ed. New York: Chapman & Hall/CRC; 2007. V. expression of thermal tolerance genes in two Drosophila species with 62. StatSoft, Inc. 2007. STATISTICA (data analysis software system), version 8.0. different acclimation capacities. J Therm Biol. 2019;84:200–7. www.statsoft.com. 36. Simões P, Pascual M. Patterns of geographic variation of thermal adapted 63. MATLAB and Statistics Toolbox Release. 2016. The MathWorks, Inc., Natick. candidate genes in Drosophila subobscura sex chromosome arrangements. BMC Evol Biol. 2018;18:60. Publisher’s Note 37. Puig M, Cáceres M, Ruiz A. silencing of a gene adjacent to the breakpoint of Springer Nature remains neutral with regard to jurisdictional claims in a widespread Drosophila inversion by a transposon-induced antisense RNA. published maps and institutional affiliations. Proc Natl Acad Sci. 2004;101:9013–8. 38. Huang YC, Dang VD, ChangNC WJ. Multiple large inversions and breakpoint rewiring of gene expression in the evolution of the fire ant social supergene. Proc R Soc B. 2018;285:20180221. 39. Lavington E, Kern AD. The effect of common inversion polymorphisms In (2L) t and In (3R) Mo on patterns of transcriptional variation in Drosophila melanogaster. G3 (Bethesda). 2017;7:3659–68. 40. Fuller ZL, Haynes GD, Richards S, Schaeffer SW. genomics of natural populations: how differentially expressed genes shape the evolution of chromosomal inversions in Drosophila pseudoobscura. Genetics. 2016;204: 287–301. 41. Karageorgiou C, Gámez-Visairas V, Tarrío R, Rodríguez-Trelles F. Long-read based assembly and synteny analysis of a reference Drosophila subobscura genome reveals signatures of structural evolution driven by inversions recombination-suppression effects. BMC Genomics. 2019;20:223. 42. Puerma E, Orengo DJ, Aguadé M. Multiple and diverse structural changes affect the breakpoint regions of polymorphic inversions across the Drosophila genus. Sci Rep. 2016;6:36248. 43. Santos M. recombination load in a chromosomal inversion polymorphism of Drosophila subobscura. Genetics. 2009;181:803–9. 44. Laayouni H, García-Franco F, Chávez-Sandoval BE, Trotta V, Beltran S, Corominas M, Santos M. Thermal evolution of gene expression profiles in Drosophila subobscura. BMC Evol Biol. 2007;10:363. 45. Ranz JM, Castillo-Davis CI, Meiklejohn CD, Hartl DL. sex-dependent gene expression and evolution of the Drosophila transcriptome. Sci. 2003;300:1742–5. 46. Michaud S, Marin R, Westwood JT, Tanguay RM. cell-specific expression and heat-shock induction of Hsps during spermatogenesis in Drosophila melanogaster. J Cell Sci. 1997;110:1989–97. 47. Lakhotia SC, Prasanth KV. tissue- and development-specific induction and turnover of hsp70 transcripts from loci 87A and 87C after heat shock and during recovery in Drosophila melanogaster. J Exp Biol. 2002;205:345–58.