Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

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

BMC Genomics volume 21, Article number: 693 (2020) Cite this article

Abstract

Background

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.

Results

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.

Conclusions

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.

Background

Among zooplankton, marine and freshwater copepods represent one of the most morphologically and functionally diverse groups [1], playing a central role in food web dynamics and biogeochemical cycles [2]. In this perspective, assessment of biotic and abiotic factors influencing copepod populations can be of primary importance to understand marine pelagic food web functioning. Phytoplankton-derived oxylipins potentially represent key factors affecting wild copepod populations [3]. These molecules are end products of well characterized enzymatic pathways activated after cell wounding, starting from lipolytic release of free fatty acids (FFAs) from complex lipids [4,5,6] and proceeding through oxygenation of FFAs by lipoxygenases (LOX) [5, 7,8,9,10,11,12,13,14,15,16].

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 oxylipin-producing diatoms has been also generated [33].

Variations in copepod egg production, hatching success and naupliar survival in response to phytoplankton abundance and composition have been investigated in several copepod species through field surveys [3, 22, 34,35,36,37,38,39,40,41,42], but information about the molecular responses of adult females from natural populations are still limited to the Northern Ariatic Sea [35].

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-per-diatom-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].

Results

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 106 cells/L) than the 30th of May (17.64 106 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.

Table 1 Abiotic variables, phytoplankton composition and oxylipins
Full size table

Similarly, oxylipin concentrations were also higher on the 23rd (202.82 ng/L) than the 30th (173.47 ng/L) of May. Also, oxylipin-per-cell production was higher on the 23rd of May (41.49 fg/diatom-cell) than the 30th (27.16 fg/diatom-cell).

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, non-significant 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).

Fig. 1
figure 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

Full size image

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, N50 = 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 non-redundant 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 unigenes were almost identical (similarity between 95 and 100%), while 76.1% of the sequences had a similarity ranging from 100 to 60% (Additional file 2: 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%).

Fig. 2
figure 2

Blast2Go Gene Ontology (GO) annotation of Temora stylifera reference transcriptome (unigenes). The number of sequences assigned to the three GO classes Biological Process (BP), Molecular Function (MF) and Cellular Component (CC) are shown

Full size image

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 up-regulated. Of the total 144 differentially expressed sequences, 31 (6 up-regulated and 25 down-regulated) received GO assignment and functional annotation (Table 2).

Table 2 Temora stylifera differentially expressed unigenes
Full size table

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.

Fig. 3
figure 3

Blast2Go Gene Ontology (GO) annotation of the differentially expressed transcript isoforms in Temora stylifera. First 27 terms of each category:Biological Process (BP), Cellular Component (CC) and Molecular Function (MF) are shown along the x-axis; the number of sequences assigned to each GO term within each GO category is displayed on the y-axis. Down-regulated sequences are indicated by blue bars and up-regulated sequences by red bars

Full size image

Interestingly, a number of specific GO terms contained isoforms that were exclusively up- or down-regulated on May 23rd in comparison to May 30th (Fig. 3). Among BP sub-categories, sequences involved in ATP metabolism (6 sequences), cell communication (4 sequences), cellular response to stimulus (4 sequences), signal transduction (4 sequences), cellular development (1 sequence), cellular localization (1 sequence) and organism interaction (1 sequence) were specifically down-regulated. In contrast, sequences involved in catabolic processes (2 sequences), cellular component organization (2 sequences), anatomical structure morphogenesis (1 sequence), microtubule-based process (1 sequence), negative regulation of metabolic process (1 sequence), pattern specification process (1 sequence) and sexual reproduction (1 sequence) were exclusively up-regulated.

Consistent differences between replicates collected on May 30th and May 23rd were supported by clustering among objects (Q-mode analysis) described by raw 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).

Fig. 4
figure 4

Heat-maps of differentially expressed sequences expressed as log-transformed raw counts. Euclidean distance was considered as distance measure and variables were clustered on the basis of the complete-linkage method. Numbers near the dendrogram of variables (R-mode) describe the main clusters identified for sequences. Cluster position of the selected differentially expressed unigenes and isoforms used for transcriptome validation is indicated. Dendrograms on the top of the heat-maps describe grouping of samples on the basis of complete-linkage clustering method (Q-mode). Data refer to the three replicates of wild Temora stylifera collected on the 30th of May 2017 and the 23rd of May. Heat-maps at the top show raw counts of differentially expressed unigenes (DE unigenes) and isoforms (DE isoforms). Heat-maps at the bottom display log-transformed raw counts of A5 Odorant binding proteins (DE A5-Obp iso) isoforms and isoforms annotated as sequences involved in reproduction and development (DE Repr & Dev iso)

Full size image

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 down-regulated 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.

Table 3 List of unigenes and isoforms tested for Temora stylifera transcriptome validation
Full size table

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, log2-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 log2-fold-changes of the three selected isoforms (Mob1b, Vasa and Pafah; log2-FC = − 9.48, 7.14 and 6.93, respectively) were much higher than RT-qPCR values. On the contrary, the expression ratio obtained after RT-qPCR analysis for the unigenes A5 Obp (log2-FC = − 4.6), Arih1 (log2-FC = − 0.85), Ppa2 (log2-FC = − 1.94), Arsb (log2-FC = − 1.78) and Crebl (log2-FC = − 1.28) resembled results obtained from the RNA-Seq analysis (log2-FC = − 4.55, − 1.76, − 3.72, − 1.29 and − 2.03, respectively).

Fig. 5
figure 5

Comparison between RNA-Seq and RT-qPCR in Temora stylifera. Relative expression ratio (log2-Fold-Change, FC) of 9 unigenes (A5, Obst, Arih1, Ppa2, Arsb, Hsp-70, Vldlr, Me31b and Crebl2) and 3 isoforms (Mob1b, Vasa and Pafah) in Temora stylifera from May 23rd vs 30th samples, measured through RNA-Seq (blue bars) or RT-qPCR (yellow bars). For RT-qPCR results, bars represented mean ± SD values and data are normalized to RNA 18S and Ubi reference genes

Full size image

Discussion

Estimating early-life history variables such as egg production, hatching success and survival rates of non-feeding (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,18,19,20,21,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 sub-categories, 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-per-diatom-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. Down-regulation of A5 Obp in T. stylifera females was confirmed by RT-qPCR results (log2-FC in RNA-Seq = − 4.6; log2-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 down-regulated. The 1-Phosphatidylinositol 3-Phosphate 5-Kinase (Fab1b) showed the strongest down-regulation (log2-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 down-regulation (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 de-phosphorylation 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 up-regulated transcripts (log2-FC of 4.45). Folate one-carbon 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, dys-regulation 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 up-regulated. 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 Ubiquitin-protein 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 (log2-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 (log2-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.

Methods

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 pre-filtered 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 Axiovert25) 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]:

$$ Hatching\ success=\frac{N\ (membranes)}{N\ \left( laid\ eggs\right)}\ 100 $$
$$ Naupliar\ survival=100-\left(\ \frac{N\ \left( dead\ nauplii\right)}{N\ (membranes)}\ 100\right) $$

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 log-transformed 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 heat-map 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 log10 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 retro-transcription 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 (Ppa2), Protein Obstructor E-Like (Obst), E3 Ubiquitin-Protein Ligase Arih1-Like (Arih1), Arylsulfatase B-Like (Arsb), Heat-Shock Protein 70 (Hsp70), Very Low-Density Lipoprotein Receptor-Like (Vldlr) and the isoforms MOB Kinase Activator 1B (Mob1b), ATP-Dependent RNA Helicase Vasa-Like (Vasa) and Platelet-Activating Factor Homolog 2 (Pafah) were selected on the basis of fold-change between samples collected on the 23rd of May 2017 with respect to samples collected on the 30th of May 2017, E-value and similarity percentage (indicating how much of the query sequence corresponds to the reference one in public databases) reported after DE analysis and functional annotation of the transcriptome.

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]:

$$ Ratio=\frac{E{(Target)}^{\varDelta Ct\ target\left( Control- Sample\right)}}{E{(Reference)}^{\varDelta Ct\ reference\left( Control- Sample\right)}} $$

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.

Availability of data and materials

Raw reads are stored into the NCBI Sequence Read Archive database under accession number PRJNA632714.

Abbreviations

HDoHE:

Hydroxy docosahexaenoic acid

RIN:

RNA integrity number

FDR:

False discovery rate

GO:

Gene ontology

Obp:

Odorant binding protein

Obst-E:

Protein obstructor E-like

Hsp70:

Heat shock protein 70

Mob1b:

Mob kinase activator 1B

Vasa:

ATP-dependent RNA helicase vasa-like

Pafah:

Platelet-activating factor homolog 2

Arih1:

Ubiquitin-protein ligase Arih1-like

Arsb:

Arylsulfatase B-like

Crebl:

cAMP responsive element binding protein-like2

Ppa2:

Serine/Threonine protein phosphatase

Mthfd:

Methylene tetrahydrofolate dehydrogenase

Tnrc18:

Trinucleotide repeat-containing gene 18

Dpt:

Dermatopontin

Mael:

Protein maelstrom

PIP:

Phosphatidylinositol phosphate

Me31b:

Putative ATP-dependent RNA helicase

Klhl12:

Kelk 12

GAG:

Glycosaminoglycans

YPP:

Yolk protein precursors

PPIn:

Phosphatidylinositide

LTER-MC:

Long-term ecological research station MareChiara

DE:

Differential expression

CMP:

Counts per million

BP:

Biological process

MF:

Molecular function

CC:

Cellular component

Ist:

Histone-3

Ubi:

Ubiquitin

Act:

Actin

References

  1. Steinberg DK, Landry MR. Zooplankton and the ocean carbon cycle. Annu Rev Mar Sci. 2017;9:413–44.

    Article  Google Scholar 

  2. Turner JT. Zooplankton fecal pellets, marine snow, phytodetritus and the ocean’s biological pump. Prog Oceanogr. 2015;130:205–48. https://0-doi-org.brum.beds.ac.uk/10.1016/j.pocean.2014.08.005.

    Article  Google Scholar 

  3. Russo E, Ianora A, Carotenuto Y. Re-shaping marine plankton communities: effects of diatom oxylipins on copepods and beyond. Mar Biol. 2018;166:9.

    Article  CAS  Google Scholar 

  4. Adelfi MG, Vitale RM, d’Ippolito G, Nuzzo G, Gallo C, Amodeo P, et al. Patatin-like lipolytic acyl hydrolases and galactolipid metabolism in marine diatoms of the genus Pseudo-nitzschia. Biochim Biophys Acta Mol Cell Biol Lipids. 2019;1864:181–90. https://0-doi-org.brum.beds.ac.uk/10.1016/j.bbalip.2018.11.008.

    Article  CAS  PubMed  Google Scholar 

  5. d’Ippolito G, Tucci S, Cutignano A, Romano G, Cimino G, Miralto A, et al. The role of complex lipids in the synthesis of bioactive aldehydes of the marine diatom Skeletonema costatum. Biochim Biophys Acta. 1686;2004:100–7.

    Google Scholar 

  6. Cutignano A, d’Ippolito G, Romano G, Lamari N, Cimino G, Febbraio F, et al. Chloroplastic glycolipids fuel aldehyde biosynthesis in the marine diatom Thalassiosira rotula. ChemBioChem. 2006;7:450–6. https://0-doi-org.brum.beds.ac.uk/10.1002/cbic.200500343.

    Article  CAS  PubMed  Google Scholar 

  7. d’Ippolito G, Lamari N, Montresor M, Romano G, Cutignano A, Gerecht A, et al. 15S-Lipoxygenase metabolism in the marine diatom Pseudo-nitzschia delicatissima. New Phytol. 2009;183:1064–71.

    Article  PubMed  CAS  Google Scholar 

  8. Lamari N, Ruggiero MV, d’Ippolito G, Kooistra WHCF, Fontana A, Montresor M. Specificity of lipoxygenase pathways supports species delineation in the marine diatom genus Pseudo-nitzschia. PLoS One. 2013;8:1–10.

    Article  CAS  Google Scholar 

  9. d’Ippolito G, Romano G, Iadicicco O, Miralto A, Ianora A, Cimino G, et al. New birth-control aldehydes from the marine diatom Skeletonema costatum: characterization and biogenesis. Tetrahedron Lett. 2002;43:6133–6.

    Article  Google Scholar 

  10. d’Ippolito G, Romano G, Caruso T, Spinella A, Cimino G, Fontana A. Production of octadienal in the marine diatom Skeletonema costatum. Org Lett. 2003;5:885–7.

    Article  PubMed  CAS  Google Scholar 

  11. d’Ippolito G, Cutignano A, Briante R, Febbraio F, Cimino G, Fontana A. New C16 fatty-acid-based oxylipin pathway in the marine diatom Thalassiosira rotula. Org Biomol Chem. 2005;3:4065–70.

    Article  PubMed  CAS  Google Scholar 

  12. d’Ippolito G, Cutignano A, Tucci S, Romano G, Cimino G, Fontana A. Biosynthetic intermediates and stereochemical aspects of aldehyde biosynthesis in the marine diatom Thalassiosira rotula. Phytochemistry. 2006;67:314–22.

    Article  PubMed  CAS  Google Scholar 

  13. Fontana A, D’Ippolito G, Cutignano A, Miralto A, Ianora A, Romano G, et al. Chemistry of oxylipin pathways in marine diatoms. Pure Appl Chem. 2007;79:481–90.

    Article  CAS  Google Scholar 

  14. d’Ippolito G, Nuzzo G, Sardo A, Manzo E, Gallo C, Fontana A. Lipoxygenases and lipoxygenase products in marine diatoms. In: Methods in Enzymology. 2018. p. 69–100.

  15. Nanjappa D, d’Ippolito G, Gallo C, Zingone A, Fontana A. Oxylipin diversity in the diatom family Leptocylindraceae reveals DHA derivatives in marine diatoms. Mar Drugs. 2014;12.

  16. Gerecht A, Carotenuto Y, Ianora A, Romano G, Fontana A, d’Ippolito G, et al. Oxylipin production during a mesocosm bloom of Skeletonema marinoi. J Exp Mar Bio Ecol. 2013;446:159–65.

    Article  CAS  Google Scholar 

  17. Fontana A, d’Ippolito G, Cutignano A, Romano G, Lamari N, Gallucci AM, et al. LOX-induced lipid peroxidation mechanism responsible for the detrimental effect of marine diatoms on zooplankton grazers. ChemBioChem. 2007;8:1810–8.

    Article  CAS  PubMed  Google Scholar 

  18. Ianora A, Miralto A, Poulet SA, Carotenuto Y, Buttino I, Romano G, et al. Aldehyde suppression of copepod recruitment in blooms of a ubiquitous planktonic diatom. Nature. 2004;429:403–7.

    Article  CAS  PubMed  Google Scholar 

  19. Barreiro A, Carotenuto Y, Lamari N, Esposito F, D’Ippolito G, Fontana A, et al. Diatom induction of reproductive failure in copepods: the effect of PUAs versus non volatile oxylipins. J Exp Mar Bio Ecol. 2011;401:13–9.

    Article  CAS  Google Scholar 

  20. Barreiro A, Guisande C, Maneiro I, Vergara AR, Riveiro I, Iglesias P. Zooplankton interactions with toxic phytoplankton: some implications for food web studies and algal defence strategies of feeding selectivity behaviour, toxin dilution and phytoplankton population diversity. Acta Oecol. 2007;32:279–90.

    Article  Google Scholar 

  21. Buttino I, De Rosa G, Carotenuto Y, Mazzella M, Ianora A, Esposito F, et al. Aldehyde-encapsulating liposomes impair marine grazer survivorship. J Exp Biol. 2008;211:1426–33.

    Article  CAS  PubMed  Google Scholar 

  22. Carotenuto Y, Ianora A, Buttino I, Romano G, Miralto A. Is postembryonic development in the copepod Temora stylifera negatively affected by diatom diets? J Exp Mar Bio Ecol. 2002;276:49–66.

    Article  Google Scholar 

  23. Ceballos S, Ianora A. Different diatoms induce contrasting effects on the reproductive success of the copepod Temora stylifera. J Exp Mar Bio Ecol. 2003;294:189–202.

    Article  Google Scholar 

  24. Carotenuto Y, Ianora A, Miralto A. Maternal and neonate diatom diets impair development and sex differentiation in the copepod Temora stylifera. J Exp Mar Bio Ecol. 2011;396:99–107.

    Article  Google Scholar 

  25. Ianora A, Romano G, Carotenuto Y, Esposito F, Roncalli V, Buttino I, et al. Impact of the diatom oxylipin 15S-HEPE on the reproductive success of the copepod Temora stylifera. Hydrobiologia. 2011;666:265–75.

    Article  CAS  Google Scholar 

  26. Lauritano C, Carotenuto Y, Miralto A, Procaccini G, Ianora A. Copepod population-specific response to a toxic diatom diet. PLoS One. 2012;7.

  27. Lauritano C, Borra M, Carotenuto Ylenia Y, Biffali E, Miralto A, Procaccini G, et al. First molecular evidence of diatom effects in the copepod Calanus helgolandicus. J Exp Mar Bio Ecol. 2011;404:79–86.

    Article  CAS  Google Scholar 

  28. Lauritano C, Borra M, Carotenuto Y, Biffali E, Miralto A, Procaccini G, et al. Molecular evidence of the toxic effects of diatom diets on gene expression patterns in copepods. PLoS One. 2011;6.

  29. Lauritano C, Procaccini G, Ianora A. Gene expression patterns and stress response in marine copepods. Mar Environ Res. 2012;76:22–31. https://0-doi-org.brum.beds.ac.uk/10.1016/j.marenvres.2011.09.015.

    Article  CAS  PubMed  Google Scholar 

  30. Lauritano C, Carotenuto Y, Procaccini G, Turner JT, Ianora A. Changes in expression of stress genes in copepods feeding upon a non-brevetoxin-producing strain of the dinoflagellate Karenia brevis. Harmful Algae. 2013;28:23–30. https://0-doi-org.brum.beds.ac.uk/10.1016/j.hal.2013.05.004.

    Article  CAS  Google Scholar 

  31. Lauritano C, Carotenuto Y, Vitiello V, Buttino I, Romano G, Hwang JS, et al. Effects of the oxylipin-producing diatom Skeletonema marinoi on gene expression levels of the calanoid copepod Calanus sinicus. Mar Genomics. 2015;24:89–94. https://0-doi-org.brum.beds.ac.uk/10.1016/j.margen.2015.01.007.

    Article  PubMed  Google Scholar 

  32. Carotenuto Y, Dattolo E, Lauritano C, Pisano F, Sanges R, Miralto A, et al. Insights into the transcriptome of the marine copepod Calanus helgolandicus feeding on the oxylipin-producing diatom Skeletonema marinoi. Harmful Algae. 2014;31:153–62.

    Article  CAS  PubMed  Google Scholar 

  33. Asai S, Sanges R, Lauritano C, Lindeque PK, Esposito F, Ianora A, et al. De novo transcriptome assembly and gene expression profiling of the copepod Calanus helgolandicus feeding on the PUA-producing diatom Skeletonema marinoi. Mar Drugs. 2020;18. https://0-doi-org.brum.beds.ac.uk/10.3390/md18080392.

  34. Halsband-Lenk C, Pierson JJ, Leising AW. Reproduction of Pseudocalanus newmani (Copepoda: Calanoida) is deleteriously affected by diatom blooms - a field study. Prog Oceanogr. 2005;67:332–48.

    Article  Google Scholar 

  35. Lauritano C, Romano G, Roncalli V, Amoresano A, Fontanarosa C, Bastianini M, et al. New oxylipins produced at the end of a diatom bloom and their effects on copepod reproductive success and gene expression levels. Harmful Algae. 2016;55:221–9.

    Article  CAS  PubMed  Google Scholar 

  36. Ianora A, Bastianini M, Carotenuto Y, Casotti R, Roncalli V, Miralto A, et al. Non-volatile oxylipins can render some diatom blooms more toxic for copepod reproduction. Harmful Algae. 2015;44:1–7.

    Article  CAS  Google Scholar 

  37. Wichard T, Poulet SA, Boulesteix AL, Ledoux JB, Lebreton B, Marchetti J, et al. Influence of diatoms on copepod reproduction. II. Uncorrelated effects of diatom-derived α,β,γ,δ-unsaturated aldehydes and polyunsaturated fatty acids on Calanus helgolandicus in the field. Prog Oceanogr. 2008;77:30–44.

    Article  Google Scholar 

  38. Poulet SA, Escribano R, Hidalgo P, Cueff A, Wichard T, Aguilera V, et al. Collapse of Calanus chilensis reproduction in a marine environment with high diatom concentration. J Exp Mar Bio Ecol. 2007;352:187–99.

    Article  Google Scholar 

  39. Poulet SA, Wichard T, Ledoux JB, Lebreton B, Marchetti J, Dancie C, et al. Influence of diatoms on copepod reproduction. I Field and laboratory observations related to Calanus helgolandicus egg production Mar Ecol Prog Ser 2006;308 2016:129–142.

  40. Carotenuto Y, Ianora A, Di Pinto M, Sarno D, Miralto A. Annual cycle of early developmental stage survival and recruitment in the copepods Temora stylifera and Centropages typicus. Mar Ecol Prog Ser. 2006;314:227–38.

    Article  Google Scholar 

  41. Halsband-Lenk C. Metridia pacifica in Dabob Bay, Washington: the diatom effect and the discrepancy between high abundance and low egg production rates. Prog Oceanogr. 2005;67:422–41.

    Article  Google Scholar 

  42. Pierson JJ, Halsband-Lenk C, Leising AW. Reproductive success of Calanus pacificus during diatom blooms in Dabob Bay, Washington. Prog Oceanogr. 2005;67:314–31.

    Article  Google Scholar 

  43. Ribera d’Alcalà M, Conversano F, Corato F, Licandro P, Mangoni O, Marino D, et al. Seasonal patterns in plankton communities in a pluriannual time series at a coastal Mediterranean site (gulf of Naples): an attempt to discern recurrences and trends. Sci. 2004;68:65–83.

    Google Scholar 

  44. Mazzocchi MG, Buffoni G, Carotenuto Y, Pasquali S, Ribera d’Alcalà M. Effects of food conditions on the development of the population of Temora stylifera: a modeling approach. J Mar Syst. 2006;62:71–84.

    Article  Google Scholar 

  45. Mazzocchi MG, Dubroca L, García-Comas C, Di Capua I, Ribera d’Alcalà M. Stability and resilience in coastal copepod assemblages: the case of the Mediterranean long-term ecological research at station MC (LTER-MC). Prog Oceanogr. 2012;97–100:135–51. https://0-doi-org.brum.beds.ac.uk/10.1016/j.pocean.2011.11.003.

    Article  Google Scholar 

  46. Mazzocchi MG, Licandro P, Dubroca L, Di Capua I, Saggiomo V. Zooplankton associations in a Mediterranean long-term time-series. J Plankton Res. 2011;33:1163–81.

    Article  Google Scholar 

  47. Russo E, d’Ippolito G, Fontana A, Sarno D, D’Alelio D, Busseni G, et al. Density-dependent oxylipin production in natural diatom communities: possible implications for plankton dynamics. ISME J. 2020;14:164–77. https://0-doi-org.brum.beds.ac.uk/10.1038/s41396-019-0518-5.

    Article  CAS  PubMed  Google Scholar 

  48. Di Capua I, Mazzocchi MG. Population structure of the copepods Centropages typicus and Temora stylifera in different environmental conditions. ICES J Mar Sci. 2004;61:632–44.

    Article  Google Scholar 

  49. Christie AE, Fontanilla TM, Roncalli V, Cieslak MC, Lenz PH. Diffusible gas transmitter signaling in the copepod crustacean Calanus finmarchicus: identification of the biosynthetic enzymes of nitric oxide (NO), carbon monoxide (CO) and hydrogen sulfide (H2S) using a de novo assembled transcriptome. Gen Comp Endocrinol. 2014;202:76–86. https://0-doi-org.brum.beds.ac.uk/10.1016/j.ygcen.2014.04.003.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Christie AE, Fontanilla TM, Roncalli V, Cieslak MC, Lenz PH. Identification and developmental expression of the enzymes responsible for dopamine, histamine, octopamine and serotonin biosynthesis in the copepod crustacean Calanus finmarchicus. Gen Comp Endocrinol. 2014;195:28–39. https://0-doi-org.brum.beds.ac.uk/10.1016/j.ygcen.2013.10.003.

    Article  CAS  PubMed  Google Scholar 

  51. Legrand E, Forget-Leray J, Duflot A, Olivier S, Thomé JP, Danger JM, et al. Transcriptome analysis of the copepod Eurytemora affinis upon exposure to endocrine disruptor pesticides: focus on reproduction and development. Aquat Toxicol. 2016;176:64–75. https://0-doi-org.brum.beds.ac.uk/10.1016/j.aquatox.2016.04.010.

    Article  CAS  PubMed  Google Scholar 

  52. Tarrant AM, Baumgartner MF, Hansen BH, Altin D, Nordtug T, Olsen AJ. Transcriptional profiling of reproductive development, lipid storage and molting throughout the last juvenile stage of the marine copepod Calanus finmarchicus. Front Zool. 2014;11:1–15.

    Article  CAS  Google Scholar 

  53. Roncalli V, Sommer SA, Cieslak MC, Clarke C, Hopcroft RR, Lenz PH. Physiological characterization of the emergence from diapause: a transcriptomics approach. Sci Rep. 2018;8:12577. https://0-doi-org.brum.beds.ac.uk/10.1038/s41598-018-30873-0.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Roncalli V, Cieslak MC, Sommer SA, Hopcroft RR, Lenz PH. De novo transcriptome assembly of the calanoid copepod Neocalanus flemingeri: a new resource for emergence from diapause. Mar Genomics. 2018;37:114–9. https://0-doi-org.brum.beds.ac.uk/10.1016/j.margen.2017.09.002.

    Article  PubMed  Google Scholar 

  55. Christie AE, Roncalli V, Lenz PH. Diversity of insulin-like peptide signaling system proteins in Calanus finmarchicus (Crustacea; Copepoda) – possible contributors to seasonal pre-adult diapause. Gen Comp Endocrinol. 2016;236:157–73. https://0-doi-org.brum.beds.ac.uk/10.1016/j.ygcen.2016.07.016.

    Article  CAS  PubMed  Google Scholar 

  56. Christie AE, Fontanilla TM, Nesbit KT, Lenz PH. Prediction of the protein components of a putative Calanus finmarchicus (Crustacea, Copepoda) circadian signaling system using a de novo assembled transcriptome. Comp Biochem Physiol. 2013;8:165–93. https://0-doi-org.brum.beds.ac.uk/10.1016/j.cbd.2013.04.002.

    Article  CAS  Google Scholar 

  57. Semmouri I, Asselman J, Van Nieuwerburgh F, Deforce D, Janssen CR, De Schamphelaere KAC. The transcriptome of the marine calanoid copepod Temora longicornis under heat stress and recovery. Mar Environ Res. 2019;143:10–23. https://0-doi-org.brum.beds.ac.uk/10.1016/j.marenvres.2018.10.017.

    Article  CAS  PubMed  Google Scholar 

  58. Zhou C, Carotenuto Y, Vitiello V, Wu C, Zhang J, Buttino I. De novo transcriptome assembly and differential gene expression analysis of the calanoid copepod Acartia tonsa exposed to nickel nanoparticles. Chemosphere. 2018;209:163–72. https://0-doi-org.brum.beds.ac.uk/10.1016/j.chemosphere.2018.06.096.

    Article  CAS  PubMed  Google Scholar 

  59. Roncalli V, Cieslak MC, Lenz PH. Transcriptomic responses of the calanoid copepod Calanus finmarchicus to the saxitoxin producing dinoflagellate Alexandrium fundyense. Sci Rep. 2016;6:25708. https://0-doi-org.brum.beds.ac.uk/10.1038/srep25708.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Xie F, Xiao P, Chen D, Xu L, Zhang B. miRDeepFinder: a miRNA analysis tool for deep sequencing of plant small RNAs. Plant Mol Biol. 2012;80:75–84. https://0-doi-org.brum.beds.ac.uk/10.1007/s11103-012-9885-2.

    Article  CAS  Google Scholar 

  61. Lenz PH, Roncalli V, Hassett RP, Wu LS, Cieslak MC, Hartline DK, et al. De novo assembly of a transcriptome for Calanus finmarchicus (crustacea, copepoda) - the dominant zooplankter of the North Atlantic Ocean. PLoS One. 2014;9.

  62. Ning J, Wang M, Li C, Sun S. Transcriptome sequencing and de novo analysis of the copepod Calanus sinicus using 454 GS FLX. PLoS One. 2013;8.

  63. Almada AA, Tarrant AM. Vibrio elicits targeted transcriptional responses from copepod hosts. FEMS Microbiol Ecol. 2016;92:1–11.

    Article  CAS  Google Scholar 

  64. Biessmann H, Le D, Walter MF. Microarray-based survey of a subset of putative olfactory genes in the Anopheles gambiae. Insect Mol Biol. 2005;14:575–89.

    Article  CAS  PubMed  Google Scholar 

  65. Su CY, Menuz K, Carlson JR. Olfactory perception: receptors, cells, and circuits. Cell. 2009;139:45–59.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Swarup S, Williams TI, Anholt RRH. Functional dissection of odorant binding protein genes in Drosophila melanogaster. Genes Brain Behav. 2011;10:648–57.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Swarup S, Morozova TV, Sridhar S, Nokes M, Anholt RRH. Modulation of feeding behavior by odorant-binding proteins in Drosophila melanogaster. Chem Senses. 2014;39:125–32.

    Article  CAS  PubMed  Google Scholar 

  68. Matsuo T. Genes for host-plant selection in Drosophila. J Neurogenet. 2008;22:195–210.

    Article  CAS  PubMed  Google Scholar 

  69. Matsuo T, Sugaya S, Yasukawa J, Aigaki T, Fuyama Y. Odorant-binding proteins OBP57d and OBP57e affect taste perception and host-plant preference in Drosophila sechellia. PLoS Biol. 2007;5:e118 https://0-doi-org.brum.beds.ac.uk/10.1371/journal.pbio.0050118.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  70. Heuschele J, Selander E. The chemical ecology of copepods. J Plankton Res. 2014;36:895–913.

    Article  Google Scholar 

  71. Giordano G, Carbone M, Ciavatta ML, Silvano E, Gavagnin M, Garson MJ, et al. Volatile secondary metabolites as aposematic olfactory signals and defensive weapons in aquatic environments. Proc Natl Acad Sci U S A. 2017;114:3451–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Mollo E, Garson MJ, Polese G, Amodeo P, Ghiselin MT. Taste and smell in aquatic and terrestrial environments. Nat Prod Rep. 2017;34:496–513. https://0-doi-org.brum.beds.ac.uk/10.1039/C7NP00008A.

    Article  CAS  PubMed  Google Scholar 

  73. Kâ S, Carotenuto Y, Romano G, Hwang JS, Buttino I, Ianora A. Impact of the diatom-derived polyunsaturated aldehyde 2-trans,4-trans decadienal on the feeding, survivorship and reproductive success of the calanoid copepod Temora stylifera. Mar Environ Res 2014;93 2014:31–37. doi:https://0-doi-org.brum.beds.ac.uk/10.1016/j.marenvres.2013.07.009.

  74. Michalec FG, Kâ S, Holzner M, Souissi S, Ianora A, Hwang JS. Changes in the swimming behavior of Pseudodiaptomus annandalei (Copepoda, Calanoida) adults exposed to the diatom toxin 2-trans, 4-trans decadienal. Harmful Algae. 2013;30:56–64.

    Article  CAS  Google Scholar 

  75. Leising AW, Pierson JJ, Halsband-Lenk C, Horner R, Postel J. Copepod grazing during spring blooms: can Pseudocalanus newmani induce trophic cascades? Prog Oceanogr. 2005;67:406–21.

    Article  Google Scholar 

  76. Dove SK, Dong K, Kobayashi T, Williams FK, Michell RH. Phosphatidylinositol 3,5-bisphosphate and Fab1p/PIKfyve underPPIn endo-lysosome function. Biochem J. 2009;419:1–13.

    Article  CAS  PubMed  Google Scholar 

  77. Ikonomov OC, Sbrissa D, Delvecchio K, Xie Y, Jin JP, Rappolee D, et al. The phosphoinositide kinase PIKfyve is vital in early embryonic development: Preimplantation lethality of PIKfyve−/− embryos but normality of PIKfyve+/− mice. J Biol Chem. 2011;286:13404–13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  78. O’Neill C. Phosphatidylinositol 3-kinase signaling in mammalian preimplantation embryo development. Reproduction. 2008;136:147–56.

    Article  PubMed  CAS  Google Scholar 

  79. Arai H. Platelet-activating factor acetylhydrolase. Prostaglandins Lipid Mediat. 2002;68–69:83–94.

    Article  Google Scholar 

  80. Kordan W, Strzeżek J, Fraser L. Functions of platelet activating factor (PAF) in mammalian reproductive processes: a review. Pol J Vet Sci. 2003;6:55–60.

    CAS  PubMed  Google Scholar 

  81. Hoornaert I, Marynen P, Baens M. CREBL2, a novel transcript from the chromosome 12 region flanked by ETV6 and CDKN1B. Genomics. 1998;51:154–7.

    Article  CAS  PubMed  Google Scholar 

  82. Kono N, Arai H. Platelet-activating factor acetylhydrolases: an overview and update. Biochim Biophys Acta. 2019, 1864:922–31. https://0-doi-org.brum.beds.ac.uk/10.1016/j.bbalip.2018.07.006.

  83. França CT, CSN LWS, Carmagnac A, Lin E, Kiniboro B, Siba P, et al. IgG antibodies to synthetic GPI are biomarkers of immune-status to both Plasmodium falciparum and Plasmodium vivax malaria in young children. Malar J. 2017;16:386. https://0-doi-org.brum.beds.ac.uk/10.1186/s12936-017-2042-2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Hazenbos W, Förster I, Clausen B, Takeda J, Kinoshita T. In: Omine M, Kinoshita T, editors. Inflammatory defects caused by GPI-anchor deficiency in macrophages BT - paroxysmal nocturnal hemoglobinuria and related disorders. Tokyo: Springer Japan; 2003. p. 247–9.

    Chapter  Google Scholar 

  85. Klumpp S, Krieglstein J. Serine/threonine protein phosphatases in apoptosis. Curr Opin Pharmacol. 2002;2:458–62. https://0-doi-org.brum.beds.ac.uk/10.1016/S1471-4892(02)00176-5.

    Article  CAS  PubMed  Google Scholar 

  86. Leslie SN, Nairn AC. cAMP regulation of protein phosphatases PP1 and PP2A in brain. Biochim Biophys Acta. 1866;2019:64–73. https://0-doi-org.brum.beds.ac.uk/10.1016/j.bbamcr.2018.09.006.

    Article  CAS  Google Scholar 

  87. Tsafriri A, Chun SY, Zhang R, Hsueh AJW, Conti M. Oocyte maturation involves campartmentalization and opposing changes of cAMP levels in follicular somatic and germ cells: studies using selective phosphodiesterase inhibitors. Dev Biol. 1996;178:393–402.

    Article  CAS  PubMed  Google Scholar 

  88. Richard FJ, Tsafriri A, Conti M. Role of phosphodiesterase type 3A in rat oocyte maturation1. Biol Reprod. 2001;65:1444–51.

    Article  CAS  PubMed  Google Scholar 

  89. Lee MS, Bonner JR, Bernard DJ, Sanchez EL, Sause ET, Prentice RR, et al. Disruption of the folate pathway in zebrafish causes developmental defects. BMC Dev Biol. 2012;12:12. https://0-doi-org.brum.beds.ac.uk/10.1186/1471-213X-12-12.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  90. Affleck JG, Neumann K, Wong L, Walker VK. The effects of methotrexate on Drosophila development, female fecundity, and gene expression. Toxicol Sci. 2005;89:495–503. https://0-doi-org.brum.beds.ac.uk/10.1093/toxsci/kfj036.

    Article  CAS  PubMed  Google Scholar 

  91. Stanisławska-Sachadyn A, Brown KS, Mitchell LE, Woodside JV, Young IS, Scott JM, et al. An insertion/deletion polymorphism of the dihydrofolate reductase (DHFR) gene is associated with serum and red blood cell folate concentrations in women. Hum Genet. 2008;123:289–95. https://0-doi-org.brum.beds.ac.uk/10.1007/s00439-008-0475-y.

    Article  CAS  PubMed  Google Scholar 

  92. Kompis IM, Islam K, Then RL. DNA and RNA synthesis: Antifolates. Chem Rev. 2005;105:593–620. https://0-doi-org.brum.beds.ac.uk/10.1021/cr0301144.

    Article  CAS  PubMed  Google Scholar 

  93. Steck SE, Keku T, Butler LM, Galanko J, Massa B, Millikan RC, et al. Polymorphisms in methionine synthase, methionine synthase reductase and serine hydroxymethyltransferase, folate and alcohol intake, and colon cancer risk. Life Genomics. 2008;1:196–204. https://0-doi-org.brum.beds.ac.uk/10.1159/000136651.

    Article  CAS  Google Scholar 

  94. Armant O, März M, Schmidt R, Ferg M, Diotel N, Ertzer R, et al. Genome-wide, whole mount in situ analysis of transcriptional regulators in zebrafish embryos. Dev Biol. 2013;380:351–62. https://0-doi-org.brum.beds.ac.uk/10.1016/j.ydbio.2013.05.006.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  95. Catherino WH, Leppert PC, Stenmark MH, Payson M, Potlog-Nahari C, Nieman LK, et al. Reduced dermatopontin expression is a molecular link between uterine leiomyomas and keloids. Genes Chromosom Cancer. 2004;40:204–17.

    Article  CAS  PubMed  Google Scholar 

  96. Mizuta S, Yoshinaka R, Sato M, Sakaguchi M. Characterization of collagen in muscle of several crustacean species. Comp Biochem Physiol. 1994;107:365–70.

    Google Scholar 

  97. Clegg NJ, Findley SD, Mahowald AP, Ruohola-Baker H. Maelstrom is required to position the MTOC in stage 2-6 Drosophila oocytes. Dev Genes Evol. 2001;211:44–8.

    Article  CAS  PubMed  Google Scholar 

  98. Siomi M, Nagao A, Sato K, Nishida K, Siomi H. Gender-specific hierarchy in nuage localization of PIWI-interacting RNA factors in Drosophila. Front Genet. 2011;2:55. https://0-doi-org.brum.beds.ac.uk/10.3389/fgene.2011.00055.

    Article  PubMed  PubMed Central  Google Scholar 

  99. Pek JW, Ng BF, Kai T. Polo-mediated phosphorylation of maelstrom regulates oocyte determination during oogenesis in Drosophila. Dev. 2012;139:4505–13.

    Article  CAS  Google Scholar 

  100. Thomson T, Liu N, Arkov A, Lehmann R, Lasko P. Isolation of new polar granule components in Drosophila reveals P body and ER associated proteins. Mech Dev. 2008;125:865–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  101. Nakamura A, Amikura R, Hanyu K, Kobayashi S. Me31B silences translation of oocyte-localizing RNAs through the formation of cytoplasmic RNP complex during Drosophila oogenesis. Development. 2001;128:3233–42.

    CAS  PubMed  Google Scholar 

  102. Khan MT, Dalvin S, Waheed Q, Nilsen F, Male R. Molecular characterization of the lipophorin receptor in the crustacean ectoparasite Lepeophtheirus salmonis. PLoS One. 2018;13:–e0195783. https://0-doi-org.brum.beds.ac.uk/10.1371/journal.pone.0195783.

  103. Styhler S, Nakamura A, Swan A, Suter B, Lasko P. Vasa is required for GURKEN accumulation in the oocyte, and is involved in oocyte differentiation and germline cyst development. Development. 1998;125:1569–78.

    CAS  PubMed  Google Scholar 

  104. Tomancak P, Guichet A, Zavorszky P, Ephrussi A. Oocyte polarity depends on regulation of gurken by Vasa. Development. 1998;125:1723 LP – 1732. http://dev.biologists.org/content/125/9/1723.abstract.

  105. Behr M, Hoch M. Identification of the novel evolutionary conserved obstructor multigene family in invertebrates. FEBS Lett. 2005;579:6827–33.

    Article  CAS  PubMed  Google Scholar 

  106. Tajiri R, Ogawa N, Fujiwara H, Kojima T. Mechanical control of whole body shape by a single cuticular protein obstructor-E in Drosophila melanogaster. PLoS Genet. 2017;13:1–26.

    Article  CAS  Google Scholar 

  107. Angers S, Thorpe CJ, Biechele TL, Goldenberg SJ, Zheng N, MacCoss MJ, et al. The KLHL12–Cullin-3 ubiquitin ligase negatively regulates the Wnt–β-catenin pathway by targeting dishevelled for degradation. Nat Cell Biol. 2006;8:348–57. https://0-doi-org.brum.beds.ac.uk/10.1038/ncb1381.

    Article  CAS  PubMed  Google Scholar 

  108. Ravagnan L, Gurbuxani S, Susin SA, Maisse C, Daugas E, Zamzami N, et al. Heat-shock protein 70 antagonizes apoptosis-inducing factor. Nat Cell Biol. 2001;3:839–43.

    Article  CAS  PubMed  Google Scholar 

  109. Aruda AM, Baumgartner MF, Reitzel AM, Tarrant AM. Heat shock protein expression during stress and diapause in the marine copepod Calanus finmarchicus. J Insect Physiol. 2011;57:665–75. https://0-doi-org.brum.beds.ac.uk/10.1016/j.jinsphys.2011.03.007.

    Article  CAS  PubMed  Google Scholar 

  110. Rhee JS, Raisuddin S, Lee KW, Seo JS, Ki JS, Kim IC, et al. Heat shock protein (Hsp) gene responses of the intertidal copepod Tigriopus japonicus to environmental toxicants. Comp Biochem Physiol. 2009;149:104–12. https://0-doi-org.brum.beds.ac.uk/10.1016/j.cbpc.2008.07.009.

    Article  CAS  Google Scholar 

  111. Voznesensky M, Lenz PH, Spanings-Pierrot C, Towle DW. Genomic approaches to detecting thermal stress in Calanus finmarchicus (Copepoda: Calanoida). J Exp Mar Bio Ecol. 2004;311:37–46.

    Article  CAS  Google Scholar 

  112. Asai S, Ianora A, Lauritano C, Lindeque PK, Carotenuto Y. High-quality RNA extraction from copepods for next generation sequencing: A comparative study. Mar Genomics. 2015;24:115–8. https://0-doi-org.brum.beds.ac.uk/10.1016/j.margen.2014.12.004.

    Article  PubMed  Google Scholar 

  113. McMahon HT, Boucrot E. Molecular mechanism and physiological functions of clathrin-mediated endocytosis. Nat Rev Mol Cell Biol. 2011;12:517–33. https://0-doi-org.brum.beds.ac.uk/10.1038/nrm3151.

    Article  CAS  PubMed  Google Scholar 

  114. Lau AW, Chou MM. The adaptor complex AP-2 regulates post-endocytic trafficking through the non-clathrin Arf6-dependent endocytic pathway. J Cell Sci. 2008;121:4008–17.

    Article  CAS  PubMed  Google Scholar 

  115. Lyu Q, Jiao W, Zhang K, Bao Z, Wang S, Liu W. Proteomic analysis of scallop hepatopancreatic extract provides insights into marine polysaccharide digestion. Sci Rep. 2016;6:1–8. https://0-doi-org.brum.beds.ac.uk/10.1038/srep34866.

    Article  CAS  Google Scholar 

  116. Karageorgos L, Brooks DA, Pollard A, Melville EL, Hein LK, Clements PR, et al. Mutational analysis of 105 mucopolysaccharidosis type VI patients. Hum Mutat. 2007;28:897–903. https://0-doi-org.brum.beds.ac.uk/10.1002/humu.20534.

    Article  CAS  PubMed  Google Scholar 

  117. Chakrabarti R, Shepardson S, Karmakar M, Trdan R, Walker J, Shandilya R, et al. Extra-mitochondrial localization and likely reproductive function of a female-transmitted cytochrome c oxidase subunit II protein. Develop Growth Differ. 2009;51:511–9.

    Article  CAS  Google Scholar 

  118. Ramalho-Santos J, Varum S, Amaral S, Mota PC, Sousa AP, Amaral A. Mitochondrial functionality in reproduction: from gonads and gametes to embryos and embryonic stem cells. Hum Reprod Update. 2009;15:553–72.

    Article  CAS  PubMed  Google Scholar 

  119. Zingone A, D’Alelio D, Mazzocchi MG, Montresor M, Sarno D. Time series and beyond: multifaceted plankton research at a marine Mediterranean LTER site. Nat Conserv. 2019;34:273–310.

    Article  Google Scholar 

  120. Ianora A. Scotto di Carlo B, Mascellaro P. reproductive biology of the planktonic copepod Temora stylifera. Mar Biol. 1989;101:187–94.

    Article  Google Scholar 

  121. Schroeder A, Mueller O, Stocker S, Salowsky R, Leiber M, Gassmann M, et al. The RIN: an RNA integrity number for assigning integrity values to RNA measurements. BMC Mol Biol. 2006;7:3. https://0-doi-org.brum.beds.ac.uk/10.1186/1471-2199-7-3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  122. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29:644–52. https://0-doi-org.brum.beds.ac.uk/10.1038/nbt.1883.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  123. Roncalli V, Christie AE, Sommer SA, Cieslak MC, Hartline DK, Lenz PH. A deep transcriptomic resource for the copepod crustacean Labidocera madurae: a potential indicator species for assessing near shore ecosystem health. PLoS One. 2017;12:1–30. https://0-doi-org.brum.beds.ac.uk/10.1371/journal.pone.0186794.

    Article  CAS  Google Scholar 

  124. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011;12:323. https://0-doi-org.brum.beds.ac.uk/10.1186/1471-2105-12-323.

    Article  CAS  Google Scholar 

  125. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2012;29:15–21. https://0-doi-org.brum.beds.ac.uk/10.1093/bioinformatics/bts635.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  126. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550. https://0-doi-org.brum.beds.ac.uk/10.1186/s13059-014-0550-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  127. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B. 1995;57:289–300 http://0-www-jstor-org.brum.beds.ac.uk/stable/2346101.

    Google Scholar 

  128. Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21:3674–6. https://0-doi-org.brum.beds.ac.uk/10.1093/bioinformatics/bti610.

    Article  CAS  PubMed  Google Scholar 

  129. Conesa A, Götz S. Blast2GO: a comprehensive suite for functional analysis in plant genomics. Int J Plant Genomics. 2008;2008:619832. https://0-doi-org.brum.beds.ac.uk/10.1155/2008/619832.

    Article  CAS  PubMed  Google Scholar 

  130. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10. https://0-doi-org.brum.beds.ac.uk/10.1016/S0022-2836(05)80360-2.

    Article  CAS  PubMed  Google Scholar 

  131. Pfaffl MW. Relative expression software tool (REST) for group-wise comparison and statistical analysis of relative expression results in real-time PCR. Nucleic Acids Res. 2002;30:36e.

    Article  Google Scholar 

  132. Pfaffl MW, Horgan GW, Dempfle L. Relative expression software tool (REST) for group-wise comparison and statistical analysis of relative expression results in real-time PCR. Nucleic Acids Res. 2002;30:e36. https://0-doi-org.brum.beds.ac.uk/10.1093/nar/30.9.e36.

    Article  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank Gianluca Zazo, Augusto Passarelli, Ferdinando Tramontano, Marco Cannavacciuolo and all members of the MEDA Unit of the Stazione Zoologica Anton Dohrn for field samples. We also thank Dr. Francesca Margiotta of the MEDA Unit for providing environmental data and the BIOINforMA group of the Stazione Zoologica Anton Dohrn (https://bioinfo.szn.it/) for technical support in sequence deposition.

Funding

Ennio Russo has been supported by a Ph. D. fellowship funded by Stazione Zoologica Anton Dohrn (Open University-Stazione Zoologica Anton Dohrn Ph. D. Program).

Author information

Authors and Affiliations

  1. Stazione Zoologica Anton Dohrn, Villa Comunale, 80121, Naples, Italy

    Ennio Russo, Chiara Lauritano, Diana Sarno, Adrianna Ianora & Ylenia Carotenuto

  2. Consiglio Nazionale delle Ricerche, Institute of Biomolecular Chemistry, Via Campi Flegrei 34, 80078, Pozzuoli, Italy

    Ennio Russo, Giuliana d’Ippolito & Angelo Fontana

  3. Universität zu Köln, Aquatic Chemical Ecology Group, Zülpicher Straβe 47b, D-50674, Cologne, Germany

    Eric von Elert

Authors
  1. Ennio Russo
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Chiara Lauritano
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Giuliana d’Ippolito
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Angelo Fontana
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Diana Sarno
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Eric von Elert
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Adrianna Ianora
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Ylenia Carotenuto
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

YC conceived the study and designed the experiment, provided reagents/materials and was a major contributor in the interpretation of the data and writing the manuscript. ER carried out the experiments, acquired, analysed and interpreted the data, wrote the manuscript and prepared the figures and the Tables. CL provided reagents/materials and contributed to the interpretation of data, GdI provided reagents/material and contributed to the analysis of the data, AF provided reagents/material and contributed to the interpretation of the data, DS contributed to the analysis of the data, EvE contributed to the design of the study and the interpretation of the data, AI provided reagents/material and contributed to the design of the study. All authors reviewed drafts of the manuscript and approved its final version.

Corresponding author

Correspondence to Ylenia Carotenuto.

Ethics declarations

Ethics approval and consent to participate

Not applicable. No approval was required because experimental work was accomplished with an unregulated marine invertebrate.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

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 N50 are also indicated.

Additional file 2 Figure S1:

Blast2Go statistics output for Temora stylifera de novo reference transcriptome assembly (unigenes). 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), log2-Fold-Change (log2-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 up-regulated (positive log2-FC) isoforms.

Additional file 4 Table S3:

List of Reference Genes (RGs) tested in Temora stylifera. Protein name, function, amplicon length (AL) in base pairs (bp), primer sequence and amplification efficiencies (%E) are shown.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Russo, E., Lauritano, C., d’Ippolito, G. et al. RNA-Seq and differential gene expression analysis in Temora stylifera copepod females with contrasting non-feeding nauplii survival rates: an environmental transcriptomics study. BMC Genomics 21, 693 (2020). https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-020-07112-w

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-020-07112-w

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us