The lack of genetic variation underlying thermal transcriptomic plasticity suggests limited adaptability of the Northern shrimp, Pandalus borealis

Introduction Genetic variation underlies the populations’ potential to adapt to and persist in a changing environment, while phenotypic plasticity can play a key role in buffering the negative impacts of such change at the individual level. Methods We investigated the role of genetic variation in the thermal response of the northern shrimp Pandalus borealis, an ectotherm species distributed in the Arctic and North Atlantic Oceans. More specifically, we estimated the proportion transcriptomic responses explained by genetic variance of female shrimp from three origins after 30 days of exposure to three temperature treatments. Results We characterized the P. borealis transcriptome (170,377 transcripts, of which 27.48% were functionally annotated) and then detected a total of 1,607 and 907 differentially expressed transcripts between temperatures and origins, respectively. Shrimp from different origins displayed high but similar level of transcriptomic plasticity in response to elevated temperatures. Differences in transcript expression among origins were not correlated to population genetic differentiation or diversity but to environmental conditions at origin during sampling. Discussion The lack of genetic variation explaining thermal plasticity suggests limited adaptability in this species’ response to future environmental changes. These results together with higher mortality observed at the highest temperature indicate that the thermal niche of P. borealis will likely be restricted to higher latitudes in the future. This prediction concurs with current decreases in abundance observed at the southern edge of this species geographical distribution, as it is for other cold-adapted crustaceans.


Introduction
Increasing temperature is one major manifestation of global changes. To cope with increasing temperatures, some species may shift their distribution tracking their thermal niche, while others may persist in their local habitat (Perry et al., 2005;Lima et al., 2007; OPEN ACCESS EDITED BY Vera Maria Fonseca Almeida-Val, National Institute of Amazonian Research (INPA), Brazil Bonebrake et al., 2018). Predicting a species' ability to cope with global changes requires attaining an in-depth critical understanding of the mechanisms underpinning acclimatization and adaptation allowing species to persist and thrive in changing environments (Román-Palacios and Wiens, 2020). Acclimatization relies on phenotypic plasticity, the capacity of a given genotype to change its phenotypes according to the environment (Levins, 1963;Bradshaw, 1965). Phenotypic plasticity can act as a buffer for the negative effect of environmental changes on the performance and fitness of organisms (Ghalambor et al., 2007). Phenotypic plasticity is adaptive in environments that change predictably, i.e., when the cues for the development of a phenotype are reliable (Gavrilets and Scheiner, 1993;Reed et al., 2010;Scheiner and Holt, 2012). It represents a key mechanism for individuals and populations to persist in changing environments, and ultimately influencing both population dynamics and extinction risk (Reed et al., 2010;Ashander et al., 2016;Rescan et al., 2020).
Global changes are characterized by increased mean temperature, amplified environmental variation, and greater frequency and intensity of extreme thermal events (Alley, 2000;Rahmstorf and Coumou, 2011;Hobday et al., 2016;IPCC, 2022). These phenomena result in increasing environmental unpredictability and new selection horizons for living organisms. Under unpredictable changes, the reliability of environmental cues required to trigger the development of a given phenotype will diminish. In those conditions, phenotypic plasticity may be maladaptive, favorizing the evolution toward alternative strategies (Gavrilets and Scheiner, 1993;Reed et al., 2010;Scheiner and Holt, 2012;Tufto, 2015;Dey et al., 2016;Leung et al., 2020). For instance, individuals from multiple populations may evolve to a single phenotype that represents the optimal compromise among environmental conditions (e.g., Scheiner and Yampolsky, 1998). The lack of variation in plasticity levels may have for consequence to reduce the potential for evolutionary responses to climate changes (Oostra et al., 2018). Characterizing a species plasticity to temperature and the genetic variation underpinning these responses may help to predict a species' vulnerability to global changes (Landy et al., 2020).
Genomic tools can help to improve prediction of organisms' response to global changes by disentangling the molecular mechanisms involved in acclimatization and adaptation (De Wit et al., 2016;Kelly, 2019;Waldvogel et al., 2020). Gene expression is considered as a molecular phenotype to study phenotypic plasticity (Aubin-Horth and Renn, 2009;Pavey et al., 2010). Changes in gene expression underly phenotypic plasticity, as most changes in mRNA expression are associated with changes in protein production and abundance. It is shaped by both environmental and genetic components and evolves in response to selective pressures (Gerken et al., 2015;Leder et al., 2015). Alternatively, genome sequencing allows the inference of the evolutionary potential of a population, by characterizing and contrasting genetic compositions. Higher levels of genetic variation across a species distribution usually enable greater levels of phenotypic variation (Hansen, 2006), and a higher capacity of a species to undergo further adaptation via natural selection (Fisher, 1930;Hughes et al., 2008). Combined, these genomics tools help shedding light on the role of phenotypic plasticity and its genetic variation in the mechanisms of species' responses to global changes.
Marine ectotherms are models of choice to understand the role of phenotypic plasticity and genetic variation in a species vulnerability to global changes (Huey et al., 2012;Paaijmans et al., 2013;Seebacher et al., 2014;Pinsky et al., 2019). In fact, ectotherms are particularly sensitive to temperature fluctuations due to their limited ability to use metabolic heat and maintain their body temperature (Angilletta, 2009). The adjustment of their body temperature to environmental conditions is therefore usually coupled with physiological and behavioral plasticity (Huey and Stevenson, 1979;Stevenson, 1985;Amarasekare and Savage, 2012;Abram et al., 2017). However, the level of plasticity varies largely among ectotherms (e.g., Somero, 2005Somero, , 2010Logan et al., 2015).
The northern shrimp Pandalus borealis (Krøyer, 1838) is a circumboreal species distributed in the colder regions of North Atlantic Ocean. Adults and juveniles experience temperatures ranging from−1.7 to 11.1°C but are found to be most abundant between 0 and 5°C in the northwest Atlantic (Allen, 1959;Haynes and Wigley, 1969;Shumway et al., 1985;Garcia, 2007). Previous studies have shown that P. borealis is sensitive to increasing temperature, resulting in reduced survival, lowered recruitment, increased growth and metabolism rates Ouellet and Chabot, 2005;Daoud et al., 2007Daoud et al., , 2010Richards, 2012;Arnberg et al., 2013;Dupont-Prinet et al., 2013;Chemel et al., 2020). As a consequence, an important reduction in abundance has been observed over the last decade in the species southern distribution (Apollonio et al., 1986;Koeller, 2000;DFO, 2021DFO, , 2022Richards and Hunter, 2021). This together with the high economic value of the species highlight the urgency of understanding the vulnerability of the northern shrimp to increasing temperatures of populations at higher latitudes.
In this study, we investigated the effect of elevated temperature to shed light on the molecular mechanisms underlying thermal plasticity in P. borealis and predict its vulnerability to ocean warming. We hypothesize a high level of transcriptomic plasticity in the northern shrimp in response to temperature, as an adaptation to the predictable variation of temperature that the species encounters on an annual or a life-cycle basis (Koeller et al., 2009). We thus predicted differences in plasticity among populations from contrasted environments and underlay by genetic variation. To test these hypotheses, we collected shrimp from three origins in a study area in the northwest Atlantic, with one origin close to a steep climatic gradient (Stanley et al., 2018). We exposed shrimp to three temperatures under laboratory conditions, mimicking current conditions, and future ocean warming scenarios. We built a reference transcriptome and identified differentially expressed transcripts (DETs) associated with temperature treatment, shrimp origin, or their interaction. The DETs were grouped in modules of expression profiles and their functions were explored to understand the molecular mechanisms underlying response to temperature. We also assessed the standing genetic variation and how it could explain temperature transcriptomic response to clarify the role of acclimatization and adaptation in thermal plasticity of P. borealis. After approximately eight weeks, female shrimp were transferred into experimental tanks kept at 2, 6 or 10°C (two replicate tanks per temperature treatment) for 30 days. As this experiment was part of a project also dealing with ocean acidification, all tanks were maintained at a slightly acidified pH, 7.75 (total scale), matching current conditions of the deep waters of the Gulf of St. Lawrence (Mucci et al., 2011(Mucci et al., , 2018. The coldest temperature, 2°C, represents the temperature where shrimp are most abundant along the Labrador and Newfoundland coast (Orr and Sullivan, 2013) and within the range where most SLE shrimp were found in 2008-2021 (1-5.5°C; DFO, 2022). The intermediate temperature represents a predicted 4°C increase by the end of this century according to the RCP 8.5 scenario (IPCC, 2014) for shrimp living close to 2°C presently, and it is also the current temperature (5-7°C) in the deep waters of the Gulf of St. Lawrence, where shrimp populations are found (Bourdages et al., 2022). The highest temperature, 10°C, represents predicted conditions at the end of the century for populations that are already exposed to 6°C (Lavoie et al., 2020). Dead and moulting individuals were counted daily. Dead individuals were frozen until waste collection according to MLI policy. Shrimp from different origins were exposed to the same experimental conditions but they could not be exposed simultaneously, due to space limitations within the experimental system and also to keep the period between capture and the beginning of the experimentation constant. Consequently, exposure to all the three experimental temperatures occurred in August 2018, April 2019, and December 2019 for shrimp from SLE, ESS, and NNC, respectively. After 30 days of exposure, six randomly chosen individuals (three from each tank) were collected from the experimental system and immediately flash-frozen in liquid nitrogen. Samples were then stored at−80°C until nucleic acid extraction was conducted. The duration of exposure was chosen to mimic northern latitude deep-water temperature fluctuation where intra-seasonal variability occurs at a period of 10-30 days (Fischer et al., 2004). At the end of the experiment, we did not find any differences in shrimp size (cephalothorax length = 23.91 mm, s.d. = 1.50) among temperature (F (2,47) = 2.097, p = 0.134) or origin (F (2,47) = 1.845, p = 0.169).

DNA/RNA extraction, sequencing and bioinformatics preprocessing
We assessed genetic variation and gene expression levels among individuals by generating ddRAD-sequencing and RNA-sequencing (seq), respectively. Genomic DNA and RNA extraction of the 54 shrimp (three temperature treatments × three origins × six shrimp, combining three from each replicate tank) was carried out on ca. 2 g of muscle tissue from shrimp abdomen, using QIAGEN DNeasy Blood & Tissue kit and Qiazol RNA extraction reagent (Qiagen), respectively, following the manufacturer's protocol. Muscle was chosen to ensure enough material was available as input for sequencing, and because this tissue displayed metabolic changes in response to temperature changes in studies of other shrimp species, in addition to be responsible for general locomotion and escape-behavior responses in decapods (e.g., Rosa et al., 2014;Madeira et al., 2020). RNA-seq methodology was selected to investigate the changes in the entire transcriptome, with no a priori and less bias manner than selecting specific genes, in addition to displaying more power in revealing gene expression changes that could be difficult to discerned through more targeted approaches.
Preparation of ddRAD-seq libraries for the genomic DNA was carried out by the Institut de Biologie Intégrative et des Systèmes (IBIS at Université Laval, Quebec City, QC, Canada) with PstI and MspI restriction enzymes, and using the procedure described by Poland et al. (2012), with the following exceptions. First, a blue Pippin (SAGE sciences) was used to size libraries before PCR amplification (elution set between 50 and 65 min, on a 2% agarose gel); second, plate barcoding was used to enable sequencing on a shared Illumina NovaSeq 6,000 S4 PE150 lane, as described in Colston-Nepali et al. For the RNA samples, library construction was performed by Génome Québec, using NEB mRNA stranded library preparation kit. In addition to the 54 libraries (one per individual), a supplemental library was constructed from a pool of 120 shrimp exposed in different environmental conditions and from various populations (Supplementary Table S2), to maximize the number of detected genes during de novo transcriptome assembly. Finally, high-throughput sequencing steps for all libraries (Illumina NovaSeq 6,000 S4 PE100) were performed by Génome Québec. Pandalus borealis collection sites. Map of the southeast of Canada with the three geographic locations where shrimp used in this study originated from. Collection conditions are described in Supplementary Table S1.
Frontiers in Ecology and Evolution 04 frontiersin.org For the resulting DNA sequencing dataset, adapters were removed with Trimmomatic v0.39 (Bolger et al., 2014), and three base pairs were removed from the R2 (the MspI restriction site) to account for lower quality in the first nucleotides. Demultiplexing and quality filtering were performed with the process_radtags module of STACKS v2.4 (Catchen et al., 2013;Rochette et al., 2019), with a truncation at 135 pb and a check for the PstI restriction site. Cleaned reads were aligned to a P. borealis reference genome (newly assembled by C3G-McGill University, Montreal, QC, Canada) with mem from BWA-MEM with default parameters. Aligned paired-end reads were assembled with the gstacks modules and only samples with a mean depth of coverage of at least 5× were kept in the analysis. SNPs were first filtered with populations module of STACKS to only include SNPs present in all three origins, for at least 75% of individuals and with a minimum minor allele frequency of 5%. Furthermore, SNPs with more than 10% missingness and in Hardy-Weinberg disequilibrium in more than one origin were removed. Only the first SNP by loci was kept to avoid linkage disequilibrium. The ddRAD-seq final dataset consisted of 22,227 SNPs over 54 individuals.
Transcriptome assembly and annotation, and isoform count matrix estimation, were performed by C3G. De novo transcriptome assembly was performed following the pipeline described by Haas et al. (2013) and based on the Trinity assembly software suite. In brief, raw reads were trimmed using Trimmomatic from the 3′ end to have a minimum Phred score of 30 and a minimum length of 50 bp. Normalization was performed using the Trinity normalization utility v2.0.4 (Haas et al., 2013). An abundance estimate was computed using RNA-Seq by Expectation Maximization (RSEM) v1.2.12 (Li and Dewey, 2011) via Trinity align_and_estimate_abundance.pl utility. Expected count values from RSEM were used for downstream analysis. Transcripts were annotated with the Trinotate suite v2.0.2 (Haas et al., 2013), based on a homology search to known sequence data (BLAST+/SwissProt/Uniref90), protein domain identification (HMMER/PFAM), protein signal peptide and transmembrane domain prediction (signalP/tmHMM), comparison to currently curated annotation databases (EMBL UniProt eggNOG/GO Pathways databases). Each transcript was aligned against a custom protein database of the Pacific white shrimp Penaeus vannamei (Boone, 1931, GenBank: QCYY00000000.1). An annotation was assigned to the longest putative coding transcript based on the best BLAST hit with an associated gene name and an Expect (E) value cut-off <1e −5 whenever possible. All of the above described analyzes were performed using the GenPipes pipeline framework [PMID: 31185495].

Statistical analyzes
All statistical analyzes were carried out using R software, version 4.1.0 (R Core Team, 2020).

Mortality and molting rates
Mortality data used in this study originate from Chemel et al. (2020) and Guscelli et al. (submitted for review). We performed an univariate two-ways ANOVA test on arcsin square-root transformed percentage values (Sokal and Rohlf, 1995) to investigate how temperature affected shrimp mortality and molting rates per day, with Origin and Temperature as fixed factor tested in isolation and combined (Origin × Temperature). In addition, a Tukey's Honestly Significant Difference (HSD) test was used to conduct post-hoc analyzes, when significant effects were detected by the ANOVA.

Differential transcript expression analysis
Variation partitioning of gene expression and differential expression analyzes were performed using the R package vegan version 2.6-2 (Oksanen et al., 2022) and the Bioconductor's package DESeq2 version 1.36 (Love et al., 2014). Expected read count values from RSEM were used to identified differentially expressed transcripts (DETs). The count matrix was first pre-filtered to remove every transcripts not or sparsely expressed according to the criterion based on Chen et al. (2016). Briefly, a transcript was considered expressed if it had a count of at least 10-15 in at least some libraries, but basing the filtering on countper-million (CPM) values so as to avoid favoring genes that are expressed in larger libraries over those expressed in smaller libraries. We used a cutoff = 10/L where L is the minimum library size in millions. From the 118,609 identified transcripts, the pre-filtering process retained a total of 39,065 transcripts. To quantify the proportion of the total gene expression variation that was explained by the main factors or their interactions, we applied partial redundancy analyzes (RDA; Borcard et al., 1992), using the normalized count matrix (i.e., applying a variance stabilizing transformation (VST) to the count data, as implemented in DESeq2) as the response variable and the Origin and Temperature as categorical explanatory factors. Note that the six randomly chosen individuals were combining individuals from two tank replicates. Prior to this variation partitioning, we have tested for the tank effect on gene expression by performing a partial RDA adding tank identity as an explanatory variable.
We also identified any transcripts that showed changes in expression across the same three terms (i.e., Origin, Temperature and Origin × Temperature interaction), by building a general linear model, as implemented in DESeq2. Significance of each explanatory factor was assessed by applying a Likelihood Ratio Test (LRT) approach, comparing the full model to reduced ones (Love et al., 2014). Transcripts with FDR < 0.001 [value of p after Benjamini-Hochberg (BH) adjustment] were considered as differentially expressed. For all identified DETs, we assessed how expression levels varied across the three temperatures and origins by identifying modules of co-expressed transcripts. Transcripts whose expression profiles were highly correlated were grouped within a given module, using the WGCNA version 1.71 (Langfelder and Horvath, 2008). Enriched gene ontology (GO) terms were then identified, using Fisher's exact test as implemented in the topGO version 2.48 (Alexa and Rahnenfuhrer, 2022), to assess the functional role of DETs among origins or according to temperature treatment. For each GO term with FDR < 0.05, we calculated the proportion of identified co-expressed module.

Origin effect on gene expression
The origin effect on the measured gene expression could arise from different sources: 1-genetic differences among sampling sites, or 2environmental conditions in shrimp natural habitat (as individuals were sampled as adults, each of them might express specific genes programmed during their development). We attempted to clarify the role Frontiers in Ecology and Evolution 05 frontiersin.org of these two potential sources of shrimp gene expression by comparing the effect of population genetic variation, and environmental conditions at origins on transcript expression. Genetic variation was assessed using 22,227 SNPs called from ddRAD-seq dataset, and comparing the same 54 individuals as transcriptomic analysis. We assessed intra-population genetic diversity by calculating the expected heterozygosity (H E ) for each origin, using poppr version 2.9.3 (Kamvar et al., 2014), and genetic differentiation among collection sites with pairwise F ST , using dartR version 2.0.4 (Gruber et al., 2018). Environmental conditions at origins were characterized by monthly mean salinity and temperature for both surface and bottom sea water, obtained from Bedford Institute of Oceanography North Atlantic Model (BNAM) that averaged values over 1990 to 2015 period (Wang et al., 2018) and selected by removing highly correlated variables (|r| > 0.90; Supplementary Figure S2). To disentangle the different sources of gene expression that could explain the origin effect, we performed a variation partitioning (Peres-Neto et al., 2006;Legendre, 2008) of the transcriptomic response. This enabled us to assess the proportion of transcript expression explained by the terms Temperature, Environmental conditions at origin (hereafter Environment) and Genetic Differentiation (hereafter Genetics).
Factors explaining shrimp transcript expression were also determined by a model selection analysis using the function ordiR2step from vegan R package. We performed a model selection to obtain a parsimonious model maximizing the adjusted R 2 and p-value of a constrained ordination (Blanchet et al., 2008). We used the gene expression normalized count matrix as the response variable, and PCoA coordinates from pairwised F ST matrix (for genetic differentiation) and Environment as the explanatory variables, in addition to the interaction terms Temperature × Genetics and Temperature × Environment, to test for standing genetic variation for plasticity and effect of past environmental memory on gene expression, respectively.
Molting rates varied between 0.00 ± 0.00 (at 2°C in NNC) and 0.022 ± 0.001 (at 10°C in ESS). Mean molting rates differed with increasing temperatures for shrimp from the three origins, as indicated by the presence of an interaction between the terms Origin and Temperature (F 4, 9 = 16.412; p < 0.001; Figure 2B). Shrimp from NNC and SLE displayed no significant changes along the temperature gradient tested, but shrimp from NNC, and not those from SLE, displayed a quasi-absence of molt during the whole experiment ( Figure 2B). For ESS, molting rate increased with temperature ( Figure 2B).

Reference transcriptome
The reference transcriptome of P. borealis was built using RNA from pooled individuals in various physiological conditions and samples from the current study experiment (Supplementary Table S2), maximizing thereby the number of detected genes. We obtained a total of 2.81 × 10 9 raw paired reads, from which 97.48% remained after the trimming step. The de novo assembly resulted in 170,377 transcriptome isoforms that could be grouped into 118,609 components loosely representing genes. Transcript contig length ranged from 224 to 28,419 bp, with a N50 value of 1820 bp.
Functional annotations were performed on homology search to known sequence data based on the protein database of another shrimp species, P. vannamei. About 27.48% of the transcripts were successfully annotated to a specific function. Filtration of assembled transcriptome based on the functional annotation resulted in high quality contigs, composed of 46,828 transcript groups into 24,122 genes. Transcripts lengths from this filter assembled transcriptome ranged from 224 bp to 28,419 bp, with a N50 value of 3,049 bp.

Differential gene expression analysis revealed extensive transcriptomic plasticity
For shrimp collected from the temperature experiment, a total of 39,065 transcripts were expressed at a count greater than 10 reads in at least one sample. Variation partitioning analysis revealed marginal effects of Temperature (adj. R 2 = 0.87%; p = 0.015), Origin (adj. R 2 = 1.75%; p = 0.001) and their interaction (adj. R 2 = 1.34%; p = 0.006) on gene expression levels. Differential expression analysis detected a higher number of DETs across the temperatures tested (n = 1,607; 4.11%), than across origins (n = 907; 2.32%), and a very low number of significantly DETs were detected for the interaction term (n = 13; 0.03%). Ordination plot of distance-based RDA ( Figure 3A) showed that the gene expression level gradually changed with increasing temperature (db-RDA 1: 22.50%; p = 0.001). In addition, the gene expression level of shrimp from NNC differed markedly from those reported for shrimp from SLE and ESS (db-RDA 2: 20.11%; p = 0.001).
The patterns of expression profiles of DETs identified through co-expression module analysis differed between Temperature and Origin effects. A total of 1,583 (98.51%) and 845 (93.16%) DETs were grouped across temperatures and origins, respectively, into modules of co-expressed transcripts ( Figures 3B,C). We identified three and seven modules for DETs across temperatures and origins, respectively ( Figures 3B,C). For temperature, gene expression levels increased (module T_2) or decreased (modules T_1 and T_3) linearly with increasing temperature for all three modules identified (Supplementary Figure S1A). For origin, six out of the seven modules identified showed SLE and ESS displaying similar expression levels (modules O_2 and O_7), but differing from those of NNC ( Figure 3C).
GO enrichment analyzes revealed that DETs associated with temperature treatments or origins involved distinct functions. For instance, Temperature affected the expression of genes belonging to the three ontology categories (i.e., biological processes, cellular components and molecular functions), whereas Origin affected the expression of genes associated only to cellular components and molecular functions (Supplementary Figure S1A) Figure S1B) compared to those involved in catabolic activities (11.53%, s.d. 2.76%; Supplementary Figure S1B). In terms of cellular components, transcripts associated to protein, catalytic and channel complexes were mostly downregulated with increasing temperature whereas those involved in cell structures, such as actin and myosin, differed significantly among shrimp from different origins. Most cellular components associated with DETs across shrimp origins were part of module O_1 with the expression of these specific transcripts being higher in shrimp from ESS and progressively lower in shrimp from SLE and from NNC (Supplementary Figure S1B), except extracellular functions. These latter functions were represented by up to five modules characterized by similar expression levels for shrimp from ESS and SLE, which were different from NNC. Some DETs with extracellular functions were associated with cytoskeleton and non-membrane-bounded organelle. In terms of molecular functions, Temperature modified the expression of genes associated with nucleic and amino acid activities, whereas Origin involved DETs associated with structural activities. For temperature, transcripts associated with threonine related activities decreased with increasing temperature.

Genetics differences among origins
Shrimp displayed genetic differences among origins. Individuals from SLE had lower genetic diversity compared to shrimp from NNC and ESS (Tukey HSD, p < 0.001; Figure 4A) but differences in gene diversities were small (H E ranging from 0.215 to 0.223). No genetic differentiation was detected between NNC and SLE (F ST = 0, p = 0.940), whereas ESS was differentiated from the two other origins (F ST > 0.003, p < 0.001; Figure 4B).

Environmental conditions at origins
Environmental conditions at origins (hereafter Environment) were contrasted. We used only the January and July bottom temperatures since these two variables characterized female shrimp habitat and were not strongly correlated (r = 0.32), unlike other variables (Supplementary Figure S2 and see method section for more details about variables selection). The NCC origin had the most different environmental conditions (Euclidean distance with SLE and ESS > 4.04) compared to SLE and ESS, which were more similar (Euclidean distance between SLE and ESS = 2.93). The origin NNC was characterized by cold temperatures compared to ESS and SLE ( Figure 5).

Temperature treatment and environmental conditions at origin affected Pandalus borealis plasticity but not genetic differentiation
We did not detect any tank effect on gene expression (partial RDA analysis; adjusted R 2 = 0.047%; p = 0.429), and consequently tank identity was not taken into account in the subsequent analyzes. A variation partitioning analysis identified marginal effects of Temperature (adj. R 2 = 2.98%; p = 0.001) and Environment (adj. R 2 = 2.22%; p = 0.001) explaining transcript expression levels ( Figure 6). However, genetic differentiation among origins did not explain transcript expression levels ( Figure 6). In addition, model selection retained Temperature, Environment and the Temperature × Environment as terms explaining transcript expression (Table 1). The terms Genetics and the interaction between Temperature and Genetics were not selected.

Discussion
Understanding the molecular mechanisms underlying ectotherms' responses to elevated temperature is paramount to improve our ability to predict their vulnerability to the ongoing global change. Genomic tools can help to improve prediction of organisms' response to changing environments. In this study, we provides a first reference transcriptome of A B

FIGURE 2
Pandalus borealis mortality and molting rates as function of Temperature and Origin. (A) Mortality rate at each treatment temperature (green, orange and red for 2, 6 and 10°C, respectively) for shrimp from different origins (shape an important marine ectotherm species, Pandalus borealis. Using this new resource, we showed that changes in transcript levels are mainly explained by temperature exposure and environmental conditions at origin. The lack of genetics or genetics-by-temperature interaction effects on transcript expression levels suggested limited evolutionary potential in the species' transcriptomic plasticity in response to temperature changes.

Temperature drove transcriptomic plasticity
Our results showed that P. borealis displays a high level of plasticity in gene expression with changing temperatures. The observed transcriptomic response supports the hypothesis that this species possesses an elevated degree of plasticity to cope with the temperature variation it could encounter throughout its life-time.
The linear transcriptomic response we reported across the three temperatures suggests that the northern shrimp can cope with increasing temperatures up to 10°C, at least within 30 days of exposure. This result confirms the suggestion that the temperature levels used for the plasticity assay were within the northern shrimp range of thermal tolerance Guscelli et al. submitted for review). However, the greater mortality observed at 10°C compared to 2°C (10%) supports the idea that the highest temperature level used in this study is close to the Gene expression response to temperature exposure of P. borealis shrimp from three different origins. (A) Variation of gene expression levels. Distancebased redundancy analysis (db-RDA) plot performed on gene expression levels according to origins (shapes) and temperature treatments (colors). (B,C) Relative expression level for each co-expression modules for differentially expressed transcripts (DETs) among temperatures (T) and origins (O), respectively. Means with their standard errors are represented as dots and error bars, respectively. Numbers in parentheses indicate the number of transcripts belonging to each co-expression module.
Frontiers in Ecology and Evolution 08 frontiersin.org upper limit of thermal tolerance of P. borealis (Guscelli et al. submitted for review). Such incongruence in the responses to temperature between the mortality rates and the linearity of gene expression may be explained by the fact that only survivors could be analyzed for gene expression. Besides, other tissues than the abdominal muscle analyzed in this study, such as the gills and heart, may show thermal transcriptomic responses more closely associated to survival given their vital physiological functions (Buckley et al., 2006;Liu et al., 2013;Valenza-Troubat et al., 2022). The analysis of such tissues could extend our understanding of mechanisms involved in P. borealis responses to global changes. We also reported changes in biological functions of differentially expressed transcripts (DETs) with increasing temperatures in P. borealis. Lipid catabolic activities have been shown to decrease with increasing temperatures in polar fishes (Pörtner, 2002). Similarly, we also observed catabolic activities decreasing with increasing temperature. For P. borealis, we identified protein catabolism at low temperatures which may be specific to this species (Supplementary Figure S1). We also observed increasing expression of genes associated to metabolic activities with increasing temperature, which supports previous observations of increased metabolic rate of P. borealis at higher temperatures Arnberg et al., 2013;Dupont-Prinet et al., 2013;Guscelli et al. submitted for review).
Our study also identified a molecular mechanism possibly involved in cold tolerance in P. borealis. Transcripts associated with  Frontiers in Ecology and Evolution 09 frontiersin.org threonine related activities decreased with increasing temperature. In Atlantic salmon Salmo salar (Linnaeus, 1758), threonine-type proteins were shown to be involved in low temperature conditioning (Silva et al., 2020). These proteins were also observed in response to cold treatment in the Pacific white shrimp P. vannamei (Huang et al., 2017), where threonine-type proteins are downregulated in the warmest temperatures, as for P. borealis.

Variation in gene expression was associated with environmental factors, but not genetic variation
The origin of northern shrimp was associated with contrasting transcriptomic patterns, with shrimp from NNC showing more distinctive gene expression patterns when compared to the SLE and ESS origins. Such difference could be due to acclimatization (environment) or local adaptation (genetics) effects. We used a model selection analysis with temperature, environmental conditions at origin, and genetic differentiation to disentangle acclimatization and adaptation effects on transcript expression levels. Our results showed that experimental temperature and environmental conditions at origin are the two factors explaining differences in transcript levels.
Local adaptation can translate into differential gene expression patterns in distinct populations. In this study, we confirmed the existence of different populations in P. borealis in northwest Atlantic by detecting significant genetic differentiation among origins. Our results also showed that shrimp from SLE display the smallest genetic diversity compared to other origins. Despite such genetic differences, we failed to detect any genetic effect on gene expression, either qualitatively with genetic diversity or quantitatively with genetic differentiation (variation partitioning and model selection analyzes). Our results suggest that local adaptation do not underly differential gene expression for shrimp from these three origins. Yet, single or a few loci might be associated with local adaptation and contribute to the geographic variation in gene expression (Whitehead and Crawford, 2006). This hypothesis is now under investigation in a distinct study.
The similarity in transcriptomic patterns and environmental conditions at origins observed in this study may suggest acclimation of shrimp to temperature at origins. However, the interpretation of these results is not straight forward. The variable developmental status of female shrimp was also highlighted with difference in molting rates across origins in this. The SLE shrimp were tested in early fall, when molting occurs in the field, whereas the ESS shrimp were tested in spring, close to hatching and before molting. In contrast, NNC shrimps were tested in winter and adult females are not molting during this season (DFO, 2012(DFO, , 2021Bourdages et al., 2022). Consequently, the effect of environmental conditions at origins also included difference in developmental status of female shrimp. Note that differences in transcript levels among shrimp from distinct origins also mirror those we report for molting rates. In addition, most gene functions involved in DETs across origins were associated with cell structures. Transcripts involved in actin and myosin production, cytoskeletal motor activity, and cuticle structural components were also upregulated in SLE and ESS. Such gene functions were previously associated to the molting process in shrimp (de Oliveira Cesar et al., 2006;Corteel et al., 2012;Zhang et al., 2019).

Limited potential of adaptation to global changes of northern shrimp
Ecologically divergent populations can evolve different transcriptomic plasticity, as observed in the redband trout Oncorhynchus mykiss gairdneri (Richardson, 1836) from desert and montane climates (Narum and Campbell, 2015), and in the threespine stickleback Gasterosteus aculeatus (Linnaeus, 1758) from marine and freshwater environments (Morris et al., 2014). However, we did not observed different thermal transcriptomic platicity in P. borealis among genetically different populations. This result indicates that this species did not evolve locally adapted transcriptomic plasticity in response to the temperature range tested. Furthermore, this lack of Factors influencing transcript expression in P. borealis. Venn diagram representing the proportion (adj. R 2 ) of gene expression level explained by exposure temperature (Temperature), environmental conditions at origin (Environment) and population genetic differentiation (Genetics, i.e., pairwise F ST between origins). Number outside circles refer to total percentage of variation explained by each underlined factor. Numbers inside circles indicate the marginal contribution to transcript expression by each factor. Percentages within intersections indicate shared contribution across different variables. Significance of total and marginal contribution to gene expression was tested by permutations test, using 999 randomizations. genetic variation underlying transcriptomic responses to changing temperature suggests a limited evolutionary potential of P. borealis thermal plasticity, at least at the transcriptional level. It is unclear to what extent this limited evolutionary potential of thermal plasticity at the transcriptomic level reflects a limited evolutionary potential for the plasticity of integrative traits such as growth and reproduction. However, plastic response of gene expression could be adaptive (Koch and Guillaume, 2020), and the evolution of plasticity for the transcriptional phenotypes could propagate across hierarchical levels of the organisms (from gene expression to integrative phenotypes), as observed in another organism (Leung et al., 2022). Our results would also suggest a limited evolutionary potential of integrative traits plasticity.
In conclusion, understanding the molecular mechanisms underpinning a species' thermal plasticity and its vulnerability to global change drivers can ultimately help developing appropriate conservation and management plans, this being paramount also for species of commercial importance that supports the socio-economic development of coastal communities, including first nation ones. Our results show that individuals of P. borealis surviving 30 days exposure to elevated temperature are able to buffer the immediate impact of ongoing ocean warming by altering their gene expression. However, the lack of standing variation underlying phenotypic plasticity in response to temperature suggests a limited potential for the evolution of its plasticity. Our results emphasise its vulnerability under future global changes scenarios, as observed in other ectotherm organisms (Gunderson and Stillman, 2015;Sirovy et al., 2021). Consequently, the higher mortality reported at the highest temperature in this study also supports that the thermal niche of P. borealis will be progressively restricted to higher latitudes with the ongoing ocean warming (Hobday et al., 2016;Kleisner et al., 2017;Morley et al., 2018;McHenry et al., 2019). This prediction concurs with current decreases in abundance already observed at lower latitudes for this species, prediction also based on modeling (Barria et al. submitted for review) and other cold-adapted crustaceans (Neumann et al., 2013;Lenoir et al., 2020;Melo-Merino et al., 2020;Simões et al., 2021).

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA909194 and https://github.com/GenomicsMLI-DFO/ PANOMICS_Transcripto_manuscript.

Ethics statement
The Canadian Council for Animal Care does not require that projects using invertebrates other than cephalopods be approved by an Animal Care Committee.

Author contributions
GP obtained the funding for the genomic analyzes. DC and PC obtained the funding for the experimental work. CL and GP conceived the tested hypothesis. EG, DC, and PC conceived, designed, and ran experiment. AB performed the bioinformatic analyzes of ddRADseq data. CL analyzed the transcriptomic data, prepared figures and tables, and wrote the original draft. All authors contributed to the article and approved the submitted version.

Funding
The experimental work was funded by: (i) an OURANOS grant (554023)  The genomic work was funded by: Genomics R&D Initiative (GRDI), Government of Canada.

Conflict of interest
As PC was a Topic Editor of this special issue, he did not participate in any form to editorial decision nor the editorial processing of this article.
The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.