De Novo Transcriptome Assembly and Gene Expression Profiling of the Copepod Calanus helgolandicus Feeding on the PUA-Producing Diatom Skeletonema marinoi

Diatoms are the dominant component of the marine phytoplankton. Several diatoms produce secondary metabolites, namely oxylipins, with teratogenic effects on their main predators, crustacean copepods. Our study reports the de novo assembled transcriptome of the calanoid copepod Calanus helgolandicus feeding on the oxylipin-producing diatom Skeletonema marinoi. Differential expression analysis was also performed between copepod females exposed to the diatom and the control flagellate Prorocentrum minimum, which does not produce oxylipins. Our results showed that transcripts involved in carbohydrate, amino acid, folate and methionine metabolism, embryogenesis, and response to stimulus were differentially expressed in the two conditions. Expression of 27 selected genes belonging to these functional categories was also analyzed by RT-qPCR in C. helgolandicus females exposed to a mixed solution of the oxylipins heptadienal and octadienal at the concentration of 10 µM, 15 µM, and 20 µM. The results confirmed differential expression analysis, with up-regulation of genes involved in stress response and down-regulation of genes associated with folate and methionine metabolism, embryogenesis, and signaling. Overall, we offer new insights on the mechanism of action of oxylipins on maternally-induced embryo abnormality. Our results may also help identify biomarker genes associated with diatom-related reproductive failure in the natural copepod population at sea.


Introduction
Diatoms are the dominant component of the marine phytoplankton, being responsible for approximately 50% of primary production in the oceans, and 20% to 25% of all organic carbon fixation on the planet [1]. It has been shown that several diatom species possess a complex infochemical system that plays an important role in allelopathy [2], phytoplankton intercellular communication [3], and phytoplankton bloom termination [4], thus shaping the structure of phytoplankton communities [5]. Several field observations have also demonstrated inhibitory effects of diatoms on the reproduction of calanoid copepods during blooms occurring in the Adriatic Sea [6][7][8], in the English Channel [9], in the Baltic Sea [10], and in the North and South Pacific Ocean [11,12]. Such harmful effects of diatoms on copepod gametogenesis, hatching success, and naupliar fitness, is due to the production BLASTx similarity search of 30,339 C. helgolandicus unigenes against non-redundant (nr) protein database resulted in 21,336 (70.3% of the total) matched unigenes and 9001 sequences (29.7%) without BLASTx hits, whereas 2 unigenes had only InterProScan exceeding BLASTx size limit of 18,000 bp in Blast2GO. Slightly lower matches were obtained against Swissprot database, with 19,386 matched unigenes (63.9%). The e-value distribution showed that 21.9% and 14.9% of the unigenes have strong homology in the nr database, with BLASTx results ranging from 0 to 1 × 10 −100 and from 1 × 10 −100 to 1 × 10 −60 , respectively; whereas 89.4% of the hits have a similarity ranging from 30% to 80% ( Figure S2).
The species distribution of the best matches (top-hit) against the nr database of C. helgolandicus transcriptome shows that the highest number of matched unigenes have similarities with the copepods Eurytemora affinis and Tigriopus californicus, followed by other crustaceans, arthropods and copepod species, thus reflecting both phylogenetic relationship and abundant genomic information for those species ( Figure S3).
Mapping the 30,339 C. helgolandicus unigenes against the KEGG pathway database resulted in 6361 unigenes with assigned KO numbers, for a total of 4710 KO terms associated with the unigenes. Using a pathway reconstruction module, these KEGG annotated unigenes were further categorized into six different functional groups. Specifically, 2663 unigenes were assigned to metabolism (41.9%), 1231 to genetic information processing (19.3%), 1542 to environmental information processing (24.2%), 1312 to cellular processes (20.6%), 2569 to organismal systems (40.4%), and 2872 to human diseases, which was the largest group (45.1%). Interestingly, signal transduction belonging to environmental information processing group was the most represented KEGG category in C. helgolandicus transcriptome (22.5%), followed by cancers (15.3%) and infectious diseases (15.2%) ( Figure 2). A summary of the statistics of the unigene annotation for C. helgolandicus is shown in Table 2. There was no difference between the C. helgolandicus transcriptome and Drosophila melanogaster genome when comparing the proportions of unigenes in the top 20 2nd level GO functional categories, suggesting that the sequenced C. helgolandicus transcriptome did not lack any major functional categories of genes ( Figure S4).
Mapping the 30,339 C. helgolandicus unigenes against the KEGG pathway database resulted in 6361 unigenes with assigned KO numbers, for a total of 4710 KO terms associated with the unigenes. Using a pathway reconstruction module, these KEGG annotated unigenes were further categorized into six different functional groups. Specifically, 2663 unigenes were assigned to metabolism (41.9%), 1231 to genetic information processing (19.3%), 1542 to environmental information processing (24.2%), 1312 to cellular processes (20.6%), 2569 to organismal systems (40.4%), and 2872 to human diseases, which was the largest group (45.1%). Interestingly, signal transduction belonging to environmental information processing group was the most represented KEGG category in C. helgolandicus transcriptome (22.5%), followed by cancers (15.3%) and infectious diseases (15.2%) (Figure 2). A summary of the statistics of the unigene annotation for C. helgolandicus is shown in Table 2.

Differential Expression Analysis
Analysis of expression levels of unigenes showed that 280 sequences were differentially expressed in C. helgolandicus, with 117 unigenes significantly up-regulated and 163 unigenes significantly down-regulated in S. marinoi-fed compared to P. minimum-fed females (FDR ≤ 0.05). All DEGs were further mapped against nr and Swissprot databases, to infer the biological function in which they were involved. Overall, 117 (52 up-regulated and 65 down-regulated) out of 280 DEGs were functionally annotated against nr and Swissprot databases and assigned one or more specific GO term (Table S1). More specifically, in the biological process category, up-regulated DEGs were allocated to a higher number of different but functionally-related GO terms with respect to down-regulated DEGs (74 vs. 47, respectively), of which 88% contained <5 sequences, with respect to 75% for down-regulated genes. This suggested that up-regulated genes encoded for proteins acting in multiple and interconnected metabolic, organismal and cellular processes, such as response to stimulus/stress, protein folding, signaling, biological regulation, lipid metabolic and glutathione metabolic processes (Table S2). In contrast, down-regulated genes showed a lower level of sequence redundancy among terms and were distributed into a more restricted GO list which mainly comprised metabolic and catabolic processes of e.g., organonitrogen compounds (amino acids and nucleotides), proteins, carbohydrates and well-defined amino acids (e.g., serine, glycine and threonine, cysteine, and methionine) (Table S3). It is possible, hence, that the response of C. helgolandicus to S. marinoi feeding is, on the one hand, the overall activation (up-regulation) of several basic organismal-related and stimulus-related cellular processes, and on the other hand, a reduction KEGG functional group classification of Calanus helgolandicus transcriptome. Number of unigenes assigned to six KEGG functional groups: metabolism (blue), genetic information processing (orange), environmental information processing (black), cellular processes (green), organismal systems (purple), and human diseases (red).

Differential Expression Analysis
Analysis of expression levels of unigenes showed that 280 sequences were differentially expressed in C. helgolandicus, with 117 unigenes significantly up-regulated and 163 unigenes significantly down-regulated in S. marinoi-fed compared to P. minimum-fed females (FDR ≤ 0.05). All DEGs were further mapped against nr and Swissprot databases, to infer the biological function in which they were involved. Overall, 117 (52 up-regulated and 65 down-regulated) out of 280 DEGs were functionally annotated against nr and Swissprot databases and assigned one or more specific GO term (Table S1). More specifically, in the biological process category, up-regulated DEGs were allocated to a higher number of different but functionally-related GO terms with respect to down-regulated DEGs (74 vs. 47, respectively), of which 88% contained <5 sequences, with respect to 75% for down-regulated genes. This suggested that up-regulated genes encoded for proteins acting in multiple and interconnected metabolic, organismal and cellular processes, such as response to stimulus/stress, protein folding, signaling, biological regulation, lipid metabolic and glutathione metabolic processes (Table S2). In contrast, down-regulated genes showed a lower level of sequence redundancy among terms and were distributed into a more restricted GO list which mainly comprised metabolic and catabolic processes of e.g., organonitrogen compounds (amino acids and nucleotides), proteins, carbohydrates and well-defined amino acids (e.g., serine, glycine and threonine, cysteine, and methionine) (Table S3). It is possible, hence, that the response of C. helgolandicus to S. marinoi feeding is, on the one hand, the overall activation (up-regulation) of several basic organismal-related and stimulus-related cellular Mar. Drugs 2020, 18, 392 6 of 26 processes, and on the other hand, a reduction (down-regulation) of specific protein-amino acid-and glucose-related metabolic pathways. This output is confirmed by the enrichment analysis on up-and down-regulated gene datasets with respect to all DEGs. The analysis showed that within the biological process category, the GO cellular process (FDR = 6.3 × 10 −3 ) and protein folding (FDR = 1.1 × 10 −2 ) were significantly enriched in the up-regulated gene set with respect to DEGs. Whereas, GO terms metabolic process (FDR = 1.0 × 10 −4 ), organic substance metabolic process (FDR = 2.0 × 10 −4 ), primary metabolic process (FDR = 1.3 × 10 −3 ), nitrogen compound metabolic process (FDR = 2.2 × 10 −2 ) and carbohydrate metabolic process (FDR = 4.3 × 10 −2 ), were enriched in the down-regulated unigenes.

Effect of PUAs on Calanus helgolandicus Gene Expression
Relative expression ratio of DEGs and other unigenes selected by manual curation, in C. helgolandicus females exposed for five days to a mixed solution of heptadienal and octadienal at a final concentration of 10 µM, 15 µM, and 20 µM, with respect to females exposed to control methanol, is shown in Figures 5-9. In particular, for unigenes involved in folate and methionine metabolism/one carbon pool by folate, there was a PUAs concentration-dependent decrease in relative expression ratios of 10-FTHFDH, TS and BHMT1, down to a minimum of −3.0-fold, and an increase for MTHFR and MS, up to a maximum of 2.6-fold; whereas THTR1, DHFR, and BHMT1 did not show significant changes ( Figure 5).
Mar. Drugs 2020, 18, x 10 of 26 metabolism/one carbon pool by folate, there was a PUAs concentration-dependent decrease in relative expression ratios of 10-FTHFDH, TS and BHMT1, down to a minimum of −3.0-fold, and an increase for MTHFR and MS, up to a maximum of 2.6-fold; whereas THTR1, DHFR, and BHMT1 did not show significant changes ( Figure 5).

Figure 5.
Relative expression ratio (log2 fold change) of GOIs related to folate and methionine metabolism in Calanus helgolandicus females exposed for five days to a mixed solution of heptadienal and octadienal at a final concentration of 10 μM, 15 μM, and 20 μM, with respect to females exposed to methanol. Values are mean ± SD.
Regarding transcripts involved in embryogenesis, all three genes VM01, AUR, and DLL were down-regulated, with a dose-dependent decrease only for the latter gene (down to −2.8-fold) ( Figure  6). For transcripts involved in signaling, PTCHD3 had a concentration-dependent increase and SMO was up-regulated at 10 μM and 20 μM, whereas CI was down-regulated at 15 μM and 20 μM and NOTUM at 15 μM. The remaining genes WNT4 and HH did not show significant gene expression Relative expression ratio (log 2 fold change) of GOIs related to folate and methionine metabolism in Calanus helgolandicus females exposed for five days to a mixed solution of heptadienal and octadienal at a final concentration of 10 µM, 15 µM, and 20 µM, with respect to females exposed to methanol. Values are mean ± SD. Letters a, b and c denoted statistically different treatments after Student's t-test or One-way analysis of variance (ANOVA) for each gene (p < 0.05).
Regarding transcripts involved in embryogenesis, all three genes VM01, AUR, and DLL were down-regulated, with a dose-dependent decrease only for the latter gene (down to −2.8-fold) ( Figure 6). For transcripts involved in signaling, PTCHD3 had a concentration-dependent increase and SMO was up-regulated at 10 µM and 20 µM, whereas CI was down-regulated at 15 µM and 20 µM and NOTUM at 15 µM. The remaining genes WNT4 and HH did not show significant gene expression changes with respect to controls and PUAs concentrations (Figure 7). Among transcripts involved in response to stimulus, only MGST3 and HSC70 were significantly down-regulated (to −1.6-fold) (Figure 8). Finally, genes involved in general metabolism, PTL, and TRET1, were significantly down-regulated, while ELOVL and PSAP did not show significant changes (Figure 9).
Mar. Drugs 2020, 18, x 11 of 26 Figure 6. Relative expression ratio (log2 fold change) of GOIs related to embryogenesis in Calanus helgolandicus females exposed for five days to a mixed solution of heptadienal and octadienal at a final concentration of 10 μM, 15 μM, and 20 μM, with respect to females exposed to methanol. Values are mean ± SD.   Relative expression ratio (log 2 fold change) of GOIs related to signaling in Calanus helgolandicus females exposed for five days to a mixed solution of heptadienal and octadienal at a final concentration of 10 µM, 15 µM, and 20 µM, with respect to females exposed to methanol. Values are mean ± SD. Letters a, b and c denoted statistically different treatments after Student's t-test or One-way analysis of variance (ANOVA) for each gene (p < 0.05).

Figure 8.
Relative expression ratio (log2 fold change) of GOIs related to response to stimulus in Calanus helgolandicus females exposed for five days to a mixed solution of heptadienal and octadienal at a final concentration of 10 μM, 15 μM, and 20 μM, with respect to females exposed to methanol.
Values are mean ± SD. . Relative expression ratio (log2 fold change) of GOIs related to metabolism in Calanus helgolandicus females exposed for five days to a mixed solution of heptadienal and octadienal at a final concentration of 10 μM, 15 μM, and 20 μM, with respect to females exposed to methanol. Values are mean ± SD.

Discussion
The present study reports the first whole-body RNA-Seq analysis of the copepod Calanus helgolandicus, generating more than 30 Gb of raw sequence data and 30,339 assembled unigenes. Former studies on C. helgolandicus transcriptome analysis were performed using EST sequencing (SSH library) and Sanger sequencing methods, which produced only a few hundred nucleotide sequences [28,30]. As for transcriptome functional annotation, our study indicated that 70.3% (21,337 out of 30,339) of unigenes had homologs in the nr database, much more than those reported for Calanus finmarchicus [32], Temora longicornis [35], Acartia tonsa [36], and Calanus sinicus [33]. This could be attributed to the quality of the assembly and/or to the application of more stringent filtering parameters.
GO and KEGG classifications revealed that the assembled C. helgolandicus unigenes have diverse molecular functions and are involved in many metabolic and cellular pathways, thus reflecting gene richness of the global landscape of the transcriptome. Overall, the highest number of annotated unigenes were involved in cellular and metabolic processes, developmental processes, response to stimulus, and signaling, similar to that which has been reported previously [30,32,35,37]. Interestingly, the KEGG pathway classification analysis showed the highest number of unigenes was classified into signal transduction and human diseases, especially infectious diseases. Aquatic organisms may often have to deal with the challenges imposed by parasitic, bacterial, and viral infections in water; thus the presence of such signal transduction pathways in copepod transcriptomes could have evolved as a defense mechanism against infection, as suggested for the white-leg shrimp Litopenaeus vannamei [42].
Our differential expression analysis of C. helgolandicus females following an oxylipin-producing diatom diet was helpful for the identification of candidate genes underlying the response of C. helgolandicus to the diatom diet and their PUAs. In particular, within the 280 differentially expressed unigenes between Skeletonema marinoi-fed and Prorocentrum minimum-fed C. helgolandicus females, Figure 9.
Relative expression ratio (log 2 fold change) of GOIs related to metabolism in Calanus helgolandicus females exposed for five days to a mixed solution of heptadienal and octadienal at a final concentration of 10 µM, 15 µM, and 20 µM, with respect to females exposed to methanol. Values are mean ± SD. Letters a and b denoted statistically different treatments after Student's t-test or One-way analysis of variance (ANOVA) for each gene (p < 0.05).

Discussion
The present study reports the first whole-body RNA-Seq analysis of the copepod Calanus helgolandicus, generating more than 30 Gb of raw sequence data and 30,339 assembled unigenes. Former studies on C. helgolandicus transcriptome analysis were performed using EST sequencing (SSH library) and Sanger sequencing methods, which produced only a few hundred nucleotide sequences [28,30]. As for transcriptome functional annotation, our study indicated that 70.3% (21,337 out of 30,339) of unigenes had homologs in the nr database, much more than those reported for Calanus finmarchicus [32], Temora longicornis [35], Acartia tonsa [36], and Calanus sinicus [33]. This could be attributed to the quality of the assembly and/or to the application of more stringent filtering parameters.
GO and KEGG classifications revealed that the assembled C. helgolandicus unigenes have diverse molecular functions and are involved in many metabolic and cellular pathways, thus reflecting gene richness of the global landscape of the transcriptome. Overall, the highest number of annotated unigenes were involved in cellular and metabolic processes, developmental processes, response to stimulus, and signaling, similar to that which has been reported previously [30,32,35,37]. Interestingly, the KEGG pathway classification analysis showed the highest number of unigenes was classified into signal transduction and human diseases, especially infectious diseases. Aquatic organisms may often have to deal with the challenges imposed by parasitic, bacterial, and viral infections in water; thus the presence of such signal transduction pathways in copepod transcriptomes could have evolved as a defense mechanism against infection, as suggested for the white-leg shrimp Litopenaeus vannamei [42].
Our differential expression analysis of C. helgolandicus females following an oxylipin-producing diatom diet was helpful for the identification of candidate genes underlying the response of C. helgolandicus to the diatom diet and their PUAs. In particular, within the 280 differentially expressed unigenes between Skeletonema marinoi-fed and Prorocentrum minimum-fed C. helgolandicus females, several key processes and/or metabolic pathways, such as response to stimulus/stress (e.g., detoxification mechanism), protein folding, signal transduction, transport, immune response, carcinogenesis, carbohydrate, lipid, protein, and aminoacid/tetrahydrofolate metabolism, were significantly altered in copepod females after five days of feeding on S. marinoi. More importantly, incubation experiments with pure PUAs supported results obtained during copepod feeding on the diatom, suggesting that the observed molecular responses were related to the ingestion of oxylipins and their direct impairment or stimulation of specific gene expression. This is the first study reporting the molecular effects of the diatom PUAs octadienal and heptadienal on copepods and will help to link molecular responses to phenotypic responses observed in copepods due to oxylipin-producing diatoms. Although the concentrations of PUAs in the incubation experiment are higher than the concentrations theoretically offered to the copepods in the feeding experiment, considering the daily algal supply of 45.000 cells mL −1 and the potential PUAs production of 2 fmol cell −1 measured in the same S. marinoi strain used in our study [43], for a total of about 0.1 µM of PUAs, it has to be pointed out that quantification of oxylipins in algal cultures is always a potential production, based on a fixed incubation time during which free fatty acids are converted into oxylipins. On the contrary, oxylipin production in the copepod is an active process occurring continuously in the copepod gut, starting from algal ingestion, crushing of cells, and release of the oxylipins. It is, therefore, possible that the actual amount of oxylipins ingested by the copepod during the feeding experiment is higher than the potential production by the algal culture. Moreover, concentrations of 10-20 µM of PUAs induced the same reduction of egg viability and increase of naupliar abnormality in C. helgolandicus embryos as observed during algal feeding experiments with S. marinoi [13].
The activation of drug (xenobiotics) metabolism and detoxification systems is one of the most common responses to exogenous chemical compounds. In this study, several sequences similar to detoxification system genes were up-regulated in C. helgolandicus following feeding of S. marinoi, including microsomal glutathione s-transferase 3 (MGST3), pyrimidodiazepine synthase and xanthine dehydrogenase/oxidase. The same was observed for several genes coding for proteins involved in generic stress/stimulus response and protein folding, such as prophenoloxidase activating enzyme (PPAE), zinc finger protein OZF-like (OZF), heat shock cognate protein 70 (HSC70), E3 ubiquitin-protein ligase SINAT3 and dnaJ protein homolog 1. These results were confirmed by the KEGG pathway reconstruction of up-regulated genes which showed high percentages in folding, sorting and degradation, xenobiotics biodegradation and chemical carcinogenesis, the latter likely associated with the shared reactive functional group of heptadienal and octadienal with malondialdehyde and 4-hydroxy-2-nonenal, well-known PUAs having carcinogenic activity [44].
Interestingly, our findings differed from those described in [28], which reported down-regulation of genes involved in stress response and defense (glutathione S-transferase and cytochrome P450 enzymes) in C. helgolandicus feeding on the same S. marinoi strain. This could be related to the shorter feeding period as compared to the present experiment (2 days vs. 5 day), thus suggesting that prolonged ingestion of the diatom might induce stronger toxic effects eliciting activation of detoxification systems. Up-regulation of genes associated with repair/degradation systems in the present study confirm findings by [30] that S. marinoi induced over-expression of genes involved in protein folding or degradation in C. helgolandicus, thus protecting the adult copepod from the direct toxic effect of the diatom diet. Cellular chaperones (i.e., heat shock proteins) can transfer irreparably damaged proteins for degradation through the ubiquitin-proteasome pathway [45]. Similar up-regulation of transcripts encoding proteins involved in the ubiquitin-proteasome pathway has also been shown in the copepod C. finmarchicus following exposure to diethanolamine [46]. Our incubation experiments with pure PUAs in general confirmed results induced by the diatom diet, with 2-3-fold up-regulation of PPAE and OZF, although opposite trends were observed for MGST3 and HSC70. It is possible that the amount of PUAs indirectly ingested by copepod females in this in vitro test was too low to induce the activation of such a detoxification system, thus resembling shorter feeding trials [28].
The gene showing the highest up-regulation in C. helgolandicus females feeding on S. marinoi was the prophenoloxidase activating enzyme PPAE (8-fold), an enzyme with endo-peptidase activity involved in proteolysis. In crustaceans, activation of prophenoloxidase, together with pattern recognition proteins, complement, and coagulation cascade, is part of the organism immune response [47]. In the current study, the sequence complement C1q tumor necrosis factor-related protein 3-like was part of the up-regulated gene set, which also showed a high number of sequences allocated to KEGG pathway category "antigen processing and presentation," related to the immune system. We also identified in the whole transcriptome transcripts belonging to pattern recognition protein lectins (c-type lectin domain family 4 member k). A similar immune response was reported for the Chinese shrimp Fennropenaeus chinensis during white spot syndrome virus (WSSV) acute infection [48] and in white shrimp Litopenaeus vannamei exposed to nitrite [49]. Hence, these immune response-related genes can be activated as part of the defensive system of the copepod to cope with the harmful effects of a diatom diet, potentially associated with lipid peroxidation production induced by PUAs and other oxylipins.
Overall, our results suggest the occurrence of a mechanism of transient impairment of early stress defense response in C. helgolandicus females, induced by short-term feeding on oxylipin-producing diatoms (2 days), followed by later activation of detoxification, protein repair, and immune systems genes during longer feeding on the diatom diet (5 days). This concerted response might explain the high C. helgolandicus female survival observed in the present study (data not shown) and also reported in previous laboratory feeding experiments over fifteen days [7].
Ingestion of S. marinoi and exposure to PUAs also led to deregulation of folate and methionine metabolic processes, as well as of signaling pathways involved in oogenesis/embryogenesis, cellular proliferation and patterning, possibly leading to reduced embryo viability, teratogenesis, and apoptosis.
Folate (vitamin B9), is an essential cofactor for the de novo synthesis of purines and pyrimidines and hence DNA synthesis, methylation of DNA and proteins, through one-carbon metabolism and methionine cycle, respectively [50]. This complex metabolic pathway is well known in vertebrates, where its alteration is associated with increased risk for neural tube defects and other transgenerational developmental defects [51,52]. In arthropods, altered folate metabolism is associated with abnormal body patterning in Drosophila [53] and Artemia larvae [54]. Additionally, maternal exposure to the folate inhibitor drug methotrexate produced leg and wing deformities in surviving progeny of Drosophila [53,55]. To date, the genes involved in folate metabolism are not characterized in crustaceans, except for thymidylate synthase (TS) in L. vannamei [56]. Therefore, this is the first report of folate metabolism genes in a planktonic crustacean.
Our data showed that C. helgolandicus females feeding on S. marinoi or incubated with pure PUAs, had decreased expression level of genes encoding for 10-formyltetrahydrofolate dehydrogenase (10-FTHFDH) and dihydrofolate reductase (DHFR), both involved in replenishment of tetrahydrofolate (THF), the active form of folate [56]. Reduced expression of these genes might lead to decreased levels of THF and, hence, of downstream substrates for the production of 5,10-methylene THF, the donor cofactor for the transfer of one carbon units for purine and pyrimidine biosynthesis catalyzed by TS [57]. Inhibition of DHFR activity in zebrafish has been linked to occurrence of abnormal embryos [51], whereas 10-FTHFDH knock-out in the embryos, lead to delayed early development and subsequent anomalies [58]. Calanus helgolandicus females feeding on S. marinoi and exposed to PUAs also showed up-regulation of the methylenetetrahydrofolate reductase (MTHFR) gene, which is responsible for the irreversible conversion of 5,10-methylene THF into 5-methyl THF. This substrate is converted back into THF by the methionine synthase (MS), via methylation of homocysteine into methionine, thus linking the one carbon folate cycle to the methionine cycle and, hence, to DNA methylation [56]. Over expression of MTHFR in C. helgolandicus, therefore, might increase shuttling of 5,10-methylene THF towards the methionine cycle at the expense of purine and pyrimidine biosynthesis. Reduced conversion of dUMP to dTMP has been associated with increased uracil misincorporation, DNA damage, and apoptosis in humans [59]. In addition, treatment of mouse embryos with the known teratogen valproic acid induced an increase in MTHFR gene expression and was associated with higher risk of congenital malformations [60]. In concert with higher expression of MTHFR gene we also observed up-regulation of MS gene, probably as a response to higher substrate availability, and down-regulation of betaine-homocysteine methyltransferase 1 (BHMT), a likely compensatory mechanism to re-equilibrate methionine levels, being the enzyme responsible for the methylation of homocysteine into methionine using the amino acid betaine [57]. The potential effect of S. marinoi and PUAs on the methylation cycle is further supported by the up-regulation of DNA methyltransferase 1 (DNMT1) gene, involved in the methylation of cytosine residues, possibly leading to increased methylated products. DNA hyper-methylation may affect epigenetic control of gene transcription and, in turn, embryogenesis and development [61].
Taken together, our results suggest that the adverse effects on C. helgolandicus oocytes and developing embryos, induced by maternal exposure to diatoms and oxylipins PUAs [7], could be due to an altered folate and methionine metabolism in the females. Folate-dependent reactions are, in fact, essential for early growth and development in arthropods [62]. Since PUAs accumulate selectively in copepod gonads [63], it is possible that diatom oxylipins may directly target reproductive tissues and, in turn, affect embryogenesis [64].
The present study also showed that several transcripts involved in oogenesis, vitellogenesis, and developmental signal transduction pathways, were modulated in C. helgolandicus females after diatom feeding and direct PUAs incubation. For example, strong down-regulation was observed in both S. marinoi-fed and PUAs-exposed females of the gene coding for the vitelline membrane outer layer protein 1 homolog (VMO1). The vitelline membrane plays a major role in preventing mixing of yolk and albumen in eggs, whose incomplete separation has been associated with the production of abnormal oocytes and delayed development [65]. It also provides positional information to the developing embryos, being the spatial repository for dorso-ventral and antero-posterior embryonic patterning determinants produced by follicle cells during oogenesis [66,67]. In crustaceans, VMO1 proteins are synthesized by extraovarian tissues and then transported via the hemolymph to the developing oocytes [65]. Its role is still not fully known, although it has been suggested that altered expression of this gene in Daphnia magna females exposed to ibuprofen, was associated with failed oogenesis, abnormal oocytes possibly being re-absorbed and blocked embryogenesis [68]. It is therefore possible that reduced expression of VMO1 genes could play a role in the arrested oocyte maturation and oocyte degradation previously observed in C. helgolandicus females feeding on oxylipin-producing diatoms [69].
In most animals, the early embryonic development depends on maternally provided mRNA, which is crucial for cell cycle progression, axis patterning, and cell fate specification events [70]. Hence, maternally encoded transcripts might play a role in diatom-derived effects on C. helgolandcus. The aurora A kinase (AUR) is a maternally supplied cell cycle regulator which plays an important role in cell-cycle progression of fertilized eggs and early embryos of mouse [71]. Silencing of AUR gene leads to asynchronous mitotic cycles and fusion of mitotic spindles leading to larval lethality [72]. Similarly, the maternal transcript distal-less (DLL) functions as a homeodomain transcription factor and plays one of the major roles in limb development throughout the animal kingdom [73]. DLL has been found in insects and crustaceans, where it specifies distal structures and promotes outgrowth of the segmented appendages [74]. Finally, the Drosophila gene dorsal (DORSAL) is a maternal transcription factor essential for the establishment of dorsal-ventral polarity in the developing embryo [75]. Overall, down-regulation of these maternally-encoded transcription factors in C. helgolandicus females feeding on S. marinoi and exposed to PUAs, might thus contribute to impaired embryonic development of early non-feeding nauplii.
As for the signal transduction Hedgehog pathway, this developmental pathway plays a key role in invertebrate development and early embryogenesis [76]. In Drosophila, hedgehog signaling controls patterning of imaginal disc-derived adult structures such as appendages [77]. In the present study, we observed up-regulation of the gene encoding for the extracellular ligand sonic hedgehog protein-like isoform X1 (HH), homologous to hedgehog in Drosophila, and this matched the increase in expression level of the gene encoding for its receptor, patched domain-containing protein 3 (PTCHD3). However, the down regulation of the downstream cytoplasmic signal transducer smoothened (SMO) and transcription factor cubitus interruptus (CI) genes, suggests possible down-regulation of Hedgehog signaling pathway in C. helgolandicus after five days feeding on S. marinoi and PUAs incubation. This is also supported by the observed down-regulation of the DLL gene, one of the targets of Hedgehog. Reduction of DLL activity causes defects of distal leg segments in arthropods, including insects [78] and the crustaceans Parhyale hawaiensis [79], Daphnia magna [80] and Daphnia pulex [81].
In conclusion, the present study indicates that maternal ingestion of S. marinoi and exposure to PUAs by C. helgolandicus for five days leads to de-regulation of folate metabolism, possibly leading to decreased DNA synthesis and altered gene methylation, as well as disturbance of oogenesis/embryogenesis signaling pathways involved in cellular proliferation and patterning, thus, causing reduced embryo viability and teratogenesis. A summary of the genes and processes affected in C. helgolandicus is depicted in Figure 10. These results add new insights to the chemically-mediated ecological interactions between diatoms and copepods, providing novel information on the molecular mechanism of action of oxylipins on maternally-mediated teratogenesis of copepod embryos. Our results will also prompt further exploration of the molecular effects of less studied oxylipins, such as hydroxyacids, epoxyalcohols, and fatty-acid hydroperoxides, on copepod reproduction. These oxylipins have been found to impair copepod reproduction in the laboratory by inducing naupliar apoptosis similar to PUAs [7,13]. Recently, reduced hatching success and increased expression of stress-related genes have also been measured in C. helgolandicus females collected during diatom blooms in the Northern Adriatic Sea, when high oxylipin concentrations were measured in the natural phytoplankton assemblage [8]. Therefore, our results will also contribute to the evaluation of the expression of new biomarker genes associated with diatom-related reproductive failure in the natural copepod population at sea.
(PTCHD3). However, the down regulation of the downstream cytoplasmic signal transducer smoothened (SMO) and transcription factor cubitus interruptus (CI) genes, suggests possible down-regulation of Hedgehog signaling pathway in C. helgolandicus after five days feeding on S. marinoi and PUAs incubation. This is also supported by the observed down-regulation of the DLL gene, one of the targets of Hedgehog. Reduction of DLL activity causes defects of distal leg segments in arthropods, including insects [78] and the crustaceans Parhyale hawaiensis [79], Daphnia magna [80] and Daphnia pulex [81].
In conclusion, the present study indicates that maternal ingestion of S. marinoi and exposure to PUAs by C. helgolandicus for five days leads to de-regulation of folate metabolism, possibly leading to decreased DNA synthesis and altered gene methylation, as well as disturbance of oogenesis/embryogenesis signaling pathways involved in cellular proliferation and patterning, thus, causing reduced embryo viability and teratogenesis. A summary of the genes and processes affected in C. helgolandicus is depicted in Figure 10. These results add new insights to the chemically-mediated ecological interactions between diatoms and copepods, providing novel information on the molecular mechanism of action of oxylipins on maternally-mediated teratogenesis of copepod embryos. Our results will also prompt further exploration of the molecular effects of less studied oxylipins, such as hydroxyacids, epoxyalcohols, and fatty-acid hydroperoxides, on copepod reproduction. These oxylipins have been found to impair copepod reproduction in the laboratory by inducing naupliar apoptosis similar to PUAs [7,13]. Recently, reduced hatching success and increased expression of stress-related genes have also been measured in C. helgolandicus females collected during diatom blooms in the Northern Adriatic Sea, when high oxylipin concentrations were measured in the natural phytoplankton assemblage [8]. Therefore, our results will also contribute to the evaluation of the expression of new biomarker genes associated with diatom-related reproductive failure in the natural copepod population at sea.

Phytoplankton Culture
The centric diatom Skeletonema marinoi (strain FE6) and the dinoflagellate Prorocentrum minimum (strain FE100), were grown in f/2 medium and K medium, respectively. Phytoplankton cultures

Phytoplankton Culture
The centric diatom Skeletonema marinoi (strain FE6) and the dinoflagellate Prorocentrum minimum (strain FE100), were grown in f/2 medium and K medium, respectively. Phytoplankton cultures were grown as semi-continuous batch cultures to late-exponential phase of growth in 2-L glass jars kept in a temperature-controlled room under 18 • C, on a 12:12 Light: Dark cycle and an irradiance of 100 µE m −2 s −1 .

Copepod Collection and Feeding Experiments
Calanus helgolandicus specimens were sorted from zooplankton samples collected in the Gulf of Naples during May 2012 by using 200 µm Nansen net in oblique tows. Adult C. helgolandicus females were isolated under a Leica stereomicroscope, transferred to 1L jars (n = 20-30 copepods) filled with 0.22 µm filtered sea water (FSW) and enriched with either S. marinoi (45,000 cells mL −1 , 1 mg C L −1 ), or P. minimum (5000 cells mL −1 , 1 mg C L −1 ). Jars were kept in a temperature-controlled room under 18 • C, on a 12:12 Light: Dark cycle. The algal medium was changed daily with addition of fresh FSW and either S. marinoi or P. minimum diet. After five days, the copepods were transferred to clean jars containing 0.22 µm FSW for 24 h in order to allow gut evacuation. Three independent experiments were performed on different occasions during May 2012, each one with both the diatom and the control treatment and were considered as biological replicates. Ten females per replicate diet were transferred into 1.5 mL Eppendorf tube containing 500 µL of RNAlater ® reagent and processed according to manufacturer's instructions.

Transcriptome Sequencing
Total RNA was extracted from C. helgolandicus females using the RNeasy Micro Kit (Qiagen, Germany) following manufacturer's procedure [82]. RNA concentration (ng µL −1 ) and purity were assessed through Nanodrop ND-1000 UV-Vis spectrophotometer (Marshall Scientific, Hampton, NH, USA), whereas RNA Integrity Number (RIN) was checked on a 6000 Nano LabChip of Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). A total of 3 µg of RNA per sample (300 ng µL −1 ) was delivered to the Genomics Core Facility of the European Molecular Biology Laboratory (EMBL, Heidelberg, Germany), for library preparation using TruSeq RNA Sample Prep Kit (Illumina, San Diego, CA, USA). Fragments >200 bp were selected, purified and subsequently PCR amplified to create the final cDNA library template for sequencing. RNA-Seq was conducted on Illumina HiSeq 2000 platform with 50 bp paired end option. Two cDNA libraries obtained from copepods fed with either S. marinoi or P. minimum from one replicate experiment, were sequenced individually in two single lanes of the Illumina flow cell to obtain a deep and high qualitative coverage of C. helgolandicus transcriptome. The other four cDNA libraries obtained from copepods fed with either S. marinoi or P. minimum during two other independent experimental replicates were used for multiplexed sequencing in one single lane of an Illumina flow cell.

De novo Transcriptome Assembly and Functional Annotation
Cleaning, trimming, quality filtering and removal of the adapters was performed with the program Trimmomatic [83] combined with custom scripts [84] with the following options: ILLUMINACLIP set at 2:40:15 with the use of standard Illumina adapters sequences used as filter, LEADING and TRAILING set at 5, SLIDINGWINDOW at 5:20 and MINLEN at 30 (Script can be found here: https://github.com/silverkey/transcriptome/blob/master/preprocess_fastq_for_trinity.pl).
The remaining high quality paired-end reads across the six samples were assembled using Trinity [85] (version 2013-02-25) with default plus the following parameters: -seqType fa -JM 200G -inchworm_cpu 20 -bflyHeapSpaceInit 20G -bflyHeapSpaceMax 200G -bflyCalculateCPU -CPU 20 &. This 'reference' C. helgolandicus transcriptome consisted of de novo assembled Trinity transcripts with unique TR#_c#_g# identifiers ('Trinity predicted genes' or unigenes) and contained either singletons (transcripts with a single isoform, 'i') as well as the longest isoform of transcripts having multiple 'Trinity predicted isoforms' (TR#_c#_g#_i#) [86]. To exclude contaminations, poorly supported transcripts and artifacts, the assembled unigenes were further filtered based on a minimum expression level of more than 1 read per million mapped reads (RPKM) in at least 2 samples.
Unigenes were further assigned a putative protein function by a sequence similarity search using the BLASTx algorithm (e-value ≤ 10 −3 ) against the NCBI non-redundant protein sequence database (nr) and Swissprot database, followed by Gene Ontology (GO) functional annotation (e-value ≤ 10 −6 ), using Blast2GO PRO version 5 [87]. The Web Gene Ontology Annotation Plot (WEGO) software [88] was then used to plot and compare the GO annotation results among related organisms. Additionally, Kyoto Encyclopedia of Genes and Genomes (KEGG) Automatic Annotation Server (KAAS) database [89] was also searched using BLASTx algorithm (e-value ≤ 10 −5 ), to provide information about metabolic pathways in the dataset.

Differential Gene Expression Analysis
To identify differentially expressed unigenes between S. marinoi-fed and P. minimum-fed C. helgolandicus cDNA libraries, normalized raw reads from each replicate feeding treatment were firstly mapped back on the reference transcriptome using BOWTIE [90] and then counted using RSEM software [91] to estimate the gene expression level (https://github.com/silverkey/transcriptome/blob/ master/launch_mapping_trinity_analysis_folder.pl). The three replicate samples from each treatment were used to generate mean expression levels using RPKM method. Statistical analysis was performed using R/Bioconductor and the EdgeR package to identify differentially expressed genes (DEGs). Significance values were obtained by performing a hypergeometric test and corrected p-value using the false discovery rate (FDR) method [92]. Genes having a FDR ≤ 0.05 were considered differentially expressed and further annotated against nr and Swissprot using Blast2GO, and against KEEG database using KAAS, according to the procedure described for the reference transcriptome. Functional enrichment analysis for up-and down-regulated genes in S. marinoi-fed versus P. minimum-fed C. helgolandicus females, with respect to all DEGs, was also performed by Blast2GO using the Fisher's exact test with multiple testing correction of false discovery rate (FDR ≤ 0.05).

RT-qPCR of Genes of Interest (GOIs)
A series of functionally annotated genes of interest (GOIs) were selected from the list of DEGs according to their FDR and log2 fold change and used for validation of RNA-Seq analysis through real time-quantitative PCR (RT-qPCR). Other functionally annotated GOIs were also selected from the de novo assembled transcriptome, according to their putative involvement in the response of copepods to diatoms and tested by RT-qPCR. Specific forward and reverse oligonucleotide primers were designed using Primer3 software (v. 0.4.0) as in [27] and synthesized by Primm Labs (Milan, Italy).
RT-qPCR reactions of reference genes as well as GOIs, were performed in MicroAmp Optical 384-Well reaction plate (Applied Biosystem, Foster City, CA, USA) with optical adhesive covers (Applied Biosystem), using a Viia7 real-time PCR system (Applied Biosystem, Foster City, CA, USA). First, cDNA template was retrotranscribed from 1 µg of total RNA remaining from RNA-Seq analysis, using iScriptTM cDNA Synthesis Kit (BIORAD) and following the manufacturer's instructions, in the GeneAmp PCR System 9700 (Applied Biosystems). cDNA template (1 µL at 1:100 dilution) was then mixed with 5 µL of Fast Start SYBR Green Master Mix (Applied Biosystem) and 0.7 pmol/µL of each primer, for a final PCR sample volume of 10 µL. RT-qPCR reactions were carried out in triplicate, including at least two no-template negative controls for each primer pair, using PCR conditions previously reported in [30]. Serial dilutions of cDNA (1, 1:5, 1:10, 1:50, 1:100 and 1.500) were used to calculate reaction efficiencies (E) for all primer pairs using the equation E = 10 −1 /slope from a standard curve between Ct values and the log10 for each dilution factor. The relative expression ratio (R) of each GOI in the experimental conditions (copepods fed on S. marinoi) versus the control condition (copepods fed on P. minimum) and its statistical significance was calculated using REST (Relative Expression Software Tool) and the Pair Wise Fixed Reallocation Randomization Test [93]. The calculation is based on the efficiencies (E) and the Ct deviation between the experimental and the control groups of the target genes normalized to the reference genes EFA, GAPDH and S20, and is expressed as log2 fold change. Differences were considered significant when Randomization Test p-value was ≤0.05.

PUAs Incubation Experiments and RT-qPCR of GOIs
Calanus helgolandicus adult females were collected in the Gulf of Naples during March-May 2013 and fed with P. minimum (5000 cells mL −1 ) before the start of the experiment. Stock solutions of the polyunsaturated aldehydes (PUAs) 2-trans,4-trans-heptadienal and 2-trans,4-trans-octadienal (Sigma-Aldrich, Italy) at 10 mM and 1 mM were prepared by diluting the proper amount of PUAs in absolute methanol (J.T. Baker, Holland). Final working solutions were obtained by diluting stock solutions in FSW and mixing carefully to ensure optimal distribution. Methanol has no negative effect on copepod females and eggs up to 1% methanol in seawater [7]. The amount of aldehyde solution in each test concentration was kept below this threshold (final concentration used: 0.2% of methanol). PUA incubation experiments were performed by incubating individual C. helgolandicus female (n = 5−7) in 100-mL crystallizing dishes containing P. minimum (5000 cells mL −1 ) and 1:1 PUA mixture of heptadienal + octadienal (MIX) at final concentration of 10-15-20 µM in FSW. Females incubated in P. minimum + methanol were used as controls. Females were kept at 18 • C in a controlled temperature chamber under 12:12 Light: Dark cycle and transferred to new crystalizing dishes with fresh PUA mixture and algal diet daily. After 5 days of incubation, females were placed in 500 µL RNAlater ® (5-7 females per treatment), frozen according to the manufacturer's instructions and stored at −80 • C until RNA extraction. A total of 3 experiments were performed during March-May 2013. Extraction of total RNA, retrotranscription of cDNA and RT-qPCR analysis was performed as previously described. Stability of reference genes was screened again in copepods exposed to PUAs by RefFinder and the most stable genes identified were S20, GAPDH and UBI. Relative expression ratio (log2 fold change) of the same GOIs selected previously were calculated between C. helgolandicus females incubated in each PUAs MIX concentration (experimental condition) and C. helgolandicus females incubated in P. minimum + methanol (control condition) using REST. Data are expressed as log2 fold change.

Statistical Analysis
Pearson's correlation test was performed between log2 fold change of DEGs from RNA-Seq analysis and log2 fold change of the same genes in RT-qPCR, using GraphPad Prism software v.6 (GraphPad Software Inc., San Diego, CA, USA). The same software was also used to perform Student's t-test or One-way analysis of variance (ANOVA), followed by Tukey's multiple comparison test, to evaluate statistically significant differences in log2 fold change among 10-15-20 µM PUAs MIX treatments for each GOIs.