Introduction
The red king crab (Paralithodes camtschaticus) support economically and socially valuable fisheries in Alaska, USA (Stevens, Reference Stevens2014), many of which crashed in the early 1980s (Dew & McConnaughey, Reference Dew and McConnaughey2005; Bechtol & Kruse, Reference Bechtol and Kruse2009), and have been sustained at low levels since then (Stevens, Reference Stevens2014). Reasons for the fishery collapses are largely unknown but are thought to be due to overfishing (Bechtol & Kruse, Reference Bechtol and Kruse2009), shifts in ocean temperatures (Loher & Armstrong, Reference Loher and Armstrong2005) that shift crab distributions (Zacher et al., Reference Zacher, Kruse and Hardy2018), and increased abundance of fish that prey upon juvenile crabs (Bechtol & Kruse, Reference Bechtol, Kruse, Kruse, Eckert, Foy, Lipcius, Sainte-Marie, Stram and Woodby2010) and otherwise negatively impact recruitment (Zheng & Kruse, Reference Zheng and Kruse2000), behaviour (Daly et al., Reference Daly, Swingle and Eckert2009) and survival of early life stages, especially glaucothoes (Stevens, Reference Stevens2003). While most red king crab fishing areas have been closed for decades, catch management has supported sustainable Bristol Bay red king crab fisheries in recent years, albeit at a lower level (Siddeek, Reference Siddeek2003; Siddeek & Zheng, Reference Siddeek and Zheng2007). However, there are concerns that continued ocean warming and acidification could negatively impact red king crab populations, resulting in a complete collapse of the fisheries (Seung et al., Reference Seung, Dalton, Punt, Poljak and Foy2015).
Decreased ocean pH caused by increased pCO2 reduces growth and survival in embryonic larval and juvenile red king crab (Long et al., Reference Long, Swiney and Foy2013a; Swiney et al., Reference Swiney, Long and Foy2017), especially at the extreme low pH of 7.5 (Long et al., Reference Long, Pruisner, Swiney and Foy2019). Exposure to the extreme low pH 7.5 as a constant environment for an extended period of time is unlikely to reflect conditions these crabs will experience in nature until well beyond the year 2100 (pH 7.5 reflects a pCO2 of ~1630 μatm; Long et al., Reference Long, Swiney, Harris, Page and Foy2013b), but does reflect seasonal and regional scale conditions occasionally reached in deep waters of the Bering Sea where red king crabs live (Mathis et al., Reference Mathis, Cross and Bates2011). Nonetheless, to examine sensitivity of crabs to low pH/high pCO2, levels of CO2 even higher than 1600 μatm are commonly used (Schiffer et al., Reference Schiffer, Harms, Portner, Lucassen, Mark and Storch2013; Zittier et al., Reference Zittier, Hirse and Portner2013; Harms et al., Reference Harms, Frickenhaus, Schiffer, Mark, Storch, Held, Pörtner and Lucassen2014).
Carry-over effects of environmental conditions during the embryonic stage on larval condition has been observed in red king crab (Long et al., Reference Long, Swiney and Foy2013a). The pH/pCO2 environment of brooded embryos shifts red king crab larval sensitivity to pH/pCO2, causing shifts in larval survival at lower pH/higher pCO2 (Long et al., Reference Long, Swiney and Foy2013a). This carry-over effect has not been observed in another boreal crab, the southern Tanner crab (Chionoecetes bairdi) (Long et al., Reference Long, Swiney and Foy2016), suggesting that red king crab embryos may shift cellular physiology in a manner that has long-lasting effects on early life stages. Whether those shifts occur, and if they do occur, whether they involve sustained changes in gene expression is unknown.
Juvenile red king crab has a maximal prolonged temperature for maximal performance of between 12–17°C (as indexed by feeding ration), and upper thermal tolerance is 24°C (Long & Daly, Reference Long and Daly2017). Because their thermal limits result in a large thermal safety margin, it is unlikely that stress response induction will be a large component of their response to increased temperature, as has been observed in intertidal crabs with small safety margins (Teranishi & Stillman, Reference Teranishi and Stillman2007; Stillman & Tagmount, Reference Stillman and Tagmount2009).
Warmer temperatures reduce the intermoult duration in juvenile red king crabs, suggesting faster growth (Swiney et al., Reference Swiney, Long and Foy2017). Interactive effects of decreased pH/increased pCO2 and increased temperature are apparent in the percentage of change in carapace length in the first growth interval of juvenile red king crab (Swiney et al., Reference Swiney, Long and Foy2017). At ambient pH/pCO2, average growth was reduced by ~5% at +4°C, but not +2°C, whereas at pH 7.8, the ~5% reduction in growth was at +2 and +4°C (Swiney et al., Reference Swiney, Long and Foy2017).
In this study, we investigated sensitivity to different levels of pH/pCO2 and temperature of red king crab across life stages using a high-resolution transcriptomic approach, RNA-seq. The central goals of this study were to determine the magnitude and patterns of gene expression responses (i.e. a transcriptomic fingerprint of response to acidification and warming) as well as to make mechanistic inferences about how red king crab respond to shifts in temperature and acidification at different life stages.
Methods
Specimen collection and acclimation
Adult experiment
Female adult ovigerous crabs (N = 9) were collected during the NOAA Fisheries eastern Bering Sea bottom trawl survey, shipped to the NOAA Fisheries laboratory in Kodiak, AK in June 2011. Three females were held in each of the three pH/pCO2 conditions of ambient (pH 8.1), intermediate acidified (pH 7.8) and highly acidified (pH 7.5) at ambient temperature and salinity of local seawater in Kodiak throughout embryological development. Reduced pH water was produced by bubbling CO2 into a primary acidification tank to reduce pH of ambient inflow water to 5.5, and that low pH water was mixed with ambient water to the set-point pH in header tanks, as described in previous publications (Long et al., Reference Long, Swiney, Harris, Page and Foy2013b; Long et al., Reference Long, Swiney and Foy2016; Swiney et al., Reference Swiney, Long and Foy2017). The 2012 winter and spring had unusually cold seawater temperatures of 3–4°C between January and April when the experiment was terminated. At termination of the experiment, following larval release between 14 March and 15 April 2012 (see below), crabs were sacrificed and tissues (heart, anterior gill, posterior gill, cuticle) dissected from each specimen were preserved in RNAlater and stored at −20°C.
Notably during the adult collection (10 months prior to sampling for gene expression), shipment and embryonic brood period there were several mishaps that caused stress to embryos and potentially adults. Some embryos were physically damaged during transport from Dutch Harbor to Kodiak and subsequently died, ovigerous females were exposed to high temperature and low dissolved oxygen one day due to a chiller malfunction, the ocean acidification system malfunctioned one night causing a shift in exposure pH, and as mentioned above the incoming seawater during the winter was much colder than typical. The extent to which those events may have influenced the results for adults and larvae used in the study are unknown, but because they happened many months prior to sampling, ample time for recovery existed.
Larval experiment
There were two types of larval samples. The first type were newly hatched (Day 0) larvae that had been brooded by N = 3 mothers held at ambient temperature and three pH/pCO2 conditions (ambient pH 8.0, intermediate pH 7.8 and low pH 7.5) throughout the final 7–8 months (depending on hatch timing) of embryonic development of the year-long brooding duration. Thus, this experiment examines transcriptional differences that may have occurred due to embryonic development under different conditions.
The second type of larval experiment was a three-way split-brood design using larvae from N = 3 mothers that had been held under ambient conditions during the embryonic development period (ambient temperature and pH, Figure 1). Between 1–2 days after hatching, the larvae from each mother were split into groups of ~N = 600 larvae each. For two mothers (female IDs 1603 and 1572), there were 9 replicate groups, and for the third mother (female ID 1594), which was less fecund, there were 6 replicate groups. Two (female 1594) or three replicate (females 1603 and 1572) groups of larvae from each female were held under each of three pH conditions (ambient pH 8.0, intermediate acidified 7.8 and low pH 7.5) for 7 days prior to preservation of whole larvae in RNAlater. Larval mortality during the 7 days was 6.0 ± 3.4% at pH 7.5, 6.6 ± 2.2% at pH 7.8 and 7.5 ± 3.8% at pH 8.0 (ambient) (mean ± 1 SD). Female 1594 and 1572 had higher mortality standard deviation variances (up to 6%) among replicate groups of larvae than female 1603 (up to 2%). No larval moulting was observed during the 7 days of the experiment.
Fig. 1. Temperature and pH acclimation conditions during the juvenile red king crab acclimation experiment. All specimens started at ambient conditions, the lowest temperature and highest pH (bottom right of plot) and over first 4 days of acclimation were adjusted to the test temperature and pH/pCO2 values shown here. Data represent means ± 1 standard deviation of measurements of pH measured daily and temperature logged at 8-hour intervals across the final 20 days of acclimation. Letters above each data point indicate statistically significant differences in temperature, and numbers above each data point indicate statistically significant differences in pH (ANOVA, Tukey HSD, adjusted P < 0.05). Temperature and pH/pCO2 conditions were generated and characterized as described in Long et al. (Reference Long, Swiney, Harris, Page and Foy2013b) and Swiney et al. (Reference Swiney, Long and Foy2017).
Juvenile experiment
Juvenile red king crab were raised from ovigerous females collected in the eastern Bering Sea, shipped to the NOAA Fisheries laboratory in Kodiak, AK in June 2011, and held in ambient local seawater conditions through larval hatching and early juvenile development. On 17 July 2012 the juveniles were weighed and randomly placed into individual holding pens with flowing water at ambient temperature of 10.5°C and ambient pH of 8.06. Initial specimen weights were 12 ± 5 mg (mean ± 1 SD), median 10.6 mg, and varied from a minimum of 7 mg to a high of 26 mg. Specimens were mainly intact and in good condition, though 9 individuals (out of 107) were missing or regenerating a limb. During the study 1 individual attempted to moult on day 16–17, but was deceased on day 18. Temperature conditions were gradually adjusted from ambient to warmer (+2 and +4°C from ambient) over the first 4 days of acclimation, and remained relatively constant over the final 20 days, yielding average treatment temperatures of 10.5, 12.5 and 14.5°C (Figure 1). During the entire temperature treatment specimens were held under one of three pH conditions (pH ~8.0, ~7.8 and ~7.5; Figure 1). pH levels were not exactly the same across treatment temperatures due to the thermal influence on pCO2 (Figure 1). The tanks were set up in a flow-through design (no recirculation of water). Salinity was 31.9 on average throughout the experiment, though increased from a starting value of 31.7 on 18 July to 32.1 by 8 August 2012. Representative carbonate chemistry parameters are in Swiney et al. (Reference Swiney, Long and Foy2017). Specimens were fed daily on Gelly Belly protein and vitamin-enriched gelatin-based marine aquaculture food (Florida Aqua Farms, Dade City, FL, USA). Gelly Belly contained minimal extractable RNA, and thus was not considered a potential source of contamination. Samples for transcriptomic analyses were taken at day 24 of the acclimation on 10 August 2012 by placing whole crabs in RNAlater and storing at −20°C.
RNA extraction and RNA-seq library preparation and sequencing
Whole crabs (10 pooled larvae, individual juveniles) or 2–3 mm3 chunks of dissected tissues (adults) were individually homogenized with stainless steel ball bearings in lysis buffer using a TissueLyser II (QIAGEN 85300) and RNA isolated using RNeasy Mini Kit (QIAGEN 74106). RNA was checked for integrity using a Bioanalyzer 2100 RNA chip (Agilent). Only samples with little to no degradation observable on Bioanalyzer spectra and adequate concentration were used in downstream steps.
Library construction
For adult tissue samples, RNA from each individual and tissue was used to make libraries; total of 36 libraries (9 adults, 4 tissues per adult). For larvae, each RNA extraction from N = 10 larvae was used as an independent replicate for construction of RNA-seq libraries. There were a total of 33 larval libraries, 9 libraries from day 0 larvae, and 24 libraries from day 7 larvae. For juvenile crabs, we made N = 5 independent RNA pools using 1 μg of total RNA from each of N = 10 individual crabs from each of the 9 acclimation conditions (Table S1); total N = 45 RNA pools. Pooling RNA from all individuals focused our analysis on the mean response of each acclimation group across all individuals. Though this approach eliminates the analysis of individual variation in RNA-seq data, targeted gene-specific analyses from individual samples using more sensitive techniques (e.g. qPCR, NanoString) confirm that RNA-seq data from pooled samples reflect the general patterns of differences among sample groups when biological replication is used (Figure S1). Hence pooled sample data reflect the magnitude and directionality of differences among sample groups even though actual variances are smaller (i.e. only represent technical variation) than across all individuals in the sample.
For all samples, sequencing libraries were prepared using the Illumina Truseq RNA v2 kit following the manufacturer's recommended protocol. For each library, DNA fragment size distribution and concentration were determined using an Agilent Bioanalyzer 2100 using the high sensitivity DNA chip. Median library size was 400–450 bp across samples. Samples were submitted to the Vincent J. Coates Genomics Sequencing Laboratory at UC Berkeley for multiplexed 100 base pair paired-end sequencing on an Illumina HiSeq2000. Libraries were quantified by qPCR to ensure equal concentrations of each sample that were multiplexed on each lane. 22–24 samples were multiplexed on each lane.
RNA-seq library filtration and raw data processing
RNA-seq data were analysed on the Pittsburgh Supercomputer Center Blacklight. We used Trim Galore! (V 0.3.0) to filter raw reads: trimming low quality ends of reads (Phred score <20), trimming contaminating adapter sequence (stringency = 1), and discarding any resulting short reads (<30 bp). Trim Galore! was chosen because it implements Cutadapt for removing adapter contamination and retains paired-end matching. FastQC was run to check quality of each library before and after filtering. FLASh (V 1.2.8) was used to join overlapping paired end reads in order to make longer reads for de novo assembly and to eliminate artificial double counting of mapped read overlap regions. We aimed to initiate the assembly with the most information as possible; by making longer reads by pairing the overlapping regions of paired-end reads, we increased the information content of the data input into Trinity. Libraries are deposited at NCBI GenBank Short Read Archive (SRA), accessible with metadata through NCBI GenBank BioProject umbrella PRJNA530190 for BioProjects PRJNA523816 for larval samples (BioSamples SAMN10988674–SAMN10988704), PRJNA525257 for adult samples (BioSamples SAMN11044278-SAMN11044294) and PRJNA525256 for juvenile samples (BioSamples SAMN11044233-SAMN11044277).
Transcriptome assembly
Because there were no available genome or transcriptome resources for red king crab at the time of the study, with only a mitochondrial genome sequence available (Kim et al., Reference Kim, Choi, Park and Min2013), we constructed de novo transcriptomes using Trinity (V r2013-11-10) (Grabherr et al., Reference Grabherr, Haas, Yassour, Levin, Thompson, Amit, Adiconis, Fan, Raychowdhury, Zeng, Chen, Mauceli, Hacohen, Gnirke, Rhind, di Palma, Birren, Nusbaum, Lindblad-Toh, Friedman and Regev2011). Two transcriptomes were generated: a ‘Larval’ transcriptome using only larval samples, and an ‘All’ transcriptome using a combination of adult, larval and juvenile samples. For each transcriptome, a subset of individuals from across experimental conditions was combined in order to produce the most diverse pool of transcripts for transcriptome assembly. The assembly used a minimum kmer coverage of 2. Both transcriptomes were annotated against the SwissProt database using Trinotate (V 1.1), with BLASTx of nucleotide sequences, and BLASTp and PFAM analyses of TransDecoder (V 2.0) produced protein sequences. Trinotate annotation was performed on the entire Larval transcriptome, but for the All transcriptome, Trinotate annotation was performed on a subset of 33,241 contigs for which there were adequate mapped reads in the juvenile experiment (see below). For both transcriptomes, functional diversity of resultant annotations was further summarized by semantic similarity grouping of unique gene ontology (GO) terms using the software REVIGO (Supek et al., Reference Supek, Bosnjak, Skunca and Smuc2011) with January 2017 Gene Ontology version and 15 March 2017 UniProt-to-GO mapping file. REVIGO default parameters used were: allowed similarity of medium (0.7), UniProt database as a reference and SimRel semantic similarity measure. The transcriptomes are available on NCBI GenBank Transcriptome Shotgun Assembly (TSA) Database under accession ID PRJNA524065 (Larval) and accession ID PRJNA525256 (All).
Differential expression analysis
To quantify differences in gene expression among samples, libraries from individual specimens were mapped to the de novo transcriptome using Bowtie2 (V 2.0.6) and counted using eXpress (V 1.5.1), a read mapper that probabilistically assigns reads that ambiguously map to multiple loci, thus minimizing issues arising from redundant putative transcripts in de novo assembled transcriptomes (Roberts & Pachter, Reference Roberts and Pachter2013). Larval samples were mapped to the Larval transcriptome and adult and juvenile crab samples were mapped to the All transcriptome. Read mapping rates averaged 81.1 ± 5.7% (mean ± 1SD) across samples, though mapping rates of larval samples (83.0 ± 4.2%) were about 6% higher than other samples (77.2 ± 6.4%). Those mapping rates indicate high quality transcriptomes (Conesa et al., Reference Conesa, Madrigal, Tarazona, Gomez-Cabrero, Cervera, McPherson, Szcześniak, Gaffney, Elo, Zhang and Mortazavi2016). Library size (# of mapped reads) was compared visually using boxplots (Figure S2). Any library that was small enough that the box did not overlap with other boxes was removed from analysis. Transcripts had to have expression of at least 2 reads/million mapped in 4 samples to be included in the analysis.
Statistical analysis of differential gene expression was performed in R (V 3.5.1) using the Bioconductor package EdgeR (V 1.2.4) (Robinson et al., Reference Robinson, McCarthy and Smyth2010). Multidimensional scaling of filtered count data using the function plotMDS was used to visualize overall differences among gene expression from each treatment. Within each group of samples (juvenile, adult, larval) we fit the data to a fully-factorial negative-binomial generalized linear model (GLM) with main and interaction effects for all factors (e.g. temperature, pH, age, tissue, female), and estimated common, trended and tagwise dispersions, as recommended in the edgeR manual. Likelihood ratio tests in edgeR were conducted to identify differentially expressed transcripts in datasets to examine the main factors in the study: larval origin, pH, developmental time, adult tissue and pH, and juvenile pH and temperature (Table 1). A subset of all possible likelihood ratio tests was made to highlight the variation in response to the main factors of the central hypotheses tested in this study (Table 1). We used likelihood ratio tests with false-discovery rate correction of P < 0.05 and a fold change of ≥2 to identify differentially expressed transcripts. For each dataset of samples, we generated a set of unique transcripts among the likelihood ratio tests (Table 1). Those unique transcripts were used to generate a data set of median centred log2 transformed FPKM data. For all of the datasets, the differentially expressed transcripts were expression filtered and clustered using the Software Cluster (V 3.0) (http://bonsai.hgc.jp/~mdehoon/software/cluster/software.htm) to remove transcripts for which the expression difference among groups was lower than 4-fold (maximumFPKM – minimumFPKM ≤2.0). Transcripts were clustered into a unique minimal set of clusters (i.e. the minimum number of clusters to organize the diversity of expression profiles across samples) using k-means clustering on transcripts starting at N = 10 clusters, and manually reducing the number of clusters if two or more clusters had the same expression profiles. The resulting clusters were visualized in Java Treeview (V 1.1.6) (http://jtreeview.sourceforge.net). Graphical presentation of data was generated in Treeview or R. All scripts used in this study are available on GitHub (https://github.com/safay/RNA_seq).
Table 1. Pairwise statistical comparisons (likelihood ratio tests in edgeR) used in the red king crab study to differentially expressed transcripts of interest
The asterisk (*) indicates sets of highly differentially expressed transcripts that were further analysed.
a Relaxing the false discovery rate correction for heart analyses did not dramatically increase the number of differentially expressed genes. The number of transcripts changed from 2 at FDR = 0.01 to 4 at FDR = 0.1.
b The juvenile interaction set was further expression filtered at a minimum fold change of 8 across samples (only transcripts with maxFPKM-minFPKM ≥3 were retained) to yield a more manageable sized data set of 1742 transcripts that had the strongest expression differences.
Results
Red king crab transcriptome
The two assembled transcriptomes had different sizes and numbers of annotated contigs. The Larval transcriptome had 145,780 contigs of which 18,049 (12.4%) were functionally annotated by BLAST, and 17,291 (11.9%) were functionally annotated with gene ontology (GO) terms (Table S2). There were a total of 138,777 GO terms in the Larval transcriptome, of which 8356 were unique (Table S3). Summarization of those unique GO terms by REVIGO yielded 5288 terms for the Biological Process ontology, 906 terms for the Cellular Component ontology and 2162 terms for the Molecular Function ontology. The N50 of the Larval transcriptome was 1324 bp, and N90 was 6204 bp. The All transcriptome had 242,779 contigs (Table S2). The N50 of the All transcriptome was 2424 bp, and N90 was 9460 bp. A subset of the All transcriptome containing 33,241 contigs for which there were significant expression levels across samples was used for annotation. In that subset, 23,132 (69.9%) were functionally annotated, and transcripts had a N50 of 5302 bp and an N90 of 13,117 bp. The All transcriptome was annotated with a total of 173,920 GO terms, 5244 of which were unique (Table S3). Summarization of those 5244 GO terms by REVIGO yielded 3315 terms for the Biological Process ontology, 1270 terms for the Molecular Function ontology, and 659 terms for the cellular component ontology. The two transcriptomes had the same semantic similarity supergroup as the most abundant term in biological process and cellular compartment ontologies (Renal phosphate ion absorption and Slit-Robo signalling complex, respectively; Figure S3, Table S3). After that, there was not strong agreement in the REVIGO summarization of the functional annotation diversity in the two transcriptomes. In sum, the All transcriptome had a larger number of annotated contigs, but diversity of annotation was higher in the Larval transcriptome, and GO term analysis suggests these transcriptomes overlap in only a few of the largest super annotations within each ontology (Table S3, Figure S3).
Adult responses to pH
Read mapping of adult samples to the All transcriptome resulted in 35,957 (14.8%) contigs that passed read mapping filtering (2 reads/million mapped in ≥4 samples). Examination of tissue and pH in the same analysis produced a large set of differentially expressed transcripts (N = 17,905, Table 1), of which 7267 passed expression variation filtering (≥4-fold expression differences among groups). Not surprisingly, tissue type had a much stronger influence on the number of differentially expressed transcripts than did acclimation pH (Table 1). Comparisons between different tissues yielded from 11–21 K differentially expressed transcripts (Table 1), except between anterior and posterior gills there were only 26 differentially expressed transcripts (10 of which passed expression variation filtering), so transcriptome data from anterior and posterior gills were pooled for subsequent tissue-specific comparisons (see below). To maximize the potential to identify differentially expressed transcripts across pH treatments, the data were reanalysed by individual tissue (Table 1). In gill tissue, there were 1659 differentially expressed transcripts (629 of which passed expression variation filtering, Table S4). The other tissues had far fewer differentially expressed transcripts: there were only 46 differentially expressed transcripts in cuticle (38 of which passed expression variation filtering and only 3 of which were functionally annotated), and two differentially expressed transcripts in heart tissue, neither of which passed expression variation filtering (Table 1).
Gill tissue differential expression profiles were categorized into six clusters, among which two expression patterns were evident (Figure 2). In one pattern, expression was lower at pH 8.1 and 7.5 than at pH 7.8, though the degree of difference varied among clusters (Figure 2, clusters 1–4). In the other pattern, expression was higher at pH 8.1 and 7.5 than at pH 7.8 (Figure 2, clusters 5 & 6). In no cases were there monotonic differences in gene expression observed across acclimation pH (i.e. positive or negative correlation of expression with pH). The largest gene expression differences were less than 8-fold, with the largest differences being between pH 7.5 and pH 7.8 in clusters 1 and 5 (Figure 2). In general, gene expression levels were similar or the same between pH 8.0 and pH 7.5 (Figure 2). There were differences among clusters in terms of cluster size and annotations (Table 2). While there were some overlapping functional annotations among clusters with similar expression profiles, it is clear that each cluster is largely functionally distinct. The most abundantly represented biological process term summary differed across all clusters (Table 2). Typically, annotations for the molecular function ontology occurring in more than one cluster were only abundantly represented in one cluster. In two cases an abundantly represented molecular function term summary was the highest represented in more than one cluster: metallocarboxypeptidase activity (clusters 1 & 2), and beta-N-acetylhexosaminidase activity (clusters 5 & 6; Table 2). Only two annotations were observed in more than two clusters: heme binding (in 4 clusters) and structural molecule activity (in 3 clusters; Table 2).
Fig. 2. Gill differential expression in adult red king crabs acclimated to different pH levels. Each point represents the mean ± 95% confidence interval. Grey bars atop each plot are cluster IDs. Plotted from data in Table S4.
Table 2. Cluster annotations using REVIGO for adult red king crab gill tissue response to pH
See Table S4 for complete annotation information. Terms are ordered by frequency of appearance; number in parentheses indicates the number of GO terms in each semantic summary. Bold terms are those that appear in more than one cluster, italicized terms are similar functional groups that appear in more than one cluster. * indicates annotation that appears in more than two clusters.
The three functionally annotated cuticle proteins that were differentially expressed (Table 1) were two contigs encoding cuticle protein AM1999a, and one contig encoding peritrophin-1, all of which were strongly repressed in crabs at pH 7.5 relative to crabs held at pH 7.8 and 8.0.
Of the differentially expressed transcripts between anterior and posterior gill (Table 1), three were annotated. Two of the transcripts had functional homology to homeobox protein abdominal-A, which were expressed at about 16-fold higher in posterior gills. The third transcript had functional homology to a whey acidic protein, which was about 128-fold higher expressed in anterior gills. The remaining transcripts that differed between anterior and posterior gills were all expressed at a higher level in anterior gills. None of the transcripts that differed between anterior and posterior gills were differentially expressed with respect to pH.
Larval responses to pH
Read mapping of larval samples to the Larval transcriptome resulted in 18,669 contigs (12.8%) that passed read mapping filtering (2 reads/million mapped in ≥4 samples). There were very few genes with different expression levels in day 0 larvae from females held across pH levels. Only N = 49 transcripts were identified as differentially expressed, of which N = 16 passed expression variation filtering (Table 1). Just two of those transcripts were functionally annotated (lectoxin and cytochrome P450) (Table S5).
There were large differences in transcriptome profiles between day 0 and day 7 larvae under the same pH conditions (i.e. mother brooded embryos at ambient pH, larvae held at ambient pH) that were consistent across broods (Table 1, Figure 3). Of the 4147 differentially expressed transcripts, 2478 passed expression filtering, of which 780 were functionally annotated (Table 1, Table S5). Female 1572 yielded larvae with differing expression profiles in clusters 1 & 4 (Figure 3), but otherwise there were no strong differences among females. These transcriptional differences are thus reflective of developmental shifts that occur in ambient conditions during the first week of development. Transcripts up-regulated weakly from day 0 to day 7 (Figure 3, cluster 1) included 969 transcripts, 249 of which were annotated. Annotations matched a diverse set of genes with functions ranging from cuticle, carbohydrate metabolism, ion transport, ribosomal and mitochondrial function, and others (Table S5). Transcripts down-regulated weakly from day 0 to day 7 (Figure 3, cluster 2) included 894 transcripts, 442 of which were annotated, and similar to cluster 1 encoded proteins for a diverse array of cellular functions (Table S5). Transcripts up-regulated strongly from day 0 to day 7 (Figure 3, cluster 3) included 421 transcripts, 55 of which were annotated and principally had functions in cuticle (Table S5). Transcripts down-regulated strongly from day 0 to day 7 (Figure 3, cluster 4) included 192 transcripts, of which 34 were annotated and mainly had cuticular functions (Table S5).
Fig. 3. Larval red king crab expression profiles differed between day 0 and day 7 in three females held at ambient pH (8.1). Each point represents the mean ± 95% confidence interval. Grey bars atop each plot are cluster IDs. Plotted from data in Table S5.
At day 7 there were 942 transcripts that differed in expression across broods, but little of that variation was related to exposure pH during the first 7 days of larval development (Figure 4). The differences in gene expression resolved into 6 clusters of differing sizes. Overall, the largest expression fold-differences among clusters were about 40-fold (i.e. female 1572 at pH 8.0 in clusters 1 vs 5), but most expression differences were smaller (Figure 4). Cluster 1, in which there were different responses to pH across broods, had 17 transcripts, only one of which was annotated (Cuticle protein AM1239; Table S6). Cluster 2 had pH sensitivity in one brood (from female 1572) whereby expression was lower at pH 8.0 than the acidified pHs; however expression levels in larvae from female 1594 were generally higher, though not pH sensitive (Figure 4). There were 201 transcripts in cluster 2, of which 51 were annotated to functions likely involved in cuticular or exoskeletal processes (e.g. chitinase, carbohydrate sulphotransferase, chitooligosaccharidolytic beta-N-acetylglucosaminidase) as well as protein and sugar modification (Table S6). Cluster 3, in which expression was slightly lower in larvae from female 1594, and slightly pH sensitive from pH 8 to 7.8 in larvae from female 1572 (Figure 4), had 390 transcripts, 133 of which were annotated. Functional annotations of cluster 3 transcripts include extracellular matrix and cuticular processes, carbohydrate and protein turnover, mitochondrial processes and gene expression regulation (Table S6). Relative to other clusters, cluster 4 had functionally identical expression levels among all sample groups with no pH sensitivity, though expression levels were higher enough from female 1594 to identify these transcripts as differentially expressed (Figure 4). Protein spaetzle was the most abundant transcript in cluster 4 (Table S6). Cluster 5, which showed the highest degree of pH-sensitivity among larval broods, had 110 transcripts, 48 of which were annotated. All but one of the transcripts in cluster 5 were cuticle proteins and pro-resilins (Table S6). Finally, cluster 6, which like cluster 1 had differences among broods with little pH sensitivity, contained 23 transcripts, 5 of which were annotated functionally for extracellular processes, and heme-binding (Table S6).
Fig. 4. Larval red king crab expression profiles after 7 days of exposure to pH levels differed across maternal origin (female) more strongly than across pH levels. Each point represents the mean ± 95% confidence interval. Grey bars atop each plot are cluster IDs. Plotted from data in Table S6.
Juvenile responses to pH and temperature interaction
Read mapping of juvenile crab samples to the All transcriptome resulted in 32,955 contigs (13.6%) that passed filtering (2 reads/million mapped in ≥4 samples). Overall characterization of transcriptomes by the first two dimensions on a multidimensional scaling plot indicated that there was a strong effect of acclimation temperature, and strong differences between expression at the lowest pH and the two higher pH (Figure 5). For each acclimation temperature, samples from pH 8 and 7.8 were similar to each other, but different from samples at pH 7.5 (Figure 5). Samples held at pH 8 and 7.8 were non-overlapping among the different acclimation temperatures, although the samples held at the intermediate and warmest temperatures were more similar to each other than those at the coolest temperatures. Samples held at the lowest pH were the most different from samples at higher pH at the coldest acclimation temperature, and the most similar at the highest acclimation temperature (Figure 5). That result suggests that elevated temperature causes an overall shift in the first dimension of gene expression that is most similar to extremely low pH.
Fig. 5. Multidimensional scaling plot of filtered count data for juvenile red king crab (Figure S2). Replicate samples generally grouped together and are covered by circles. Values inside of circles indicate pH level, and temperature level is given next to the circles. Colours also indicate increasing temperature (blue (coolest), green, red (warmest)) and decreasing pH (light to dark shades of each colour). In general, each temperature had a distinct graphical placement but the placement of pH 8.0 and 7.8 samples within a temperature were generally similar. Arrows indicate the movement of multidimensional data from pH 8 and 7.8 to pH 7.5 for each treatment temperature.
There were 2673 transcripts identified as differentially expressed across pH levels, 1594 of which passed filtering, and 286 of the filtered transcripts were functionally annotated (Table 1). There were 3039 transcripts that were identified as differentially expressed across temperatures, 1765 of which passed filtering, and 442 of which had some functional annotation (Table 1). There were 7512 transcripts identified as differentially expressed by analyses that identified interactive effects of temperature and pH (i.e. transcripts in which response to pH was temperature dependent) (Table 1). Filtration of the differentially expressed transcripts from interactive effects by at least one sample having a value of ≥3.0 yielded a set of N = 1742 transcripts, N = 321 of which were annotated. It was clear that there were N = 6 distinct clusters of gene expression patterns made up of N = 772 transcripts (Figure 6; Table S7), and one large group of transcripts that, although identified as statistically differentially expressed, had an order of magnitude smaller difference among groups (Figure S4). The filtration stringency for that large group of transcripts was relaxed for re-clustering (Table S8), but that clustering did not yield any additional obvious clusters or patterns in the data (i.e. strongly different expression trends with temperature or pH). Thus, those data were not further analysed.
Fig. 6. Juvenile red king crab clustered expression data. See Table S7 for cluster data. Temperature and pH values are represented as categories in this plot, actual pH and temperature data are in Figure 1. Symbols represent mean ± 1 standard deviation. Error bars for 95% confidence intervals around each point are smaller than the point symbol in every case. Grey bars atop each plot are cluster IDs. Plotted from data in Table S7.
Of the N = 6 clusters, two clusters (1 and 6) were strongly temperature responsive, one cluster (2) was strongly pH responsive, two clusters (4 and 5) demonstrated a temperature*pH interaction, and one cluster (3) had a pattern of gene expression that was not strongly responsive to either factor in the study (Figure 6). In the temperature responsive clusters, expression was graded from low expression in the cold acclimation treatment to high expression in the warm acclimation treatment (Figure 6, clusters 1 and 6). In the pH responsive cluster, expression levels were similar across acclimation temperatures, and increased with decreasing pH (Figure 6, cluster 2). In the temperature*pH interaction clusters expression levels were more responsive in samples from low and middle acclimation temperatures, in which expression was at a relatively lower level at high and intermediate pH, but expression increased dramatically at low pH. In contrast, in those interaction clusters, expression was relatively high for the warm acclimation samples at high pH, but expression decreased at low pH (Figure 6, clusters 4 and 5).
All of the clusters with environmentally responsive gene expression patterns (i.e. all but cluster 3) had relatively low levels of annotation at 25–35% (Table S7). Cluster 3 had much lower annotation of 9%, nearly all of which were identified as unknown proteins. Overall, 76% of the annotated differentially expressed transcripts were identified as functionally involved in formation and regulation of the cuticle.
The most common functional annotation across all clusters was for cuticular matrix proteins (Tables 3 & 4, Table S9). These are proteins that form the physical matrix of the cuticle upon which deposition of chitin and calcification takes place. Annotations involved with cuticle proteins were 84% of the annotated transcripts in pH-responsive cluster 2 (Table 4), and were more abundant in the interaction clusters 4 and 5 at 36% and 56% of annotated transcripts, respectively (Table 4). Clusters 1 and 6 had 22% and 17%, respectively, of the annotated transcripts as cuticle proteins (Table 4). There were a total of 23 different cuticle protein transcripts identified. The most abundant, cuticle protein CP1158, was only observed in cluster 2 (Table 4). The greatest diversity of cuticle proteins was observed in cluster 5, which contained 15 different cuticle protein transcripts (Table 4).
Table 3. Summary of the juvenile red king crab gene expression data for each of the clusters from Figure 6
Each transcript was assigned to only one of the functional annotation categories. See Tables S7 and S9 for complete annotation detail.
Table 4. Cuticle matrix proteins that were found to be differentially expressed in juvenile red king crab
Annotations by BlastP, BlastX, Pfam and EggNog (Tables S7 and S9).
The next most common functional annotations were for cuticular proteins that regulate cuticular elasticity, chitin binding and mineralization (Table 5). Those were more highly represented in the temperature-sensitive clusters 1 and 6, and made up 51% and 57% of the annotated transcripts in clusters 1 and 6, respectively (Table 5). Clusters 4 and 5 also contained cuticular regulation transcripts, but to a lesser extent at 21% and 9% of the annotated transcripts, respectively (Table 5). In contrast to the diversity of cuticular matrix proteins observed, there was low diversity in annotation of transcripts involved with cuticular regulation (Table 5, Table S6). All but two of the transcripts included the annotation of pro-resilin, though some transcripts also contained other annotations including chitin_bind_4, larval cuticle protein A3A, and salivary glue protein Sgs-3 (Table 5). Pro-resilin suggests multiple potential functions, including modulation of cuticle elasticity, calcium ion binding, chitin binding (Table S9). The two transcripts that had calcium ion binding as a functional annotation but were not a pro-resilin were otopetrin-1 and chondroitin proteoglycan 2 (Table 5).
Table 5. Proteins involved in cuticle plasticity, chitin and calcium binding that were differentially expressed in juvenile red king crab
Annotations by BlastP, BlastX, Pfam or EggNog included ‘calcium ion binding’ in each of these transcripts (Tables S7 and S9).
a These transcripts are often also annotated as pro-resilin.
The next most abundant classes of functional annotations were for protein regulation, including proteases, protease inhibitors and protein modification, the latter of which were only observed in the temperature-response cluster (Table 3). Transcripts with a functional annotation for immune function appeared in the temperature*pH interaction clusters only. Surprisingly, there was only one transcript with a functional annotation suggesting it could be involved with metabolism, a ferric-chelate reductase 1 (Table S9). Also surprisingly, there was only one differentially expressed transcript with a functional annotation that potentially involves ionic or pH regulation, a sodium-driven chloride-bicarbonate transporter (Table S9).
Discussion
In this study we present an initial analysis of the transcriptomic response of multiple life history stages of red king crab (Paralithodes camtschaticus) to environmental changes, namely pH/pCO2 and temperature (in juveniles). Transcriptomic responses of crabs to shifts in pH/pCO2 and temperature have been previously demonstrated to be informative in gill tissue of the boreal spider crab Hyas araneus following exposure to reductions in pH similar to our study (Harms et al., Reference Harms, Frickenhaus, Schiffer, Mark, Storch, Held, Pörtner and Lucassen2014). We present the first transcriptome assemblies for red king crab, and an analysis of shifts in transcriptome-wide gene expression that occur across pH levels in adult tissues, newly hatched larvae, one-week old larvae, and in juveniles where interactive effects of temperature were examined alongside pH. We chose to focus on transcripts with expression differences across samples of 4-fold or greater in order to present the most robust gene expression differences; in other words, a conservative analysis that is likely to minimize false positives. The main findings from our study are that responses to pH are observed at the level of gene expression shifts in adult gill tissue, but not in heart or cuticle, that there are large transcriptional shifts with larval age that are pH-independent, that there are large effects of maternal origin on transcriptome profiles in larvae, and that juvenile gene expression responses to elevated temperature are stronger than responses to reduced pH, except when pH levels are extremely low (pH 7.5). While the mechanistic conclusions drawn from this study are only as good as the transcriptome assemblies and annotations of transcripts, variation in gene expression patterns among treatments help to illuminate environmental and biological factors to which the transcriptome is sensitive in red king crab. When significantly improved genomics resources are available for red king crab, re-mapping the RNA-seq libraries generated in this study may yield improved inferences about the functional genomic responses of red king crab to shifts in the pH and thermal environment.
Adults
The largest transcriptomic differences observed in adult tissues were among tissues (Table 1). While not surprising that there were tissue-specific transcripts expressed, this observation is important to highlight as a positive control demonstrating that our transcriptome is complete enough to find large-scale differences where they should exist. The differences across tissues were at least an order of magnitude larger than the largest differences across pH seen in gill tissue (the most responsive tissue to pH) (Table 1). We did not report a systematic analysis of the tissue-specific gene expression profiles as that was outside the scope of the pH-sensitivity hypotheses tested in this study. But, such an analysis, as well as construction of tissue-specific transcriptomes, is a logical extension of what could be accomplished with the data collected in this study.
We expected to observe a large number of differentially expressed transcripts across pH levels in adult cuticle tissue, because calcification and the ability to calcify the new exoskeleton has been shown to be compromised at pH 7.7 (Long et al., Reference Long, Swiney and Foy2013a). However, we did not observe a significant response of cuticle at the transcriptome level, as there were only 3 transcripts that changed between pH 8.1 and 7.8 (Table 1). It is possible that there is a threshold effect on cuticle response to low pH whereas cuticle in pH 7.7-exposed adults has a transcriptomic response more like what we observed at pH 7.5 (i.e. more differentially expressed transcripts at extreme low pH). Importantly, the overall transcriptomic response of cuticle to reduced pH was small, suggesting that this tissue may not have a plastic physiology that can compensate for environmental shifts through changing the expression of genes.
A potential confounding factor in gene expression for any tissue, but especially cuticle, is variation in moult-stage across the females in the study. All crabs were in intermoult, which for red king crab is a period from 3 months (Dvoretsky & Dvoretsky, Reference Dvoretsky and Dvoretsky2012) to over 9 months (Stevens, Reference Stevens2009). However, when during the intermoult our crabs were sampled is unknown as we did not stage crabs more specifically (e.g. by examination of mouthpart setae morphology; Stevens, Reference Stevens2009), nor would we have had enough replicates to assess effects of moult stage with only N = 3 specimens per treatment, N = 9 total.
Heart tissue had only two differentially expressed transcripts between the most extreme pH levels, both of which were only mildly differentially expressed (more than 2-fold but less than 4-fold) and neither of which were annotated. While there are no studies to address functional differences in cardiac activity with pH in red king crab, boreal spider crabs are known to shift heart rate with pH as adults (Walther et al., Reference Walther, Sartoris, Bock and Portner2009). Reduced heart rate has also been observed in embryonic and larval temperate intertidal zone porcelain crabs at reduced pH (Carter et al., Reference Carter, Ceballos-Osuna, Miller and Stillman2013; Ceballos-Osuna et al., Reference Ceballos-Osuna, Carter, Miller and Stillman2013). Crab cardiac transcriptomes have been shown to be highly responsive to shifts in temperature (Stillman & Tagmount, Reference Stillman and Tagmount2009; Ronges et al., Reference Ronges, Walsh, Sinclair and Stillman2012), so our results suggest that if there are changes in cardiac performance with pH in red king crab, those changes are likely regulated post-transcriptionally.
Anterior and posterior gill tissue transcriptomes did not differ to a large extent, suggesting that there is no functional differentiation among gills in red king crab, as has been shown in other crabs (Havird et al., Reference Havird, Mitchell, Henry and Santos2016). The ability to detect more differentially expressed transcripts across pH levels in gill, compared with other tissues, may have been because after pooling anterior and posterior gills, the number of replicates in pairwise comparisons was effectively larger increasing the power of the statistical analyses to detect differences between sample groups. In other words, each adult animal was represented by two gill samples instead of one. Because of the fact that we employed a stringent expression difference cut-off of 4-fold among samples, it is likely that we pruned transcripts that were identified as differentially expressed because of the increased sample size (i.e. fold-change differences were small and only detectable with a higher-powered analysis facilitated by increased sample size). Of the nearly 1700 differentially expressed transcripts at 2-fold expression difference, about one third were represented by transcripts that differed by at least 4-fold across samples (Table 1), so the pruning was indeed conservative.
The N = 180 differentially expressed transcripts in gill tissue with respect to pH (Table 1) represented a diverse array of cellular functions (Table 2) indicating that cellular shifts in gill tissue in response to pH shifts occurred in a number of different protein, lipid and carbohydrate-based cellular components. While the functional impacts of those transcripts have not been empirically determined, based on the functional annotations, we discuss possible cellular responses to changing environmental pH. Shifts in the expression of transcripts encoding beta-N-acetylhexosaminidase, UDP-glucose 4-epimerase, cuticular structural constituents, and extracellular matrix constituents (Table 2) suggest that pH induced structural differences in gill tissue chitin and cuticle. Shifts in gill cuticle structure could facilitate shifts in ion transport, as has been shown in lobster gill epithelia (Lucu & Towle, Reference Lucu and Towle2010). Shifts in the expression of transcripts encoding RNA polymerase II domains along with transcription factors and DNA binding transcripts indicate shifts in gill tissue nuclear activity (Table 2), suggesting that pH induces the expression of genes that themselves cause shifts in gene regulation. Shifts in the expression of transcripts involved with fatty-acid processes (i.e. long-chain-enoyl-CoA hydratase, long-chain fatty-acid-CoA ligase, beta,beta-carotene 15,15′-dioxygenase) (Table 2) indicates potential shifts in cellular metabolism (e.g. mitochondrial processes) or signalling (i.e. hormonal processes). The shift in steroid hormone receptor activity with pH along with shifts in lipid modification activity (Table 2) could be interpreted as a shift in paracellular signalling in gill tissue. Transcripts encoding metal ion binding (e.g. zinc, heme binding) and metal ion regulated processes (e.g. metallocarboxypeptidase activity) were highly responsive to pH in gill tissue. Shifts in the expression of those genes are not specifically linked to specific biological processes, but do indicate that the fundamental biology of gill tissue cells (e.g. proliferation, metabolism) may have differed in response to pH acclimation.
Larvae
Reduced pH has been suggested to negatively impact the survival of red king crab embryos and larvae (Long et al., Reference Long, Swiney and Foy2013a) as well as embryos and larvae of other crabs (Carter et al., Reference Carter, Ceballos-Osuna, Miller and Stillman2013; Ceballos-Osuna et al., Reference Ceballos-Osuna, Carter, Miller and Stillman2013; Schiffer et al., Reference Schiffer, Harms, Portner, Mark and Storch2014b). Results of transcriptomic analyses presented in this study strongly suggest that wide variation in environmental pH largely does not influence regulation of gene expression in late embryos/newly hatched larvae or in the expression of genes in larvae following 7–8 days of exposure to reduced environmental pH (Table 1, Figures 3 & 4). There were strong differences in gene expression between newly hatched and one-week old larvae (Table 1, Figure 3), and among broods from different females (Table 1, Figure 4). The large differences in the transcriptome between day 0 and day 7 larvae likely reflect shifts in cellular processes important for larval development and activity. The exoskeleton calcification of whole red king crab larvae was insensitive to reduced pH in a 7-day exposure, though 7-day old larvae have cation content that is different from newly hatched larvae (Long et al., Reference Long, Swiney and Foy2013a). Thus, at the organismal and the transcriptomic levels it appears that the red king crab embryonic and larval stages are relatively unable to modify their biology in response to environmental variation in pH, in contrast to what has been observed in H. araneus (Schiffer et al., Reference Schiffer, Harms, Portner, Lucassen, Mark and Storch2013, 2014a, 2014b). The transcriptomic inferences from our study could have been influenced by the handling stress of the ovigerous females during shipment and holding (see Methods), as larval survival through later larval stages was lower in this batch of larvae than typical (K. Swiney pers. obs.). However, we do not have any direct evidence to link the conditions during embryonic development with the lowered larval survival, which could have occurred for many unrelated reasons (e.g. pathogens present later in larval development).
Female-to-female variation in offspring condition has been observed in red king crabs, though causal factors are unknown (Swiney et al., Reference Swiney, Eckert and Kruse2013). Factors that could explain the differences in transcriptome profiles among broods from different females could include female condition at the time of oogenesis (i.e. egg provisioning with specific proteins, mRNAs and regulatory RNAs), embryonic exposure in early development, and influence of maternal and paternal genotype. Our results highlight the importance of increasing numbers of biological replicates in split-brood studies of field-caught animals as the variation in environmental response among broods can be larger than the mean response in a population, as has been previously demonstrated in porcelain crab embryos and larvae (Carter et al., Reference Carter, Ceballos-Osuna, Miller and Stillman2013; Ceballos-Osuna et al., Reference Ceballos-Osuna, Carter, Miller and Stillman2013).
Juveniles
The temperatures used in this study were not close to the thermal limits of juvenile red king crab, where 1-day exposure did not cause any mortality until temperatures above 23°C, and there was no increase in mortality after a 45-day exposure to 15.5°C (Long & Daly, Reference Long and Daly2017). Thus, transcriptomic responses to temperature observed in this study are not likely to reflect extreme cellular stress responses that occur at temperatures near lethal limits. However, energetic responses to temperature were observed to change at temperatures between 10–15°C, including daily food intake increasing from 20–40% body mass between 10–15°C (Long & Daly, Reference Long and Daly2017), and 35% reduction in daily growth rate between 12 and 16°C, despite the food intake increase (Long & Daly, Reference Long and Daly2017). Thus, it is likely that transcriptomic differences observed in this study are reflective of a physiological response that reflects an incomplete compensatory change in organismal energetics. A reduction in carapace length (a proxy for body size) at the first moult was observed between +2 and +4°C at ambient pH, and between ambient temperature and +2°C at pH 7.8 (Swiney et al., Reference Swiney, Long and Foy2017), suggesting that the incomplete compensation in energetics had an effect on integrated growth at the time of moult. Additionally, condition (body mass normalized to carapace size) of juvenile red king crabs was overall reduced by 25% at pH 7.8 relative to pH 8.0 (ambient) (Long et al., Reference Long, Swiney, Harris, Page and Foy2013b), and mass increase (growth) with moulting was negatively correlated with decreasing pH (Long et al., Reference Long, Pruisner, Swiney and Foy2019), suggesting an energetic deficiency resulting from rearing at lowered pH.
In this study we attempted to avoid including specimens in the final pre-moult stage, which is characterized by large physiological shifts in water volume, ion transport and gene expression (Lee & Mykles, Reference Lee and Mykles2006; Chang & Mykles, Reference Chang and Mykles2011; Wittmann et al., Reference Wittmann, Benrabaa, Lopez-Ceron, Chang and Mykles2018). We did not observe shifts in gene expression known to be regulated as the immediate pre-moult occurs (e.g. mTOR pathway genes; Wittmann et al., Reference Wittmann, Benrabaa, Lopez-Ceron, Chang and Mykles2018), even though those genes were represented in our transcriptome (Table S2). The intermoult duration for first instar crabs under ambient pH and pH 7.8 and temperature conditions of ambient, +2 and +4°C used in this study ranged from 28–42 days (times were shorter at increased temperature, but not different across pH conditions) (Swiney et al., Reference Swiney, Long and Foy2017) and survival was near or at 100% through 50 days across all conditions (Swiney et al., Reference Swiney, Long and Foy2017), with 100% of specimens undergoing the first moult except for in the combined effects of pH 7.8 and +4°C (Swiney et al., Reference Swiney, Long and Foy2017). From those data, it is unlikely that our specimens held at ambient pH or pH 7.8 would have been nearing the immediate pre-moult stage prior to sampling, and thus the transcriptomic differences we have observed would be more generally characterizing the intermoult period across different samples, though it is likely that specimens in warmer conditions were closer to moulting (~4–6 days) than specimens in cooler conditions (~16–18 days) at the time of sample preservation.
The largest pH and temperature-dependent transcriptomic responses observed in this study for juveniles were for genes involved in the formation and biomineralization of cuticle (Table 3). Those changes strongly suggest that shifting the types and quantities of those gene products is a central component of compensatory shifts made by juvenile red king crabs to maintain size and major cation content of the exoskeleton across conditions. This inference comes from the fact that prior studies of response of juvenile red king crabs did not find large differences in size, mass, magnesium or calcium content after the first moult (Long et al., Reference Long, Swiney, Harris, Page and Foy2013b).
The expression diversity of transcripts functionally annotated as cuticle proteins (Table 3) suggests that there may be important functional specialization with respect to the cuticular formation or biomineralization process occurring when the animals are in a different energetic state, or when faced with altered environmental conditions. The fact that there was high specificity to the expression profiles of specific cuticle protein transcripts (Table 3) would suggest that there is functional diversification across that group of proteins. While the extent of functional diversification across these proteins is largely not characterized, prior studies have characterized some relevant variation in the cuticle protein transcripts we observed. For example, ACP20 (adult specific cuticle protein 20) is a hard, but flexible, glycine-rich cuticle protein found in adult arthropods that is negatively regulated by juvenile hormone and positively regulated by ecdysteroid (20E) such that it is chiefly up-regulated in pre-moult cuticle (Charles et al., Reference Charles, Bouhin, Quennedey, Courrent and Delachambre1992; Braquart et al., Reference Braquart, Bouhin, Quennedey and Delachambre1996; Lemoine et al., Reference Lemoine, Mathelin, Braquart-Varnier, Everaerts and Delachambre2004). ACP20 was only expressed in clusters 1 and 5 (Table 3), where there was a strong temperature-dependence of expression (Figure 6) suggesting that there may have been specimens nearing immediate pre-moult in our study. Cuticle protein CP1158, which was solely differentially expressed in a pH-sensitive cluster (cluster 2; Figure 6) is known to be down-regulated in pre-moult and post-moult in the blue crab (Portunus pelagicus) (Kuballa et al., Reference Kuballa, Merritt and Elizur2007), but its functional role or why it should be specifically pH sensitive is unknown. Glycoprotein SgAbd-8 is also expressed in a pH sensitive manner, though this protein has previously been shown to be differentially expressed in response to temperature in crab tissues (Ronges et al., Reference Ronges, Walsh, Sinclair and Stillman2012). The present study thus yields further information regarding potential functional differentiation among the cuticular proteins as well as targets for future mechanistic research in compensatory changes in the cuticle structure and function under environmental change.
Conclusions
In this study we have presented the first transcriptome assemblies for red king crab, Paralithodes camtschaticus, and used those assemblies to analyse differential gene expression across the transcriptome in response to acidification and warming in adult, larval and juvenile crabs. We present genomics resources to the scientific community engaged in transcriptomic analyses of king crabs. Overall, our inferences about transcriptome responses to reduced pH and temperature suggest that juvenile red king crabs are more sensitive to warming than acidification, and that responses to acidification at the transcriptomic level occur at different levels of pH across life-stages, with juveniles being less pH sensitive than adults (at least in gill tissue of adults). While inferences regarding the specific biological responses associated with changes in gene expression are likely to change as resources for red king crab genomics enabled studies continue to improve (i.e. better genomic assemblies and annotation), our inferences about general sensitivities to temperature and pH across life stages of red king crab are robust and unlikely to shift.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S002531541900119X.
Financial support
Funding through NOAA NMFS Grant OAPFY12.02.AFSC.003 to RF, and computing resources through XSEDE Allocation DEB140012 to JS. The findings and conclusions in the article are those of the authors and do not necessarily represent the views of the National Marine Fisheries Service, NOAA. Reference to trade names does not imply endorsement by the National Marine Fisheries Service, NOAA.
Author contributions
Designed and conducted live animal experiments: RF, KS. Designed molecular and bioinformatic approach: SF, JS. Conducted molecular benchwork: MA, SF. Bioinformatic Analyses: SF, JS. Manuscript writing: JS. Manuscript editing: RF, KS.
