Environmental Stress Responses and Experimental Handling Artifacts of a Model Organism, the Copepod Acartia tonsa (Dana)

Handling animals during experiments potentially affects the differential expression of genes chosen as biomarkers of sub-lethal stress. RNA sequencing was used to examine whole-transcriptome responses caused by laboratory handling of the calanoid copepod, Acartia tonsa. Salinity shock (S=35 to S=5) was used as positive stress control; individuals not exposed to handling or other stressors served as negative stress control. All copepods were grown from eggs to adults without being handled or exposed to any stressors prior the experiment. Survival of nauplii and adults was estimated for up to 10 min of exposure to handling stress and salinity shock. Only adults exhibited decreased survival (44±7% with 10 min of exposure) in response to handling stress and were selected for definitive experiments for RNA sequencing. After 10 min of experimental exposures to handling stress or salinity shock, adults were incubated for 15 min or 24 h at normal culture conditions. A small number of significantly differentially expressed genes (DEGs) were observed 15 min after exposure to handling stress (2 DEGs) or salinity shock (7 DEGs). However, 24 h after exposure, handling stress resulted in 276 DEGs and salinity shock resulted in 573 DEGs, of which 174 DEGs were overlapping between the treatments. Among the DEGs observed 24 h after exposure to handling stress or salinity shock, some commonly-used stress biomarkers appeared at low levels. This suggests that a stress-response was induced at the transcriptional level for these genes between 15 min and 24 h following exposure. Since handling stress clearly affects transcriptional patterns, it is important to consider handling when designing experiments, by either including additional controls or avoiding focus on impacted genes. Not considering handling in gene expression studies can lead to inaccurate conclusions. The present study provides a baseline for studying handling stress in future studies using this model organism and others.


INTRODUCTION
Copepods provide a principal link in the transfer of energy from phytoplankton to higher trophic levels in the marine food webs, and are preferred prey for predators, such as juvenile fishes and shrimps (Turner, 2004). Given the high natural abundance of copepods and their importance for marine ecosystems, understanding how stressors affect copepods is a concern for estuarine and marine ecology. Copepods are widely used in environmental monitoring as indicators of ecosystem health (Beaugrand, 2009). Hence, stress responses of copepods used for diagnostic (e.g., ecotoxicology testing) or experimental purposes might result in wrong conclusions when interpreting or extrapolating results from procedures that entail experimental handling.
Transcriptional biomarkers that are commonly used to indicate sub-lethal effects of stress include detoxification enzymes (i.e., cytochrome P450 and Glutathionine-S-transferase), as well as stress-related proteins, or chaperones, that protect macromolecules from damage (Davies and Vethaak, 2012;Amiard-Triquet and Berthet, 2015). In general, it is an overlooked issue that some of these biomarkers may respond to stress associated with handling, capture, collection, and other events in the experimental setup of both laboratory and field studies. Failure to consider the effects of experimental handling on gene expression during studies of environmental stress could cause erroneous conclusions about the data, by either increasing the risk of false positive results or by masking treatmentspecific effects. Experimentally-induced and handling-related stress has been extensively studied in larger crustaceans (Fotedar and Evans, 2011). To our knowledge, only a few studies, e.g., Aruda et al. (2011) and Rahlff et al. (2017), have examined handling stress in copepods with targeted methods, which entail the evaluation of specific transcriptional biomarkers selected to evaluate certain stressors, typically by real-time quantitative PCR.
The aim of this study is to examine the transcriptome-wide effects of handling stress on the calanoid copepod, Acartia tonsa. Because of the growing interest in A. tonsa as a model species for experimental studies, as well as an indicator species for environmental monitoring, and a valuable live-feed species for aquaculture industries, there are important issues for selecting appropriate biomarkers and establishing an accurate baseline description of a non-stressed copepod (Kwok et al., 2015). In both laboratory and field studies, plankton nets are commonly used for collection and size separation of copepod life stages (Uye and Kuwata, 1983;Rahlff et al., 2017). Because of this, the use of plankton nets was used to represent handling stress in this study.
A. tonsa is a robust species that, when acclimated, can persist in salinities ranging from 1 to 72 S, with an optimal salinity around 15 to 22 S . Even though A. tonsa is more tolerant to salinity variation than other Acartia species, abrupt change in salinity has been documented as a significant stressor (Lance, 1964;Chinnery and Williams, 2004;Calliari et al., 2006). Abrupt changes in salinity that exceed 10-15 S relative to the ambient level of A. tonsa have been shown to decrease survival more than 50% (Cervetto et al., 1999). For a positivestress control, we used an extreme salinity shock from S = 35 to S = 5 to provoke a response at both the transcriptional and physiological levels.
In addition to the overall lack of data in relation to handling stress, there is no information regarding the extent to which handling will cause changes at the transcriptional level. Since this is a first approach to examine transcriptomewide handling stress, the treatments used may be considered somewhat "extreme, " and are designed to ensure a response at the transcriptional level. The intention of the present study is to establish a foundation for future studies for this model species and others, in which handling stress can be described in greater detail.

Experimental Cultures
The experimental design comprised three treatments: control (no stress) (Figure 1A), handling stress ( Figure 1B) and salinity shock ( Figure 1C). Embryos of mixed age and stage were harvested from three stock cultures and transferred to 2 L Nalgene© polycarbonate bottles (Thermo Fisher Scientific, USA) containing 0.2 µm filtered seawater (S = 35), where they were incubated at 17 • C for hatching. For each of the treatments triplicate cultivation flasks were set up as biological replicates. During cultivation, as well as during all experimental treatments, the flasks were kept at 17 ± 1 • C, with dim light and gentle aeration. R. salina were fed to the copepods in excess (>800 µg C L −1 ; Berggreen et al., 1988) daily. Oxygen content was measured daily with a hand-held oxygen-probe (Handy Polaris 2, OxyGuard International A/S, Denmark) and exhibited values ranged from 6.9 to 7.5 mg O 2 L −1 .
Nauplii (body length: 130 ± 16 µm, 4 days of development n = 45) and adults (prosome length: 730 ± 54 µm, 15 days of development) were monitored for survival. Adults (764 ± 42 µm, 15 days of development n = 48) used for RNA sequencing were grown from eggs to the desired life stage without being handled. Prosome lengths were measured by photographing the copepods with a Nikon SM218 microscope with 13.5x magnification, mounted with a Nikon Digital sight DS-U3 camera (20×magnification) and subsequently analyzing the images using the software package NIS-Elements BR 4.40 (Nikon Instruments Europe, B.V., The Netherlands). The control consisted of Acartia tonsa grown from eggs to adults prior to the experiment without being handled, followed by incubation at regular culture conditions for 15 min or 24 h, prior to fixation with RNAlater. (B) In the handling stress treatment, adult individuals were placed on Nitex plankton net material with mesh size of 250 µm, and then kept at stock culture conditions for 15 min or 24 h before fixation with RNA later. (C) The salinity shock treatment consisted of animals exposed to salinity S = 5 for 10 min, followed by incubation at stock culture conditions for 15 min or 24 h. Each treatment was performed in triplicate, with each sample containing 10 adult individuals of A. tonsa.

Experimental Design
The experimental design was repeated twice, once to estimate survival and a second time for RNA sequencing (Figure 1). The initial determination of survival was used to identify both the life-stage and the exposure time for the definitive experiment for RNA sequencing.
For the control, nauplii or adults were removed from the triplicate cultures in Nalgene© bottles using a 25 mL automatic pipettor (NS 29.2/32,Witeg,Germany) and transferred to 1 L glass beakers containing 975 mL seawater, with the same conditions as described for the stock culture. The copepods were held in the control treatment for 25 s or 10 min before being transferred to new glass beakers. The copepods exposed to salinity shock were transferred to seawater with a salinity of S = 5; copepods exposed to handling stress were transferred to a Nitex plankton net (mesh size 54 µm for nauplii, 250 µm for adults) that was not submerged in water. For salinity shock and handling stress, the copepods were exposed for 25 sec, 1, 5, or 10 min before being transferred to glass beakers containing seawater with the same conditions as the stock culture. To distinguish between alive and dead, neutral red stain was added (15 mg/L seawater; Elliott and Tang, 2009). After rinsing with distilled water, dead and alive copepods were counted and survival was estimated. All three treatments were performed with 4 replicates.
Based on the results of the initial survival experiment, the definitive experiment for RNA sequencing was designed to include adult copepods and an exposure time of 10 min for handling stress or salinity shock, after which copepods were incubated for 15 min or 24 h at stock culture conditions. Seawater was then gently removed by inverted suction and the copepods were preserved in 20 mL RNAlater. Ten individuals from each triplicate experimental treatment were immediately transferred to 1 mL fresh RNAlater and stored at −20 • C (Figures 1A-C).
Additionally, survival was estimated for copepods 24 h after exposure to 10 min salinity shock or handling stress with a control of non-handled individuals, as described for the initial survival estimation.
Statistical analysis and preparation of graphics were done using R (www. R-project.org, ver. 3.4.0). From the counts of dead and alive copepods, survival (%) was calculated for each exposure time and analyzed by linear regression. Differences between treatments were analyzed based on two-tailed t-tests of the regression coefficients (Sokal and Rohlf, 1995).

RNA Extraction, Library Preparation, and Sequencing
Total RNA was extracted from A. tonsa using the RNeasy Mini Kit (Qiagen GmbH, Hilden, Germany). For each treatment (control, handling stress, salinity shock), three biological replicates were analyzed, each consisting of 10 pooled individuals (prosome length = 758 ± 67 µm).
After removal of excess RNAlater (Sigma Aldrich, St. Louis, MO, USA), the copepods were homogenized in 50 µL RLT buffer using disposable micro-pestles, after which 550 µL RLT buffer was added. The samples were vortexed for 1-2 s, and then processed according to the manufacturer's protocol, with a final elution volume of 30 µL in RNase-free water.
RNA quality was assessed using an Agilent Tapestation 2200 with RNA High Sensitivity Assay (Agilent Technologies, Santa Clara, CA, USA). DNase treatment was not done, since previous extractions treated with the Turbo DNA-free kit (Ambion, Life Technologies, Carlsbad, CA, USA) exhibited significant degradation. The RNA quality profiling of A. tonsa exhibited a merged peak of 18S rRNA and 28S rRNA, presumably resulting from "hidden break" in 28S rRNA typical of many arthropods, which causes 28S rRNA to run at about the same size as 18S rRNA (McCarthy et al., 2015).
Library preparation was done using 360 ng total RNA from each sample following the manufacturer's protocol for the Illumina Stranded mRNA Library Preparation Kit (Illumina, Inc., San Diego, CA, USA). The libraries constructed from the 18 samples (3 treatments × 3 replicates × 2 incubation times) were multiplexed and sequenced in 2 runs across 4 lanes on the NextSeq500 platform (Illumina, San Diego, CA, USA), with a mid-output 150 cycle kit (FC-404-2001, Illumina, Inc., San Diego, CA, USA) with 75 bp paired-end reads and a sequencing depth of 25 million reads per sample. Library preparation and sequencing were carried out at the Center for Genome Innovation at the University of Connecticut (Storrs, CT, USA).
A reference transcriptome was determined from RNA extracted from a single individual of A. tonsa (female, prosome length 722 µm) selected at random from the control with 24 h of incubation. Total RNA (140 ng) was sequenced in 4 lanes on a NextSeq500 platform using a mid-output 300 cycles kit (FC-404-2003, Illumina, Inc., San Diego, CA, USA), with 150 bp paired-end reads resulting in ∼350 million reads.

De Novo Transcriptome Assembly and Differential Gene Expression Analysis
FastQC (ver. 0.7; Andrews, 2010) was used to validate the quality of the raw sequence reads. Illumina adapter sequences and low-quality reads (Phred score < 20) were removed using Trimmomatic (ver. 3; Bolger et al., 2014) in paired-end mode, with a sliding window across an average of 4 bases. Initial read biases, introduced by random hexamer priming under cDNA synthesis, were corrected by removing the first 12 bp of each read (Hansen et al., 2010). Reads > 50 bp after quality trimming were retained, resulting in a total of ∼225 million reads. The reference transcriptome was assembled de novo with Trinity (ver. 2.3.2; Grabherr et al., 2011) using default parameters for pairedend reads, with normalization to decrease run time and memory requirements.
The completeness of the reference transcriptome was evaluated using the Benchmarking Universal Single-Copy Orthologs (BUSCO, ver. 2; Simão et al., 2015), which defines a set of eukaryotic core genes to test the proportion and completeness of these genes in the transcriptome assembly. Bowtie2 (ver. 2.2.6; Langmead and Salzberg, 2012) was used to examine the RNA sequencing read representation of the assembly by realigning the input reads to the de novo transcriptome. Contig N50 and E90N50 statistics were computed based on the scripts Frontiers in Marine Science | www.frontiersin.org included in the Trinity software package, as well as transcript abundance estimation using Kallisto (ver. 0.43.0;Bray et al., 2016). The A. tonsa Transcriptome Shotgun Assembly project has been deposited at DDBJ/EMBL/GenBank under the accession GFWY00000000. The version described in this paper is the first version, GFWY01000000.
The reads from the experimental samples were pseudoaligned (i.e., rapid determination of compatibility between reads and targets, without the need for alignment) to the reference transcriptome and quantified using Kallisto (ver. 0.43;Bray et al., 2016), with 100 rounds of bootstrapping. The bootstrapping option in Kallisto accounts for technical variability and is used to estimate the probability of correct assignment to a transcript. Differential gene (and transcript) analysis was performed with Sleuth (ver. 0.29; Pimentel et al., 2017), using the likelihood ratio test (LRT) and the Wald test in R Core Team (2017) to estimate significant results (ver. 3.4.0; R Core Team, 2017). Statistically significant (q-values < 0.05) differential gene expression is reported as beta values, which are bias estimators of the foldchange that accounts for the technical variability of transcripts and are reported as natural log values (Pimentel et al., 2017; see Supplementary Material S4 for results of Sleuth analysis of non-annotated transcripts).
The Trinity transcript identifications for non-annotated transcripts were added into the gene-level analysis in Sleuth. Transcripts showing statistically significant differential expression (q-value < 0.05) were isolated and annotation was attempted using Blast2Go and the RefSeq database with the arthropod taxonomy filter, in order to maximize the proportion of identified genes in the analysis (Götz et al., 2008) (Supplementary Material S2). Based on the BlastX results (E ≤ 10 −3 ), gene symbols for identified transcripts were added to the gene-level analysis in Sleuth. Transcripts that could not be identified were excluded from the analysis of differential gene expression.
Functional enrichment of differentially expressed genes was performed using clusterProfiler (ver. 3.6.0; Yu et al., 2012) for gene ontologies (GO) of cell compartments (CC), biological processes (BP) and, molecular functions (MF). The genes annotated by Trinotate from the reference transcriptome were used as background gene list. Differentially up-and downregulated genes for handling stress and salinity shock (15 min and 24 h after exposure) were used as input data. The functional enrichment was performed with a Fisher Exact test (p-value < 0.05, FDR < 0.1). Graphs and subsequent data handling was done in R Core Team (2017) using the ggplot2 package (ver. 2.2.1).

Survival
Survival of A. tonsa nauplii was statistically significantly (p < 0.001) affected by salinity shock (i.e., exposure to S = 5 for up to 10 min), in comparison to the control ( Figure 2E, Table 1E). Naupliar survival in the control was 91 ± 4% (mean ± SD), while survival after salinity shock after 25 sec exposure was 95 ± 1%, declining to 22 ± 4% after 10 min exposure. Naupliar survival after handling stress was 97 ± 2% for exposure up to 10 min, which did not differ from the control, with an average survival of 97 ± 1% ( Figure 2D, Table 1D). Survival for handling stress and salinity shock differed significantly by regression analysis (p < 0.001) ( Figure 2F, Table 1F).
Adult individuals of A. tonsa were significantly affected both by salinity shock (p < 0.001) and handling stress (p < 0.001) for exposure up to 10 min (Figures 1A,B, Tables 1A,B) in comparison to the control, which exhibited survival of 98 ± 1% up to 10 min exposure time. Survival declined from 91 ± 2% after 25 s to 56 ± 7% after 10 min of exposure to handling stress. Survival declined in the salinity shock from 92 ± 1% after 25 s to 70 ± 5% after 10 min of exposure. Survival was higher when exposed to salinity shock than when exposed to handling stress, and a significant relation was found by regression analysis (p < 0.05, Figure 2C, Table 1C).
Survival of salinity shock and handling stress was estimated 24 h post 10 min of exposure. The control exhibited a survival of 99 ± 1%, salinity shock 47 ± 1% and handling stress 54 ± 2%. In the samples, with individuals exposed to handling stress, 34 ± 1% exhibited physical damage of the antennae, setae and antennules. Of the damaged individuals, 74 ± 3% were categorized as dead during staining. For the salinity shock samples, only 4 ± 3% exhibited damage, of which 53 ± 1% was categorized as dead.

De Novo Reference Transcriptome and Annotation
A total of ∼225 million reads >50 bp in length (after quality trimming) was retained for reference transcriptome assembly. The de novo assembled transcriptome consisted of 60,688 contiguous consensus sequences (contigs) grouped into 27,171 Trinity components ("genes") with a GC content of 38.49%. Statistics based on all transcript contigs had an N50 value of 1,874 bp, with an average contig length of 1,222.45 bp from for a total of 74,188,026 assembled bases ( Table 2).
The quality of different Trinity transcriptome assemblies was evaluated using Bowtie2 for realignment of the reads to the reference, BUSCO evaluation of completeness, and E90N50 profiles of contig length ( Table 2).
Of the reference input RNA sequencing reads realigned with Bowtie2, 90.35% were represented in the assembly of the chosen reference transcriptome ( Table 2). The remaining unassembled reads, likely corresponded to low-expressed transcripts with insufficient coverage to enable assembly, was of low quality or resulted from aberrant reads.
The Ex90N50 transcript contig length of 2,731 bp was computed by combining Kallisto (ver. 0.43.0; Bray et al., 2016) FIGURE 2 | Linear regression analysis of survival (%) vs. time elapsed for various treatments compared to the control treatment. The handling stress treatment consisted of placing copepods on a Nitex plankton net screen (adult mesh size: 250 µm; nauplii mesh size: 54 µm). The salinity shock treatment consisted of exposure of copepods to S = 5. (A) A. tonsa adults exposed to handling stress (black) vs. control (gray). The two treatments differed statistically significantly from (Continued) FIGURE 2 | each other ( Table 1, A). (B) A. tonsa adults exposed to a salinity shock (black) vs. the control (gray). Survival in the two treatments differed statistically significantly ( Table 1, B). (C) A. tonsa adults exposed to salinity shock (black) vs. individuals exposed to handling stress (gray). Survival in the two treatments differed statistically significantly from each other ( Table 1, C). (D) A. tonsa nauplii exposed to handling stress (black) vs. control (gray). Survival in the two treatments did not differ statistically significantly from each other (see Table 1, D). (E) A. tonsa nauplii exposed the salinity shock (black), vs. control (gray). The two treatments differed statistically significantly from each other ( Table 1, E). (F) A. tonsa nauplii exposed to the salinity shock (black) vs. individuals exposed to handling stress (gray). The two treatments differed statistically significantly from each other ( Table 1, F).  Figure 2. SE slope is the standard error of the given regression coefficients. N is the sample size. SSE is the error sum of squares. T is the calculated "T-test" value and df is the degree of freedom where p is the probability value. and the ExN50 statistic script included in the Trinity package. Since N50 statistics discard read coverage, E90N50 gave an indication of whether deeper sequencing would result in higher quality assembly. The ExN50 profile peaked at N69, with a contig length of 3,075 bp ( Table 2).
Considering the overall realignment, BUSCO profile, and ExN50 profile, we evaluated the reference transcriptome to be of acceptable quality for the differential gene expression analysis. The Trinotate annotation pipeline resulted in identification of 45% of the assembled Trinity transcripts ( Table 2). The remaining 55% of unidentified transcripts were excluded from the differential gene expression analysis, after ensuring that the significantly differentially expressed transcripts could not be identified in any way (see Supplementary Material S4 for Sleuth analysis including unidentified transcripts.). The average pseudoalignment of the experimental sample reads to the reference transcriptome using Kallisto was 82.2 ± 2.6% (mean ± SD) ( Table 2).
From the Trinotate annotation, 80,047 Gene Ontologies (GO) and 24,463 Kyoto Encyclopedia of Genes and Genomes (KEGG) ontology (KO) terms were associated with the genes. It should be noted that multiple GO-terms could be assigned to the same gene. The GO-terms consisted (p < 0.05, FDR < 0.2) of 48% Biological Processes (BP), 30% Cellular Compartments (CC) and 22% Molecular Functions (MF) (see Supplementary Material S5 for all enriched GO terms).

Differentially Expressed Genes
The patterns of differential gene expression for copepods from the handling stress, salinity shock and control 15 min after exposure were not distinct in the Principal Component Analysis (PCA, Figure 3A). The control exhibited more distinct clustering from handling stress and salinity shock treatments 24 h after exposure ( Figure 3B), suggesting that time following exposure has a significant effect on gene expression.
Handling stress resulted in 276 DEGs 24 h after exposure, of which 177 were up-and 99 down-regulated (Figures 4, 6, Table 3). None of the same DEGs were overlapping 15 min and  24 h after exposure. The salinity shock resulted in 7 DEGs (6 upand 1 down-regulated) 15 min after exposure, which increased to 396 DEGs (221 up-and 175 down-regulated) 24 h after exposure (Figure 4, Table 3).
Among the up-regulated DEGs 15 min after exposure to salinity shock, the majority of enriched GO-terms were transport mechanisms, especially related to ER homeostasis and proteins (see Supplementary Material S5 for the full list, Table 4 for top 10 of the GO-terms with most involved genes). The remainder of the enriched GO-terms were related to metabolic, homeostatic, and developmental processes ( FIGURE 4 | Comparative analysis of differentially expressed genes (DEGs) in Acartia tonsa exposed to 10 min handling stress or 10 min salinity shock after 15 min or 24 h. The Venn diagram was constructed from overlapping transcripts in R with the library VennDiagram (ver. 1.6.17).
The overlap between handling stress and salinity shock 24 h after exposure included 174 DEGs. For the overlapping 112 up-regulated DEGs, 223 GO-terms were enriched (see Supplementary Material S5 for the full list, Table 5 for top 10 of the GO-terms with most involved genes). No KO-terms were enriched. For the 62 down-regulated DEGs, 176 GO-, and 7 KEGG terms were enriched (Supplementary Material Table 5, Table 5).
The remaining non-overlapping 102 and 202 DEGs for handling stress and salinity shock 24 h post exposure, respectively, may be stressor-specific see Supplementary Material S5 for the full list, Table 6

DISCUSSION
The results from our examination of transcriptome-wide responses and survival in relation to handling stress and salinity shock in A. tonsa clearly indicate that handling stress is a significant factor for experimental manipulation of this species. Failure to consider this parameter in, for instance, biomarkerrelated transcriptional studies is likely to lead to inaccurate interpretation of the stressors' impact.
The survival rates of the exposed nauplii were not affected by handling stress. This may be because of their small size and lack of fragile appendages, allowing a boundary layer of seawater to form around them as protection from physical interaction with the plankton net. Since handling stress did not negatively impact naupliar survival, we used only adult individuals of A. tonsa for the definitive experiments for RNA sequencing analysis. The use of adult females is most usual for copepod incubation experiments to obtain physiological rates (e.g., oxygen consumption, grazing, growth, and specific egg production) based on published results.
Adult copepod survival decreased by 44 ± 7% (mean ± SD) after 10 min exposure to handling stress. The average mortality reported by Jepsen et al. (2007) was 17 ± 1% per day, which was considered to result from handling stress in adult individuals of A. tonsa every 12 h during their study. In comparison, other studies have reported constant daily mortality ranging from 5 to 10% (Medina and Barata, 2004;Drillet et al., 2014;Nilsson et al., 2017). The observed handling stress effect on survival in the present study is considered very high and indicative of impacts of handling stress for adult individuals of A. tonsa. Handling of adult copepods should therefore be done carefullyand quickly-in order to minimize the imposed stress, in contrast to nauplii, which do not seem to be affected on survival.
In comparison to handling stress, salinity shock resulted in a slightly higher adult survival (29 ± 5%) after 10 min of exposure. Calliari et al. (2008) found that an abrupt reduction in salinity from S = 35 to S = 4 resulted in 31% decrease in survival, which is comparable to the decrease of 22-35% observed in the present study.
Survival 24 h after exposure to 10 min handling stress or salinity shock declined 45 ± 1.7% and 51 ± 1.6%, respectively. For handling stress, this was an additional decline of ∼7% compared to 15 min post exposure. For the salinity shock, the post-exposure period of 24 h resulted in higher mortality than for handling stress, with an additional decline in survival of ∼22%. This suggests that handling stress results in higher instant mortality than salinity shock, while salinity shock is more lethal over an extended post-exposure period. Part of the instant mortality may be due to the physical damages (i.e., broken antennae) observed for about one-third of the individuals, most of which were categorized as dead. In addition to the elevated mortality, the higher number of DEGs suggests that salinity shock is more stressful than handling.
In order to replicate the impact of "collection stress" and laboratory handling stress on A. tonsa, it was important to establish a baseline for copepods that were not exposed to handling stress. Since this is the first study examining impacts of handling stress at the transcriptome-wide level, there is no available information on the time course of the copepods' responses after handling stress. To avoid this as a bias in our experiments, we cultured A. tonsa from eggs to the desired life-stages in ideal culture conditions without disturbing the copepods.
For field studies, it is often impossible to know how long the copepods may have been in contact with a plankton net during collection and before preservation (e.g., Mack et al., 2012). According to the field guide by Goswami (2004), it is FIGURE 5 | Heat map of differentially expressed genes (q < 0.05) for the handling stress and salinity shock treatments 15 min post exposure relative to the control. The color scale of transcript expressions is normalized as z-scores; up-regulated transcripts are red and down-regulated are blue. The heat map was generated in R using the library gplots (version 3.0.1). recommended to collect zooplankton by towing a plankton net at slow speed (1.5-2.0 knots) for 5-10 min. Additional collection time is used in recovering the plankton net and handling the copepods for different experimental purposes. The mortality of field-collected copepods, however, has been shown to range from 0 to 90% based on direct observation and 13-37% based on neutral red staining (Elliott and Tang, 2009). This is within the range of the survival decrease observed in the present study (44 ± 7%) and suggests that observed mortality in field-collected copepods may be the result of handling stress. For laboratory studies, copepods are usually exposed to plankton net screens for shorter periods. But often other tools used for transferring copepods, such as pipettes and tweezers, may also result in physical damage and stress.
Based on the observed mortalities, we chose the exposure time of 10 min for each stressor, which is within the range of handling time described in published studies, although sufficient to ensure that we would induce a transcriptional response (Elliott and Tang, 2009;Mack et al., 2012).
An additional source of uncertainty regarding responses at the transcriptional level was the impact of the length of time post-exposure before preservation for analysis.
We found only two (1 up-and 1 down-regulated) and seven (6 up-and 1 down-regulated) DEGs 15 min after exposure to handling stress and salinity shock, respectively (Figure 5, Table 3). None of these DEGs were in common between the treatments or have previously been used as transcriptional biomarkers. Handling stress resulted in 276 DEGs and salinity FIGURE 6 | Heat map of overlapping differentially expressed genes (q < 0.05) for the handling stress and salinity shock treatments 24 h post exposure relative to the control. Similarities between samples are shown as a dendrogram, with hierarchical clustering based on Pearson correlation with complete distance determination. The color scale of transcript expressions is normalized as z-scores; up-regulated transcripts are red and down-regulated are blue. The heat map was generated in R using the library gplots (version 3.0.1). shock in 396 DEGs 24 h after exposure. Of these, 174 DEGs (112 up-and 62 down-regulated) were overlapping. The up-regulated expression of these genes may provide general protection against multiple stressors. However, the time period following exposure clearly had significant impact on whether there was a detectable and measureable stress response. The one up-regulated DEG 15 min after exposure to handling, Inositol-Pentakisphosphate 2-Kinase (IPPK), was assigned two BPs (melanosome transport and determination of left/right symmetry). IPPK in yeast is involved in the transcriptional regulation of responses to environmental and nutritional changes; in plants, it is involved in stress signaling, and in mouse embryonic development (e.g., Tsui and York, 2010). The role of IPPK is thus very diverse among species and therefore the enriched BPs seems difficult to explain. In A. tonsa, the early up-regulation of IPPK may indicate its role in initiating a transcriptional response to handling stress.
The gene product of the down-regulated SID1 Transmembrane Family Member 1 (SIDT1) is involved in RNA-interference (RNAi), by transporting dsRNA across cellular membranes (e.g., Whangbo et al., 2017). The down-regulation suggests that gene silencing of RNAi inhibited genes are being removed, which allows transcription.
The majority of enriched GO terms for up-regulated DEGs 15 min post exposure to salinity were related to protein transport, Endoplasmic Reticulum (ER), and protein homeostasis. This is an indication of ER stress, which typically induces the unfolded protein response (UPR) in order to improve the imbalance between protein load and folding capacity of the ER (Hori et al., 2006;Hetz and Papa, 2017).
Enrichment of GO terms related to metabolic processes indicated the need for cellular energy and is in agreement with the observation by Calliari et al. (2006) that A. tonsa modulates its energy balance in relation to salinity stress. The enrichment of vitamin digestion and absorption (KO , Table 4) for MAP kinaseinteracting serine/threonine protein kinase 2 (MAPK) could reflect a need for energy. However, MAPK is also induced in response to The top 10 (or fewer in some cases) GO-terms with most genes involved are shown for up-and down regulated DEGs for each treatment. The DEGs were sorted by the number of involved genes, and then by p-values for the most significant enrichments (p). Cat., ontology category; can be BP, Biological Process; CC, Cellular Compartment; MF, Molecular Function; or KEGG. p, false discovery rate corrected p-value (FDR < 0.2). Change, imply if gene expression were up-or down regulated in relation to the negative control. ID, ontology ID. Genes, abbreviations for DEGs; full names can be found in Supplementary Material S3. No DEGs 15 min post exposure were overlapping between handling-stress and salinity shock. Enrichment of GO and KEGG terms were done using the R-package ClusterProfiler (ver. 3.6.0, Yu et al., 2012) with the Trinotate annotated transcriptome as background gene list. environmental stress as a part of a signaling cascade, hence the two signaling pathways in Table 4 (Waskiewicz, 1997). The KEGG enrichment was done using general KEGG Ontology (KO) terms, which are related to the usual model organisms, human and mouse. The enriched terms and the actual functions may therefore differ in relation to copepods.
Five of the overlapping DEGs 24 h post exposure to handling stress and salinity shock were enriched for the BP, neurogenesis. These include Histone deacetylase 1 (HIDAC1),  serine/threonine-protein phosphatase 4 regulatory subunit 2 (PP4R2), helicase domino (DOM), longitudinals lacking protein (LOLA3) and selenide water dikinase (SPS1). Even though they are enriched in the GO-term, neurogenesis, it is noteworthy that only LOLA3 is directly linked to neurogenesis (e.g., Goeke et al., 2003). The gene product resulting from HIDAC1 is mainly a regulator of gene expression for other genes responsible for histone de-acetylation (e.g., Kelly and Cowley, 2013). DOM is, like HIDAC1, also responsible for transcriptional regulation by chromatin remodeling (Sif, 2004). Especially HIDAC1 is also enriched for terms related to chromatin remodeling, like chromatin (CC) and transcription corepressor activity (MF). PP4R2 has functional roles in cell development, differentiation, apoptosis, tumor progression and DNA-repair (e.g., Shui et al., 2007;Nakada et al., 2008;Liu et al., 2012).
In addition to being enriched in GO-terms related to transcriptional regulation, oogenesis is enriched for HDAC1, DOM, Ecdysone-induced protein 74EF isoform B (E74EB) and beta-1,3-galactosyltransferase BRN (BRN). Both E74EF and BRN are involved in oogenesis (Goode et al., 1996;Paul et al., 2005). When exposed to stressful conditions, energy of an individual tends to be reallocated from fecundity and growth to survival mechanisms (López-Maury et al., 2008;de Nadal et al., 2011). Thus, egg production of copepods decreases when the surrounding environment is sub-optimal (Calliari et al., 2006;Peck and Holste, 2006). The up-regulation of oogenesis may indicate that homeostasis in A. tonsa has been restored to such extent that there is energy for egg production.
Even though fecundity-related mechanisms were enriched for up-regulated overlapping DEGs 24 h post exposure to handling stress and salinity shock, stress related mechanisms were also present among the enriched GO-terms (Table 5, Supplementary Material S5). This includes cellular response to DNA damage stimulus, ubiquitin protein ligase binding and chaperone-mediated protein folding (Table 5).
Heat shock proteins (hsps), especially heat shock protein 70 kDa (hsp70), have been frequently used as a transcriptional indicator of stress in copepods (Voznesensky et al., 2004;Tartarotti and Torres, 2009;Lauritano et al., 2011Lauritano et al., , 2016Nilsson et al., 2013;Chan et al., 2014;Petkeviciute et al., 2015;Aguilera et al., 2016;Smolina et al., 2016;Rahlff et al., 2017). We found significant, but small, down-regulation of hsp70 in response to salinity shock 24 h after exposure (Supplementary Material S1; https://figshare. com/articles/S1a_Trinotate_Annotation_A_tonsa_xls/5928799/ 1). Heat shock cognate 71 kDa (hsc70), which is a member of the hsp70 family, was slightly up-regulated for handling stress and salinity shock 24 h after exposure (Supplementary Material S1; https://figshare.com/articles/S1a_Trinotate_Annotation_ A_tonsa_xls/5928799/1). Aruda et al. (2011) found that the heat shock proteins, hsp70a, hsp21, and hsp22, had significantly higher expression 3 h after handling with a plankton net, but did not find any significant differences 2 h after exposure. Rahlff et al. (2017) found significant changes in expression of hsp70 in response to handling, which was reduced to negligible levels after 24 h. These prior studies suggest that the expression of hsp70 peaks within 24 h after stress exposure, and may explain why we did not observe an increase in expression level of this gene 24 h after handling stress. Hsp70 responses in relation to other stressors (e.g., temperature) seem to be in agreement with this explanation. Petkeviciute et al. (2015) found a 63.8-fold increase in hsp70 transcripts for A. tonsa after 45 min exposure to 30 • C. This corresponds with the findings of Rahlff et al. (2017), where a heat shock of 28 • C for 3 h resulted in significant up-regulation of hsp70, which was measurable after 30 min and peaked with a 185-fold increase after 1.5 h. A smaller peak of 60.4-fold increase remained 4 h after the heat shock (Rahlff et al., 2017). The peak in hsp70 expression that occurred a few hours after exposure, which subsequently declined, could explain our findings. In general, the expression levels observed here were low, implying that hsp70 and a number of other genes may show peak up-regulation within 24 h after exposure to stressors.
The down-regulated hsp70 24 h post exposure to salinity was enriched for the nucleus (CC , Table 6), precatalytic spliceosome (CC , Table 6), presynapse (Supplementary Material S5), late endosomal microautophagy (BP, Supplementary Material S5), cellular response to topologically incorrect protein (BP, Supplementary Material S5), and perichromatin fibrils CC, Supplementary Material S5). Many of these terms are stress-related, and could have included additional GO-terms that are relevant in relation to hsp70 "chaperone mediated protein folding" and "de novo protein folding" (Supplementary Material S5).
In summary, handling stress clearly affects both biomarkers and transcriptome-wide patterns of differential gene expression of A. tonsa, and these stress responses probably take place within 24 h after exposure to a stressor. Some of the differentially expressed genes were in common between the handling stress and salinity shock treatments, suggesting that these may play an important role in protection against multiple stressors. Due to the small, but significant differences in expression levels of some of the commonly used biomarkers, these genes should be used with caution in stress-related studies, since they potentially peak within 24 h after exposure.
The limited response at the transcriptional level 15 min following exposure to handling stress suggests that organisms collected in plankton tows of short duration and immediately used in incubation experiments or being preserved are likely to exhibit transcriptional profiles that represent their in situ physiological state (e.g., Häfker et al., 2017). It is, however, important to be aware of that handling has the potential to affect gene expression regardless of animal species. Thus, handling should therefore be considered as a factor when examining stress at the transcriptional level.

DATA AVAILABILITY
RNA-Seq data available in the NCBI Sequence Read Archive with the bio-project accession number PRJNA407266.

AUTHOR CONTRIBUTIONS
BN, PJ, AB, and BH conceived the study; BN, PJ, AB, and BH designed the experiments; BN, PJ, and AB collected and analyzed the data; BN wrote the paper; BN, PJ, AB, and BH contributed substantially to interpreting the data and developing the manuscript, and take full responsibility for the content of the paper.

FUNDING
This work was supported by the Villum Foundation [Project AMPHICOP no. 8960,to BH].