Transcriptomics of differences in thermal plasticity associated with selection for an exaggerated male sexual trait

Abstract


Abstract Background
The information about the magnitude of between-individual differences in thermal plasticity and identi cation of the underlying molecular mechanisms are key to understand the evolution of thermal plasticity.In particular, genes underlying variation in the physiological response to temperature can provide raw material for selection acting on plastic traits.Using RNAseq, we investigate the transcriptional response to temperature in males and females from the bulb mite populations selected for the increased frequency of one of two discrete male morphs ( ghter-and scrambler-selected populations) that differ in relative tness depending on temperature.

Results
At decreased temperature, males from ghter-selected populations showed higher transcriptomic plasticity, as indicated by a signi cant selection treatment by temperature (decreased vs. control) interaction effect on the expression of 40 genes, 38 of which were overexpressed in ghter-selected populations in response to temperature decrease.These genes were mostly associated with carbohydrate metabolism.At increased temperature, no selection-by-temperature interaction in gene expression was detected.Hence, between-morph differences in response to increased temperature are most likely determined by genes consistently differing in expression between morphs.These genes were associated with protein metabolism, ion transport, lipid metabolism and oxidoreductase activity among others.In females, we did not nd genes with selection-treatment-speci c response to temperature decrease or increase, but both sexes differed between selection treatments in expression of 79 genes.They can be treated as candidates for genes underlying temperature sensitivity of tness differences between females from ghter-and scramblerselected populations.

Conclusion
Different mechanisms underly the divergence in thermal response between populations differing in sexually selected traits prevalence at decreased vs. increased temperature.While temperature decrease was associated with higher transcriptomic plasticity of males with more elaborate armaments, differential response to temperature increase likely depended on genes associated with their distinct thermal tolerance.Selection on males drove gene expression patterns in females.These patterns could be associated with temperature-dependent tness differences between females from ghter-vs.scrambler-selected populations reported in previous studies.Our study shows that selection for divergent male sexually selected morphologies and behaviors has a potential to drive divergence in metabolic pathways underlying plastic response to temperature in both sexes.

Background
Temperature is a key environmental factor affecting various aspects of life history, physiology and behavior, with its effect being particularly important in ectotherms [1,2].Plastic response to change in thermal conditions requires substantial shifts in metabolism and, as such, is often associated with trade-offs in resource allocation (e.g., a trade-off between body maintenance and reproduction or immunity [1,3,4].Importantly, the intensity of these trade-offs may vary between individuals due to intraspeci c phenotypic differentiation in tness-maximizing strategies [5].Such phenotypic differentiation is often, at least to some extent, determined genetically [5,6,7,8,9].In consequence, within a single population/species different genotypes may vary in the physiological response to temperature change (e.g.[5,6,10]), which produces non-parallel reaction norms of different genotypes, termed also genotype-by-environment interaction.Genes responsible for this non-parallelism can be important drivers of evolution, providing raw material for selection acting on plasticity and plastic traits (reviewed in [11]).One of the factors determining the response of individuals to changing temperature is their thermal preference, i.e, the range of thermal tolerance and the thermal optimum which denotes the highest performance [12].Physiological heterogeneity in this aspect can cause complex dynamics of thermal responses within a population, driven by the strategy-speci c metabolic trade-offs.Identifying the metabolic pathways underlying these trade-offs across different temperatures is necessary to understand their prevalence and role in evolution.
Sex is one of the factors that may affect the way individuals deal with a thermal challenge.A recent transcriptomic study found that, although males and females of Drosophila montana use mainly the same mechanisms to adjust to cold, there are also important differences between the sexes in their response to temperature decrease [13].In addition to the between-sex variation, substantial differentiation of thermal plasticity may occur between individuals of the same sex.This is particularly likely in males because traits involved in sexual selection show high variation associated with condition-dependence [14,15].For example, thermal adjustment might be less costly for males in better condition as they can redirect more resources to maintenance when needed.On the other hand, high condition males often develop disproportionally large sexual ornaments/armaments, which might be particularly costly in harsh conditions (reviewed in [16]).Metabolic costs of sexual behaviors (e.g.reproductive performances or ghts) are also affected by temperature [16,17].Thus, males differing in the elaboration of the sexually selected traits might experience different stress-related trade-offs and might be selected to respond to environmental challenges (including temperature) using different physiological pathways.
In some species, variation in male phenotypes is discrete rather than continuous, with different male morphs utilizing different reproductive strategies and exhibiting morphological polymorphism (revied by [18]).Such discrete male polymorphisms provide excellent models for studying to what extent variation in male phenotype is associated with the differences in mechanisms of thermal response.Transcriptomic approaches are particularly useful here as they offer an insight into the thermal response at the scale of the entire genome.At the level of gene expression, intra-sexual differences in thermal response might manifest in two, nonexclusive, ways.First, different male phenotypes (morphs) might adjust the expression of different sets of genes at different temperatures and the scale of this adjustment might differ (leading to different levels of transcriptomic plasticity).This would result in a signi cant morph-by-temperature interaction in the expression level of these genes.Second, some genes that differ in expression between morphs but do not change expression between temperatures (or change expression concordantly in both morphs) may underlie the differences in the cost of temperature adjustment and/or alter thermal tolerance.In this case, we expect gene expression differences between morphs, but not the temperature-by-morph interaction.
Here, we use the male-dimorphic bulb mite, Rhizoglyphus robini, to investigate the transcriptional response to temperature in two discrete sexually selected phenotypes.The bulb mite is characterized by two heritable [7,19,20] male morphs differing in the morphology of a sexually selected trait [21], reproductive behavior [22] and gene expression pro les [23].Fighters possess an elaborated sexually selected trait -enlarged third pair of legs that are used to attack (and sometimes kill) other males, while scramblers have unmodi ed legs, are less aggressive and unable to kill rivals [22].Generally, ghters develop from larger juveniles [24] and are of better condition than scramblers (but there is also a signi cant parental effect on morph expression), achieving higher reproductive success and winning reproductive competition in mixed-morph populations [25].However, the bene ts of the ghter strategy seem to depend on the thermal environment, as the frequency of ghters in populations decreases with increasing temperature [26,27].This effect can be partly attributed to developmental effects of environment (likely mediated by the effect of temperature on body size), but it has also a genetic component.Indeed, the reproductive advantage of ghters over scramblers is higher at lower temperature [27], which may be associated with their higher transcriptomic plasticity (ability to adjust gene expression) at that temperature (see [28]).In addition, the metabolic costs of production and/or maintenance of ghters' weapon and performing aggressive behavior likely increase with temperature, contributing to the observed thermal pattern of ghter male tness.Indeed, genes involved in general metabolism are overrepreseted among those differentially expressed between morphs [23], supporting this hypothesis.These metabolic differences should be mirrored by females from ghter and scrambler selected populations.Previous studies found that genetic variants underlying ghter morph have negative consequences for female tness being a classic example of intralocus sexual con ict [29, see also 30].Fecundity of females from ghter-selected populations is reduced [29], especially at increased temeperature [31].However, at lower tempearures the opposite patterns is observed -females from scrambler-selected populations have lower tness than that of females from ghter-selected populations [31], indicating that pleiotropic effects of genes associated with morph expression are environment-sensitive.Thus, as females are not morphologically differentiated between each other, there must be a substantial difference in how the two male morphs respond to temperature that extends beyond the sexual weapon maintenance itself and involves some metabolic pathways active in both sexes.On the other hand, recent studies suggest that sexual antagonism might be compensated by selection for increased condition and/or more e cient pury ng selection in gther-selected poulations [21].This additional factor should, however, not affect the above predctions, which should then hold even if the apparent differences between females from gher-and scrambler-selecred populations are absent.
Here, we investigate trascriptomic basis of the differences in thermal responses driven by selection for male morph expression.We assess gene expression in males and females from populations selected for increased/decreased frequency of ghters ( ghterand scrambler-selected populations, respectively) at the temperature at which they have evolved (23°C) as well as at decreased and increased temperature,18°C and 28°C, the latter being highly stressful [32].In our analyses, we address the four speci c hypotheses. 1) As ghters do better at lower temperature, we hypothesise that they exhibit higher transcriptomic plasticity than scramblers when temperature decreases (i.e., ghters show higher change in expression of some genesthan scramblers after they are subjected to decreased temperature).2) Such ghter-speci c response should be also visible in gene expression of females from ghter-selected populations, which could explain their tness advantage compared to females from scrambler-selected populations at low temperature.3) Lower tness of ghters at increased temperature might be associated with lower transcriptomic plasticity compared to scramblers at higher temperature, and 4) this effect should also be observed in females.The lack of suport for these hypotheses would suggest that inter-morph differences in response to temperature result from the differential expression of genes that in uence thermal tolerance and/or costs of thermal adjustment rather than differences in thermal plasticity levels.Such genes should be included among those differentially expressed between morphs.To test the above hypotheses, we identify genes with morph-speci c expression response to temperature (which is indicated by a signi cant selection treatment by temperature interaction at 18°C vs. 23°C and 28°C vs. 23°C (hypotheses 1 and 3, respectively) and check whether the same genes show selection treatment-speci c thermal response in females from the same populations (hypotheses 2 and 4).We also discuss which of the genes differing in expression between selection treatments/morphs regardless of temperature (the effect of selection treatment, no interaction) could be responsible for between-morph differences in thermal tolerance.

Results
We sequenced a total of 79 male and female RNA samples, each consisting a pool of ca.50-60 individuals, originating from 4 ghter-selected and 4 scrambler-selected populations described in [21].For each population, gene expression was investigated in samples of individuals that had developed and been maintained at 18°C, 23°C or 28°C.We obtained a total of 2.3 billion reads (2x150bp), which mapped to 62,708 gene models (henceforth called genes).After ltering out the genes with the total number of mapped reads across all the samples of a given sex lower than 10, our dataset consisted of 58,672 and 59,823 genes in the male and female dataset, respectively.The per sample number of reads ranged from 16 M to 40.8 M, with an average of 29.8 M (SD 5.0 M).

Expression differences between selection treatments
In males, there were 2,632 genes (ca.4% of all assayed genes) that differed in expression between ghter-and scrambler-selected populations regardless of temperature.Of these genes, 1553 (59%) were overexpressed in ghter-selected populations, whereas 1079 (41%) in scrambler-selected populations.The expression level differed between the ghter-and scrambler-selected populations more than two-fold for 587 (22.3%) and more than four-fold for 129 (4.9%) differentially expressed genes (Fig. 1A).
Genes differentially expressed between selection treatments in males were enriched for GO categories (biological processes and metabolic functions) associated with protein metabolism, ion transport, metabolism of lipids and oxidoreductase activity.The full list of signi cantly enriched categories is given in Table 1.Among genes differing in expression more than two-fold, we found an overrepresentation of the following categories: proteolysis, cysteine-type peptidase activity, serine-type endopeptidase activity, lipid metabolic process, oxidoreductase activity (acting on paired donors).Gene expression in females was much less divergent across selection treatments (Fig. 1A, see also Supplementary material Fig. S2).There were 202 genes (ca.0.3% of all assayed genes) differentially expressed between selection treatments in females, 112 of them (55.4%)being overexpressed in ghter-selected populations and 90 (44.6%) in scrambler-selected populations.Expression of 77 (38.1%) genes differed more than two-fold and 19 (9.4%) more than four-fold.No GO categories were signi cantly enriched for genes that differed between selection treatments in females.Genes with known GO terms for which the change in expression was the highest included such categories of biological processes and molecular functions as nucleic acid binding, carbohydrate metabolic process and hydrolase activity (hydrolyzing O-glycosyl compounds), ATP binding and ATP hydrolysis activity (3 of the genes overexpressed in ghter selected populations) and DNA-binding transcription factor activity, zinc ion binding, sequencespeci c DNA binding and regulation of transcription (DNA-templated), iron ion binding, oxidoreductase activity (acting on paired donors, with incorporation or reduction of molecular oxygen), heme binding and obsolete oxidation-reduction process (two of the genes overexpressed in scrambler-selected populations).
We detected 79 genes that differed in expression between selection treatments for both males and females (expression of 30 and 11 differed more than two-fold and more than four-fold, respectively) and in all cases the direction of change was the same in both sexes (Fig. 2A).These genes were signi cantly enriched in N,N-dimethylaniline monooxygenase activity, obsolete oxidationreduction process, hydrolase activity (hydrolyzing O-glycosyl compounds), oxidoreductase activity (acting on paired donors) and carbohydrate metabolic process.Genes with the highest overexpression in the ghter-selected populations were involved in hydrolysis, carbohydrate metabolism and the usage of ATP, while one of the two genes with the highest overexpression in the scrambler-selected populations was cytochrome P450 3A9-like isoform involved in oxidation-reduction process (the other gene could not be annotated).

Expression differences between temperatures
In males, differential expression analyses identi ed as many as 16,830 genes (29% of all assayed genes) differentially expressed between the decreased (18°C) and control (23°C) temperature, 9,127 (54.2%) of which showed higher expression at 18°C.In contrast, only 5,263 (9%) genes were differentially expressed between the increased (28°C) and control temperature, with similar numbers of genes increasing and decreasing expression at 28°C (2,328 and 2,881, respectively).In females, we found 15,599 (26%) genes that differed in expression between the decreased and control temperature, 6,971 (44.7%) of which increased and 8,628 (55.3%) decreased their expression at 18°C.We found 9,375 (16%) genes that differed in expression between the increased and control temperature, with 4,719 of them (50.3%)overexpressed and 4,657 (49.7%) underexpressed at 28°C.The number of signi cantly differentiated genes that changed expression between temperatures more than two-or four-fold are given in Table 2.The functional categories associated with expression differences between temperatures were diverse and to a large extent depended on whether the control temperature was contrasted with 18°C or 28°C (Supplementary material Fig. S3 and S4).The fraction of genes that responded to temperature in both males and females was higher at 18°C than at 28°C (Table 2).
Selection treatment speci c response to temperature In males, 40 genes (Table 3) showed signi cant temperature by selection treatment interaction at the decreased (18°C vs. 23°C) temperature.For most of them (38 genes), the difference between selection treatments was higher at 18°C than at control temperature due to increased expression of ghter-selected populations (Fig. 1B).No genes with signi cant temperature by selection treatment interaction were found when comparing the increased (28°C) to control (23°C) temperature.Genes that showed a selection-treatment-speci c response to temperature were enriched in hydrolase activity (hydrolyzing O-glycosyl compounds) and carbohydrate metabolic process.For none of the genes with signi cant temperature by selection treatment interaction in males this interaction was signi cant in females. .However, intra-speci c variation in response to thermal and other environmental stimuli has received surprisingly little attention (with exception of some studies of variation in response to nutritional conditions; e.g.[40,41,42].This lack of research is unfortunate as such variation is likely to provide the raw material for adaptation [43,11].Our study attempts to ll this gap by analyzing transcriptomic responses to thermal challenge in bulb mites from populations experimentally evolved with speci c regard to morph expression .

Transcriptomics of differences in thermal plasticity between the male morphs
We show that most transcriptional response to temperature is shared between both male morphs, which differ by a sexually selected trait and associated reproductive tactics.Importantly, however, there are also signi differences in how male morphs respond to thermal challenge.This agrees with the observations that relative tness bene ts of expressing a given morph-speci c tactic are temperature-sensitive [27], potentially driving variation in morph frequencies and associated metabolic variation in habitats where temperature varies over space and/or time as well as between habitats differing in thermal conditions.
Notably, selection-treatment-speci c response to temperature in gene expression was observed only when temperature was decreased.We found forty such genes, which likely underly between-morph differences in plasticity levels.These genes were enriched in GO terms associated with carbohydrate metabolism and were involved in hydrolyzing O-glycosyl compounds which bind carbohydrates to other molecules.This strongly suggests that the glycosidic bonds synthesis can play important role in shaping differences in cold adjustment between ghters and scramblers.Interestingly, a vast majority of these forty genes (38) were overexpressed in ghter-selected populations, but not in scrambler-selected ones in response to decreased temperature.This, together with the fact that ghter advantage in direct competition with scramblers over the access to females is higher at lower temperature [27] suggests that ghters are better in plastic adjustment of their carbohydrate metabolism to temperature decrease (consistently with our hypothesis 1).It is possible, for example, that ghters are able to switch to the more e cient usage of carbohydrates at lower temperature and, hence, to transfer these resources more e ciently into tness-related traits (e.g., to acquire more energy to use in ghts with rivals).Indeed, the e cacy of carbohydrate absorption and utilization in ectotherms depends on temperature and the usage of carbohydrates is a plastic trait [e.g.44].Our study provides evidence that the levels of such plasticity might vary within a species, translating to temperature-dependent tness differences between phenotypes.
Some of the genes that consistently differ in expression between ghter-and scrambler-selected populations (without selection treatment by temperature interaction) might also be responsible for inter-morph differences in thermal adjustment costs, thermal tolerance and/or thermal optima.Such genes might contribute to tness differences between morphs at both decreased and increased temperature.They may be especially important at elevated temperature, in which we did not nd genes with signi cant selection treatment by temperature interaction, despite clear evidence that ghter tness advantage decreases at elevated temperatures [see 26].Among genes differing in expression between selection treatments we found signi cant enrichment of some metabolic-related pathways, in particular genes associated with different aspects of amino acid and protein metabolism.These could contribute to between-morph differences in thermal tolerance as temperature affects the energetic costs of protein conformational changes and in uences the rate of protein unfolding [2].Protein synthesis and degradation is the main contributor to the basal metabolism in humans and rats [2].Importantly, protein metabolism also plays a role in thermal plasticity of ectotherms as temperature has been shown to in uence the abundance of protein and RNA in taxonomically diverse species, indicating thermally induced changes in protein synthesis and turnover [45].Similarly, ion transport can play a role in between morph differences in thermal performance.This is because costs of ion pumping decrease with decreasing temperature [2].
Our experiment compared gene expression patterns between males from morph-selected populations rather than between male morphs per se.Most of expression differences between them are likely associated with the expression of a given morph or tness bene t to a given morph.However, some differences might result from the evolution in populations differing in morph frequency, i.e. from differences in the social context.For example, as the proportion of ghters was reduced in scrambler-selected populations, the risk of attack was lower and the intensity of sperm competition was likely higher in these populations.As in our design the effects of morph expression and tness bene ts to a morph cannot be separated from the effect of a social context, this might be perceived as a limitation of the study.However, as temperature affects morph proportions in populations [26,27], the effects of morph expression and indirect selection resulting from varying morph frequencies are likely to interplay in natural mite populations and thus investigating them in concert is fully justi ed.
Transcriptomic differentiation and its thermal plasticity in females from scrambler-and ghter-selected populations Transcriptomic pro les of females were far less divergent between selection treatments, as there were 10 times fewer genes differentially expressed between females from scrambler-and ghter-selected populations than between males from these treatments.This is consistent with previous ndings in the same species, which identi ed only 11 genes differentially expressed in females from ghter-and scrambler-selected lines compared to 532 genes differentially expressed between males from the two treatments [46].The numbers of differentially expressed genes are comparable between [46].and our study when considering genes with at least two-fold expression difference (Fig. 2A).A limited divergence of gene expression between females is also consistent with the lack of morphological differentiation between females from the two selection treatments.However, 202 genes differing in expression between selection treatments in females (and 19 of them differing more than 4-fold) indicate that selection for male morph prevalence is re ected in female expression pro les.While some expression differences between females from the two selection treatments migh be a correlated response to selection on males, others may be driven by selection pressures resulting from different social context or by pleiotropic effects of genes associated with male morph expression or bene cial to a given morph.
If expression differences are driven by the social context-dependent selection the genes and pathways affected by selection may differ between females and males.In principle, the pleiotropy-effect candidates can be identi ed among genes that diverged in expression between morph-selected populations in both sexes or the genes with a signi cant treatment by temperature interaction in both sexes.Interestingly, however, no genes were identi ed in the latter category.Notably, among 79 genes differentially expressed between ghter and scrambler selected populations in both sexes, the expression always changed in the same direction in both males and females, and many such genes exhibited large expression differences (see Fig. 2A, Table 4).Even more intriguingly, 13 of these genes showed also signi cant change in response to temperature decrease in both males and females and eight out of ten with known GO categories were associated with nucleic acid binding, suggesting their role in transcription regulation or related processes.Thus, pleiotropic effects of genes associated with these processes might play an important role in plastic response to low temperature in females.
Table 4 Overview of genes with more than 2 fold difference in expression between selection treatments in both sexes.Genes with four fold or higher difference in expression and in bold.A majority (two thirds) of the 79 genes differentially expressed between ghter-and scrambler-selected populations in both sexes had higher expression in ghter-selected populations, suggesting the ghter expression is associated with large pleiotropic effects on females.Some of these effects may be related to sexual con ict [29].However, as the current study did not nd evidence for increased sexual con ict in the ghter-selected populations studied here (possibly due to large population sizes), most of these effects must be associated with indirect selection for better condition and increased purifying selection in ghter selected populations [21].The association of these genes with oxidation-reduction processes and carbohydrate metabolism suggests that selection for an aggressive and highly competitive phenotype in males is re ected in metabolism of females.

Name
Among all genes differing in female expression between selection treatments no signi cant enrichment of any GO category was found.This suggests that selection imposed by the treatments acts on diverse aspects of female physiology.This view is further supported by annotations of genes with the highest between-treatment differences in expression.They were associated with such general metabolic categories as nucleic acid binding, regulation of transcription, carbohydrate metabolism, oxido-reductase processes, and ATP hydrolysis and binding.Genes from these categories play a role in general cell metabolism, affecting many physiological processes, and are likely expressed in multiple tissues.Consistent with this, morph frequencies in a population affect diverse aspects of social environment, such as relative levels of precopulatory and postcopulatory sexual selection, copulation duration [47] or chemical communication in a population [48] hat are expected to drive complex responses in female physiology and behavior.
Some of the expression differences between females from the two selection treatments can be associated with thermal speci city of the effect that morph selection was shown to exert on female tness in [31].The exact mechanism behind this pattern is not known.It has been hypothesized to be related to metabolic costs: ghters may harbor genes increasing metabolism that have the same or similar effect in both sexes.These costs might be especially pronounced at higher temperatures, when energy expenditure is high.The fact that genes overexpressed in females from ghter-selected populations are related to key metabolic pathways supports this intuition.Moreover, the processes in which differentiated genes are involved are either energetically demanding and hence sensitive to temperature (e.g.transcription processes) and/or are the main drivers of energy acquisition and use (e.g.carbohydrate metabolism, ATP-related processes).However, the generality of gene ontologies makes it di cult to identify precise pathways, beyond their involvement in processes associated with transcription, carbohydrate metabolism and oxidation reduction processes.

Conclusion
The understanding of the magnitude of between-individual differences in thermal plasticity and identi cation of the underlying physiological mechanisms are key to understand the evolution of plasticity.We found that in the bulb mite transcriptional response to temperature is to a large extent consistent in populations where males utilize different reproductive tactics, but the populations were clearly differentiated in their expression of many genes involved in diverse metabolic processes.This differentiation can play an important role in shaping intra-sexual differences in thermal tolerance.
Our results indicate that the differences in thermal response of two male morphs at decreased vs. increased temperature are caused by different molecular mechanisms.At lower temperature, males from ghter-selected lines show higher transcriptomic plasticity, likely leading to increased tness of ghter vs. scrambler males (in accordance with hypothesis 1).Between-morph tness differences at increased temperature are, on the contrary, shaped by a set of genes differing in expression between morphs (contrary to hypothesis 3) or by the deferences in protein-coding sequences; some of these genes probably underly thermal tolerance and cost of thermal adjustment.No selection-treatment-speci c thermal adjustment, predicted by hypotheses 2 and 4, was detected in females.Therefore, our results suggest that temperature-dependent tness differences in females from ghter-and scrambler-selected populations are caused by genes that show similar expression differences between morphs regardless of temparature.
In general, we show that selection for divergent male morphologies and behaviors can drive divergence in metabolic pathways underlying the plastic response to temperature.Such morphological and behavioral differentiation in uences diverse aspects of physiology and metabolism, which could account for intra-species differences in thermal tolerance.Selection acting on males can also affect female physiology, shaping expression of genes associated with different aspects of metabolism that can potentially affect female tness at different temperatures.

Experimental populations
The experiment was performed on the bulb mite populations selected for the increased prevalence of ghters or scramblers [21].
Four independent populations have been maintained in each selection treatment (eight in total).Each generation, 500 females and 500 males of a given morph were placed together in 200 ml plastic containers (110x80x40mm) with a plaster of Paris soaked with water and yeast as the food source.Individuals were left to interact for 6 days, after which they were transferred to fresh containers, where females could lay eggs.The transfer ensured that offspring hatched from eggs laid during the rst ve days of the procedure did not contribute to the next generation.After 24h of oviposition, the adults were removed and the eggs were left to develop.Two weeks later, adult individuals were used to establish the next generation.
The populations have been maintained at 23°C under constant darkness and high humidity at the Faculty of Biology, Adam Mickiewicz University in Poznań, Poland.After 13 generations the ghter-selected populations contained ca.99% ghters, while scrambler-selected populations contained ca.83% scramblers.At generation 13, a subset of each population was transferred to the Institute of Environmental Sciences, Jagiellonian University, Kraków, Poland, where the populations were maintained without any selection for ca. 5 weeks (ca. 2 generations), until the start of the experiment.The populations were maintained at 23°C in identical containers as in Poznań and could expand freely.Then, population densities were reduced by transferring a part of each population to a new container to prevent overcrowding.This procedure was done after the rst pair of ghter-and scrambler-selected populations (block 1; see below) were subjected to the experiment, ca. a week before the assay on the second pair of populations (block 2).

The experiment
For each ghter and scrambler-selected population, gene expression was measured at three temperatures: 18°C, 23°C (the temperature in which the population evolved) and 28°C (Fig. 2).While 28°C is highly stressful for lab-adapted mites [32] as it is associated with both decreased egg-laying rate and shorter life, 18°C imposes only a mild, if any, stress -egg-laying rate is lower than at 23°C, but it is compensated by an increased lifespan (see [31] for more discussion).
The experiment was performed in 4 blocks, each block consisted of a single ghter-and a single scrambler-selected population (Fig. 2A); the blocks started ca.one week each after another.In each block, previously mated females were placed in 6 containers (30 mm in diameter, 100 females per container) at 23°C to lay eggs for two days, after which they were removed.The containers with eggs were placed at different experimental temperatures to develop.Two containers were put at each of the temperatures, with random assignement of the containers to thermal treatments.
Gene expression was measured in adult males and females at each temperature (Fig. 2B).When most males reached maturity, ca.50 individuals from each container were transferred to RNAzol®RT and frozen at -80°C (two samples per temperature, one per each of the containers).For females, an additional procedure was performed before collecting the samples for RNA extraction to minimize the effect of possible interaction with males from their own populations that differed in morph proportions depending on selection treatment and temperature.One to two days old females were transferred to a container with ca.70 scrambler males from the stock population acclimated to the experimental temperature for at least 24h.The females were transferred on consecutive days until we collected ca.110-120 females per container.We then waited for another 3 days, during which individuals could interact freely.Then a sample of ca.60 females from each container was placed in RNAzol®RT and frozen at -80°C.

RNA extraction and sequencing
RNA extractions were performed using the RNAzol®RT protocol [49].Samples in RNAzol were homogenized when still frozen, which increase the RNA yield.The quality of extracted RNA was checked with the Agilent 2100 Bioanalyzer System and samples showing signs of RNA degradation were excluded.As a result, 79 samples were retained.We prepared indexed stranded RNAseq libraries using NEBNext® Ultra™ II Directional RNA Library Prep Kit for Illumina (New England BioLabs) and combined the libraries into a single pool for sequencing.Sequencing was performed on six lanes of the MGISEQ-2000 platform [see 50] for comparison between MGISEQ-2000 and Illumina HiSeq 4000).

Read trimming and mapping
The quality of sequencing results was assessed using FastQC software.Reads were then trimmed using trimmomatic (version 0.39) using TrueSeq3 adapters and maximum mismatch of 2, palindrome clip threshold set to 30 and simple clip threshold set to 10. Sliding windows trimming was performed with window size of 4 and the required quality set to 15.We also removed low quality bases from the beginning and the end of reads, with the minimum quality required to keep a base set to 3. Finally, reads shorter than 36 bp were discarded.Filtered reads were then mapped to the reference bulb mite genome [21] using STAR (version 2.7.6a; [51]), keeping only those junction-containing reads that passed ltering steps (--outFilterType BySJout), setting minimum overhang for junctions to 8 (--alignSJoverhangMin 8, --alignSJDBoverhangMin 8), minimum and maximum intron size to 20 and 100000 respectively (--alignIntronMin 20 --alignIntronMax 100000), maximum gap between two mates to 100000 (--alignMatesGapMax 100000) and the threshold for the ratio of mismatches to read length to 0.04 (--outFilterMismatchNoverReadLmax 0.04).

Differential gene expression and gene ontology enrichment analyses
Bam les were read into R with Rsamtools (1.34.1).Exon locations were extracted from the genome annotation le with GenomicFeatures package (1.34.8).Read counts for each sample for each gene were obtained using summarizeOverlaps function from Bioconductor package GenomicAlignments (1.18.1), using "union" mode of counting and ignoring read strands.
Differential expression analysis was done with Bioconductor package DESeq2 (1.30.1; [52] for males and females separately.First, differentiation between samples was analyzed with generalized Principal Component Analyses, which allows for the reduction of dimensionality of non-normally distributed data [53], using glmpca function from Glmpca package [53].Variance stabilizing transformation of raw read counts [54] with vst function was performed prior to generalized PCA.Then, differential expression analysis was performed using DESeq function.Reads from all the lanes were collapsed within each sex, population and temperature combination.The model included selection treatment, temperature and their interaction.We de ned a set of contrasts to test for differential expression between selection treatments across temperatures, between increased and control temperature and decreased and control temperature.Additionally, to look for the genes in which temperature-related expression patterns differed between selection treatments, we tested for a temperature by selection treatment interaction at both decreased (18°C vs 23°C) and increased (28°C vs 23°C) temperature.In all these comparisons, we used the false discovery rate (FDR) ≤ 0.05 to identify genes signi cantly differing in expression.The analyses were then rerun with Log2Fold change threshold set to 1/-1 (two-fold change in expression) and 2/-2 (four-fold change in expression) to distinguish genes with subtle expression differences from those for which the change in expression is more substantial.
To search for signi cantly enriched gene ontology terms, we used gene ontology terms downloaded from the functional annotation of the bulb mite genome [21] and Bioconductor package goseq (1.42.0; [55], which enables gene ontology analyses for non-model organisms and accounts for gene length bias in the data.Null distribution for GO category membership was generated using the Wallenius distribution and the FDR cutoff was set at 0.02.

Declarations
Ethics approval and consent to participate Not applicable Consent for publication Not applicable Figures Figure 1 Overview of genes differentially expressed between selection treatments in males and females (panel A) and genes showing signi cant selection treatment by temperature interaction (panel B).Black rectangles represent genes with more than 4-fold difference in expression, orange -genes with difference in expression higher than two-fold and lower than four-fold, bluesigni cantly differentiated genes that did not exceed a two fold difference in expression threshold.

Table 1
Gene ontology categories signi cantly enriched in genes differentially expresses between morph selection treatments in males.

Table 2
Numbers of differentially expressed genes for different temperatures.

Table 3
Overview of genes with signi cant selection treatment by temperature (cold vs. control) interaction in males.Genes with the highest changes in expression (more than four-fold change) are marked in bold.Log2 fold changes > 0 denote larger difference between temperatures in ghter selected lines, Log2 fold changes < 0 denote larger difference between temperatures in scrambler selected populations.
35,34]37,ion39]nscriptomics enables studying mechanisms of environmental adjustment at the scale of entire genomes [see33,34], providing a unique opportunity to investigate the interplay between different physiological pathways during the response to changing environment.Such approaches have allowed better understanding of the mechanisms underlying thermal plasticity [e.g.35,36,37,38,39]