The dataset comprised of RNA-Seq data generated for individual female N. flemingeri is publicly available through the National Center for Biotechnology Information (BioProject: PRJNA324453). The time series included weekly or semi-weekly sampling of deep-collected females followed by laboratory incubations starting at collection (Wk0) to 9.5 weeks post-collection (Wk9.5). In a previous study, the data from Wk0 to Wk7 were analyzed to characterize the reproductive program (Roncalli et al., 2018b).

As described previously (Roncalli et al., 2018b), Neocalanus flemingeri adult females were collected from depth on September 20, 2015 in Prince William Sound (station “PWS2”; depth:∼740 m; Latitude 60° 32.1′N; Longitude 147° 48.2′W) during the fall cruise of the Seward Line Long-term Observation Program (LTOP)¹. Briefly, a multiple opening and closing plankton net (0.25 m² cross-sectional area; 150 μm mesh nets; Multinet-Midi, Hydro-Bios) was towed vertically from 700 m. Collections from 400–700 m were immediately diluted, and adult females of N. flemingeri were live sorted and either preserved within 35 min of the retrieval of the net [“Wk0”, n = 6; RNAlater Stabilization Reagent (QIAGEN)] or sorted into 20L carboys for the incubation experiment. Carboys were stored in large insulated coolers and kept in darkness at 5–6°C for transportation to the University of Alaska Fairbanks (UAF).

At UAF, laboratory incubations were set up on September 24, 2015 with 30 sets of four adult females Falcon tissue-culture flask with 600 ml of seawater collected from depth (400–600 m). Flasks were incubated in the dark at 5°C. Between weeks 1 and 7, three tissue-culture flasks were randomly chosen at weekly intervals and harvested. Thereafter, flasks were harvested at semi-weekly intervals to Wk9.5. Females in each harvested flask were assessed for survival and any visible morphological changes associated with oocyte development. In addition, surviving females were preserved in RNALater Stabilization Reagent (QIAGEN) and stored at −80°C. Flasks were also checked for the presence of eggs and/or nauplii, which were first observed at 7.5 weeks with two more clutches released per female between 7.5 and 10 weeks. Overall RNA-Seq was performed on 61 females with six individuals per time point between Wk0 and Wk7 (n = 48), and 13 individuals after females started to release clutches of eggs (Wk7.5: n = 3, Wk8.5: n = 2, Wk9: n = 2 and Wk9.5: n = 6; Supplementary Table S3).

Total RNA extraction and RNA-Seq followed previously described protocols (Roncalli et al., 2018a, b). Briefly, extracted total RNA was extracted from individual females using QIAGEN RNeasy Plus Mini Kit (catalog # 74134)] in combination with a Qiashredder column (catalog # 79654) and shipped on dry ice to the Georgia Genomics Bioinformatics Core². Double-stranded cDNA libraries [Kapa Stranded mRNA-seq kit (KK8420)] were sequenced using a High Output Flow Cell (150 bp, paired end Illumina Next-Seq 500 instrument). Reads with a length ≥50 bp and a Phred score ≥30 were kept for downstream analysis (FASTX Toolkit v.2.0.0; Illumina Basespace Labs) leading to the removal of ca. 7% of reads per sample, and resulting in an average of 20.7 million reads per sample (range: 20 to 32 million reads) (Supplementary Table S1). Quality-filtered libraries (n = 61) were mapped to an existing adult female N. flemingeri reference transcriptome (transcripts = 140,841) (Roncalli et al., 2018a) using kallisto software (default settings; v.0.43.1) (Bray et al., 2016). The BioConductor package edgeR (Robinson et al., 2010) was used to remove transcripts with very low expression levels. This filter removed all transcripts that were expressed at less than 1 cpm in a set of replicate samples. Application of this filter reduced the total number of reference transcripts to 85,505.

The output from the kallisto mapping step was used to compute normalized relative expression per individual female and per transcript as reads per kilobase per million (RPKM) (Mortazavi et al., 2008). Relative expression in RPKM were log-transformed [Log₁₀(RPKM + 1)] and used in a hierarchical clustering step to identify relationships among all females (Wk0 to Wk9.5). The analysis included all 61 individuals and was performed with the function hclust in the base R package applying the average linkage method (UPGMA) and all other default settings (Müllner, 2013). The clusters (CL) resulting in this analysis, which group females with similar expression profiles, were compared and tested for statistical differences in gene expression (see below).

Three separate statistical comparisons were completed to identify differentially expressed genes (DEGs): (1) pairwise comparisons between diapausing females (CL1) and post-diapause females (all other clusters) to establish how relative expression changed over time from emergence to spawning; (2) generalized linear model followed by pairwise comparisons between three clusters of post-diapause females (CL2-CL4; Wk1 to Wk7) to define transcriptional transitions during the reproductive program; and (3) generalized linear model followed by pairwise comparisons between clusters of spawning and post-reproductive (“end of life”) females (Wk7.5 to Wk9.5) to reveal changes in gene expression as females approach end of life.

Statistical analyses were performed on the transcripts mapped to the reference transcriptome using edgeR. Prior to statistical testing, RNA-Seq libraries were normalized using the TMM methods (trimmed means of M values) as implemented by edgeR. Differentially expressed genes between clusters of females in pairwise comparisons were identified using the Fisher’s exact test with p-values corrected for False Discovery Rate (Benjamini–Hochberg procedure) (FDR; p ≤ 0.05). In multiple cluster comparisons DEGs were identified using the generalized linear model (glmFit function) with adjusted p-values (Benjamini–Hochberg procedure) (FDR; p ≤ 0.05) followed by pairwise comparisons between clusters using likelihood ratio test (glmLRT).

For the functional analysis, the annotated N. flemingeri adult female reference transcriptome (transcripts = 59,544) was used to identify DEGs with significant annotation against the Gene Ontology database (GO terms) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (enzyme commission number, EC) (Roncalli et al., 2018a). To identify biological processes (GO terms) that were over-represented among the DEGs, enrichment analyses were performed using topGO software (Alexa and Rahnenfuhrer, 2010). Up- and down-regulated genes with GO annotation identified from each pairwise comparison were independently compared against the 59,544 annotated transcripts in the reference transcriptome (Roncalli et al., 2018a) using the Fisher exact test with a Benjamini–Hochberg correction (p-values ≤ 0.05 (v. 2.88.0, default algorithm weight01). Visualization of the enrichment results was done using ReviGO software (similarity setting to medium = 0.7) (Supek et al., 2011). The software summarizes lists of GO terms removing the redundancy using a simple clustering algorithm that relies on semantic similarity measures.

Additional enrichment analyses of the pathways involved in energy metabolism were performed on the DEGs identified in diapausing vs. post-diapause comparisons (1) using ClusterProfiler software (R package; v. 3.8.1) (Yu et al., 2012). For this analysis, enrichKEGG function (p-values ≤ 0.05 corrected for FDR) was used to compare DEGs (combined up- and down-regulated) annotated with EC numbers against the KEGG database (toType = ncbi-proteinid) using the fruit fly D. melanogaster as reference (organism = dme).

Results from enrichment analyses (GO and KEGG) were used to generate heat maps of DEGs involved in specific metabolic pathways (glycolysis, citric acid cycle (TCA), proteolysis and ubiquitination), cellular homeostasis, immune response, programmed cell death and chromatin silencing. Genes involved in glycolysis, glycerophospholipid metabolism, TCA and ubiquitination were obtained from the KEGG database using the insect D. melanogaster as the reference organism (map00010, map00564, map00020, map04120). Genes annotated with the GO term “cellular homeostasis” [GO:0019725], “immune response” [GO:0006955], “programmed cell death” [GO:0010942, GO:0012501] and “chromatin silencing” [GO:0006342] were searched using AMIGO software GOOSE (May, 2019; Ashburner et al., 2000). As implemented by the software, the selected GO terms and their descendants were searched using the LEAD SQL wiki called find descendants of the node “nucleus” with “nucleus” replaced with the specific GO term name (e.g., “cellular homeostasis”). Heat maps of the genes involved in each pathway/process of interest were generated by averaging relative expression [Log₁₀(RPKM + 1)] of all individuals within a cluster and comparing across clusters. Statistical significance for differential expression of the genes among the clusters was confirmed by searching the genes in the list of DEGs.

Survival rate in the harvested flasks averaged 87% (range: 70 to 100%) during the 9.5-week incubation period with the lowest survival recorded in flasks harvested at Wk7.5. The first morphological changes associated with oogenesis were noted in a few females harvested at Wk3 (Roncalli et al., 2018b). Visible changes included the appearance of light coloration in the dorsal region indicative of early oogenesis. By Wk8.5 all females were in a late stage of oogenesis: dark-colored eggs were clearly visible in the oviducts along the mid-dorsal and posterio-lateral regions of the prosome. Furthermore, by this time eggs and/or nauplii were present in the flasks suggesting that most females had released at least one clutch of eggs. All individuals harvested in Wk9 had either many mature eggs in their oviducts, or they appeared to be “partially spent” as indicated by a smaller number of eggs in the oviducts. At the last sampling time point (Wk9.5), 75% of the females in the flasks were still alive, but most were inactive, and several individuals were “spent” with no or very few mature eggs remaining in the oviducts.

Agnostic cluster analysis of individuals by normalized expression profiles (transcripts >1 cpm = 85,505) segregated the 61 females into two major branches (Figure 1). The “early” branch included all females that had not started to release eggs (n = 48, Wk0 to Wk7), while females in the “late” branch had started to spawn (n = 13; Wk7.5 to Wk9.5). Within each branch females were grouped into clusters, which were named based on the stage in the reproductive program (Roncalli et al., 2018b). In the “early” branch, females aggregated into three clusters (CL1, CL2, and CL3). “Diapause” (CL1) included the six Wk0 individuals, which were clearly separated from all other females within the early branch. The remaining 42 individuals separated into CL2 “early to mid-oogenesis” with 29 individuals from the Wk1 to Wk5 harvests, CL3 “pre-spawning” with 12 individuals from the Wk6 and Wk7 harvests (Figure 1). One outlier, a Wk5 individual was excluded from all downstream analyses.

The 13 females in the “late” branch separated into clusters as well. One cluster was composed of two sub-clusters (CL4a and CL4b), which included females with eggs in the oviduct. Females from CL4 “spawning” separated in CL4a with eight individuals (Wk7.5 to Wk9.5) and CL4b with three Wk9.5 individuals. Individuals from these two sub-clusters (CL4a and CL4b) were combined in comparisons with diapausing females (CL1) and in the characterization of transitions between major clusters. They were only considered as separate groups for the analysis focused on “end-of-life” processes. The second cluster (CL5) included two spent (= post-reproductive) Wk9.5 females. The grouping of only two females in CL5 named “end-of-life” is consistent with microscopic observations prior to preservation for RNA-Seq, which identified these two females as “spent.”

Paired comparisons between diapause (CL1) and post-diapause females (CL2-CL4) underscored the large transcriptional differences associated with emergence from diapause (Table 1 and Supplementary Table S2). A high number of DEGs (n = 7,224) were identified in the comparison between “diapause” (CL1) and “early to mid-oogenesis” (CL2) females, marking the termination of diapause and the beginning of oogenesis. The number of DEGs increased to more than 10,000 between “diapause” (CL1) and “pre-spawning” (CL3) and “spawning” (CL4) females (pairwise comparisons, p ≤ 0.05 after FDR correction) (Supplementary Table S2). While ‘reproductive process’ was common among the DEGs in all comparisons, other processes represented over 50% of the GO terms among the annotated DEGs (Supplementary Figure S1). These other GO terms included “single organism process,” “multicellular organismal process,” “biological regulation,” and “localization”.

The first transition between “diapause” (CL1) and “early to mid-oogenesis” (CL2) females was characterized by a general activation of cellular processes with four being enriched among the up-regulated genes with no enriched down-regulated processes (Table 1 and Figure 2A). Genes involved in “reproductive process,” “cellular redox homeostasis,” “response to stimulus,” and “intracellular signal transduction” were all up-regulated, consistent with a physiologically more active state and the initiation of the reproductive program in “early to mid-oogenesis” (CL2) females (Table 1 and Figure 2A).

Three GO terms remained enriched in the paired comparisons between “diapause” (CL1) and “pre-spawning” (CL3) and “spawning” (CL4), however, some enriched processes switched from up-regulated to down-regulated (Table 1). Notable among these processes was “reproductive process,” which switched to down-regulated in “spawning” (CL4) females. Other processes enriched among down-regulated genes were “immune response,” and “muscle contraction” (Table 1). Processes that were enriched among the up-regulated genes in “pre-spawning” (CL3) and “spawning” (CL4) females included “proteolysis,” “chromatin silencing,” and “protein ubiquitination” (Table 1).

An analysis of the enzymatic pathways represented among the DEGs showed a temporal progression of enriched pathways in pairwise comparisons between diapause and post-diapause clusters of females. Significant over-representation of DEGs involved in glycolysis and glycerophospholipid metabolic pathway characterized the “early to mid-oogenesis” (CL2) females in comparison with the “diapause” (CL1) females. Consistent with the enrichment of the “proteolysis” GO term, “pre-spawning” (CL3) females showed an enrichment of DEGs involved in enzymatic pathways of protein processing. In the “spawning” (CL4) females, DEGs involved in the ubiquitin pathway were enriched among the up-regulated genes (Table 1).

Several genes that are often used as indicators of diapause termination in insects were also significantly regulated in N. flemingeri in the transition between “diapause” (CL1) and “early to mid-oogenesis” (CL2) females (Figure 3A). “Diapause termination” genes such as proliferating cell nuclear antigen (PCNA), cyclin-dependent kinase 1, cyclin-B, tribbles, sestrin, and 5′-AMP-activated protein kinase (AMPK) were differentially expressed between “diapause” (CL1) and “early to mid-oogenesis” (CL2) females (Figure 3A). The expression of these genes was also lower in females in the ‘pre-spawning’ and “spawning” clusters (CL3, CL4), but not significantly different compared with “early to mid-oogenesis” (CL2) females (Figure 3A). Low expression was also observed for several anti-longevity genes (e.g., edem 1, sirtuin 6) in CL1 females (Figure 3A). In contrast, relative expression of several transcripts annotated to phosphoenolpyruvate carboxy-kinase (PEPCK), a metabolic enzyme involved in gluconeogenesis, and pro-longevity genes (e.g., embryonic lethal abnormal vision, cheerio) was low in all post-diapause females (CL2, CL3, and CL4) compared with “diapause” (CL1) females (Figure 3A).

The three clusters of post-diapause females (CL2 to CL4) represented major transitions in transcriptional profiles associated with “cellular redox homeostasis” and “metabolic process” in addition to “reproductive process.” A total of 13,788 DEGs were identified among the clusters of post-diapause females (CL2-CL4) with ca. 7,000 DEGs between “early to mid-oogenesis” (CL2) and “pre-spawning” (CL3), and ca. 12,000 DEGs between ‘pre-spawning’ (CL3) and ‘spawning’ (CL4) females (Supplementary Table S2). “Reproductive process” was the top GO term among the DEGs (37 to 45%) between adjoining clusters of post-diapause females, consistent with sequential changes in expression related to the progression of the reproductive program described by Roncalli et al. (2018b) (Supplementary Figure S1C). The remaining DEGs (>55%) were involved in conserved processes such as “single organism process,” “multicellular organismal process,” “biological regulation,” and “localization” (Supplementary Figure S1C).

Females grouped in “pre-spawning” (CL3) cluster were on average 2 weeks further along the reproductive program than females within the “early to mid-oogenesis” (CL2) cluster (Table 1). Eight GO terms belonging to five functional groups were enriched among the DEGs between “early to mid-oogenesis” (CL2) and “pre-spawning” (CL3) females, as summarized by reviGO analysis (Supplementary Table S3 and Figure 2B). The transition was characterized by an overall down-regulation in the “pre-spawning” (CL3) females compared with those in “early to mid-oogenesis” (CL2). Six out of eight process were enriched among the down-regulated genes in the “pre-spawning” (CL3) females and this included “reproductive process” as well as multiple GO terms linked to organismal and cellular maintenance, such as “response to stress,” “cellular redox homeostasis” and “immune response” (Supplementary Table S3 and Figure 2B). As part of cellular redox homeostasis, expression of genes like glutaredoxin and sulfhydryl oxidase was lower in “pre-spawning” (CL3) and “spawning” (CL4) females compared with females from “diapause” (CL1) and “early to mid-oogenesis” (CL2) (Figure 3B). Expression of two peroxidasins were high in diapause females and continued to be up-regulated in “early to mid-oogenesis” (CL2) compared with “pre-spawning” (CL3) and “spawning” (CL4) females. Relative expression of genes involved in immune response (c-type lectins, a single protein Toll and embryonic polarity dorsal protein) was high in “early to mid-oogenesis” (CL2) compared with “pre-spawning” (CL3) females, and remained low in “spawning” (CL4) females (Figure 3B).

The transition from “early to mid-oogenesis” (CL2) to “pre-spawning” (CL3) included the up-regulation of genes involved in chromatin modification, as shown by the enrichment of two GO terms associated with this process (Supplementary Table S3 and Figure 2B). Significant up-regulation of seven longitudinal-lacking proteins, involved in chromatin condensation was found in “pre-spawning” (CL3) females compared with females from the previous two clusters (Figure 3B). High expression in “pre-spawning” (CL3) females was found also for a single protein MRG15 and a longitudinal-lacking protein although their expression was already high in “early to mid-oogenesis” (CL2) females compared with “diapause” (CL1) ones (Figure 3B). Expression of all these genes continued to be high also in “spawning” (CL4) females (Figure 2C).

Females grouped in CL4 were actively spawning and had been incubated for an average of 8.8 weeks in the lab. Despite the high number of DEGs (>12,000), only two GO terms were enriched in the comparison between “pre-spawning” (CL3) and “spawning” (CL4) females: “regulation of programmed cell death” and “protein ubiquitination” (Figure 2C and Supplementary Table S3). Both processes were enriched among the up-regulated genes. Relative expression of seven genes involved in the “regulation of programmed cell death,” which included several transmembrane BAX proteins, two metacaspases, a single caspase 1, was significantly higher in “spawning” (CL4) females compared with all other clusters (Figure 3B). In contrast, relative expression for the apoptosis inducing factor 3 was high in “early to mid-oogenesis” (CL2) and “pre-spawning” (CL3) females but not in “spawning” (CL4) (Figure 3B).

Transitions between clusters were characterized by large changes in the relative expression of genes involved in energy metabolism. The progressive change in expression is shown in a heat map that visualizes how sets of genes associated with specific pathways were up-regulated in one or two clusters (Figure 4). With the exception of PEPCK (Figure 3A), expression of genes involved in energy metabolism was low in “diapause” (CL1) females (Figure 4). The up-regulation of these genes is consistent with an increase in basal metabolism as females emerged from diapause. In addition, the data suggest that females differentially regulated energy sources as they progressed through the reproductive program.

The increase in transcriptional activity between “diapause” (CL1) and “early to mid-oogenesis” (CL2) females was marked by high expression of all nine glycolytic enzymes (hexokinase, phosphoglucose isomerase, phosphofructokinase, fructose-bisphosphate aldolase, glyceraldehyde 3-phosphate dehydrogenase, phosphoglycerate kinase, phosphoglycerate mutase, enolase and pyruvate kinase), which were significantly higher in CL2 females compared with those in all other clusters (Figure 4). A similar pattern was observed in the expression of genes involved in glycerophospholipid metabolism, an enriched pathway. Relative expression of 12 phospholipases (class A1, A2, B) and a single monoacylglycerol lipase (MGL) was significantly higher in “early to mid-oogenesis” (CL2) females compared with all other clusters (Figure 4). Females within the ‘early to mid-oogenesis’ (CL2) cluster were also characterized by the up-regulation of genes involved in the citric cycle (TCA), which regulates metabolic rate (Figure 4). The eight TCA enzymes, citrate synthase, aconitase, isocitrate dehydrogenase, α-ketoglutarate dehydrogenase, succinyl-CoA synthetase, succinate dehydrogenase, fumarase, malate dehydrogenase and citrate synthase, were up-regulated in “early to mid-oogenesis” (CL2) and “pre-spawning” (CL3) females in comparison with both “diapause” (CL1) and “spawning” (CL4) ones (Figure 4).

The significant high expression of 17 proteolytic genes such as carboxypeptidases, metallo-peptidases and amino peptidases in “pre-spawning” (CL3) females compared with all other clusters, is consistent with the proteolysis signal identified by enrichment analysis. Eleven genes were exclusively differentially expressed in ‘pre-spawning’ (CL3) females, and these included carboxypeptidases (6 of 7), three proteases (ADAM, cysteine, serine easter), a single laminin, and matrix metallo-peptidase 9 (Figure 4). Two aminopeptidases two metallo-peptidases and a single carboxypeptidase D2 were up-regulated in “pre-spawning” (CL3) and “spawning” (CL4) females, although relative expression of these genes tended to be lower in CL4 females compared with CL3 ones. A single disintegrin was highly expressed in “pre-spawning” (CL3) females as well as “early to mid-oogenesis” (CL2) and “spawning” (CL4) females (Figure 4).

The expression pattern of ‘spawning’ (CL4) females was distinct from all others: relative expression of most of the genes involved in energy metabolism was either low or intermediate with the exception of the high expression found for genes involved in ubiquitination (Figure 4). Relative expression of E1, E2, and E3 ubiquitins was significantly higher in “spawning” (CL4) females compared with expression in “diapause” (CL1), “early to mid-oogenesis” (CL2) and “pre-spawning” (CL3) (Figure 4). An exception to this pattern, the E3 ubiquitin-protein ligase ARIH1, which catalyzes ubiquitination of target proteins together with ubiquitin-conjugating enzymes, was highly expressed in both pre- and spawning females (CL3 and CL4) (Figure 4).

The “spawning” and “end of life” phases comprising of females in CL4a (n = 8), CL4b (n = 3) and CL5 (n = 2) were analyzed separately to elucidate transcriptional transitions that occurred during the release of multiple clutches of eggs (CL4a, CL4b) and post-reproduction (CL5) (Figure 1). A total of 2,023 DEGs were identified among the females in the three clusters (GLM) with more than 700 DEGs identified between CL4a and CL4b and CL4a and CL5 and 583 DEGs between CL4b and Cl5 (Supplementary Table S2). Four GO terms were enriched among the DEGs (Figure 5A). Two of the GO terms, “positive regulation of programmed cell death” and “protein ubiquitination,” were also enriched in the transition between “pre-spawning” (CL3) and “spawning” (CL4) females (Figure 3B), while the other two GO terms were specific to transcriptional changes among clusters in the late branch (“phosphorylation” and “phosphate-containing compound metabolic process”; Figure 5A). Relative expression of genes involved in these processes is shown in Figure 5B. For comparison, this figure also includes relative expression of these genes in “pre-spawning” (CL3) females (first column on left). The transition between CL4a and CL4b females was characterized by an overall decrease in the expression of genes involved in phosphorylation (e.g., casein kinase, CaM kinase, aurora kinase), cellular integrity [collagen alpha-2, phosphatidylinositol-3,5-bisphosphate 3-phosphatases (ptdlns)] and cellular redox homeostasis (GSH synthetase and DnaJ homolog C16). With the exception of a single flotillin (structural integrity), all genes shared a similar pattern of expression: low expression in “pre-spawning” (CL3), significant up-regulation in CL4a, followed by lower expression in CL4b, which decreased even further in CL5.

The transition between CL4a and CL4b was characterized by the up-regulation of genes involved in “programmed cell death” (Figure 5B). Two programmed cell-death proteins (2, 5), and a growth hormone-inducible transmembrane protein were significantly up-regulated in CL4b compared with CL4a and CL5 females (relative expression: CL4b > CL5 > CL4a) (Figure 5B). Relative expression of four other genes involved in this process were lower in CL4a compared with Cl4b and Cl5, and this difference was significant (Figure 5B). Genes involved in “protein ubiquitination,” which characterized the transition between “pre-spawning” (CL3) and “spawning” (CL4) females (Figure 3B) became increasingly up-regulated with cluster number (Figure 5B). Relative expression of all four genes, E3 ubiquitin ligase, E2 conjugating enzyme and ubiquitin carboxyl-terminal hydrolase calypso was low in CL4a females in comparison with CL4b ones, which was followed by another increase in expression in CL5 females (Figure 5B).