Cellular stress responses to chronic heat shock and shell damage in temperate Mya truncata

Acclimation, via phenotypic flexibility, is a potential means for a fast response to climate change. Understanding the molecular mechanisms underpinning phenotypic flexibility can provide a fine-scale cellular understanding of how organisms acclimate. In the last 30 years, Mya truncata populations around the UK have faced an average increase in sea surface temperature of 0.7 °C and further warming of between 1.5 and 4 °C, in all marine regions adjacent to the UK, is predicted by the end of the century. Hence, data are required on the ability of M. truncata to acclimate to physiological stresses, and most notably, chronic increases in temperature. Animals in the present study were exposed to chronic heat-stress for 2 months prior to shell damage and subsequently, only 3, out of 20 damaged individuals, were able to repair their shells within 2 weeks. Differentially expressed genes (between control and damaged animals) were functionally enriched with processes relating to cellular stress, the immune response and biomineralisation. Comparative transcriptomics highlighted genes, and more broadly molecular mechanisms, that are likely to be pivotal in this lack of acclimation. This study demonstrates that discovery-led transcriptomic profiling of animals during stress-response experiments can shed light on the complexity of biological processes and changes within organisms that can be more difficult to detect at higher levels of biological organisation. Electronic supplementary material The online version of this article (10.1007/s12192-018-0910-5) contains supplementary material, which is available to authorized users.


Introduction
Acclimation via phenotypic flexibility is the process in which an individual adjusts to a change in environmental conditions allowing maintained performance and fitness. It has been proposed as the fastest, and potentially most important, response to climate change, particularly for long-lived species with extended generation times and the consequent limited capacity for genetic adaptation (Somero 2012;Bay and Palumbi 2015). Understanding the physiological stress that results from a range of global change-related abiotic and biotic variables, and the associated acclimation mechanisms, is pivotal to deciding an organism's fate in a given local population (Somero 2012). Many of the previous studies on the acclimation of marine invertebrates have concentrated on responses to warming and used upper lethal temperatures as a measure of successful acclimation to altered conditions (Somero 2010;Peck et al. 2014). More recently, the rise in next-generation sequencing technologies applied to non-model, environmental, species has provided more detailed data on the molecular mechanisms underpinning physiological stress, and hence delivers a more fine-scale cellular understanding of how organisms acclimate, or in some cases fail to acclimate, to changed environmental conditions (Somero 2012;Bay and Palumbi 2015;Moya et al. 2015;Clark et al. 2017).
The subject of this study, Mya truncata (Superorder = Imparidentia, Order = Myida, Superfamily = Myoidea, Family = Myidae), is a large sediment-burying bivalve. It is a long-lived species [reported ages up to 40 years (Hewitt and Dale 1984)] that plays an important functional role in benthopelagic coupling, sediment stabilisation, and bioturbation (Queiros et al. 2013). In terms of biotic interactions, M. truncata in temperate habitats are an important food source for many species, such as the sea star Asterias rubens and the gastropod Buccinum undatum (Himmelman and Dutil 1991;Morissette and Himmelman 2000;Gaymer et al. 2001a;Gaymer et al. 2001b). Despite the large biomass and important functional role of M. truncata within ecosystems, very little is known about the acclimation potential of this species, particularly at the molecular level.
Mya truncata has a boreal-arctic distribution spanning a wide latitudinal range. In marine databases (www. marinespecies.org, www.marlin.ac.uk), it is reported as occurring in the Arctic down to the Bay of Biscay in Europe (~79°N to 45°N) (Ballerstedt 2002;Gofas 2004;Oliver et al. 2016), but populations have also been identified in the Mediterranean: Monaco and the southern (Atlantic) shores of Spain [specimens in the Muséum National d'Histoire Naturelle, Paris (E. Harper, pers. comm., Checa, (1993)]. Ecological observations of this species have concentrated on the Southern North Sea as part of a long-term benthic macrofauna monitoring programme Witbaard et al. 2005). These studies have shown that numbers of M. truncata are drastically reduced with decreasing latitude, with up to 100-fold less individuals at latitudes around the Frisian front [< 1 individual per m 2 (Amaro et al. 2003)] compared to Arctic latitudes [> 100 individuals per m 2 (Welch et al. 1992)]. In addition, M. truncata frequently failed to recruit at the Frisian Front. From 1987 to 2001, the species was absent from survey records at the Frisian Front (Amaro et al. 2003), and individuals had low numbers of ripe oocytes indicating reduced reproductive fitness . Based on these data, it would appear that distributions of M. truncata are impacted by temperature, but that the species, in general, has a moderate window of physiological capability. Populations in the Arctic (circa. 79°N) can experience 9 months at sub-zero temperatures with 1.5 m of ice cover (Hewitt and Dale 1984), whereas temperate populations experience no ice cover and summer temperatures in excess of 15°C (in the Mediterranean for example). To date however, apart from the environmental information described above, there are no experimental data on either upper lethal temperatures or the thermal tolerances of M. truncata.
All across M. truncata's boreal-arctic distribution populations are experiencing warming as a result of anthropogenic climate change. In the last 30 years, M. truncata populations around the UK have faced an increase in sea surface temperature of 0.7°C (Frost et al. 2016). Using a medium emissions scenario, the UK Climate Projections 2009 (UKCP09) predict further future warming, in all marine regions adjacent to the UK, of between 1.5 and 4°C. Given the previous surveys suggesting a decrease in population densities with latitude, data are required on the ability of M. truncata to tolerate multiple physiological stresses, and most notably, chronic increases in temperature, in order to predict the likely ecological fate of this species in the future.
The aim of the present study was to evaluate the ability of Mya truncata to perform damage-repair under a future climate change scenario. Temperate subtidal M. truncata were acclimated to chronic warming (2 months at 16 ± 1°C; 1°C above current maxima) which was followed by artificially induced shell damage. This experiment represented the multi-faceted nature of survival in a changing world with warming and the ecological interaction of simulated predation. Transcriptomic profiling of individual animals over a 2-week time course after shell damage was used to provide insight into the molecular mechanisms underpinning responses to these two combined physiological stresses and the cellular impacts associated with life at a higher temperature.

Animal collection, habitat parameters and animal husbandry
Specimens of Mya truncata were hand-collected by SCUBA divers at the UK National Facility for Scientific Diving from Dunstaffnage Bay, North West Scotland (56°27′ 06.5′′ N 5°2 6′ 02.2′′ W) between April and May 2016 (Fig. 1a). All specimens were sexually mature adults (mean shell length = 65 ± 1 mm S.E). All collections were from a depth of 10-15 m. Data on habitat temperatures were obtained from Martin Sayer (UK National Facility for Scientific Diving). A temperature logger, continuously running on a mooring, on the edge of Dunstaffnage Bay at 10 m depth (Saulmore Point, 56°27′ 14.6′′ N 5°24′ 48.6′′ W) was used to measure temperature every 10 min-which was then converted into a daily average-from January 2000 to October 2016 (Fig. 1). These data were used to select an appropriate elevated temperature (16°C) for the heat-stress experiment. Heat-stress in this study was defined as 1°C higher than any ambient condition animals from the sampled population would experience in their natural habitat (Fig. 1b). In addition, an increase in temperature of just 1°C is a very conservative prediction of warming in the region sampled (using a medium emissions scenario, the UK Climate Projections 2009 (UKCP09) predict warming, in all marine regions adjacent to the UK, of between 1.5 and 4°C) and allows sub-lethal effects to be measured (such as acclimation mechanisms).
After collection, clams were immediately transferred to the Scottish Association of Marine Science (SAMS) flow-through aquarium (salinity 35-38 ppt) where they were randomly distributed into three tanks. An elastic band was placed around each specimen to prevent gaping and animals were laid on top of approximately 20 cm of sand, which they could bury into. The animals were maintained under a 12:12 h simulated light:dark cycle and fed a pre-mixed microalgal blend twice weekly throughout acclimation and experimentation (Shellfish Diet, Varicon Aqua solution Ltd., UK).

Experimental design
The experiment entailed an initial chronic, long-term exposure to elevated temperature followed by the response to shell-damage, which was assessed at the molecular level over a 2-week time course. For the long-term heat exposure, 40 individuals were kept at 16 ± 1°C (1°C higher than they would experience in their natural habitat, as per the habitat temperature data, Fig.  1b) for 2 months. Twenty of these heat-exposed animals were then subjected to shell damage. Shells were damaged by drilling three holes (approximately 2.5 mm in diameter) through the shell, just inside the pallial line close to the ventral edge, within reach of the mantle margin, using a 10.8 V Lithium-Ion Dremel cordless modelling drill (at a speed of 5000-35,000 rpm) fitted with a round-tipped bit to minimise any trauma to the underlying soft tissue as in Sleight et al. (2015). To inflict damage, animals were lifted out of the water, placed on crushed ice on a bench in the aquarium where their shells were drilled and returned immediately to the aquarium tank. Control animals were similarly removed from the tank and placed on ice, before returning to the tank, to control for handling stress. In order to assess repair after damage, shells were removed from the animal and examined for evidence of repair on both internal and external surfaces (Fig. 2). The remaining 20 animals were left as undamaged for controls. All animals were then maintained at 16°C before sampling of both control and damage individuals over a time course at 1, 3, 5, 7, and 14 days post-treatment (n = 4 per treatment, per time point).

Haemocyte counts and shell damage-repair
At each time point of the experiment (see above), animals were removed from each tank and were photographed to assess the degree of repair in the damaged shells (n = 4). The stage of repair was assigned to one of four qualitative repair categories (Fig. 3). Given the shell damage method caused a breach in the external barrier of the animal, it was likely to cause an immune challenge, and hence it was predicted the immune system would be involved in the response. As a proxy for immune activity, haemolymph (fluid equivalent to blood in most invertebrates [250 μL]) was extracted from the adductor muscle of each individual into a syringe containing 250 μL of 4% formaldehyde in 0.22 μm filtered sterile seawater and fitted with a 25 gauge needle. The total number of haemocytes (invertebrate immune cells) per millilitre (THC) of haemolymph was counted using an improved Neubauer haemocytometer under a light microscope. Mantle tissue samples were then taken across the three mantle folds from a region under the middle hole as described in Sleight et al. (2015).

Sequencing and bioinformatics
Total RNA was extracted from each mantle tissue sample on ice using Tri-Reagent (Bioline, UK), and purified using RNeasy columns including a DNase treatment (QIAGEN, UK), all performed according to manufacturer's instructions. The RNA samples were analysed for concentration and quality by spectrophotometer (NanoDrop, ND-1000) and tape station analyses (Agilent 2200 TapeStation). The best three RNA samples from the four animals extracted at each time point was sent to the sequencing centre. cDNA libraries were made for each individual (n = 30; 15 treated animals from 5 time points (n = 3 at each time point), with a corresponding number of controls) and library preparation was conducted by the Earlham Institute, Norwich, UK (formerly The Genome Analysis Centre). Stranded libraries were prepared using the NEXTflex™ Rapid Illumina Directional RNA-Seq Library Prep Kit and sequencing was carried out on a Hi-Seq 2000, generating 100 base paired-end reads. Reads went through an initial quality control process conducted by the Earlham Institute that removed Illumina adaptor sequences and ribosomal RNA reads. Reads were then further cleaned for quality (Phred score 30) and minimum read length (80 bp) using the eautils (v1.1.2) tool fastq-mcf (https://github.com/ ExpressionAnalysis/ea-utils/blob/wiki/FastqMcf.md). The cleaned reads were normalised using Trinity's (v2.2.0) In silico Read Normalisation tool (Haas et al. 2013), with default parameters. The left and right reads for each library were normalised; all of the left and right reads were then concatenated and the resulting file was normalised a second time with default parameters. The concatenated, normalised reads were de novo assembled using Trinity (v2.2.0) with default parameters (Grabherr et al. 2011). For a quantitative assessment of orthology completeness, the assembled transcriptome was subject to BUSCO analysis (Benchmarking Universal Single-Copy Orthologs, v3) using the entire Metazoa dataset (metazo_odb9 downloaded January 2018). Transcript abundance was estimated by alignment-based quantification using Trinity (v2.2.0) utilities (Grabherr et al. 2011;Haas et al. 2013). Transcripts were aligned to the de novo transcriptome using bowtie with default parameters and transcript abundance estimates were calculated using RNA-Seq by Expectation-Maximisation (RSEM). The gene- Fig. 3 Photographs of shell repair categories. a Category 1, thin clear film viewed from outside of shell. b Category 1 viewed from inside of shell. c Category 2, translucent brown proteinacious film viewed from outside of shell. d Category 2, viewed from inside of shell. e Category 3, partially calcified brown film viewed from outside of shell. f Category 3 viewed from inside of shell. g Category 4, fully calcified layer, viewed from outside of shell. h Category 4 viewed from inside of shell. 1 mm scale down the right-hand side. N.B. For illustrative purposes only, these are not M. truncata shells from this study as animals in this study did not heal past category 3 level abundance estimates (raw counts) for each of the libraries were constructed into a matrix for downstream expression analyses (using the Trinity abundance_estimates_to_matrix.pl script). Transcript abundance estimation was quality checked using the Trinity Perl-to-R 'PtR' toolkit and principle component analysis was used to check for batch effect and outliers (Haas et al. 2013). Outliers were identified (one sample from the 1-day time point and two samples from the 5-day time point) and removed, which significantly unbalanced the experimental design and hence the 5-day time point could not be included in downstream analysis ( Supplementary  Information, Fig. S1, for more details about outlier were removal). Differentially expressed genes, between control and damaged animals at each time point, were identified using the Bioconductor (v3.4) edgeR package in R (v3.1.1) with a false discovery rate (FDR) of 5% and a twofold change cut-off (Robinson et al. 2010;McCarthy et al. 2012). In the same analysis, genes that showed significant temporal changes in response to damage, termed here, time-dependant damage-response genes, were identified. The full R script for all downstream analyses and the count data matrix are available (Supplementary Files S1 and S2).
The longest isoform of each gene was extracted from the transcriptome for annotation as in An et al. (2014), using the Trinity utility script 'get_longest_isoform_seq_per'. The longest isoforms of each gene were compared to a local NCBI non-redundant (NR) database (updated 01 June 2016) using Basic Local Alignment Search Tool (blastx, cut-off <1e −10 ) to search for sequence similarity and putative gene annotation (Altschul et al. 1990).
To aid interpretation of the hundreds of up-and down-regulated genes, and to provide a visual qualitative assessment of the different biological processes at each time point, STRING (v10.0) was used to produce protein-protein interaction networks (Franceschini et al. 2013;Szklarczyk et al. 2015). All of the differentially expressed genes at each time point were compared to a local UniProt/SwissProt human database (updated 05 October 2016) using Basic Local Alignment Search Tool (blastx, cut-off <1e −10 , Table S1) and these identifiers used as input into STRING. STRING only accepts identifiers from UniProt/SwissProt accession numbers and human data are the most abundant and best annotated dataset available, providing the most comprehensive coverage compared with other model species, in spite of the evolutionary distance (Clark, pers. comm.). Protein interaction scores in STRING were set to 'high confidence 0.7', only the query proteins were included in the network and all non-interacting proteins were removed from the analysis. Functional enrichment analysis of each network was performed using the whole human genome as a statistical background (as above, STRING uses human annotations and therefore compared to human genome as a background). In addition to the STRING protein-protein interactions, all annotations (blastx against NR and UniProt/SwissProt human) of the differentially expressed genes were manually interrogated and assigned a putative function.

Availability of data and material
All raw sequence data has been submitted to NCBI SRA and is publically available (accession number: SRP115712). The mantle transcriptome assembly is available upon request and the code and data matrix used in downstream analysis are available in the electronic supplementary material.

Results
Only 3, out of 20 damaged individuals, reached healing category 3 'holes covered by a brown film which shows signs of partial calcification' after 14 days (Figs. 3 and 4). There was no difference in THC between control and damaged individuals at any time point ( Supplementary Information, Fig. S2).
All of the RNA libraries (n = 30) were cleaned, normalised and assembled into a de novo transcriptome (Table 1). The assembly and downstream transcript abundance estimates were quality-checked (Table 1 and Supplementary  Information, Fig. S1) and deemed suitable for further analysis, 16% of Trinity genes were assigned putative annotation using blastx sequence similarity searching against the NR database (below an e value of 1e −10 ).
Protein-protein interaction networks were constructed out of the annotated, differentially expressed genes between damaged and control animals at each time point. The networks were used to provide a visual qualitative assessment of the different biological processes at each sample time. One day after damage treatment, transcription and protein transport were the differential processes between damage and control animals. Three days after damage, transcripts associated with RNA processing, cell cycle regulation, wound healing, immunity, ion transport, biomineralisation, stress responses, protein turnover, lipid metabolism, amino acid metabolism and muscle contraction were differentially expressed. Seven days after damage, differences in transcriptomic profiles between experimental and control animals were those involved in RNA processing, cell cycle regulation, stress response, amino acid metabolism and apoptosis. Finally, 14 days after damage, the differences in the transcripts were those governing: cell cycle regulation, stress responses, lipid metabolism, apoptosis, immunity, cell differentiation and development. Each network was analysed for functional enrichment (FDR < 0.05) using GO terms, against the whole genome background. One day after damage, there were no functional enrichments whereas, after 3 days, 12 biological processes (including immune and stress responses) and 9 cellular processes (mainly extracellular and vesicle associated processes) were enriched. Subsequently, 7 days after damage, 14 biological processes (including biomineralisation and immune processes) and 8 cellular processes (again mainly extracellular and vesicle related processes) were enriched. By 14 days post-treatment however, transcripts for only two biological processes (stress and immune) and no cellular processes were enriched (Table 2).
In addition to investigating differential expression within each time point, genes that showed significant temporal changes across time in response to damage (termed here, time-dependant damage-response genes) were also identified. The top 50 most significantly different (smallest FDRs), annotated genes, were selected from 6103 identified timedependant damage-response genes (Supplementary File S4). Genes that were increasingly up-regulated over the time series (in damaged compared to control animals) included putative candidates in biomineralisation such as collagen, calumenin, matrilin and pif; the immune response such as calumenin (n.b. multi-functional), theromacin and L-rhamnose-binding lectin, and also the stress response, such as the detoxification genes cytochrome P450 3A24 and the chaperone heat shock protein 60.

Discussion
The present study shows that M. truncata displays evidence of a lack of acclimation under a predicted end-of-century climate change scenario with a consequential inability to successfully repair shell damage. Transcriptomic profiling revealed consistent evidence for stress responses to be the key players in the molecular mechanisms underlying this lack of acclimation and ability to repair shell damage. At each time point, genes  were differentially expressed between control and damaged animals (Fig. 5) and, those differentially expressed genes were functionally enriched with known cellular stress response candidates (Table 2, Fig. 6, Supplementary Information, Supplementary file 1). In addition, genes that showed the most significant temporal changes in expression between damaged and control animals across the time series, i.e. time-dependant damage-response genes, also included classic stress response genes such as the detoxification genes cytochrome P450 3A24 and the chaperone heat shock protein 60 (Supplementary file 4). Previous work has shown that M. truncata frequently have low fecundity and a failure to recruit at the Frisian Front (a region of relatively similar latitude and temperature range to the specimens used in the present study). Furthermore, previous work has also shown that molluscs (Margaritifera margaritifera) living at high latitudes with a longer lifespans have better wound and shell-healing capacity than those living at lower latitudes (Ziuganov et al. 2000). The difference in repair capacity between high and low latitude populations of M. margaritifera was hypothesised to be due to the reduced energy expenditure for growth, and therefore greater energy availability for shell and tissue regeneration to sustain homeostasis. Interestingly, neither the Frisian Front nor Dunstaffnage Bay are close to the southern range boundaries of this species (E. Harper, pers.comm., Checa 1993) and therefore do not represent populations living in the warmest sea temperatures this species is known to inhabit. M. truncata are broadcast spawners with planktonic larvae (Amaro et al. 2003), but to date, very little is known about the population genetics of this species or the dispersal capabilities of the larvae. The most likely explanation therefore is that there is a lack of connectivity between these northern populations and those in the Mediterranean, which are more adapted to warmer waters. Unpublished population genetics data on a sister species, Mya arenaria, support this explanation (specifically, distinct populations in Scotland, northern Europe, Mediterranean and Spain/Portugal; M. De Noia & D. Vendrami, pers. comm.). Previous studies have shown that vulnerability to warming in a species is context dependent and does not follow a strict latitudinal pattern (Sagarin and Somero 2006;Osovitz and Hofmann 2007). Hence, combining the experimental data presented in the present study, with the recent trend of oceanic warming around the UK of 0.7°C in the last 30 years (Frost et al. 2016), suggests that some temperate M. truncata populations are being pushed past their viable thermal tolerance limit and possibly impacting the availability of cellular energy levels, reflected by an inability to repair shell damage.
The transcriptional profiling results show activation of cellular stress responses. Functional enrichment analysis of the differentially expressed genes between control and damaged animals demonstrated that cellular stress responses were enriched at 3, 7 and 14 days post-damage. These responses could be split into two broad categories, the classical heat shock response and antioxidant production.
Heat shock proteins (hsps) are involved in protein folding and chaperoning. They are either constitutively expressed as part of normal cellular function or induced in response to stress (Hartl 1996) Fig. 5 Smear plots of the up-and down-regulated genes between damaged and control animals at each time point a 1 day after damage, b 3 days after damage, c 7 days after damage and d 14 days after damage. Grey dashed lines indicate log 2-fold cutoff and genes which were found to be significantly different by edgeR are indicated in red, the number of significantly up-and downregulated genes indicated on the left of the smear proteins and also target degraded proteins by regulating their removal (Fink 1999). In an earlier study, a mantle transcriptome was produced from animals sampled within days of collection and revealed that two hsp genes were highly expressed, heat shock protein 90-alpha 1 and heat shock protein 70 . It was suggested that one possible explanation for the high background expression of hsps, particularly the inducible hsp70, was that the animals were exhibiting low levels of thermal stress. In the present study, multiple hsps were significantly up-regulated in damaged animals, compared to undamaged controls, at the 3-, 7-and 14day time point (Supplementary File 3) and in addition, two hsps were detected in the top 50 most significant timedependant damage-repair genes ( Fig. 6 and Supplementary File 4). These up-regulated hsps included two members of the hsp70 superfamily, one of which was clearly hsp7012A, an hsp60 and two small heat shock proteins (Table 3).
Detailed comparison of the other HSP70 family member described in this study with the hsp70 identified in Sleight et al. (2016) showed clear differences with 71.8% sequence identity at the amino acid level ( Supplementary Information, Fig. S3). Supplementary Information,  Fig. S3). BLAST sequence similarity searching discounted the likelihood of either sequence as being GRP78 and therefore they were both designated as hsp70 (Table 3). It is likely that since the hsp70 in Sleight et al. (2016) was identified under 'normal' environmental conditions, this is putatively the constitutive form (hsc70), whilst the transcript identified in the present study is the inducible form, (hsp70) upregulated in response to chronic heat shock and shell damage, future functional analyses will be required to confirm this.  Fig. 6 STRING Database predicted protein-protein interactions built from differentially expressed genes between control and damaged treatments. a 1 day after damage, b 3 days after damage, c 7 days after damage and d 14 days after damage. Line colour and type depict molecular actions and effects as per the network key and manually added italicised labels indicate most likely associated biological processes or pathways. Please see Table S1 in supplementary information for all abbreviations The induction of heat shock protein genes in response to various types of stress has been documented in many marine invertebrates and although these can be viewed as a generic response, the up-regulation of particular genes is often context-dependent (Schill et al. 2003;Chang 2005;Clark and Peck 2009;Huang et al. 2013). The ability to predict the vulnerability of a species to environmental change is best achieved at the molecular level, as sub-lethal effects across a range of functions can be more easily detected and evaluated (Clark and Peck 2009). In addition to identifying up-regulated hsps within and across time points in the present study, protein-protein interaction networks indicated that upregulated hsps could be acting as molecular hubs, interacting with a variety of functional pathways including immune response, ion transport and ubiquitination/protein turnover (Fig.  6) (Csermely 2004;Korcsmaros et al. 2007). More specifically, protein-protein interactions within functional pathways revealed that up-regulated hsps could be regulating stressinduced cell death. The hsp cluster in the 3-day post-damage protein-protein interaction network included three connections to a mesencephalic astrocyte-derived neurotrophic factor gene, which has been shown to mediate stressinduced cell death . This is in line with experiments showing that failure to repair stress-induced unfolded proteins and to remove aggregated proteins leads to the activation of apoptotic cell death pathways (Powell et al. 2005). The present study provides further evidence that hsps underpin cellular responses to stress and their high expression in damaged animals, who were unable to repair shell damage after chronic exposure to elevated temperature (+ 1°C), points to a molecular mechanism which is likely to be critically important in stress-tolerance this species.

The sequences of both transcripts included the classic motifs for members of the HSP70 superfamily, in particular the R-A-[RK]-F-E-[ED]-[LM] motif characteristic of HSP70, HSC70 and GRP78 proteins (
The other major functional category of up-regulated genes was antioxidant activity. Oxidative stress can be described as an imbalance between the production of reactive oxygen species (ROS) and the ability of an organism to counteract or detoxify their harmful effects via antioxidants (Jones and Go 2010). Transcriptomic profiling revealed that oxidative stress was most prominent 14 days after damage where antioxidant genes such as glutathione synthase-like (Noctor and Foyer 1998;Wu et al. 2004) and spermine oxidase-like (Bouchereau et al. 1999) were significantly up-regulated Table 3 BLAST sequence similarity searching results for the up-regulated heat shock proteins identified in both Sleight et al. (2016) File 3). In addition, cytochrome P450-one of many genes which encodes an enzyme that most commonly catalyses a monooxygenase reaction to metabolise xenobiotics-was detected as one of the 50 most significant timedependent damage-response genes (Supplementary file 4). As well as generating potent ROS during xenobiotic metabolism, cytochrome P450s have also been shown to have protective effects against oxidative stress, and therefore, overall, must contribute to an organism's oxidative balance (Caro and Cederbaum 2004;Lingappan et al. 2014). Accurately assessing oxidative balance in an organism, and interpreting what different oxidative states mean for ecology and evolution, is an incredibly complex problem (Cohen et al. 2010;Haussmann and Marchetto 2010). The combination of chaperone production and the response to ROS, as seen here however, is indicative of a failing physiological state.
The differentially expressed genes between control and damaged animals were functionally enriched with immune processes both 3 and 7 days post-damage (Table 2). Like all protostome invertebrates, molluscs rely entirely on the innate defences for host protection. These defences comprise two main parts: recognition and direct effector mechanisms. Regarding the recognition phase, 3 days post-damage, differentially expressed genes were functionally enriched with five receptor-based GO terms (for example, Toll-like receptor 3 signalling pathway, FDR = 0.02). More detailed interrogation of the gene sequences that were significantly up-regulated in the damaged animals includes one encoding a fibrinogen-like domain (Supplementary File S3). Rather than this sequence being the full fibrinogen molecule itself (which is primarily involved on clotting in vertebrates), it is more likely that it is part of a gene that encodes one of the hypervariable family of fibrinogen-related domain (FReD) proteins (Hanington and Zhang 2011). FReDs are pathogen recognition receptors that are responsible for binding non-self-moieties. The presence of FReDs in the differentially expressed genes is therefore highly indicative of an immune response by M. truncata. Other genes that point to recognition events in the present study include those that encode for various lectins (most notably ficolin) as well as C1q and leucine-rich repeat-containing proteins (LLRs). Lectins are glycoproteins that have two sugarspecific binding sites. Their binding properties facilitate recognition and agglutination of non-self-materials with ficolins binding n-acetyl glucosamine, a molecule present in the cell walls of bacteria. C1q is expressed in a wide variety of tissues, including the bivalve mantle (Gerdol et al. 2015) and, although it is regarded as a pattern recognition receptor in these animals, in mammals it is part of the classical complement system that binds antibody. Bivalves do not express immunoglobulin antibodies but their C1q protein binds lipopolysaccharide, peptidoglycan, and β-glucan (Zhang et al. 2008;Wang et al. 2012), all of which are pathogen-associated molecular patterns (PAMPs) on microbial surfaces. Lastly, among the differentially expressed genes found in M. truncata that represent immune recognition are the LLRs. In many species, including vertebrates, proteins with leucine-rich repeats are well represented as hypervariable receptors and are best known as Toll or Toll-like molecules (Zheng et al. 2005) and the variable lymphocyte receptors of lamprey (Pancer et al. 2004). All of these genes, except C1q, were upregulated at day 3 of the study, which fits with the roles of the encoded proteins being part of an early immune response. Fibrinogen, ficolin and LLR genes were still highly expressed at 7 days, and ficolin, C1q and LLR were still up-regulated by day 14, possibly demonstrating immune recognition occurring in response to ingress of potential pathogens from the surrounding water through the wounds during the slow healing process.
Genes differentially expressed between control and damaged animals were also functionally enriched with processes relating to immune effector mechanisms, for example the GO cellular process melanosome 3 days after damage (FDR = 0.042). In addition, genes up-regulated in damaged animals include those encoding theromacin, a thioster-containing protein and some apoptosis-related proteins. Melanosomes are 'lysosome-related' organelles in which melanins are synthesised and stored. It is well known that melanins are a crucial component of the invertebrate immune response (Grimaldi et al. 2012), and parallels have been made between melanosomes in mammalian melanocycte and melanin-rich granules in invertebrate immune cells. Mydlarz et al. (2008), for example, describe melanosomes detected in amoebocytes adjacent to protective melanin bands in infected sea fans, which were absent in uninfected controls. Of particular relevance to the present study, melanin production has repeatedly been demonstrated as an effector mechanism in bivalve immunity (Luna-Acosta et al. 2017). Theromacin is an antimicrobial protein, originally identified in a leech, which is active against Gram-positive bacteria (Tasiemski et al. 2004). A similar gene has been reported for bivalves, which presumably has the same function (Xu et al. 2010). Thioester-containing proteins are a large and diverse group that are present in both vertebrates and invertebrates. In terms of immunity, they are best known as complement factors, C3 and α 2 -macroglobulin. In higher chordates, C3 acts as an opsonin that is deposited on pathogen surfaces, while α 2 -macroglobulin acts as a caging and anti-proteinase molecule (Williams and Baxter 2014). It is likely that the thioester-containing gene in M. truncata encodes a protein with similar effects. The up-regulated apoptosis-related genes are likely to encode proteins involved in programmed cell death, and would be expected to be upregulated during recovery from drilling, as apoptosis helps to remove redundant and degenerating cells, as well as acting as a defence strategy against viral threats. Taken together, the results discussed above demonstrate that the immune system plays an important role in response to temperature and shelldamage stress in M. truncata.
It was surprising that there were no significant differences in the total haemocyte counts between the control and damaged animals ( Supplementary Fig. S2, File 3). Previous studies have found that, under stressful conditions (e.g. pollution, acidification and temperature), haemocyte numbers can increase or decrease (Dyrynda et al. 2000;Matozzo et al. 2012). It is possible that any rise in haemogram in the present study may have been masked in the damaged animals by the loss of haemolymph following drilling (through the hole). An alternative explanation could be that the haemocytes migrated to the damage site during the repair process, which has been documented in other bivalves (Mount et al. 2004;Kadar 2008;Li et al. 2016) and hence an active haemocyte production process could have been masked. The exact role of haemocyctes in both normal biomineralisation and shell repair is still unclear; in both scenarios, these cells are likely to be playing a dual role in the immunoprotection of the animal at its external barrier, and in shell production. In addition to total haemocyte numbers, the different haemocyte types (here, DCC) is also an important consideration. Bivalves possess diverse haemocyte types, primarily granular and agranular subpopulations, which perform different immune functions (e.g. Wootton et al. 2003). The DCC under varying conditions, therefore, should ideally be recorded. The low total number of haemocytes in the present study made it impossible to obtain reliable DCC values. Thus, it remains a possibility that there were variations in this parameter, especially as there are always large inter-animal variations in invertebrate animals taken from the wild (Wootton et al. 2003).
Although the primary purpose of the present study was to understand general bivalve cellular stress responses, some of t h o s e r e s p o n s e s w e r e i n e v i t a b l y s h e l l r epair and biomineralisation mechanisms given the specific experimental shell-damage challenge imposed. Inflicting shell damage stimulates calcification pathways to repair the shell and hence, this experimental manipulation has previously been used to study biomineralisation mechanisms (Mount et al. 2004;Fleury et al. 2008;Sleight et al. 2015;Huning et al. 2016). Only three animals produced a brown film with signs of partial calcification to occlude the drill hole (typical of healing category 3) after 14 days. Despite this lack of physiological repair, transcriptomic profiling revealed up-regulated biomineralisation genes. At the 3-and 7-day time points, the protein-protein interaction networks were significantly functionally enriched with cellular processes linked to biomineralisation, such as extracellular regions, membranebound vesicles and extracellular exomes ) ( Table 2). In addition, at the 7-day time point, osteoblast differentiation and ossification terms were also functionally enriched, albeit with terminology specific to human biomineralisation as human annotations were used for input into STRING (Table 2). The annotated differentially expressed genes also included putative biomineralisation candidates as previously identified in molluscs (and therefore not in the human-centric STRING database, Supplementary files 3 and 4). Pif is a well-characterised shell matrix protein (Suzuki et al. 2009;Suzuki et al. 2013) that is present in the M. truncata shell proteome ) and has a mantle-specific expression pattern when compared to other tissues in this species ). In the present study, pif was highlighted as a one of the most significant timedependant damage-response genes as it was progressively up-regulated in damaged animals throughout the course of the experiment (LogFC values 1 day = 1.4, 3 days = 3.8, 7 days = 4.1 and 14 days = 4.2, Supplementary File 4), and hence, the data presented here provides further evidence for the important role of this gene in molluscan shell biomineralisation. In summary, these analyses show that discovery-led transcriptomic profiling of animals during stress-response experiments can shed light on the complexity of important biological processes and changes within organisms that cannot be detected at higher levels (we could detect biomineralisation signals in the transcriptomic data despite the lack of progress in physiological shell repair).

Conclusions
Here, we present the first experimental assessment of the cellular response to environmental stressors in an important temperate bivalve species, M. truncata. Failure to acclimate to increased temperature, as evidenced by the inability to rapidly repair shell damage, suggests that some temperate M. truncata populations are being pushed past their viable thermal tolerance limit for future ecological sustainability. At the cellular level, this lack of acclimation was clearly associated with key molecular mechanisms, and most specifically the up-regulation of classic stress response genes such as hsps, antioxidants and immune genes.