RNAseq Reveals Sensitive, Concentration-Dependent Transcriptional Markers of Copper in Mytilus californianus Larvae and Adults

Copper contamination of coastal waters is a long-standing problem in many regions. Copper water quality criteria in marine waters is often determined with bivalve embryo-larval toxicity tests which measure survival and normal development as endpoints. Gene expression data is increasingly incorporated into such assays as a complementary and sensitive marker of contaminant exposure or toxicity. Here, we measured the impacts of copper on Mytilus californianus larval transcriptional profiles, and identified sensitive biomarkers of copper exposure. Key functional categories that were identified among these genes include biomineralization/shell formation, metal binding, and development. Finally, we compared the transcriptional response of larvae to that of adult gill tissue, and show in both datasets that patterns of declining transcript expression occur at lower copper concentrations than those required to induce increases in transcript expression, suggesting that down-regulated genes serve as the most sensitive marker of copper exposure.


INTRODUCTION
In recent decades, improvements in coastal water quality have been achieved by regulation of discharges of contaminants from point sources. These regulations have been based in part on toxicity testing of marine species and on documenting the accumulation of contaminants in their tissues (EPA, 1985). Mussels of the genus Mytilus have been regularly used for such testing because they occur in almost every coastal habitat and are filter feeders. As larvae, Mytilus are highly sensitive to their chemical environment, so short-term bivalve embryo-larval toxicity assays are standard to determine the toxicity of metals, effluent, and numerous other chemicals (EPA, 1995;E50 Committee, 2013), and LC50 data arising from laboratory testing of these effects serve as the criteria for setting water quality regulations (Rivera-Duarte et al., 2005). The bay mussel species, Mytilus edulis and M. galloprovincialis, have been most widely used in contaminant testing, but their deployment presents challenges because together with M. trossulus, they are almost impossible to distinguish morphologically and can hybridize with one another (Bierne et al., 2003;Riginos and Cunningham, 2005), making it very difficult to standardize testing protocols without genetic analysis of the animals (Arnold et al., 2009).
Bivalve embryo-larval toxicity assays currently use abnormal development as the most sensitive endpoint, and rarely consider other detrimental effects that may occur at low contaminant concentrations. Over the past two decades, advances in large-scale gene expression analysis have created the potential to provide more sensitive indicators of toxicity, while simultaneously characterizing biochemical and physiological regulation in response to a toxin (Ankley et al., 2006;Calzolai et al., 2007;Schirmer et al., 2010;Hahn, 2011). Transcriptional characterization may also reveal unexpected pathways involved in the toxicity of a chemical (Poynton et al., 2007). Ultimately from these datasets, biomarker genes can be identified that are highly correlated with toxin exposure (Nordberg, 2010), or with some negative outcome at the whole-organism or population level (Ankley et al., 2010;Connon et al., 2010).
Copper is a common environmental pollutant in coastal waters globally (Sadiq, 1992), but Southern California has been a hotspot for copper contamination (Schiff et al., 2003(Schiff et al., , 2007Rivera-Duarte et al., 2005). Mytilus mussels are a key component of ecosystems in this region, and thus provide an environmentally relevant model for copper research. The effects of copper on larval Mytilus spp. survival and development are relatively well studied (Johnson, 1988;Hoare et al., 1995a,b;Arnold et al., 2009), yet there has been little investigation into the biochemical effects of copper on this sensitive stage, nor on the potential molecular drivers of developmental abnormality. All investigations of Mytilus spp. global transcriptome responses to toxins have been conducted in adults (Venier et al., 2006;Negri et al., 2013;Varotto et al., 2013;Xu et al., 2016), so we know little about how early life history stages respond to toxins at the transcriptional level, especially at low, sublethal concentrations.
Traditional endpoints of toxicity assays (i.e., developmental abnormality) are often measured in concentration-response experiments, which allow for modeling of the response to a toxin and for consistent determination of biologically relevant concentration thresholds, such as the lowest observed effect concentration (Walker et al., 2000). The transcriptional response to a drug or toxin can also be analyzed in a dose-responsive fashion (Daston, 2008;Ji et al., 2009Ji et al., , 2011. Concentration-responsive transcriptional analysis has been used to some extent in ecotoxicological investigations (Poynton et al., 2008;Brulle et al., 2010;Whitehead et al., 2010), yet most transcriptomic studies assess only a few concentrations, if more than one is even used. Transcriptional dose-response or concentration-response experiments allow for a more complete characterization of the discrete physiological responses to different levels of a toxin . More specifically, concentration-specific biomarkers can be identified, which could provide a highresolution metric of toxin-exposure level in environmental monitoring. Biomarkers that respond to low concentrations of a contaminant can be especially useful for identifying sublethal physiological changes that may not be evident at the whole-organism level.
This work presents an analysis of concentration-responsive transcription in Mytilus californianus larvae and adults after short-term (48 h or 24 h, respectively) exposures to copper. This species was chosen because M. californianus is morphologically distinct from bay mussels, possesses a ribbed shell, and cannot hybridize with other Mytilus species. Reports indicate that incorrect identification of bay mussel species can confound efforts to standardize regulatory criteria and have called for data collected in M. californianus to be added to toxicity databases (Arnold et al., 2009). We aimed to identify low-concentration transcriptional biomarkers of short-term copper exposure, and to characterize functional pathways associated with a range of copper exposure concentrations. We also compared the concentration-dependent copper response of adults and larvae, and identified transcriptional biomarkers of copper exposure and toxicity. The biomarkers identified in this study represent robust indicators of copper exposure that could be incorporated into environmental monitoring and toxicity testing.

Preparation of Mussel Larvae
Adult M. californianus were collected from Santa Monica, California and transported to the Wrigley Marine Science Center (WMSC) on Catalina Island, where they were held in a subtidal cage hanging from the dock for 1.5 years prior to spawning. Although M. californianus does not appear to exhibit significant population differences throughout its range (Addison et al., 2008), we further sought to avoid potential population-level sources of genetic variation by collecting all the mussels from the same site in Santa Monica. Animals were collected for two trial spawns during May-June of 2013. During both experiments, 15-20 animals were removed from the subtidal cage, and placed directly on ice for 1 h. To induce spawning, mussels were then rinsed in fresh DI water, scraped clean of any fouling organisms, and transferred to a tank of filtered seawater at 23-25 • C. Mussels were sometimes fed a small quantity of Shellfish Diet 1800 (Reed Mariculture) to induce spawning if thermal shock alone was not sufficient.
Once an adult mussel started spawning, it was removed from the main tank, rinsed thoroughly with filtered seawater, and transferred to a separate beaker containing 0.2 µm filtered seawater. Each individual that spawned was held alone in a beaker to prevent cross-contamination of gametes. When spawning was complete, the adult mussels were removed from the beakers, and the appearance of eggs and motility of sperm were examined under low power on a compound microscope. Once eggs had transformed from club-shaped to round, sperm from a single male was added to eggs of a single female to reach an average density of ∼5 sperm per egg. Fertilization was generally completed within ∼30 min and the number of fertilized embryos, evidenced by the formation of a polar body and first division of the zygote, was counted.

Embryo-Larval Toxicity Assay
Copper solutions were prepared in 0.2 µm filtered natural seawater (34 parts per thousand salinity) using a 1 mM stock solution of copper sulfate (CuSO 4 ). Each copper exposure experiment consisted of one control container and nine containers with increasing copper concentrations-2, 3.1, 4, 6, 8, 10, 15, 20, and 25 µg/L Cu. Larval exposure concentrations were selected based on EC50 and LC50 values in the literature (e.g., Martin et al., 1981;Arnold et al., 2009), as well as pilot work conducted by our lab. For the purposes of our experiment, a greater number of exposure concentrations was more important than technical replicates of each concentration. The experiment relied on the notion that more concentrations are important to effectively model concentration response of transcripts, as opposed to fewer concentrations with more replicates, which would be favorable for use of ANOVA or differential expression analysis ("regression" design vs. "ANOVA" design, described by Graney, 1993). As our study primarily sought to model transcriptional dose responses, we opted for higher resolution in the concentration, with counts from each larger container serving as replicates. This is similar to the approach taken in Ji et al. (2009), the study that introduced and first applied the Sigmoidal Dose Response Search program (described in more detail below).
After spiking copper into each container, jars were mixed thoroughly and allowed to equilibrate for 1 h. Embryos were transferred into the 1L treatment containers at a density of 10 per mL (total count of 10,000 per container), and incubated at 20 • C with a 16:8 h L:D cycle. After a 48-h incubation, 5 replicate 10 mL samples were taken from each culture and preserved with EtOH for counting and morphological analysis. The 48-h larval exposure time was chosen to match the industry standard for 48-h bivalve toxicity assays (EPA, 1995), which is based on the standard time it takes for bivalve embryos to reach the D-hinge larval stage. Microscope counts of larvae were used to calculate the number of larvae surviving and normally developing in each container. To calculate the proportion of control-normalized survival and normal development, each estimated number was divided by the mean estimated number of larvae surviving in the no-copper control. These normalized proportions were analyzed in the r package 'drc' (Ritz et al., 2015), with replicate counts from each concentration serving as replicates for the model. Four-parameter log-logistic curves were fitted to each dataset to calculate 50% lethal concentration (LC50) and 50% normal development effective concentration (EC50) values. The remainder of each culture was collected on a 20 µm sieve, concentrated by centrifugation at 2,000 × g for 4 min, and larval pellets were frozen at −80 • C for subsequent RNA extraction.

RNA Sequencing of Larval Samples
A total of 20 samples were processed for RNAseq (9 copper concentrations and 1 control seawater sample × 2 copper exposure experiments). Total RNA was extracted using Trizol (Thermo Fisher) in MaxTract tubes (Qiagen). RNAseq libraries were prepared for next generation sequencing using a modified version of the Illumina TruSeq protocol, and the 20 libraries were multiplexed and sequenced across two lanes on the Illumina HiSeq 2500 platform as 50 bp single end reads. In addition to the libraries prepared for experimental samples, two longer libraries were prepared for transcriptome assembly. The first was prepared using larvae exposed to 8 µg/L Cu for 48 h, and was sequenced on one lane of an Illumina HiSeq 2500 as 150 bp single end reads. The second library was prepared from larvae reared under control conditions for 48 h, processed using the NEBNext Ultra Directional RNA Library Prep Kit for Illumina, and run over 10% of an Illumina NextSeq lane as 75 bp paired-end reads. The resulting raw sequence data is deposited in the Sequence Read Archive (SRA) database on NCBI under BioProject ID PRJNA639354.
Peptides were predicted for the resulting assembly (Supplementary File S1), and the contigs and predicted peptides were annotated using Trinotate v 2.0.2 with default parameters (Haas et al., 2013). Full taxonomic paths for hits from the nt NCBI BLAST database 1 were attained by parsing the NCBI taxonomic files. UniProt annotations from Trinotate were filtered to only include hits with E-values lower than 1e-15, and hits to the nr NCBI BLAST database (see text footnote 1) were filtered to only include best hits with E-values lower than 1e-20. All annotations were searched and only those containing the term 'Metazoan' in the taxonomic path were retained. Contigs with a 'Metazoan' annotation from either the UniProt or nt database were retained. The final assembly consisted of 23,025 contigs with an average length of 1905.38 bp (Supplementary File S2).

Larval RNAseq Mapping
For mapping, RNAseq reads were quality trimmed and contaminating adapter sequence was removed using Trimmomatic v0.33 (Bolger et al., 2014) with default parameter settings. Reads arising from mitochondria were removed by mapping to the M. californianus mitochondrial genome using BBMap v34 (minid = 0.95 ambiguous = all sssr = 1.0) (Bushnell, n.d.). The remaining reads were mapped to our M. californianus transcriptome assembly detailed above with BBMap v34 (minid = 0.95, ambiguous = all, sssr = 1.0, nhtag = t). Read counts were obtained using FeatureCounts v1.5.0-p3 (Liao et al., 2014), counting on the feature level and only allowing for uniquely mapping reads to be counted (-T 8 -O -F). The two samples from the highest copper concentrations (25 µg/L) were removed at this point, as the majority of read count values were 0, and the samples had consisted of very few living larvae. The resulting count table was further processed in EdgeR (Robinson et al., 2010;McCarthy et al., 2012) to generate normalized trimmed mean of M-values (TMM) counts. Methods for larval transcriptome assembly and mapping are summarized in Supplementary Figure S1.

Adult Mussel Copper Exposure and Microarray-Based Gene Expression Analysis
Similar experiments were conducted with adult animals to compare adult and larval transcriptional responses and look at commonalities and differences across these life history stages. For adults, copper exposure was performed at 15 • C in individual acid-washed polycarbonate containers containing 1L of seawater that was spiked with either 3, 6, 9, 15, 30, 60, 90 or 120 µg/L copper by dilution of a CuSO 4 stock solution. Concentrations were selected to span a range where no toxicity was expected for adult animals (3 ppb is lower than dissolved copper water quality criteria) through to near lethal concentrations for the adult animals (determined based on pilot work by our lab). The concentration range for adult exposure was different than for larvae because for animals of each life history stage, we wanted to elicit a full range of transcriptional responses, from nontoxic exposure through a toxic response that was associated with mortality. Because adult mussels have a much higher tolerance for copper than larval mussels, it was necessary to expose adult animals to a wider range of copper concentrations to elicit a similar response. Importantly, there was overlap in the doses that each life-history stage experienced, so larvae experienced 2, 3.1, 4, 6, 8, 10, 15, 20, and 25 µg/L and adults were exposed to 3, 6, 9, 15, 30, 60, 90 or 120 µg/L. The concentrations were extended to higher levels in the adults to reflect their higher resistance to copper, with the goal of inducing a maximal transcriptional response. Three individual mussels of identical size (∼4 cm long) were placed in each jar with gentle aeration to maintain saturated levels of dissolved oxygen and to promote water movement. Every 4 h the water in each jar was completely replaced with a liter of freshly prepared copper-spiked seawater to compensate for copper draw-down. Two replicate containers of natural seawater without addition of copper served as controls. The exposure was continued for 24 h, and then the mussels were quickly dissected and samples of gill tissue were flash frozen and archived at −80 • C. Adult mussels were exposed for 24 h because there was no standardized incubation time for adult mussel exposures, and out pilot data showed that markers of cell stress (namely HSP70 mRNA induction) were maximally expressed at 24 h of exposure. The 24-h timepoint was also intended to capture an acute response which compares more directly to the goals of the EPA's standard mussel embryo-larval toxicity assay. Gill was analyzed because in adult mussels, the gill and mantle represent the largest tissues, with gill tissue being deployed previously in many physiological studies because of its central role in feeding and respiration and ease of dissection. More importantly, gill tissue is known to be the key site of dissolved metal uptake in aquatic mollusks (reviewed in Marigómez et al., 2002), and has thus been used regularly in studies of molluscan response to metals. Understanding tissue-specific physiological responses is ideal, but in early larvae gills are still developing (Cannuel et al., 2009), so direct tissue comparison would not be possible.
Total RNA was isolated from the 30 samples of gill tissue using Trizol (Invitrogen) and the resulting RNA was purified further across RNAeasy filter columns (Qiagen). Amplified RNA (aRNA) was prepared for each individual gill sample as previously described (Connor and Gracey, 2011). The 30 fluorescently labeled aRNA samples were hybridized to 40 M. californianus arrays (Connor and Gracey, 2011) using an interwoven loop design which yielded a balanced design in which aliquots of each labeled RNA sample were hybridized to either two or four arrays with fluorescent dye reversal. In total, expression data were generated for six replicate mussels sampled under control conditions, and three replicate mussels for each of the 8 concentrations of copper. The resulting expression data were deposited in the ArrayExpress public database as accession number E-MTAB-810. Genes whose expression varied across the series of samples were identified using analysis of variance (ANOVA) in MAANOVA Version 0.98-7 (Kerr et al., 2000) in a model in which dye was treated as a fixed effect and array was treated as a random effect, testing for differences in expression between mussels sampled after different copper concentration exposure. The p-values were adjusted for false discovery rate (Benjamini and Yekutieli, 2001) and a statistical threshold of p < 0.05 defined a list of genes that were statistically significantly differentially expressed across treatments. These differentially expressed genes were analyzed with respect to copper exposure as described below. The microarray probe sequences were aligned to the RNAseq-based transcriptome assembly using a reciprocal BLASTn search so that sequence annotation across the microarray and RNAseq datasets could be compared. Methods for adult Mytilus microarray printing, hybridization and analysis are summarized in Supplementary Figure S1.

Data Analysis
For both larval and adult data, concentration-responsive gene expression was analyzed with Sigmoidal Dose Response Search (SDRS) v0.04 (Ji et al., 2011) to identify genes with sigmoidal expression, and to parameterize the expression dose-response curve of each gene. For genes with significant (p < 0.05) sigmoidal expression, the program calculates an expression EC50, representing the concentration at which the gene's expression was induced (or repressed) by 50% of its maximum expression value (Supplementary Figure S2). The low and high dose for the 95% confidence interval of the expression EC50 are also calculated. The lower threshold of the 95% confidence interval for the EC50 Frontiers in Marine Science | www.frontiersin.org value, which is the lowest concentration of significant induction or suppression of the gene, was considered to be the lowest observed transcriptional effect concentration (LOTEC) for genes (Supplementary Figure S2).
For larval expression data, Weighted Gene Co-expression Network Analysis (WGCNA) (Langfelder and Horvath, 2007) was used to create gene co-expression networks and identify modules of genes with similar expression patterns. WGCNA was run using all genes with a summed TMM across all samples > 300, a module-merging threshold of 0.15, and distinct modules being assigned different colors. Functional enrichment of gene groupings obtained with SDRS and WGCNA were conducted using Gene Ontology (GO) (Ashburner et al., 2000) enrichment analysis in the Cytoscape (Shannon et al., 2003) plugin program, BINGO (Maere et al., 2005). Enrichment was calculated using a hypergeometric test for overrepresentation (p < 0.05), with a Benjamini and Hochberg False Discovery rate correction.

Embryo-Larval Development and Survival
The proportion of larvae surviving (survival) and the proportion of larvae reaching a normal D-hinge shape (normal development) in each copper treatment at the end of the 48-h period exhibited significant sigmoidal concentration-responsive patterns in each trial (Figure 1). Based on the modeled curves, the EC50 (corresponding to the effect of copper on normal development) in exposure trial 1 was calculated to be 6.04 µg/L Cu, and the LC50 (corresponding to the effect of copper on larval survival) was 6.5 µg/L Cu. The larvae in trial 2 exhibited greater sensitivity to copper, with an EC50 of 3.78 µg/L, and an LC50 of 3.96 µg/L. Thus, despite differences in the apparent sensitivity of the different larval cultures, LC50 and EC50 values were relatively similar to one another in a given trial.

The Transcriptional Concentration-Response to Copper in Larvae
In larvae, 1,483 transcripts exhibited sigmoidal expression in at least one trial (Supplementary Figure S3), and 317 transcripts exhibited sigmoidal expression in both trials. Of these common transcripts, 248 were upregulated in response to increasing copper (Figure 2A), and 61 were downregulated ( Figure 2B).
Calculation of expression EC50 and LOTEC values provided a quantitative measurement of each gene's responsiveness to copper allowing genes to be organized with respect to their transcriptional sensitivity to copper. Plotting the expression LOTEC and EC50 values with the corresponding phenotypic data showed that nearly all downregulated genes exhibit an expression LOTEC and EC50 that are lower than the EC50 that was calculated based on the gross phenotypic marker of normal development (Figure 3 and Supplementary Figure S4). For example in Trial 1, the mode LOTEC and mode expression FIGURE 1 | Control-normalized survival (A) and normal development (B) were plotted against copper concentration for trial 1 (blue line) and trial 2 (dashed black line). Data indicate that larvae were slightly more sensitive to copper in trial 2.
EC50 were 3.46 µg/L and 3.29 µg/L, respectively, while the normal development EC50 was 6.04 µg/L Cu, and in Trial 2, the mode LOTEC and mode expression EC50 were of 0.001 µg/L and 2.12 µg/L, respectively, while the normal development EC50 was 3.78 µg/L. In contrast, almost all upregulated genes had expression LOTECs and EC50s that exceeded the normal development EC50 in each respective trial (Figure 3 and Supplementary Figures S4B,D), with a mode LOTEC of 9.63 µg/L and mode expression EC50 of 12.29 µg/L compared to a normal development EC50 of 6.04 µg/L Cu in Trial 1, and a mode LOTEC of 9.63 µg/L and mode expression EC50 of 14.23 µg/L compared to a normal development EC50 of 3.78 µg/L in Trial 2. Thus, overall transcript repression occurs at lower copper concentrations than transcript induction, and in many cases the repression of specific genes can be detected before the deleterious effects of copper exposure are phenotypically evident as abnormal development and mortality. GO enrichment analysis revealed that the genes induced by copper were enriched in functional categories related to the proteasome core complex, proteasome regulation, the ESCRT III complex, ribonucleotide binding, oxidoreductase activity, hydrolase activity, peptidase activity, and copper ion binding. Copper-repressed genes were enriched in just two categories, polysaccharide binding and metabolism.

Modularity of Dose-Responsive Gene Sets
Visual inspection of the data indicated that SDRS did not capture the full-extent of the transcriptional response because a proportion of genes did not show a significant sigmoidal expression response. Clustering co-expressed genes into modules has proven a good way to classify gene sets that respond in the same way to external stimuli, and often share involvement in the same functional pathways (Langfelder and Horvath, 2007). Therefore, WGCNA was used to construct 23 modules of co-expressed genes, of which five modules were significantly correlated (p < 0.05), with the experimental parameter (copper concentration) or the observed phenotypic response (normal development or survival) (Supplementary Figure S5). Three of the modules (colors white, yellow, and steel blue), which collectively comprised 326 transcripts, captured down-regulated expression in response to copper (Supplementary Figure S5). Indeed, 85% (52 transcripts) of the sigmoidally down-regulated genes identified by SDRS were represented in these three modules, indicating that WGCNA effectively captured an expanded set of copper-repressed transcripts. Examination of the expression EC50 and LOTEC values for the SDRS genes associated with each module indicated that the steel blue module contained transcripts that were most sensitive to copper, while transcripts in the yellow module were the least sensitive ( Table 1). The most prominent functional category of down-regulated genes was related to development (Table 1 and  Supplementary Table S1) and included genes involved in cell differentiation and organogenesis (Supplementary Figure S6A), neural and muscular development (Supplementary Figure S6B), and muscle formation (Supplementary Figure S6E). Another functional grouping was genes associated with biomineralization and shell matrix proteins, which were primarily members of the white and yellow modules ( Table 1; Supplementary  Table S1, and Supplementary Figure S6C). Finally, genes coding for divalent metal-containing enzymes were well represented in the white module with genes primarily coding for proteins involved in the uptake or transport of zinc, iron, and molybdenum ( Table 1, Supplementary Table S1, and Supplementary Figure S6D).
Weighted Gene Co-expression Network Analysis also revealed two modules of upregulated transcripts (magenta and saddle brown) that were significantly correlated with copper, and significantly anti-correlated with normal development and survival (Supplementary Figure S5). The magenta module had the most significant correlation with copper concentration of any module, and consisted of genes whose induction was first detected at 8-10 µg/L, and whose expression was further elevated FIGURE 3 | The frequency of sigmoidal genes' lowest observed transcriptional effect concentration (LOTEC) was plotted as a histogram for downregulated genes in blue -trials 1 (A) and 2 (C), and upregulated genes in gray -trials 1 (B) and 2 (D). These values were compared to normal development (black line) in each trial. LOTEC thresholds are more sensitive indicators of copper-induced expression changes than EC50 values, which is evidenced by the shift toward lower copper concentrations compared to frequency histograms of EC50 values (Supplementary Figure S4). at progressively high copper concentrations. The saddle brown module was characterized by substantial upregulation between 8 and 15 µg/L copper. While neither module was enriched for any specific GO terms, the saddle brown module contained numerous cysteine peptidase genes, as well as several additional genes with activity related to superoxide metabolism, lysosome processes, ubiquitination, and autophagy. Cysteine peptidases are proteases that contain an active site cysteine. Cysteine-containing proteins are redox sensitive, and can play a role in modulating metal ions that cause oxidative stress (Giles et al., 2003). The magenta module also contained several genes involved in the oxidoreductase activity or mediating the response to oxidative stress, including DBH-like monooxygenase protein 1 homolog, glutathione S-transferase 1, peroxisome biogenesis factor 10, putative ferric-chelate reductase 1 homolog, and sulfotransferase 1C4. Three subunits of cytochrome P450, an enzyme involved in responding to toxic compounds, were also upregulated as part of the magenta module.

The Transcriptional Concentration-Response to Copper in Adult Gill Tissue
ANOVA was used to identify a subset of 1,012 transcripts that were significantly differentially expressed across the adult gill tissue dataset, and of these 88% exhibited a sigmoidal concentration-response to copper according to SDRS, with 572 induced transcripts and 323 repressed transcripts (Figure 4). As was observed in larvae, transcript repression occurs at lower concentrations of copper than transcript induction in adult gill tissue, with the mode EC50 value for downregulated transcripts in adults at 12.5 µg/L, while the mode EC50 for upregulated transcripts was 22.5 µg/L (Supplementary Figure S7). Examination of copper-induced transcripts revealed evidence of a cell stress response, with many of the induced genes being associated with protein folding, apoptosis, and oxidative stress/redox chemistry (Figure 4, Supplementary Table S2, and Supplementary Figure S8). For example, several genes that participate in the detoxification of reactive electrophilic compounds by catalyzing their conjugation to glutathione were induced, including three isoforms of the glutathione-S-transferase family of proteins (A2, Mu4, omega-1), and glutaredoxin 1. Consistent with the induction of enzymes that utilize glutathione, data indicated that expression of both the regulatory and catalytic subunits of glutamate-cysteine ligase, the enzyme which catalyzes the first step in the biosynthesis of glutathione (Krejsa et al., 2010), were induced at high concentrations of copper. Other induced transcripts with roles in the detoxification of electrophilic compounds included thioredoxin reductase 3 which catalyzes the reduction of thioredoxin (Supplementary Figure S8) (Conterato et al., 2013), and peroxiredoxin 1 which controls the level of peroxides (Ishii and Yanagawa, 2007). Also induced were three transcription factors that are activated by changes in levels of cyclic AMP and which are reported to play a pivotal role in the cellular response to oxidative stress, cyclic AMP-dependent transcription factor ATF-3, cAMP-responsive element-binding protein 1, and cAMP-responsive element-binding protein-like 2 (Hai et al., 1999;Lee et al., 2009).
The list of repressed transcripts was less coherent with respect to the functional roles of the members (Supplementary Table S2). At lower concentrations of copper we identified a down-regulation of transcripts that are associated with the cell cycle including Centrosome associated protein CEP25, centrosomal protein of 135 KDa, centriolin and G1/S-specific cyclin-D2 (Figure 4 and Supplementary Figure S8). At higher copper concentrations we detected the reduced expression transcripts associated with adenylate metabolism, such as adenylate kinase and adenosine deaminase, while at the highest concentrations we detected a decline in the expression of a number of subunits of Tubulin.

Larval Versus Adult Gill Transcriptional Response
Sigmoidal Dose Response Search analysis was more effective at capturing the expression dynamics of the adult response that that of the larvae. This likely reflects the toxicity of copper on the different life stages, with adult mussels showing no mortality even after exposure to 120 µg/L whereas exposure of larvae to just 25 µg/L resulted in close to 100% mortality. Thus, while adult mussels were capable of mounting an adaptive transcriptional response throughout the concentration range, the larvae are highly impacted at the high concentrations and die. Only 48 genes were identified as significantly sigmoidal by SDRS in both life history stages (Supplementary Table S3). This included 5 downregulated genes and 43 upregulated genes. The low number FIGURE 4 | Gene expression dose response profiles detected in gill tissue of adult Mytilus californianus exposed to increasing concentrations of copper. Heatmaps depict induced (A) and repressed (B) genes that exhibited a statistically significant sigmoidal transcriptional dose response. Each row corresponds to a single gene and each column corresponds to a particular concentration of copper. The average expression of each gene at each concentration is represented by a color, with yellow or blue indicating that the gene is up-regulated or down-regulated relative to the expression observed in control mussels. The genes are ranked and grouped according to their observed EC50 value, with more responsive transcripts located at the top of each heatmap. Selected genes whose expression may be pertinent to the response to copper are indicated next to each grouping.
of shared genes indicates that the larval and adult responses to copper largely consist of different sets of genes. The genes that were upregulated in both datasets were generally involved in the cell stress response. This list included seven HSP-70 isoforms and HSP-binding proteins, stress-induced phosphoprotein 1 and two subunits of the T-complex chaperonin protein that assists in protein folding (Supplementary Table S3). In both datasets, we detected a particularly robust induction of sequestosome 1, a ubiquitin-binding protein involved in the transfer of ubiquitinated proteins to the proteasome, and a marker of autophagy (Myeku and Figueiredo-Pereira, 2011). The five genes that were commonly downregulated in adults and larvae did not exhibit a consistent functional theme (Supplementary Table S3).

DISCUSSION
The M. californianus Larval Response to Copper Normal development EC50s observed in this experiment were consistent with EC50s of Mytilus galloprovincialis, Mytilus edulis, and Mytilus californianus from previous studies (Martin et al., 1981;Arnold et al., 2009;Funk, 2014;EPA, 2016). The measured EC50s of 6.04 and 3.78 µg/L copper in the first and second trial, respectively, were higher than the U.S. Environmental Protection Agency recommended chronic water quality criteria of 3.1 µg/L copper, and the trial 1 EC50 was higher than the recommended acute criteria of 4.8 µg/L (EPA, U. S, n.d.). The observed differences in LC50 and EC50 between trials is similar to previous reports showing that mussel larvae-based water quality testing conducted by different entities can yield different results, and that copper tolerance varies between populations (Hoare et al., 1995a). In this study, trials were conducted with progeny from the same population, but from different parental crosses. Bivalves exhibit natural high-standing genetic variation, which can result in a range of phenotypes from different crosses of bivalves from the same population (Curole and Hedgecock, 2007). The differences in toxicity responses to copper reported here are expected to be linked to genetic variation between the parents used in each cross. This is manifest as left or right-shifted differences in copper concentration-response curves in each of the crosses.
Despite differences in the toxicity results of the two trials, the transcriptional response was qualitatively similar between the trials, with progeny from the more sensitive cross exhibiting a similar transcriptional response at lower concentrations than the more copper-resistant cross. Sigmoidal concentration response modeling and co-expression analysis revealed two broad patterns of modulation-transcriptional repression was a marker of exposure to lower copper concentrations, while transcript induction occurred at higher copper concentrations (Figures 2, 3). Generally, transcriptional repression, as measured by the expression LOTEC and EC50 values, occurred at copper concentrations that were below the observed normal development EC50, indicating that repressed genes will probably serve as the most sensitive biomarkers of copper exposure (Figure 3 and Supplementary Figure S3). Indeed, a number of transcripts appear to be repressed at copper concentrations at which toxicity is not apparent in the larvae with respect to development or mortality, so these genes appear to indicate copper exposure, and provide an early warning that copper concentrations are approaching toxic levels. On the other hand, transcript induction was generally initiated at higher copper concentrations, at which toxicity was evident in the larvae. Thus, transcripts that are elevated by copper can serve as markers of toxicity and the induction of a stress response, but they aren't as useful as predictors of impending stress.

Functional Characterization of Sensitive Larval Biomarkers
Sensitive downregulated transcripts were characterized by several broad functional categories. First, genes involved in development were well-represented among the transcripts that were repressed at low copper concentrations ( Table 1,  Supplementary Table S1, and Supplementary Figure S6A). Downregulation of development genes has been observed before in two studies of post-larval red abalone and scallops exposed to low copper concentrations (Zapata et al., 2009;Silva-Aciares et al., 2011). Specifically, developmental genes related to the nervous and muscle system, including numerous calcium binding proteins, were notable among downregulated genes ( Table 1, Supplementary Table S1, and Supplementary Figures S5B, S6E). While metals such as mercury and lead are well known for inhibiting nervous system function (Weis, 2014), data reflecting the effects of copper on neural development is sparse, although exposure of zebrafish embryos to copper reduced the formation of sensory neuroblasts that are used in swimming (Johnson et al., 2007). Several genes with regulatory roles in muscle development, such as calponin, were also down-regulated in response to copper, which is plausible given that muscle and neuronal development are closely integrated in early larval stages (Dyachuk and Odintsova, 2009). Although developmental pathways are not commonly recognized biochemical targets of copper toxicity in invertebrates, altered expression of these genes is probably not that surprising given that copper exposure results in morphological abnormalities in larva. Fundamental neuronal and muscular developmental pathways are activated early in development, generally during or immediately after the trochophore stage (Voronezhskaya et al., 2008;Dyachuk and Odintsova, 2009), and continue to progress through the early veliger stage at 48 h post-fertilization. Thus, drastic structural malformation and erratic swimming behavior of abnormal larvae could be explained by the systemic down-regulation of developmental programs and neurological function. Here, we have discovered some potential drivers of this phenomenon which precede developmental abnormality, and may be indicative of the early traces of failure. The exact mechanism by which copper causes developmental and neurological malfunction, and whether copper generally affects these genes systematically or targets primarily a few regulators, is still unclear.
A second key functional group that comprised many sensitive, down-regulated biomarkers was biomineralization and shell matrix protein genes. This group included several known and purported shell formation genes (Supplementary Figure S6C): protein PIF (Suzuki et al., 2009), chitin synthase (Schönitzer and Weiss, 2007), several tyrosine and tyrosinase genes (Aguilera et al., 2014), and components of the tyrosine metabolic pathway (Liu et al., 2015). Copper suppression of shellformation pathways in larvae has not been observed, yet this offers additional potential explanation of the eventual failure of abnormal larvae to develop a normal D-shaped shell in the early veliger stage. The suppression of calcium binding proteins, which was mentioned as a pathway involved in nervous system development, could also be associated with failure in shell development.
Finally, genes coding for divalent metal-containing enzymes were down-regulated across a range of low copper concentrations ( Table 1, Supplementary Table S1, and Supplementary Figure S6D). Interestingly, most of these genes code for enzymes that bind to divalent cations other than copper, including Zn 2+ , Mo 2+ , and Fe 2+ . Copperinduced downregulation of iron and zinc binding stress-protein transcripts was observed previously in juvenile abalone (Silva-Aciares et al., 2011). Metals often bind organic ligands without much specificity (Williams, 1981(Williams, , 1984, and copper is the most likely to replace almost any other divalent metal at a binding site (Irving and Williams, 1953). Copper is already known to replace zinc in some proteins (Viarengo, 1985), so as copper concentration increases relative to other metal ions, enzymes may exhibit even more binding promiscuity. Thus, downregulation of metal-binding enzymes that have become functionally disrupted by erroneous copper binding may be a necessary and protective response mechanism. Similarly, downregulation of metal-binding enzymes may serve as a mechanism to maintain homeostasis by limiting metal requirements and thus limiting metal uptake from the environment.
Metallothioneins in our study showed either no response in adults, or induction and down-regulation, depending on the isoform, in larvae. Previous research also confirms that metallothionein gene expression is not a consistent marker of copper exposure. In zebra mussel larvae, metallothionein exhibited an inconsistent response to copper (Navarro et al., 2011), and adult gill in the same study also showed no response of metallothionein to copper. In our study of larvae, MT 10-III and IV were downregulated in response to copper, while MT 20-IV was upregulated. Previous research shows that MT-10 is a constitutive isoform, while MT-20 is a metal-inducible isoform, and that each isoform exhibits differential copper binding (Vergani et al., 2007). These differences could explain the opposing expression patterns of the two MT isoforms observed in our study.

The Toxic Response Across Life History Stages
We compared adult and larval transcriptional responses to look at commonalities and differences across life history stages. The functional themes associated with the sets of up-regulated genes were somewhat consistent between larvae and adults. Genes involved in oxidative stress, damaged protein turnover, and signaling were up-regulated in both life history stages. As many of these are generic cellular stress response pathways, which are highly conserved across taxa as well (Kültz, 2005), it is not surprising that these pathways are common in larvae and adults. The specific genes in each of these categories were quite different, however (Supplementary Tables S2-S4), which indicates that larvae and adults may be using different means to reach the same end of cellular defense and protein turnover. The copper-induced expression of certain genes such as heat-shock protein 70, sequestosome-1, and glutathione-S-transferase have been reported in other studies of both adult and larval mollusks suggesting that these genes likely represent universal markers of copper stress (Navarro et al., 2011;Negri et al., 2013;e.g., Zapata et al., 2009;Varotto et al., 2013;Xu et al., 2016).
The most striking expression response observed in both developing larvae and a mature adult tissue was that transcript repression began at relatively low copper concentrations, and preceded any evidence of transcript induction. This conserved pattern is consistent with what we might expect from an energetic model of stress tolerance. In such a framework, the earliest phases of stress response involve limitation of growth, development, and reproduction to divert energy to maintenance and homeostasis (Sokolova, 2013). As these processes all require metals to some extent, limitation of these processes could also reduce metal requirements and thus reduce copper uptake from the environment. As stress intensifies, pathways involved in oxidative stress and protein chaperoning are activated to reduce further harm to vital metabolic and respiratory processes.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/, PRJNA639354; https://www.ebi.ac.uk/arrayexpress/, accession number E-MTAB-810.

AUTHOR CONTRIBUTIONS
MH contributed to the experimental design, executed the experiments, analyzed the data, and authored the manuscript. JM contributed to the research idea and experimental design. AG contributed to the research idea and experimental design, data analysis, and authored the manuscript. All authors contributed to the article and approved the submitted version.