Enjoying the warming Mediterranean: Transcriptomic responses to temperature changes of a thermophilous keystone species in benthic communities

Information about the genomic processes underlying responses to temperature changes is still limited in non‐model marine invertebrates. In this sense, transcriptomic analyses can help to identify genes potentially related to thermal responses. We here investigated, via RNA‐seq, whole‐transcriptomic responses to increased and decreased temperatures in a thermophilous keystone sea urchin, Arbacia lixula, whose populations are increasing in the Mediterranean. This species is a key driver of benthic communities’ structure due to its grazing activity. We found a strong response to experimentally induced cold temperature (7°C), with 1,181 differentially expressed transcripts relative to the control condition (13°C), compared to only 179 in the warm (22°C) treatment. A total of 84 (cold treatment) and three (warm treatment) gene ontology terms were linked to the differentially expressed transcripts. At 7°C the expression of genes encoding different heat shock proteins (HSPs) was upregulated, together with apoptotic suppressor genes (e.g., Bcl2), genes involved in the infection response and/or pathogen‐recognition (e.g., echinoidin) and ATP‐associated genes, while protein biosynthesis and DNA replication pathways were downregulated. At 22°C neither HSPs induction nor activation of the previously mentioned pathways were detected, with the exception of some apoptotic‐related activities that were upregulated. Our results suggest a strong transcriptional response associated with low temperatures, and support the idea of low water temperature being a major limitation for A. lixula expansion across deep Mediterranean and northern Atlantic waters.

, to more specific effects. On the other hand, mean temperature increases and heat waves also result in a number of lethal and sublethal effects on particular species and/ or populations, including coral reef bleaching (Hughes et al., 2017), alteration of animal migration behaviour, and shifts of marine taxa distribution patterns, among many others (e.g., Deutsch, Ferrel, Seibel, Pörtner, & Huey, 2015;Hoegh-Guldberg & Bruno, 2010;King et al., 2018).
Although the function of HSPs is well-known in some marine invertebrates, additional knowledge on the involvement of other molecular pathways, such as antioxidant genes, apoptosis-associated and immune-associated genes, is needed to uncover other relevant mechanisms involved in thermal stress responses in ecologically relevant species (Gleason & Burton, 2015;Zhu et al., 2016). One relatively recent approach to investigate rapid organismal responses to environmental perturbations, to identify potential physiological networks, and to discover candidate genes and isoform variants involved in their responses, is to explore the whole transcriptional profiles using RNA-seq techniques (e.g., Evans, Pespeni, Hofmann, Palumbi, & Sanford, 2017;Xu et al., 2018;Zhu et al., 2016). Although the relationship between mRNA transcript abundance and protein abundance is still not clear (Feder & Walser, 2005), some studies have shown a correlation between these two variables (Maier, Güell, & Serrano, 2009). Changes in gene expression are considered to be sensitive indicators of stress and potential predictors of organismal physiology under experimental conditions (Buckley, Gracey, & Somero, 2006;Feder & Walser, 2005;Schoville, Barreto, Moy, Wolff, & Burton, 2012).
The black sea urchin Arbacia lixula (Linnaeus 1758) has tropical affinities (Tortonese, 1965) and an amphi-Atlantic distribution across shallow rocky ecosystems, being the Moroccan coast its northern-most distribution limit in the east Atlantic. This sea urchin entered the Mediterranean basin during the last Pleistocene interglacial period (Pérez-Portela et al., 2019;Wangensteen et al., 2012), and it is now a common species across the whole Mediterranean (Palacín, Turon, Ballesteros, Giribet, & López, 1998;Tortonese, 1965). Densities of this species significantly increased in some Mediterranean areas during the recent decades (Francour et al., 1994;Harmelin et al., 1995;Hereu et al., 2012), and it is among the key drivers structuring littoral communities due to its grazing activity (Bonaviri, Fernández, Fanelli, Badalamenti, & Gianguzza, 2011). The species is capable of shifting complex littoral macroalgal beds into "barren grounds"-areas of high densities of sea urchins deprived of erect seaweeds and dominated by crustose coralline algae - Gianguzza et al., 2011).
Several authors have predicted that the foreseen global warming might have a positive effect on its reproduction output and larval survival (Francour et al., 1994;Gianguzza et al., 2014;Visconti et al., 2017;Wangensteen, Dupont, Casties, Turon, & Palacín, 2013;Wangensteen, Turon, Casso, & Palacín, 2013). This potential effect, if real, will represent a worrisome increase of the impact of this sea urchin on littoral ecosystems in a near future Wangensteen, Dupont, et al., 2013;. On the other hand, it seems that the distribution of A. lixula is constrained by low temperatures, like the low sea surface temperature provoked by the southward Portugal Current (Martins, Hamann, & Fiúza, 2002), which might be the cause of its absence along the Atlantic coast of Europe (Wangensteen et al., 2012). In this sense, experiments to investigate the potential of A. lixula to invade deep waters, analysing the combined effect of pressure (from 1 to 250 atm) and temperature (from 5 to 15°C) on the survival of embryos and larvae, showed that the combination of high temperatures and pressures, rather than temperature per se, might be the major factor limiting the distribution of the species at depth (Young, Tyler, & Fenaux, 1997). In contrast, more recent studies have demonstrated higher mortality rates, larval growth abnormalities and significant delays in settlement at the lowest experimental temperatures tested on this species (experimental temperatures from 18 to 22°C: Privitera, Noli, Falugi, & Chiantore, 2011;and from 16 to 19°C: Wangensteen, Dupont, et al., 2013). According to these studies, the abundance of A. lixula in the Mediterranean might be constrained by the low winter temperature of colder years, when mean temperatures can drop to 11°C, because gonad maturation is then considerably impaired (Lejeusne et al., 2010;Wangensteen, Dupont, et al., 2013). However, whereas the mentioned studies provided insights on the effects of thermal variation on the early development stages of A. lixula, almost nothing is known about its effects on the general performance of adult individuals, which can have different thermal sensitivity (Buckley & Huey, 2016). The capability of adult individuals to acclimatize and endure thermal changes is highly relevant from an evolutionary perspective. It not only affects their own physiological performance and/or the quality of their gametes, but it can also result in negative transgenerational carry-over effects on hatchability and larval size of the next generation, which have been shown after prolonged periods of parental exposure to elevated temperatures in some sea urchins (Zhao et al., 2018). In sea urchins, transcriptomes from different tissue types and larval thermal stress responses have been characterized (e.g., Clark et al., 2019;Gaitán-Espitia, Sánchez, Bruning, & Cárdenas, 2016;Gillard, Garama, & Brown, 2014;Jia et al., 2017;Pérez-Portela, Turon, & Riesgo, 2016;Runcie et al., 2012). But, to our knowledge, transcriptome-wide screenings have never been used for measuring responses to thermal variation in adult individuals of this animal group.
The aim of this study was to explore the short-term transcriptional response to thermal changes among individuals of the subtropical sea urchin A. lixula. We set three specific objectives for our study: (a) To quantify and compare transcriptional responses to both high and low temperature treatments in A. lixula under experimental conditions; (b) To identify some of the most important candidate genes involved in rapid thermal responses in sea urchins; and (c) To determine the existence of common genes involved in responses to increasing and decreasing temperatures.
Many studies on global warming focus on the negative effect of rising temperatures, but in this study, we worked under the hypothesis that A. lixula will experience higher stress when subjected to low rather than to high temperatures. Based on previous transcriptional information from marine invertebrates under thermal stress (e.g., Gleason & Burton, 2015;Zhu et al., 2016), we also expect changes of expression patterns in different gene pathways during our temperature treatments, including genes encoding HSPs, apoptosis and antiapoptosis mechanisms, ATP-associated genes due to an increase of energy demand to maintain cell homeostasis, antioxidant genes since extreme temperatures can increase cells' oxidative stress, and immune-associated genes (Xu et al., 2018). The information obtained here will be relevant to understand the ecophysiological patterns of sea urchins exposed to thermal changes. We also discuss the significance of our findings for the foreseeable ecological spread of this keystone species in the Mediterranean.

| Experimental design
To quantify rapid transcriptomic responses of A. lixula under thermal assays, we exposed adult sea urchins (test diameter 40-50 mm) to three different treatments under controlled laboratory conditions for 20 hr: control (CT) with sea water at 13 ± 1°C, sea water temperature at 7 ± 0.5°C (T7), and sea water temperature at 22 ± 0.5°C (T22). We set the temperature exposure time to 20 hr because previous experiments of thermal stress responses in other marine invertebrates demonstrated maximum peaks of expression between the first 6-24 hr, depending on the genes (e.g., Kim et al., 2017;Zhu et al., 2016).
It is important to note that our goal was to submit the test organisms to an acute thermal change to measure their responses, not to mimic highest or lowest seasonal temperatures in the area. The treatment temperatures were chosen to represent an important shift with respect to the controls (13°C, the surface water temperature at this location when sea urchins were collected in wintertime) while remaining within realistic values for our area of study, the NW Mediterranean. Thermal sensitivity and resistance of organisms are not constant over time and often shift in response to seasonal conditions (Buckley & Huey, 2016 (Marbà et al., 2015;Pastor, 2012). The global average for the coldest month of the year (February) in the Mediterranean is 14.5°C, with a lower average value (12°-13°C) found at the NW Mediterranean (Pastor, 2012) (see Appendix S1, Figure S1). Since the species' thermal history can determine the thresholds of stress response (Osovitz & Hofmann, 2005) and thermal sensitivity can change over the seasons, we made a preliminary assessment of the tolerance limits of our NW Mediterranean population at that time of the year (so-called here "trials"), with several temperatures assayed over a 20 hr period and visual inspection of the state and activity level of 10 sea urchins per temperature treatment. Specimens used for the trials were not used for further experiments and were returned to the sea after experimentation, nor were samples collected for transcriptomic analysis during the trials. For the trials, we used 22°, 24°, and 26°C as upper thermal limits, and 12°, 9° and 7°C as lower limits. 7° and 22°C marked the lower and upper thresholds, respectively, at which all individuals used for the trials remained alive, visually healthy (intact skin, no algae or microorganism colonies growing up over the animal surface and no massive spine lost) and active (feet and spines movement). For the cold treatment, 7°C (a decrease of 6°C relative to the control) was the limit temperature achievable in winter in shallow embayments in the NW Mediterranean (e.g., Ordóñez et al., 2015), while for the warm treatment we increased temperature by 9°C (relative to the control), being 22°-23°C the conditions encountered in mid-summer in the study area (e.g., De Caralt, González, Turon, & Uriz, 2018;Marbà et al., 2015;Pastor, 2012). Over 22°C, experimental animals either died or presented clear signs of infection with microorganism colonies over the skin and/or massive loss of spines. We emphasize that, while sea urchins thrive at this temperature and higher in summer, we were performing an acute exposure treatment during wintertime, so we had to adjust our treatments accordingly.
Our experimental design for transcriptomic analysis consisted of two different experiments: A "Low temperature" experiment comparing the control condition at 13° ± 1°C and experimental condition at 7° ± 0.5°C, hereafter named as "Control versus T7", and a "High temperature" experiment comparing the control condition at 13° ± 1°C and experimental condition at 22° ± 0.5°C, hereafter named as "Control versus T22" (see Figure 1). Samples used as control condition were the same for both experiments, since all treatments were run at the same time and laboratory. After the acclimation period of 48 hr, each sea urchin was placed in an independent aquarium to avoid interactions among specimens. Each aquarium had constant airflow and the seawater temperature was set at the required temperature (13°, 7° or 22°C) prior to adding the sea urchins. Temperature of the aquaria was controlled with HOBO loggers (one per aquarium).
Aquaria with different treatments were randomly allocated across the wet-lab space to avoid any bias related to their spatial distribution. Spain embedded in paraffin, cut in 5 μm sections using a Microm HM325 Microtome, and stained in haematoxylin-eosin as described in  and Garcia-Cisneros et al. (2018).
Sex was then determined under the optical microscope.

| Coelomocytes collection and RNA sequencing
Coelomocytes consist of several cell types contained in the coelomic fluid and are immune effectors in echinoderms (Matranga et al., 2000;Smith et al., 2018). They have been used as biomarkers of stress due to their prompt response to changing environmental conditions (Matranga, Bonaventura, & Di, 2002;Matranga et al., 2000Matranga et al., , 2005Pinsino et al., 2008) that can reduce the protective capacity of these cells and rapidly induce activation of the heat shock proteins expression (Matranga et al., 2000;Pinsino et al., 2008). Additionally, these cells showed higher thermal response capacity than other tis-

| Sequence processing and de novo assembly
The software FASTQC v. 0.10.0 (www.bioin forma tics.babra ham. ac.uk) was used to visualize and measure the quality of the raw reads generated in the HiSeq2000. Adapters and bases with low quality (phred scores < 33) were trimmed off, and a length filter was applied to keep only sequences of > 25 bases using TrimGalore v. 0.2.6 (www.bioin forma tics.babra ham.ac.uk). High-quality reads were rescreened in FASTQC to ensure a good quality of the samples after trimming. A basic scheme of the most important steps of our pipeline is presented in Figure 2.

| Differential expression analyses and annotation
Reads from all replicates in each experiment were aligned against the corresponding "reference" transcriptome as per experiment (see Figure 2). Paired reads after trimming were mapped using Bowtie2 v. 2.2.1 (Langmead & Salzberg, 2012)

as implemented in
Trinity (Grabherr et al., 2011). RSEM v. 1.2.11 (Li & Dewey, 2011) was then run to generate a table with read counts, and unmapped reads were discarded. In the "reference" transcriptomes, transcripts of the same trinity component were treated as different isoforms.
We retained information of differential expression of all isoforms detected for a given gene ( were observed between males and females (p-adjusted > .01), and "sex" was not considered as a variable in further analyses.
For differential gene expression analyses, read counts were first normalized in DESeq2, and then a negative binomial model was fit to accurately estimate differential expression. The significance value for multiple comparisons was adjusted to 0.01 with the function "padj" (Benjamini-Hochberg adjustment) as implemented in DESeq2.
Transcripts with significantly different expression values relative to the controls will be hereafter called "DE" transcripts. Principal com-  Using the GO annotation results from the de novo assemblies of the two experiments, we obtained the GO terms associated to the differentially expressed isoforms, which were then input (together with their associated log 2 fold-change) to the REVIGO web server (Supek, Bošnjak, Škunca, & Šmuc, 2011) to obtain summaries of GO terms. Results were graphically represented with the "treemap" R package. Size of the rectangles was adjusted to reflect the log-2 fold-change in REVIGO. Differentially expressed isoforms without blast hit, unknown function and/or without annotation for each experiment were further assessed with the InterProScan 5 software

| Data filtering and de novo assembly
A total of 18 RNA-seq data sets were generated in this study. For de novo assembling of each reference (see Figure 2) we used 12 RNAseq data sets. For quantifying transcriptomic responses, we used 11 data sets per experiment since one sample from treatment T7 and another from T22 were discarded for gene expression analyses because of their low quality (see Figure 1). Data sets have been deposited in Mendeley Data (https://doi.org/10.17632/ 5673n 552yj.1) and the NCBI (BioProject n° PRJNA642021). The number of trimmed reads used for de novo assembly, as per sample replicate and treatment, are detailed in Appendix S1 (Tables S1 and S2). All replicates had over 26 million reads.
The de novo assembly "CT + T7", used as a reference for the "Control versus T7" experiment, included 141.5 Megabases that rendered 211,650 transcripts (including both genes and their different isoforms), and 19.6% of them had blast hit with known proteins of metazoans (see species blast hit distribution in Appendix S1, Figure   S2). The reference assembly "CT + T22" for the "Control versus T22" experiment included 147.4 Megabases, and rendered 219,655 different transcripts, from which 17.9% had blast hit (see species blast hit distribution in Appendix S1, Figure S2). Both de novo assemblies were very comparable (and had 99.5% transcripts in common), presenting relatively high N 50 values, between 1,102 and 1,114, meaning that over 50% of the transcripts were longer than 1,100 bases. Details of the de novo assemblies for the two different experiments are presented in Appendix S1 (Table S1). Both, "CT + T7" and "CT + T22", showed high completeness when compared with BUSCO conserved orthologue databases of eukaryotes and metazoans (see Appendix S1, Table S3). For the reference assemblies, "CT + T7" and "CT + T22", 194 and 4,293 transcripts, respectively, had blast hits against proteins and/ or nucleotides of bacteria and were removed from subsequent analyses. In fact, most differences between the reference assemblies "CT + T7" and "CT + T22" were due to the amount of bacterial transcripts.

| General results of differential expression analyses
The differential expression analyses revealed changes in gene expression between controls and temperature treatments in both experiments, "Control versus T7" and "Control versus T22". Additionally, we detected a remarkable difference in the magnitude of the transcriptomic responses between experiments, which was over six-fold greater in number of DE transcripts in the "Control versus T7" experiment, as explained below. We also observed differences in gene expression among different isoforms of the same genes.
In the "Control versus T7" experiment, we detected  Figure 2). Only 35 transcripts (19.7% of the total DE transcripts) were assigned annotation and known function (Appendix S1,   A total of 84 and three GO terms were found associated to DE genes in the "Control versus T7" and "Control versus T22" experiments, respectively (Appendix S1, Table S4

| D ISCUSS I ON
The response of marine organisms to thermal shifts is probably different across the species' range of distribution (Donelson et al., 2019). In our study, we investigated transcriptional responses of a keystone species, the black sea urchin, in the northern part of its range of distribution (NW Mediterranean). We found contrasting responses to low (7°C) and high (22°C) temperatures, with the former eliciting a much stronger reaction. Such differences were related to both the magnitude of the transcriptional response (e.g., number of up-and downregulated transcripts and gene expression fold-change) and the diversity of genes and pathways involved in these responses.
The capacity of ectotherm species to thrive across wide temperature ranges is, in part, based on their ability to modulate the expression of genes encoding proteins involved in the physiological, metabolic and cellular stress responses (Kim et al., 2017;Runcie et al., 2012;Stillman, 2003;Tomanek, 2010). Resistance to acute sublethal temperatures is an adaptive trait that varies among species of the same genus from different latitudes and habitats (Stillman, 2003;Yao & Somero, 2012). In general, marine tropical species are more heat tolerant than their temperate and cold counterparts (Somero, 2010). Paradoxically, analyses of both marine and terrestrial ectotherms suggest that tropical, or the warm-adapted species, may be more threatened by global warming because they live closer to their upper physiological thermal limit, and have higher metabolic rates that accelerate quicker than in colder species under rising thermal conditions (e.g., Somero, 2010;Stillman, 2003).
According to this expectation, A. lixula, a heat tolerant species with subtropical affinities (Tortonese, 1965;Wangensteen et al., 2012), could be threatened by global warming across the warmest areas of its geographical distribution (Elmasry, Razek, El-Sayed, Omar, & El Sayed, 2015;Rilov, 2016), where it might be closer to its thermal physiological limits. However, in the northwestern Mediterranean this species is in the coldest part of its range of distribution, which encompasses both sides of the tropical and subtropical Atlantic (Wangensteen et al., 2012), and thus it could be more limited by cold temperatures. Current Mediterranean sea warming may be removing thermal limitations for this species (Francour et al., 1994;Gianguzza et al., 2014;Visconti et al., 2017;Wangensteen, Dupont, et al., 2013; allowing an increase in its abundance in the Mediterranean. In general, it is difficult to determine whether changes of expression in particular genes have important functional consequences, because for each gene the threshold for metabolic and physiological downstream effects can be different, and relatively small changes in gene expression of only a few genes can be as functionally important F I G U R E 6 Gene ontology treemaps for annotated differentially expressed genes in Control versus temperature 7°C. Down-and upregulated categories at 7°C are presented as separated figures for "Biological processes" and "Cellular components". Molecular functions are presented in Appendix S1 ( Figure S4). The size of the rectangles reflects the log 2 fold-change associated to the differentially regulated categories. [1, multicellular organism development; 2, positive regulation of GTPase activity; 3, regulation of ARF protein signal transduction; 4, single-organism cellular process; 5, DNA replication; 6, polyamine biosynthetic process; 7, histone H3 deacetylation; 8, DNA-templated transcription initiation; 9, microtubule-based movement; 10, cation transport; 11, transmembrane transport; 12, ER to Golgi vesicle-mediated transport;  to Hsp70s (Ngosuwan, Wang, Fung, & Chirico, 2003). In general, these HSP chaperones are involved in the strong and mild thermal stress responses, including protein-folding reactions to avoid protein denaturation of adults, early development stages and eggs of sea urchins (e.g., González et al., 2016;Matranga et al., 2002;Matranga et al., 2000;Runcie et al., 2012). Their presence might be involved in the wide thermal distribution of some particular marine species (see Zhu et al., 2016, and references herein), and the HSP family seems to be a mechanism to cope with the stress associated with cold in A. lixula (e.g., NW Mediterranean). On the other hand, no upregulation of genes encoding HSPs was detected at 22°C relative to the control condition in A. lixula.
Under conditions of thermal stress, protein refolding by HSPs may not be efficient enough, and misfolded protein degradation can be necessary to restore cell homeostasis (Mosser et al., 2000).
Therefore, other mechanisms such as proteolysis to eliminate disfunctional proteins via the Ubiquitin proteosome pathway and, finally, apoptosis to eliminate damaged cells, can be activated (Logan & Somero, 2011;Somero, 2010;Zhu et al., 2016). We only detected signs of Ubiquitin proteosome pathway activation (upregulation of some genes involved in this pathway) in the 7°C treatment, with the upregulation of the gene sequestosome 1 (Appendix S2), which is an autophagosome cargo that detects proteins for autophagy and has been previously identified in echinoderms (Bitto et al., 2014), and the e3 ubiquitin-protein ligase, which targets damaged proteins for transport and degradation by the proteasome (Ardley & Robinson, 2005).
In addition, we observed differential expression of several apoptosis-associated genes in both treatments, 7 and 22°C. Several studies demonstrated that sea urchins hold a complex apoptotic system (Agnello & Roccheri, 2010;Lesser, 2012). We found transcriptional changes at 7°C in apoptosis suppressor genes such as the Bcl2 (upregulated, Appendix S2), widely distributed in different marine invertebrates (see Lesser, 2012), and in genes containing death domains (downregulated: fas-associating death domain-containing protein and death ligand signal enhancer, Appendix S2) that induce cell apoptosis through the regulation of caspase activation (Agnello & Roccheri, 2010;Zhu et al., 2016). These findings suggest the activation of some particular pathways to control the programmed cell death at low temperatures. The upregulation at 22°C of the gene immediate early response 3-interacting protein 1-like (Appendix S2), which is a molecule involved in protein transport between the sarcoplasmic reticulum and Golgi apparatus and that mediates apoptosis in human cells (https://www.unipr ot.org), suggests that apoptosis is also occurring as a response to increased experimental temperatures.
Additionally, a Serine threonine-protein kinase pim3, an enzyme involved in the regulation of cell transport and survival, which prevents apoptosis by inducing the release of the anti-apoptotic Bcl2 mentioned before (Cross et al., 2000) was also overexpressed at 7°C, whereas a Serine threonine-protein phosphatase 6, with opposite function to the kinase enzyme (Cross et al., 2000), was downregulated at 22°C. Another interesting finding is the opposite pattern of gene expression found between experiments for the Wsc domain-containing protein 1 (downregulated at 7°C and upregulated at 22°C) (Appendix S2). Different members of the Wsc family are identified as putative receptors of stress and required for the heat shock response and the maintenance of cell wall integrity in yeasts (Lodder, Lee, & Ballester, 1999). The Wsc members are upstream regulators of other serine-threonine kinases, the protein kinase C1 (PKC1) and mitogen-activated protein kinase (MAPK), which can promote apoptosis (Cross et al., 2000;Lodder et al., 1999 Previous experiments on echinoderms demonstrated the effect of thermal stress on the immune capacity of coelomocytes; this effect is greater at higher rather than lower temperatures in the sea cucumber Apostichopus japonicas (Wang, Yang, Gao, & Liu, 2008).
However, in A. lixula, it was the lowest temperature the one that triggered a higher immune response in terms of gene expression. The echinoidin, senescence associated-gene, cytohesin-like and tripartite motif-containing protein 3 (TRIM) (Appendix S2) involved in the infection response and/or pathogen-recognition process against bacteria, fungi and viruses (Ozato, Shin, Chang, & Morse, 2008;Smith et al., 2006) were upregulated at 7°C. In addition, the gene interleukin-17 (Appendix S2) which is a cytokine inducing and mediates pro-inflammatory responses in metazoans and stimulates phagocytosis in echinoderms (Beck et al., 1993), was also upregulated at 7°C. None of these immune genes were, however, upregulated (or, when detected, were downregulated) at the highest experimental temperature (e.g., TRIM), suggesting no immune response at 22°C.
The differentially expressed genes for the low and high temperature experiment were associated to different GO categories that provide additional information. These GO categories summarize the most significant biological processes, cellular components, and molecular functions that were up-and downregulated during the experimental response in A. lixula. For the high temperature experiment, we could only recover GO terms of three transcripts, and therefore, there is limited information to reach conclusions on the GO categories for this experiment. However, we detected the down-regulation of two interesting GO terms, the "notch signaling pathway" with the associated gene neurogenic locus notch (Notch1), and the "integral component of membrane" with the associated gene encoding a notch ligand, the delta protein. Notch is a calcium-dependent cell signalling system involved in different functions, including cell differentiation, proliferation and apoptosis. In general, notch inhibits apoptosis and induces cell proliferation but, in vitro studies using different cell lineages, it was shown that hyperthermia reduced Notch1 expression and apoptosis in some cell lineages, whereas an opposite pattern was obtained in other cell lineages (Basile, Biziato, Sherbet, Comi, & Cajone, 2008). Although the effect of the notch downregulation at high temperatures in coelomocytes is not completely clear, it suggests the existence of an alternative pathway of apoptosis under thermal stress.
Among the GO terms upregulated during cold exposure that add further information was the "tyrosine metabolism" term, which is related to cell protection against stress, including the upregulation of HSPs, cytoskeletal stabilization and apoptosis decrease (Baird, Niederlechner, Beck, Kallweit, & Wischmeyer, 2013). This major GO term also includes the subordinate "Positive regulation of apoptotic process", which can induce apoptosis when protein refolding by HSPs is not efficient enough. The induction of HSPs during thermal stress can considerably increase the energy demand in cells (Dong, Yu, Wang, & Dong, 2011;Tomanek, 2010). This increased energy demand is reflected in the overrepresentation of the GO category "ATP hydrolysis", a catabolic process that releases energy previously stored in the form of ATP, and the upregulation of the V-type proton ATPase gene (see Appendix S2), a proton pump found within the "proton-transporting V-type ATPase, V0 domain" term.
Likewise, the terms "protein folding" and "protein transport", the last one subordinate to the "ATP hydrolysis" category, are linked to protein transport to the Sarcoplasmic reticulum for folding reaction to avoid protein denaturalization by HSPs. Hence, the "Sarcoplasmic Reticulum" category, corresponding to a key organelle involved in the thermal stress response that ensures that misfolded proteins are directed towards a degradative pathway to the central cytoplasmic proteolytic machinery (Malhotra & Kaufman, 2007), was also overrepresented at 7°C. Actually, the induction of expression of Hsp70s has been directly associated to the accumulation of unfolded proteins in the sarcoplasmic reticulum (Rachel, Tyson, & Stirling, 1997;Rao et al., 2002), which are later eliminated if refolding fails by retrograde transport across the reticulum membrane (Kostova & Wolf, 2003). Other minor upregulated GO terms, at the biological process and molecular function, were "oxidation-reduction processes" (1 GO term) and "oxidoreductase activitiy" (2 GO terms).
These terms suggest that low temperature affects the intracellular redox state in coelomocytes.
Among the downregulated GO terms at the 7°C treatment we found "Neurotransmitter transport" with the associated differentially expressed genes Creatine transporter and Trafficking protein particle complex subunit 2 protein. The Creatine transporter is essential for normal brain function in humans and tissues with high energy demands because, together with other molecules, maintains ATP levels (Christie, 2007). The downregulation of these genes and pathways could be a potential response to energy competition with the induction of HSPs during thermal stress. The 7°C treatment also seemed to inhibit nuclear replication, as represented by the downregulation of the "Nuclear origin of replication recognition complex" and "DNA replication" terms, among others. The origin recognition complex is an ATP-dependent system that, together with other factors, enables the initiation of DNA replication in eukaryotic cells (Li & Stillman, 2012). Cells under stressful conditions must prevent cell division in favour of protective functions (Jonas, Liu, Chien, & Laub, 2013), as well as to avoid entering in a new DNA replication cycle if there is DNA damage (Lee et al., 2009). We also found downregulation of the "intracellular signal transduction" term, with the subordinate "cell redox homeostasis" and "smoothened signalling pathway" terms, and the "protein O-linked glycosylation" term. Smoothened is a key transmembrane protein involved in a critical cell-to-cell communication system for tissue homeostasis. Glycosylation, on the other hand, is one of the most common post-transcriptional modifications during protein biosynthesis, which contributes to increase protein solubility and stability against proteolysis, and can also be involved in their correct folding (Shental-Bechor & Levy, 2008). Hence, the downregulation of these last two terms may reflect a negative effect of low temperatures on protein biosynthesis and stabilization, and homeostasis control in coelomocyte cells.
Another remarkable result from our work is the large number of gene isoforms found in A. lixula transcriptome. Different gene isoforms are mostly generated by alternative splicing, whose function is to increase the diversity of mRNA expressed by the genome (e.g., Kelemen et al., 2013;Stamm et al., 2005). It has been demonstrated that alternative splicing can promote from neutral or subtle effects on transcripts, and finally proteins functioning, to drastic physiological changes (see a review in Kelemen et al., 2013). Therefore, the presence of different gene isoforms is relevant when studying gene expression and physiological pathways. Although our objective was not to analyse differential splicing in A. lixula, our results show different levels of expression of some gene isoforms under thermal treatments.
In summary, and despite the limited proportion of annotated with activation of many genes whose functions could be related to stress responses in the form of chaperone production, apoptosis regulation, ATP-associated genes, enhancement of the immune system and redox processes, and downregulation of gene pathways related to protein biosynthesis and DNA replication. Nevertheless, contrary to that found in other studies (e.g., Gleason & Burton, 2015;Zhu et al., 2016) no activation of genes encoding antioxidant enzymes was detected in our experiments. As we initially expected, a markedly lower response is found in the warm treatment, with no activation or deactivation of the previously mentioned pathways, with the exception of the apoptosis regulation. Although some caution is needed, as we have characterized transcriptional changes and not protein levels, the differential patterns found in these genes strongly indicated that sea urchins are more stressed under lowered experimental temperatures.
We acknowledge that we have tested only acute thermal conditions, without any progressive acclimation. This is an unrealistic scenario but was chosen to elicit a short-term measurable response.
This response was much more marked against lower than higher temperatures, which indicates potential to compensate for cold stress.
However, our results indicate that A. lixula might require energy expenditure to withstand the stress associated with low temperatures, while it does not undergo relevant transcriptional changes when exposed to warm temperatures. This is coherent with the notion of a thermophilous species living near the colder limit of its physiological tolerance, as found also when analysing reproductive and larval features (Wangensteen, Dupont, et al., 2013;. Future research should consider a wider panoply of temperature regimes and populations, combined with acclimation periods, to explore the potential effect of global warming and heat waves across warmer areas where the species inhabits. Additionally, considering different reproductive seasons would be also interesting since organisms may display different thermal sensitivity depending on their gonad maturation stage, and may experience changes in energy allocation. Increasing gene annotation quantity and quality for this species is also desirable to obtain meaningful physiological conclusions in further studies.
It has been suggested that the tropicalization of NW Mediterranean can lead to a shift in dominance between the temperate common sea urchin Paracentrotus lividus, which will suffer from warming temperatures, and the thermophilous black sea urchin A. lixula (Carreras et al., 2020;Gianguzza et al., 2011;Wangensteen, Dupont, et al., 2013;. Such a shift can have drastic ecological impacts, as both species are conspicuous engineer species shaping benthic communities (Bulleri, Benedetti-Cecchi, & Cinelli, 1999;Bonaviri et al., 2011). Specific biological and genomic studies are needed to understand the adaptive capabilities of A. lixula to ongoing warming, but our results add to the available evidence that colder rather than warmer temperatures may be a limiting factor for A. lixula. The absence of clear signs of stress at warm temperatures in adults of A. lixula, together with information on larvae development and gonad maturation (Wangensteen, Dupont, et al., 2013;, support the hypothesis of a positive effect of winter warming on the species' reproduction output and larval survival. The ongoing expansion of the species across the littoral coast of the Mediterranean, with the concomitant impacts of its grazing activity on littoral communities, may be exacerbated in the near future by rising winter temperatures in the NW Mediterranean.

AUTH O R CO NTR I B UTI O N S
All authors contributed to the design of this study and were involved in the aquarium experiments. R.P.P., A.R., X.T., and O.S.W.
analysed the data. R.P.P. wrote the first draft of the manuscript and created figures and tables. X.T., A.R., and C.P. contributed to improve the first draft, and all authors revised the final version of the manuscript.