Global warming is causing plastic and evolutionary changes in the phenotypes of ectotherms. Yet, we have limited knowledge on how the interplay between plasticity and evolution shapes thermal responses and underlying gene expression patterns. We assessed thermal reaction norm patterns across the transcriptome and identified associated molecular pathways in northern and southern populations of the damselfly Ischnura elegans. Larvae were reared in a common garden experiment at the mean summer water temperatures experienced at the northern (20 °C) and southern (24 °C) latitudes. This allowed a space-for-time substitution where the current gene expression levels at 24 °C in southern larvae are a proxy for the expected responses of northern larvae under gradual thermal evolution to the predicted 4 °C warming. Most differentially expressed genes showed fixed differences across temperatures between latitudes, suggesting that thermal genetic adaptation will mainly evolve through changes in constitutive gene expression. Northern populations also frequently showed plastic responses in gene expression to mild warming, while southern populations were much less responsive to temperature. Thermal responsive genes in northern populations showed to a large extent a pattern of genetic compensation, i.e. gene expression that was induced at 24 °C in northern populations remained at a lower constant level in southern populations, and were associated with metabolic and translation pathways. There was instead little evidence for genetic assimilation of an initial plastic response to mild warming. Our data therefore suggest that genetic compensation rather than genetic assimilation may drive the evolution of plasticity in response to mild warming in this damselfly species.

Pooled RNA samples (n=24) were sent to BGI, Hong Kong where 160 bp libraries were constructed using the TruSeq technology. Next, the libraries were 101 bp Paired-End sequenced on a Hiseq4000 platform. Raw reads were visually inspected with FastQC version 0.11.5 (Andrews, 2010) and trimmed using Trimmomatic version 0.36 (Bolger et al. 2014) for adapter sequences. The trimmed reads were de novo assembled with Trinity version 2.2.0 (Grabherr et al., 2011) using default parameters (kmer length = 25). Duplicate transcripts were removed by clustering the assembly using CD-HIT-EST version 4.6.6 (Li & Godzik, 2006) at 95% sequence similarity resulting in an assembly containing only the representative sequences of the detected clusters. Furthermore, Repeatmasker version 4.0.6 was used to remove transcripts containing ribosomal RNA. The completeness of the assembly was assessed using Benchmarking Universal Single-Copy Orthologs (BUSCO) version 1.22 (Simão, Waterhouse, Ioannidis, Kriventseva, & Zdobnov, 2015) which measures the number of near-universal single-copy orthologs in the assembly that are stored in the Arthropods database from OrthoDB. Only transcripts with at least one read per million in at least two samples were kept.