RNA-Seq and differential gene expression analysis in Temora stylifera copepod females with contrasting non-feeding nauplii survival rates: an environmental transcriptomics study

Copepods are fundamental components of pelagic food webs, but reports on how molecular responses link to reproductive success in natural populations are still scarce. We present a de novo transcriptome assembly and differential expression (DE) analysis in Temora stylifera females collected in the Gulf of Naples, Mediterranean Sea, where this copepod dominates the zooplankton community. High-Throughput RNA-Sequencing and DE analysis were performed from adult females collected on consecutive weeks (May 23rd and 30th 2017), because opposite naupliar survival rates were observed. We aimed at detecting key genes that may have influenced copepod reproductive potential in natural populations and whose expression was potentially affected by phytoplankton-derived oxylipins, lipoxygenase-derived products strongly impacting copepod naupliar survival. On the two sampling dates, temperature, salinity, pH and oxygen remained stable, while variations in phytoplankton cell concentration, oxylipin concentration and oxylipin-per-diatom-cell production were observed. T. stylifera naupliar survival was 25% on May 23rd and 93% on May 30th. De novo assembly generated 268,665 transcripts (isoforms) and 120,749 unique ‘Trinity predicted genes’ (unigenes), of which 50% were functionally annotated. Out of the 331 transcript isoforms differentially expressed between the two sampling dates, 119 sequences were functionally annotated (58 up- and 61 down-regulated). Among predicted genes (unigenes), 144 sequences were differentially expressed and 31 (6 up-regulated and 25 down-regulated) were functionally annotated. Most of the significantly down-regulated unigenes and isoforms were A5 Putative Odorant Binding Protein (Obp). Other differentially expressed sequences (isoforms and unigenes) related to developmental metabolic processes, protein ubiquitination, response to stress, oxidation-reduction reactions and hydrolase activities. DE analysis was validated through Real Time-quantitative PCR of 9 unigenes and 3 isoforms. Differential expression of sequences involved in signal detection and transduction, cell differentiation and development offered a functional interpretation to the maternally-mediated low naupliar survival rates observed in samples collected on May 23rd. Down-regulation of A5 Obp along with higher quantities of oxylipins-per-litre and oxylipins-per-diatom-cell observed on May 23rd could suggest oxylipin-mediated impairment of naupliar survival in natural populations of T. stylifera. Our results may help identify biomarker genes explaining variations in copepod reproductive responses at a molecular level.

In the last two decades, extensive evidence was reported about impaired reproductive success in copepod females fed oxylipin producing diatoms, which led to detrimental effects on egg production rates, egg hatching and survival of non-feeding nauplii (NI/NII) through a maternal effect [9,[17][18][19][20][21][22][23][24][25]. Since 2011, a number of studies have started to inspect the effects of oxylipin producing diatoms on the molecular responses of copepod females, evaluating variations in the quantitative expression of selected genes of interest [26][27][28][29][30][31] and applying a suppression subtractive library approach to gain insight into copepod responses at a transcriptomic level [32]. Very recently, the de novo assembled transcriptome of copepod females feeding on oxylipinproducing diatoms has been also generated [33].
In the present survey, we investigated the molecular responses of adult females of the calanoid copepod Temora stylifera from the Gulf of Naples (GoN), where it dominates the autumnal copepod community [40,[43][44][45][46]. The GoN has been traditionally described as an oligotrophic basin showing low phytoplankton densities and consequent low oxylipin concentrations. However, we recently showed that high oxylipin-per-litre concentration and oxylipin-perdiatom-cell productions seasonally occur in this area [47]. Several studies have already investigated the population dynamics of T. stylifera in th GoN, exploring whether abiotic factors and life-history traits could explain the marked seasonality of this copepod in the area [40,44,48]. However, no genomic and transcriptomic information are available for this species.
Through a High-Throughput Sequencing approach, we generated a de novo assembled transcriptome of adult T. stylifera females. We also performed a Differential Expression (DE) analysis between specimens collected on two consecutive weeks (the 23rd and the 30th of May 2017), when early non-feeding nauplii with opposite survival rates (25% vs 93%, respectively) were laid.
Analyses of de novo assembled transcriptomes were reported to explore the biosynthetic pathways of gaseous signals [49], the enzymatic processes leading to hormone biosynthesis [50], reproductive processes [51,52], including diapause [53][54][55][56] as well as responses to stress [57,58] and phycotoxins [59] in several pelagic copepod species. Our results offer the opportunity to understand if molecular responses of T. stylifera females from natural populations can help to better explain different naupliar survival rates in relation to environmental (temperature, salinity, pH and oxygen), biological (phytoplankton abundance and composition) and biochemical (phytoplankton-derived oxylipins) variables [47].

Environmental, chemical and biological variables
Information about abiotic (temperature, salinity, pH and oxygen), phytoplankton and oxylipin variations on the two selected dates are reported in Table 1. Abiotic variables did not show wide variations between the two sampling dates. In contrast, more pronounced variations were detected in phytoplankton community abundance and composition, when considering major phytoplankton groups (i.e. coccolithophores, dinoflagellates and phytoflagellates < 10 μm) as well as the most abundant diatom genera (i.e. Chaetoceros, Skeletonema, Leptocylindrus, Pseudo-nitzschia, Thalassiosira and the mixed group "other diatoms"). In general, phytoplankton was less concentrated on the 23rd of May (14.74 10 6 cells/L) than the 30th of May (17.64 10 6 cells/L). In particular, coccolithophores, dinoflagellates and the diatom genus Chaetoceros occurred at higher concentrations on the 23rd of May than the 30th, while higher densities of phytoflagellates < 10 μm and of the "other diatoms" group were observed on the 30th of May than the 23rd.
T-test results demonstrated that early-life history traits estimated for T. stylifera on the two sampling dates (May 23rd and 30th) differed significantly in terms of survival rates of NI nauplii (25 and 93% of survival, respectively, p < 0.001, N = 15, Fig. 1). By contrast, nonsignificant differences (p > 0.01, N = 15) were observed in the number of faecal pellets (an indirect measure of feeding rates) (61.6 ± 3.39 and 72.35 ± 4.13 pellets per female per day, respectively), the number of spawned eggs (62.8 ± 11 and 75.36 ± 11.19 eggs per female per day, respectively) and the percentage of egg hatching success (63.4 ± 12.3% and 89.95 ± 4.05%, respectively).

De novo assembly and functional annotation of Temora stylifera transcriptome
Total RNA extracted from pools of T. stylifera females collected on May 23rd and 30th had an average concentration of 232.7 ng/μl, with RIN = 10 and 260/280 as well as 260/230 ratios~2. Illumina-based RNA-Seq generated a total of~132 million reads, after quality cleaning. The same number of reads was achieved for both the forward and the reverse cDNA filaments, supporting consistency in the sequencing output. Raw reads are stored into the NCBI Sequence Read Archive database under accession numbers PRJNA632714. The de novo assembly made with Trinity on high quality reads generated 268,665 transcripts (isoforms) (average length of 517.6 bp, N 50 = 665), and contained 120,749 'Trinity predicted genes' (unigenes), i.e. non-redundant transcripts with unique TR#_c#_g# identifiers (Additional file 1: Table S1). Both the full (transcript isoforms) and the reference (unigenes) transcriptome, the latter consisting of the longest transcript isoform of each predicted gene, were processed for functional annotation. However, detailed description of annotation results is here provided only for the reference transcriptome.
Blast2Go mapping outputs indicated that almost 10% of the matching unigenes showed very high homology (0 < E-value< 10 − 100 ) to similar sequences in the nonredundant protein database. Overall, more than 42% of the sequences showed high probability of homology (0 < E-value< 10 − 30 ). Similarity values, which express the similarity percentage between the de novo assembled sequence and proteins in the non-redundant database, highlighted that a low fraction (1.7%) of the total  Fig. S1). The species distribution of the best matches (top-hit) against the non-redundant protein database indicated that the largest fraction of matching unigenes showed similarities with sequences of the copepod Eurytemora affinis, followed by the copepod Acartia pacifica, the cladocerans Daphnia pulex and Daphnia magna and the copepod Pseudodiaptomus poplesia. The other top-hit species were mainly crustaceans or arthropods, while three molluscs and one brachiopod were among the other first 20 top-hit species (Additional file 2: Fig. S1).
Blast2Go annotation outputs showed that 31,346 unigenes, out of 62,648 that received significant matching in BLASTx, were functionally annotated (50.04%). In total, 126,358 GO annotation terms were assigned and distributed among GO categories for Biological Process (BP, 36.77%), Molecular Function (MF, 35.57%) and Cellular Component (CC, 27.66%) (Fig. 2). The majority of recognized unigenes were assigned to metabolic and cellular processes (29%), binding and catalytic activity (49.59 and 32.55%, respectively) and cell or cell part (both 20%).

Differential expression analysis and transcriptome validation
Analysis of expression levels of T. stylifera unigenes between samples collected on May 30th and May 23rd showed that a total of 144 unigenes were differentially expressed (FDR p ≤ 0.05). In particular, on May 23rd 108 genes were down-regulated, while 36 genes were upregulated. Of the total 144 differentially expressed sequences, 31 (6 up-regulated and 25 down-regulated) received GO assignment and functional annotation ( Table 2).
In order to have a wider spectrum of gene functions and to allow a more detailed description of the molecular responses of T. stylifera females on the two sampling dates, differential expression analysis was performed also on transcript isoforms. Among isoforms, 331 sequences were differentially expressed (FDR p ≤ 0.05), 199 were down-regulated and 132 were up-regulated. In total, 119 differentially expressed isoforms received GO assignment and were functionally annotated (58 up-and 61 down-regulated) (Additional file 3: Table S2). In total, 563 GO terms were associated to the differentially expressed isoforms and were assigned to the three main GO categories, which were almost equally divided among BP (37.18%) and MF (34.66%), while a smaller fraction described the CC category (28.16%). Analysis of GO distribution among the three main categories was also repeated dividing up-and down-regulated isoforms (Fig. 3). Results showed similar number of up-regulated and down-regulated sequences in the different GO in terms of BP, MF and CC categories.
Consistent differences between replicates collected on May 30th and May 23rd were supported by clustering among objects (Q-mode analysis) described by raw Fig. 1 Temora stylifera responses. Average daily faecal pellet and egg production (N per female per day) measured in adult females as well as average egg hatching success and NI naupliar survival (%) for the two sampling dates (May 23rd and May 30th 2017). Differences in production or percentage were analysed through t-test (95% confidence interval). Significance level: *** < 0.001 counts of both the differentially expressed unigenes and isoforms (Fig. 4). Also, such clustering was confirmed when raw counts of isoforms involved in specific molecular pathways were selected and analysed separately (Fig. 4).
Most of the significantly down-regulated unigenes and isoforms on May 23rd were described as A5 Putative Odorant Binding Protein (Obp, annotated as sequences involved in 'response to stimulus'). Other downregulated unigenes were annotated as sequences related to developmental metabolic processes (involving chitin and collagen), protein ubiquitination, response to stress (mainly Heat Shock Protein 70), oxidation-reduction reactions and hydrolase activities. Similarly, additional down-regulated isoforms were involved in respiration, protein binding, transmembrane transport and cellular development. The significantly up-regulated unigenes and isoforms were mainly involved in reproduction, cell development and proliferation (e.g. Vitellogenin-like unigenes, RNA Helicase and Lipoprotein Receptor unigenes), transmembrane transport and reception activity.
Based on these results, 9 unigenes and 3 isoforms were selected as Gene of Interests (GOIs) for transcriptome validation through RT-qPCR analysis depending on function, fold-change, significance (adjusted p-values), sequence length, E-value and sequence similarity percentage. Although unigenes offered a narrower array of functions in comparison to isoforms, most primers were selected from unigenes to reduce redundancy due to multiple transcript isoforms within the same Trinity gene cluster. Amplicons were all in the range of 111-228 bp and showed primer amplification efficiencies between 1.9 and 2.1. The full list of primer sequences for these selected sequences of interest is shown in Table 3.
For RT-qPCR analysis, the expression of GOIs was normalized considering 18S ribosomal RNA (18S) and Ubiquitin (Ubi) as reference genes (Additional file 4: Table S3). These two genes were indicated as the most stable ones among the five candidates selected as potential reference sequences according to results provided by RefFinder [60] (mean of ranking values: 1 and 1.68, respectively).       Normalized expression of the selected GOIs for transcriptome validation supported RNA-Seq results, because 9 sequences out of 12 reflected the same up-or down-regulation patterns as in the differential expression analysis (Fig. 5). However, the genes Protein Obstructor E (Obst) and Heat-Shock Protein 70 (Hsp70) as well as the isoform MOB Kinase Activator 1B (Mob1b) showed an opposite trend in comparison to RNA-Seq results. In general, log 2 -fold change indicating differential expression in T. stylifera females collected on May 23rd in comparison to those collected on May 30th was larger in the RNA-Seq output than RT-qPCR results. In particular, RNA-Seq log 2 -fold-changes of the three selected isoforms (Mob1b, Vasa and Pafah; log 2 -FC = − 9.48, 7.14 and 6.93, respectively) were much higher than

Discussion
Estimating early-life history variables such as egg production, hatching success and survival rates of nonfeeding (NI/NII) nauplii in field and laboratory studies has been traditionally considered to infer recruitment potential in copepod populations [18,34,40,42]. In this perspective, understanding which molecular mechanisms contribute to define naupliar viability may help predict reproductive responses of natural copepod populations.
In particular, our DE analysis may allow elucidating molecular mechanisms affecting naupliar viability in T. stylifera females collected from natural populations in the Gulf of Naples (GoN). Information about environmental variables, phytoplankton community composition and abundance, oxylipin-per-litre concentrations and oxylipin-per-diatom-cell production can elucidate the influence of these variables on molecular responses of copepod females and final naupliar viability. While on a side environmental variables were proposed as major factors affecting T. stylifera population in the GoN [44,48], on the other side wide evidence has been provided about the detrimental effect of phytoplankton-derived oxylipins on the reproductive success of grazer copepods [9, 17-22, 24, 36], possibly in relation to differential expression of key genes involved in reproduction and stress responses [3,[26][27][28][29][30][31][32][33]35].
The present transcriptome analysis resulted in 268,665 transcripts (isoforms) and 120,749 predicted genes (unigenes). These numbers were sensibly higher than results reported in other calanoid copepods, such as Calanus helgolandicus [33], C. finmarchicus [61], Calanus sinicus [62], Temora longicornis [57] and Acartia tonsa [58]. Rather, our de novo transcriptome assembly for T. stylifera was much closer to results presented by [52] for C. finmarchicus, because these authors reported 241,140 transcripts and 124,157 components (corresponding to our Trinity predicted genes). The vast majority of the 20 top-hits BLAST species for the present transcriptome was represented by crustaceans and the first two top-hit species were the copepods Eurytemora affinis and Acartia pacifica. High affinity of the blasted copepod sequences to crustaceans and arthropods is increasing in the last years [33,57,63] and this indicates that publicly available copepod genomic resources are improving. Nonetheless, only nearly 50% of T. stylifera transcripts and predicted genes found hits in BLASTx and around 25% of these sequences received GO annotation. These low percentages demonstrate that much effort is still needed to improve genomic resources for copepods, even if the present results suggest improvements from previous de novo transcriptomic analyses [61,62]. Annotation output of T. stylifera transcriptome showed that most unigenes were annotated in metabolic and cellular processes BP sub-categories and in cell or cell part CC subcategories, in good agreement with previous transcriptomic studies on calanoid copepods [33,49,57,58,62].
After an overview of the de novo transcriptomic assembly and annotation, we inspected if our differential analysis could help infer maternally-related molecular processes involved in naupliar survival. We compared expression rates in T. stylifera adult females generating highly viable nauplii (on May 30th) with those of females   generating nauplii with low survival rates (on May 23rd) and found that differentially expressed sequences were mostly involved in signal perception and transduction as well as in reproductive and developmental processes. More in particular, most down-regulated transcripts on May 23rd were Putative Odorant-Binding Proteins. Most of them were A5-Binding Proteins (belonging to the Phosphatidylethanolamine-Binding Protein class), whose expression has been found in the palps and antenna of the mosquito Anopheles gambiae [64]. Odorant binding proteins (Obps) are part of the arthropod olfactory system, which has been extensively characterized in the model fly Drosophila melanogaster [65]. In this organism, suppression of Obp expression mediated an altered behavioural response [66] and modulation in the ingestion of bitter tastants [67]. In particular, inhibition of Obp A5 in D. melanogaster females led to differential responses to benzaldehyde and citral, to increased ingestion of N-phenylthiourea and papaverine as well as to reduced consumption of caffeine and denatonium benzoate. These results support that Obps in Drosophila are involved in detection of harmful chemicals. In fact, Obps were found to mediate a negative response (i.e. avoidance) of Drosophila sechiella to phyco-toxins, because the Obps 57d and 57e allowed recognition of the plant-derived hexanoic acid and octanoic acid [68,69]. One of the most important features highlighted in the Obp modulation as response to olfactory stimuli in D. melanogaster is that Obp gene expression is combinatorial [66,67]. This means that the response to an odorant stimulus is regulated by the expression of multiple Obps. Copepods base their feeding behaviour on mechanosensory setae, chemosensory sensilla (i.e. aesthetascs) or bimodal sensilla, present on antennules A1, antennas and the maxilliped. Such sensory organs allow copepods to distinguish prey morphology, to detect foraging stimuli and avoid food items on the basis of chemicals [70]. In spite of combinatorial regulation as in D. melanogaster, our differential analysis revealed down-regulation of only A5 Obps along with those sequences described as OV-16. This result could suggest specific involvement of A5 Obp in the response to algal-derived toxic compounds in T. stylifera females. Despite major differences in signal perception between aquatic and non-aquatic environments [71,72], Obps in arthropods likely encode for proteins displaying similar functions. Several experiments have reported different copepod behaviours in response to specific preys or compound classes [73][74][75]. In particular, dissolved polyunsaturated aldehydes (PUAs), a class of volatile oxilipins, induced attraction of T. stylifera females in odour-choice experiments, suggesting that the copepods showed a behavioural preference for oxylipins [73]. Total oxylipin-per-litre concentration and oxylipin-perdiatom-cell production were both higher on May 23rd than 30th. Down-regulation of A5 Obps in T. stylifera could potentially suggest that detection capacity of phytoplankton-derived harmful molecules was reduced in copepod females collected on May 23rd, leading to unselective feeding, higher ingestion of oxylipin producing diatoms and negative effects on larval survival.
In contrast to wide functional information of Obp regulation in D. melanogaster, to the best of our knowledge this is the first time that differential expression of these sequences is reported in a copepod. Downregulation of A5 Obp in T. stylifera females was confirmed by RT-qPCR results (log 2 -FC in RNA-Seq = − 4.6; log 2 -FC in RT-qPCR = − 4.55) and could have played a relevant role in determining ingestion of harmful preys and low reproductive output on May 23rd. These results suggest that Obps can represent important target genes and proteins to investigate how copepods respond to phytoplankton-derived chemicals.
In addition to sequences involved in signal perception, differential analysis supported altered functions in signal transduction, because four sequences classified in the signal transduction BP sub-category were downregulated. The 1-Phosphatidylinositol 3-Phosphate 5-Kinase (Fab1b) showed the strongest down-regulation (log 2 -FC of − 10.34). This enzyme catalyses phosphorylation of the Phosphatidylinositol 3-Phosphate to synthesize Phosphatidylinositol 3,5-Bisphosphate, PI (3,5) P2. This molecule is one of the seven regulatory Polyphosphoinositides (PPIn) that are ubiquitous in eukaryotes, where they regulate membrane trafficking in endosomes and lysosomes [76]. Ablation of the Phosphoinositide kinase Pikfyve (the FAB1B counterpart in mammals) and consequent depletions in intracellular levels of PI (3,5) P2 have been demonstrated to induce embryonic and neonatal lethality in mammals [77]. In fact, PIP2 acts precursor for Phosphatidylinositol triphosphate (PIP3), a membrane lipid acting as a docking site for proteins involved in cell survival, proliferation and differentiation [78].
In support to the down-regulation of Fab1b, the significant up-regulation of the Platelet-Activating Factor Acetyl Hydrolase (Pafah) suggests repression of the PPIn pathway. In fact, PAFAH is a particular type of Phospholipase A2 that deacetylates the Platelet-Activating Factor (PAF), inducing loss in its activity [79,80]. PAF is one of the most potent lipid mediators and occurs in membrane phospholipids. In mammals, PAF is supposed to regulate reproductive cycle and pregnancy [80] by initiating the enzymatic reactions leading to activation of PIP3, which determines high survival of offspring [78]. PIP3 also induces calcium transients that activate the cAMP-responsive element-binding protein (CREB) through phosphorylation. CREB is known to be a fundamental regulator of the signalling mechanisms in the cell and is fundamental during embryonic development [78]. In line with down-regulation of Fab1b and up-regulation of Pafah, Crebl2 (cAMP-Responsive Element-Binding Protein-Like 2) was down-regulated in copepod females collected on May 23rd, as also validated by RT-qPCR results. Therefore, differential expression of this gene could have key implications for cell cycle and cell differentiation [78,81]. Additionally, PAFAH can target phospholipids and participate to oxidative-stress responses as well as to regulation of fertility and apoptosis [82]. In humans, the glycosylphosphatidylinosotol (GPI) is known to induce pro-inflammatory responses [83,84], thus differential expression of Fab1b and Pafah could potentially suggest also dysregulation of the inflammatory response as a consequence of an oxidative stress potentially mediated by oxylipins. Also in this perspective, a synergic response of PAFAH, FAB1B and CREBL in T. stylifera females could have contributed to determine low naupliar survival on the 23rd of May.
Low naupliar survival rates could also relate to downregulation (validated also by RT-qPCR) of the Serine/ Threonine Protein Phosphatase 2A (Ppa2), because the important role of this enzyme in regulation of apoptosis has been documented [85]. Serine/Threonine phosphatases, such as PP1, PP2A and PP2B, mediate dephosphorylation of Bcl-2/Bcl-X-associated death promoter (BAD), a pro-apoptotic protein, thus leading to cell death. PP2A activity can be regulated by Phosphodiesterases (PDEs), which repress PP2A activity through regulation of cAMP levels [86], influencing oocyte maturation in pre-ovulatory follicles and meiosis in oocytes [87,88]. As highlighted by our DE analysis, Pde was down-regulated and this response could have led to higher cAMP levels and stronger stimulation of Ppa2 as compensatory mechanism to re-equilibrate dysregulated apoptosis.
Differentially expressed sequences directly involved in reproduction and development processes may also explain strong differences in naupliar survival rates at the molecular level. The enzyme Methylenetetrahydrofolate Dehydrogenase (Mthfd) was listed among the upregulated transcripts (log 2 -FC of 4.45). Folate onecarbon metabolism plays a fundamental role in development and disease response of animals [89,90]. The folate metabolism is complex, because it involves the activity of several enzymes. The active folate form is the tetrahydrofolate (THF) [91], which is a cofactor involved in the biosynthesis of thymidylate, purines, glycine, serine and homocysteine [92]. The enzyme methylenetetrahydrofolate dehydrogenase (MTHFD) converts tetrahydrofolate (THF) into 10-Formyl, 5,10-Methenyl and 5,10-Methylene derivatives which are key cofactors in the de novo synthesis of DNA [93]. For example, altered folate expression has been shown to drastically reduce fecundity of D. melanogaster females and to induce abnormal development in offspring [90]. Also, dysregulation of the gene expression pathway associated to folate metabolism was recently reported in Calanus helgolandicus females feeding on the oxylipin-producing diatom Skeletonema marinoi [33].
In this perspective, altered expression of Mthfd could potentially contribute to explain impairment of the final reproductive success of T. stylifera in our survey.
Also, the Trinucleotide repeat-containing gene 18 (Tnrc18) was significantly down-regulated. This protein is involved in cell differentiation and its contribution to development of copepod larvae may be important, because this gene is known to have a chromatin binding activity and has been proposed as a new marker gene for development in zebrafish embryos [94]. Additionally, Dermatopontin (Dpt) was also down-regulated. This gene encodes for a protein involved in collagen fibril orientation [95]; therefore its down-regulation and consequent impaired function could have led to negative effects on correct tissue formation in T. stylifera nauplii [96].
Previous laboratory experiments have reported variations in the expression of αand β-Tubulin in copepods fed oxylipin-producing diatoms [3] and similar responses were quantified from natural populations of the copepod C. helgolandicus [35]. β-TUBULIN affects spindle formation and cell division; the up-regulation of this transcript could suggest an active response of T. stylifera females to high oxylipin concentrations potentially impairing reproductive success. In line with β-Tubulin, also the Protein Maelstrom (Mael) was upregulated. MAEL has been reported as a required factor for the correct positioning of the microtubule-organizing centre (MTOC), in D. melanogaster [97]. Moreover, not only does MAEL act as a repressor for micro-RNA-7 in the nucleus, guaranteeing proper differentiation of germline cells [98], but also contributes to oocyte determination [99].
Other sequences relevant for cell development and reproduction were up-regulated in T. stylifera collected on May 23rd, such as me31b (ATP-Dependent RNA Helicase), Vldlr (Very Low-Density Lipoprotein Receptor) and Vasa (ATP-Dependent RNA Helicase Vasa). These genes are known to play key roles in germ cell formation [100,101] as well as in oocyte development and differentiation [102][103][104].
Down-regulation was observed also for other sequences encoding for key functions in development and reproduction. Among them, Obst-E (Protein Obstructor E) is known to encode for a chitin-binding protein [105] and is crucial for the correct cuticle development in arthropods [106]. In line with down-regulation of Obst-E, the unigenes Kelk 12 (Klhl12) and Arih1 (E3 Ubiquitin Protein Ligase, validated through RT-qPCR) were also down-regulated. Klhl12 mediates ubiquitination and regulates the Wnt signalling pathway, thus playing a key role in collagen export and embryonic cell development [107]. Similarly, Arih1 codifies for an E3 Ubiquitinprotein ligase and, together with Obst-E, was included by [51] in the list of the 31 genes involved in reproduction, growth and development in the calanoid copepod Eurytemora affinis, where ubiquitin genes were reported to affect gametogenesis. The down-regulated transcript Hsp-70 (Heat-Shock Protein 70) has close affinity with E3 Ubiquitin-protein ligase and it has been demonstrated to have a critical role in preventing apoptosis [108]. Hsp-70 expression is also involved in thermo-tolerance mechanisms as well as in protection to xenobiotic exposure and to general stress in copepods [29,[109][110][111]. Previous laboratory studies have reported down-regulation of Hsp-70 in the copepod C. helgolandiscus fed oxylipin-producing diatoms [26,28,112], suggesting that down-regulation of Hsp-70 in T. stylifera could also have occurred in response to these chemicals.
Differential analysis also suggested altered functions in protein endocytosis and polysaccharide degradation, which we hypothesize could have contributed to the low naupliar viability observed on May 23rd. Among these sequences, the Adaptor Protein Complex 2 (Ap-2) showed a significant down-regulation (log 2 -FC of − 7.41). Ap-2 is involved in clathrin-dependent endocytosis in which proteins are incorporated into vesicles surrounded by clathrin (CVV) [113]. Ap-2 is therefore strictly related to both the endosomal and the lysosomal systems [114] and is essential to fundamental cellular processes [113]. Also, Arsb (Arylsulfatase B), mainly related to digestion of polysaccharides [115], was a an additional down-regulated unigene and such differential expression was supported by RT-qPCR results. ARSB plays a central role in degradation of glycosaminoglycans (GAG) and altered regulations of this gene can result in lysosomal excretion of polysaccharides [116].
In final instance, several differentially expressed transcripts involved in the oxidative phosphorylation chain, such as Cytochrome oxidase subunit I and Cytochrome oxidase subunit II, were also down-regulated (log 2 -FC ranging from − 2.73 to − 5.75). Considering that oocyte mitochondrial dysfunction in terms of respiration and electron transport can alter fertility and embryo development [117,118], this response is in line with the low naupliar viability observed in samples collected on May 23rd. The potential effects of phytoplankton-derived oxylipins on these transcripts is plausible, because significant down-regulation of the Cytochrome P450-4 (Cyp) was also reported in wild C. helgolandicus females from the Northern Adriatic Sea, when collected at sampling sites showing high diatom and oxylipin concentrations [35]. These results were confirmed in the laboratory after feeding C. helgolandicus females from the Northern Adriatic Sea and the English Channel with the oxylipin-producing diatom Skeletonema marinoi [26].

Conclusions
In general, information collected at the transcriptomic level revealed that T. stylifera females generating nauplii with low viability showed down-or up-regulation of sequences encoding for proteins involved in metabolic processes crucial for final reproductive success. In particular, key genes involved in ATP synthesis, cell differentiation and signal transduction could explain at a molecular level the maternally-mediated low naupliar survival observed on May 23rd. Differentially expressed sequences can also indicate molecular pathways regulated by T. stylifera females to counteract negative reproductive potential. In this perspective, some differentially expressed genes are known to encode for interconnected proteins participating in cellular biochemical patterns finally regulating apoptosis, cell cycle, yolk protein precursors (YPP), gene expression and autophagy. These processes could partially disentangle coordinated molecular pathways leading to altered naupliar viability. Interestingly, key genes involved in detoxification such as aldehyde dehydrogenases, superoxide dismutase or glutathione synthase were not differentially expressed in T. stylifera from the GoN.
Although we cannot exclude that biotic and abiotic factors not considered in this study may have mediated molecular responses of copepods, our data suggest that temperature, oxygen, pH and salinity played a marginal role, because these factors remained constant in the two sampling dates. Also, T. stylifera responses did not seem to be mediated by food deficiency, because differences in faecal pellet production were not significant despite the higher phytoplankton concentration on the 30th of May than the 23rd. Interestingly, both oxylipin-per lire concentrations and oxylipin-per-diatom-cell production were higher on the 23rd of May. Involvement of these chemicals in differential molecular responses and impaired naupliar survival in copepods from natural populations is plausible. Down-regulation of A5 Obp could suggest impaired detection of harmful preys in females collected on May 23rd, higher ingestion of oxylipins and lower naupliar survival. Indeed, our interpretation is tempting, but needs to be supported by larger datasets exploring wider temporal variations in the reproductive success and the molecular responses of natural T. stylifera populations. Our transcriptomic analysis could offer a guideline to identify biomarker genes of interest explaining at a molecular level copepod reproductive responses to biotic and abiotic variables, both in laboratory experiments and in the natural environment.

Field sampling and physiological responses of Temora stylifera
Environmental variables (temperature, salinity, pH and oxgen), phytoplankton community composition and phytoplankton-derived oxylipins were measured from surface water sampled at the Long-Term Ecological Research station-MareChiara (LTER-MC) [119] on the 23rd and the 30th of May 2017 as described by Russo et al. [47]. Briefly, abiotic variables were measured by deploying a SBE911 CTD. For phytoplankton and oxylipin analysis, surface water was collected with a bucket deployed from the boat. Phytoplankton cells were left to sediment in a Utermöhl chamber and counted under a Zeiss Axiovert200 optic microscope (Zeiss, Oberkochen, Germany). For oxylipin analysis, phytoplankton cells were accumulated on polycarbonate filters (2 μm mesh size). Cells were re-suspended with water and sonicated. Oxylipins were finally extracted with dichloromethane and methylated before the LC-MS/MS analysis, performed with a Q Exactive Hybrid Quadrupole Orbitrap (Thermo Scientific, Waltham, USA) using a methanol: water gradient [14].
Zooplankton samples were collected the same days at LTER-MC by the SZN Unit through oblique towing of a 200 μm Nansen net (113 cm mouth diameter) equipped with a 200 μm filtering cod-end, on-board of the R/V "Vettoria". A maximum of 15 females of T. stylifera were sorted from the zooplankton community. They were individually incubated (20°C, 12 h,12 h dark,light cycle) in crystallizing dishes filled with 100 ml of 50 μm prefiltered surface seawater collected the same day at LTER-MC. After 24 h, the females were removed from the crystallizing dishes and the number of spawned and chewed eggs as well as the number of faecal pellets were counted under an inverted microscope (Zeiss Axio-vert25) at 25x magnification [120].
The crystallizing dishes were incubated (20°C, 12 h,12 h dark:light cycle) for additional 48 h to allow eggs to hatch. Subsequently, the number of dead nauplii (stage NI) laying on the bottom of the crystallizing dishes was counted under the inverted microscope, in order to assess survival of first non-feeding nauplii (NI), and thus, maternally-related larval survivorship. The content of the crystallizing dish was then fixed by adding 15 ml of ethanol (96%), and the number of hatched membranes, the non-viable eggs and total nauplii present in the container were counted. Hatching success of the eggs as well as naupliar survival rates were calculated following [40]: Differences in faecal pellet and egg productions (N per female per day), egg hatching (%) and NI naupliar survival (%) between samples collected on May 23rd and 30th were analysed on the basis of t-tests considering Welch's correction for unequal variances (N = 15). To avoid type I error, statistical significance was considered at 99%.

RNA extraction and de novo transcriptome assembly
A number of 30-60 adult females of T. stylifera were sorted from the zooplankton sample collected on the 23rd and the 30th of May 2017 at LTER-MC and transferred to three replicate 1.5 ml Eppendorf tubes per date. In particular, two replicates from May 23rd consisted of 8 individuals, while one replicate contained 10 individuals. The three replicates from May 30th all consisted of 10 specimens. Samples were immediately frozen in liquid nitrogen after removing excess water and finally stored at − 80°C for later total RNA extraction. These dates were selected because of the opposite naupliar survival rates measured in these two consecutive samples (25 and 93%, respectively; t-test considering unequal variances; N = 15; see results).
Total RNA from pools of T. stylifera females (8-10 individuals per pool), was extracted using the RNeasy Micro Kit (Qiagen, Hilden, Germany) [112], following manufacturer's specification and performing on-column DNase-I treatment. The RNA was finally eluted in 15 μl of RNase-free water. RNA concentration (ng/μl) and purity were assessed through Nanodrop ND-1000 UV-Vis spectrophotometer (Marshall Scientific, Hampton, USA). Overall RNA integrity and DNA contamination were checked by electrophoresis of at least 200 ng of RNA on a 1% agarose gel in 0.5x Tris Borate EDTA buffer (TBE) and by analysing 150-200 ng of RNA in a 6000 Nano LabChip of an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, USA), which defines RNA quality as RNA Integrity Number (RIN) [121]. RIN values > 8 were considered suitable for NGS analysis.
At least 2 μg of extracted RNA per sample (200 ng/μL) were delivered to Genomix4Life S.r.l. (Laboratory of Molecular Medicine and Genomics of the University of Salerno, Salerno, Italy) for library preparation, sequencing and de novo transcriptome assembly. Six cDNA libraries, each one generated from a single RNA sample extracted from groups of 8-10 individuals, were prepared using TruSeq RNA Sample Prep Kit (Illumina) according to manufacturer's recommendations and pooled such that each index-tagged sample was present in equal-molar amounts. The pooled samples were subjected to cluster generation and multiplexed sequencing using an Illumina HiSeq 2500 platform (Illumina, San Diego, USA) in a 2 × 100 paired-end format. Raw reads were cleaned, trimmed and clipped with BBDuk (https://jgi. doe.gov/data-and-tools/bbtools/) setting a minimum phred score (Q) of 20 (base call accuracy of 99%), and a minimum length of 35 nucleotides. The quality of the reads before and after trimming was checked with the software FASTQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). High quality paired-end reads from all samples were used as input for transcriptome assembly using Trinity [122]. A filter for contaminants was performed by BLASTing the transcripts against the NCBI nr database, discarding all the sequences having a significant hit (E-value ≤0.0001) against bacteria or vegetal cells. Two different transcriptomes were generated, a 'full' transcriptome comprising all Trinity assembled transcripts (isoforms) and a 'reference' transcriptome of 'Trinity predicted genes', consisting of assembled transcripts with unique TR#_c#_g# identifiers (unigenes). This latter transcriptome 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#) [123].

Differential expression analysis and functional annotation of Temora stylifera transcriptome
Quantification of de novo assembled transcripts (isoforms) and unigenes abundance for each sample was assessed using RSEM software provided by the Trinity package [124], after mapping back the reads on the 'full' and the 'reference' transcriptome, respectively, using STARS [125]. Mean expression levels of isoforms and unigenes read counts, expressed as Counts Per Million (CPM), were used as input to perform Differential Expression (DE) analysis between samples collected on May 23rd with respect to samples collected on May 30th, using the Trinity DeSeq2 package [126]. Statistical significance was obtained by performing a hypergeometric test and corrected p-value using the False Discovery Rate (FDR) method [127] and isoforms or unigenes having a FDR ≤0.05 were considered differentially expressed.
Functional annotation of T. stylifera 'full' (isoforms) and 'reference' (unigenes) transcriptomes was performed using the comprehensive bioinformatics tool Blast2Go [128,129]. Sequence similarity was found using BLASTx function (E-value cut-off set to 1 − 3 ) [130], which compares a nucleotide query sequence translated in all reading frames against a non-redundant protein sequence database (nr) used to find potential translation products of an unknown nucleotide sequence. Default parameters were selected, with a number of HITs equal to 1. Subsequently, Gene Ontology annotation was associated to mapped sequences and divided by the three GO terms: Biological Process (BP), Molecular Function (MF) and Cellular Component (CC).
Differences between samples from May 23rd and May 30th were highlighted through heat-maps based on logtransformed raw counts of differentially expressed unigenes and isoforms. Moreover, separate heat-maps were represented to highlight differences in A5 Odorant Binding Protein isoforms and in isoforms annotated as sequences involved in reproduction and cell development processes. Euclidean distance was calculated among variables, which were clustered on the basis of the 'complete-linkage' method. Cluster analysis and heatmap representation was performed using R 3.6.3 implemented in R-Studio.

Transcriptome validation through RT-qPCR
To validate Illumina sequencing and differential expression results, six reference genes (RGs), previously optimized in Calanus helgolandicus [27,28], were tested in T. stylifera through RT-qPCR to identify the most stable genes in samples used for transcriptome assembly. RGs were: Actin (Act), Histone 3 (Ist), two ribosomal units (18S RNA and the ribosomal protein S20) and Ubiquitin (Ubi). Firstly, primer specificity and amplification efficiency of the six reference genes were assessed. Complementary DNA (cDNA) needed as template for RT-qPCR analysis was retro-transcribed from 1 μg of T. stylifera total RNA, extracted as described before, in a final volume of 20 μl using iScriptTM cDNA Synthesis Kit (Bio-Rad, Hercules, USA) following the manufacturer's instructions. PCRs were performed on a GeneAmp PCR System 9700 (Applied Biosystems, Foster City, USA), with 2 μl of 10× PCR reaction buffer (Roche), 2 μl of 10 × 2 mM dNTPs (Roche), 0.8 μl of 5 U/μl Taq polymerase (Roche), 1 μl of 20 pmol/μl of each primer (forward and reverse), 1 μl of T. stylifera template cDNA and nuclease-free water to 20 μl of final volume. The PCR program consisted of a denaturation step at 95°C for 3 min, 40 cycles at 95°C for 30 s, 60°C for 1 min and 72°C for 30 s, and a final extension step at 72°C for 7 min. Amplified PCR products were analysed by electrophoresis on 1.5% agarose gel in 0.5x TBE buffer to check proper length specificity of the cDNA products.
RT-qPCR reactions for primer amplification efficiency were then performed in a MicroAmp Optical 384-Well reaction plates (Applied Biosystem, Foster City, USA) with optical adhesive covers (Applied Biosystem, Foster City, USA), using a Viia7 Real Time PCR system (Applied Biosystem, Foster City, USA). The final PCR volume for each sample was 10 μl, with 5 μl of SensiFAST SYBR Green Master Mix (Meridian Inc., Cincinnati, USA), 1 μl of cDNA template and 4 μl (concentration of 0.7 pmol/μl) of each primer pair. All RT-qPCR reactions were carried out in triplicate to capture intra-assay variability. Three negative controls (consisting of 1 μl of water instead of the cDNA template) were considered for each primer pair. PCR conditions for all samples analysed were set as follows: 95°C for 20 s, 40 cycles of 95°C for 1 s, and 60°C for 20 s. Dissociation protocol with a gradient (0.5°C every 30 s) from 65°C to 95°C was also used to investigate the specificity of the primers and presence of primer dimers. Gene-specific amplification was confirmed by a single peak in the melting curve analysis. To quantify gene expression, primer amplification efficiencies were calculated through six serial dilutions of cDNA (1, 1:5, 1:10, 1:50, 1:100 and 1:500) for all primer pairs. The reference equation for efficiency calculation is E = 10 -1/slope , where the slope is obtained from a standard curve between Ct values and the log 10 of each dilution factor. Stability of the six reference genes in T. stylifera samples used for transcriptome analysis (i.e. females collected on the 23rd and the 30th of May) was finally identified through RT-qPCR analyses as described above, using 1 μl of cDNA template (dilution 1:5) obtained by retrotranscription of the same total RNA employed for the Illumina sequencing. The most stable reference genes were then evaluated using RefFinder (http://leonxie.esy.es/ RefFinder/) [60], a user-friendly web-based tool for evaluating and screening reference genes from extensive experimental datasets. Based on the rankings from each program, RefFinder assigns an appropriate weight to an individual gene and calculates the geometric mean of their weights for the overall final ranking. RGs 18S and Ubi were used as reference genes for the analyses, because they were indicated as the most stable ones.
Once the most stable reference genes were identified, 9 unigenes and 3 isoforms were selected from the list of Differentially Expressed Genes (DEGs) generated by Illumina sequencing and subsequent DE analysis and used to validate differential expression results through RT-qPCR.
In particular, the unigenes Putative Odorant Binding Protein A5 (A5), cAMP-Responsive Element-Binding Protein-Like2 (Crebl), Putative ATP-Dependent RNA Helicase (Me31b), Serine/Threonine Protein Phosphatase These selected sequences were then used to design specific forward and reverse primers through the online available platform Primer3web 4.1.0 (http://bioinfo.ut.ee/ primer3/) and synthesized by Sigma-Aldrich (Merk, Germany). Desired primer size was set in the range of 18-24 bp (optimum 20) and melting temperature (Tm) between 59 and 61°C (optimum 60°C). Amplification products ranged between 110 and 240 bp. To verify primer specificity, PCRs were performed as described for RGs. Length of PCR product was verified through agarose gel electrophoresis. Subsequently, specificity and amplification efficiencies of the primers were further tested in RT-qPCR, using serial dilutions (1, 1:5, 1:10, 1: 50, 1:100) of T. stylifera cDNA template. Settings of qPCR reactions were the same as described before. As for RGs, gene-specific amplification was confirmed by a single peak in the melting curve analysis.
To validate differential expression analysis results, cDNA obtained from T. stylifera samples used for transcriptome analysis was used as template (1:10 dilution) for RT-qPCR reactions. Relative expression levels of each target gene in the samples collected on the 23rd of May was compared to samples collected on the 30th of May through the Pfaffl equation [131], using the tool REST (Relative Expression Software Tool) [132]. Results were analysed based on the following equation, where the relative expression of a target gene is compared in a 'test' sample (copepod collected on the 23rd of May) versus a 'control' (copepod collected on the 30th of May), normalized to each reference gene [131]: The E (target) is the RT-qPCR efficiency of target gene primers; E (ref) is the RT-qPCR efficiency of the reference gene primers; ΔCt (target) is the Ct deviation of control-sample of the target gene transcript; ΔCt (reference) is the Ct deviation of control-sample of reference gene transcript.
Additional file 1 Table S1: De novo transcriptome assembly of Temora stylifera. Number of reads counted in the forward and the reverse filament are shown along with the number of assembled transcripts and Trinity predicted genes (unigenes), transcripts with unique TR#_c#_g# identifiers. The latter list includes singletons as well as the longest isoform of each predicted gene. Average transcript length, median and N 50 are also indicated.
Percentage distribution of E-value (0 < E < 1 − 3 ) and sequence similarity percentage (30-100%) are displayed on the top of the figure. Bottom panel describes top 20 blast hit taxon groups; a subplot of the total hits is shown for clarity.
Additional file 3 Table S2: Temora stylifera differentially expressed isoforms that received functional annotation in Blast2Go. Trinity ID number with predicted gene and isoform identifiers, length (bp), log 2 -Fold-Change (log 2 -FC), adjusted p-value (p-adj) of statistical analysis (FDR) for each predicted genes, sequence description and functional annotation as provided by Blast2Go are shown. Sequences are ordered by p-adj values within each down-regulated (negative log2-FC) and upregulated (positive log2-FC) isoforms.
Additional file 4 Table S3: List of Reference Genes (RGs) tested in Temora stylifera. Protein name, function, amplicon length (A L ) in base pairs (bp), primer sequence and amplification efficiencies (%E) are shown.