Effects of liquid cultivation on gene expression and phenotype of C. elegans

Background Liquid cultures have been commonly used in space, toxicology, and pharmacology studies of Caenorhabditis elegans. However, the knowledge about transcriptomic alterations caused by liquid cultivation remains limited. Moreover, the impact of different genotypes in rapid adaptive responses to environmental changes (e.g., liquid cultivation) is often overlooked. Here, we report the transcriptomic and phenotypic responses of laboratory N2 and the wild-isolate AB1 strains after culturing P0 worms on agar plates, F1 in liquid cultures, and F2 back on agar plates. Results Significant variations were found in the gene expressions between the N2 and AB1 strains in response to liquid cultivation. The results demonstrated that 8–34% of the environmental change-induced transcriptional responses are transmitted to the subsequent generation. By categorizing the gene expressions for genotype, environment, and genotype-environment interactions, we identified that the genotype has a substantial impact on the adaptive responses. Functional analysis of the transcriptome showed correlation with phenotypical changes. For example, the N2 strain exhibited alterations in both phenotype and gene expressions for germline and cuticle in axenic liquid cultivation. We found transcript evidence to approximately 21% of the computationally predicted genes in C. elegans by exposing the worms to environmental changes. Conclusions The presented study reveals substantial differences between N2 and AB1 strains for transcriptomic and phenotypical responses to rapid environmental changes. Our data can provide standard controls for future studies for the liquid cultivation of C. elegans and enable the discovery of condition-specific genes. Electronic supplementary material The online version of this article (10.1186/s12864-018-4948-7) contains supplementary material, which is available to authorized users.


Background
All living systems possess the fundamental property to regulate physical or genetic states in response to environmental stimuli. Caenorhabditis elegans offers one of the unique models to study the effects of environmental changes as it can reproduce in large numbers and has a natural ability to live in both solid and liquid conditions. This eukaryotic organism provides the opportunity to study genetic mechanisms in the whole animal rather than a cell culture. In addition, it can consume a bacterial or an axenic food source thereby enabling the identification of diet-related biological processes. Placing the worms from bacteria seeded agar plates (OP50 NGM) to liquid axenic media poses a drastic change for the animals. For example, the animals display locomotive differences in the liquid by swimming via a rapid C-shape thrashing motion as opposed to crawling via a sinusoidal configuration on solid [1]. Behavior and physiological differences arise when the worms are grown in liquid axenic media (such as CeMM; C. elegans Maintenance Medium or CeHR; C. elegans Habitation and Reproduction) [2,3]. In axenic media, the worms exhibit phenotypical alterations such as delayed development, reduced fecundity, increased longevity, and altered body morphology [4][5][6]. Environment changes can trigger different gene expression profiles, and these profiles can be transmitted to the subsequent generations [7][8][9]. Changing the environment of the animals from agar to axenic media, therefore, can cause transcriptome-level alterations which may have a sustained effect on the progenies. However, the knowledge about transcriptomic responses to axenic media is limited. For example, in their microarray analysis, Szewczyk et al. could manage to obtain reliable data for only 1202 cDNAs on all three arrays consistently for CeMM-grown worms [5]. The common usage of axenic media in space, toxicology, and pharmacology studies highlights the necessity of a complete transcriptomic analysis, preferably with RNA-sequencing (RNA-seq) [6,10].
In its natural habitat, C. elegans occupies microbiota-rich decaying plant material and experiences a continuous environment shift between liquid and solid conditions along with the variations in its diet [11]. The wild-type laboratory strain of C. elegans (N2) has been cultured in both liquid and solid conditions [12], but bacteria-seeded agar plates are used more commonly [13]. It has been reported that the domestication of the N2 strain has created a selective pressure on the animals [14]. Some of the critical questions that await answers are whether: 1.) the adaptive responses depend on the genotype or the environment change; 2.) the environment change-induced gene expression profiles are transmitted to the successive generation in equal levels for different genotypes; 3.) the environment change-induced gene expression profiles are transmitted to the next generation equally for different environmental conditions.
To study the maternal effects on transcriptomic responses of C. elegans from solid to liquid environments, we maintained the P 0 worms on bacteria-seeded agar plates and F 1 in liquid cultures. By culturing F 2 generation back on agar, we identified the maintained maternal transcriptional responses. The N2 strain and a wild isolate AB1 strain were used to examine if a different genotype attributes to the transcriptional dynamics. It is important to note that even though the N2 and the AB1 strains are relatively close to each other in the phylogenetic tree [15], AB1 is not the ancestral strain of N2. Therefore, the observed differences in norm of reaction can stem from the natural genetic differences between N2 and AB1 as much as the impact of laboratory domestication on adaptive responses. We found substantial variances in gene expression profiles resulting from different genotypes, environments, and genotype-environment interactions. Our results from RNA-seq data analysis and microscopy showed consistency for the predicted and observed phenotype. Furthermore, we demonstrated that a large number of genes are expressed only in specific growth conditions, and many of these genes previously lacked experimental expressed sequence tag (EST) or cDNA evidence.

Results
To induce an adaptive response and observe its sustained effects on the following generation, the P 0 generation C. elegans was grown on agar, the F 1 generation in liquid, and the F 2 generation back on agar (Additional file 1: Figure  S1a). We carried out three sets of experiments with different food sources and animal strains. For the first experiment, the F 1 generation N2 strain was grown in axenic CeHR Medium. To examine the effect of dietary changes in the liquid cultures, we cultured the F 1 generation N2 strain in bacterial S-Medium for the second experiment. Finally, the transcriptomic response differences were investigated by growing F 1 generation AB1 strain in CeHR Medium. For each generation, RNA-seq was performed on young adult animals with three replicates to determine the transcriptional responses (Additional file 1: Table S1).
Distinct gene groups are expressed exclusively in specific environmental conditions The transcriptome-based clustering of the experiments in the dendrogram of Jensen-Shannon distances demonstrated that the genotypes are separated from each other (Fig. 1a). Exposure to different environmental conditions causes the N2 transcriptome to be separated even further than the AB1 transcriptome as illustrated by the dendrogram and principal component analysis ( Fig. 1a, b, Additional file 1: Figure S1b). These results suggest that the genotype plays a key role in adaptive response, and the different strains of the animals can show highly distinct transcriptome profiles under environmental changes. We categorized the genes by using Venn diagrams to identify if their expressions are condition specific (Fig. 1c). We considered the genes expressed if their Fragments Per Kilobase of transcript per Million mapped reads (FPKM) values are greater than one as this is a commonly used criterion for expression [16][17][18][19]. The majority of the expressed genes are conserved among the P 0 , F 1 , and F 2 generations. Nevertheless, hundreds of genes are expressed only under particular conditions. The number of the uniquely expressed genes showed a significant difference between CeHR grown AB1 and N2 strains, but not within the N2 strain experiments (Chi-squared test with Yates correction; P-value = 0.01 and P-value = 0.9, respectively). In addition, the N2 strain showed an overall decrease in the transcription levels when exposed to the changes in growth conditions (Additional file 1: Figure S2). AB1 strain, however, demonstrated an increase in the total FPKM but a decrease in the average FPKM under CeHR (Additional file 1: Figure S2). Interestingly, the number of gene expression patterns transmitted to the next generation after exposure to an environment change (F 1 F 2 ) varies among the experiments (Chi-squared test with Yates correction; P-value < 0.001). Together, these results seem to suggest that the transmission of the gene expression is not exclusively dependent on the genotype and the environment separately-but dependent on their interaction, whereas the condition-specific expression of genes is affected by the genotype.
Further analyses revealed that non-coding RNA (ncRNA) molecules are enriched (16-78%) in our pool of condition-specific transcripts indicating their function in rapid adaptive responses (Fig. 1c). To note, only the polyadenylated ncRNA molecules were included in the analysis, and piRNA and miRNA molecules were excluded from the analysis due to their short lengths. The ncRNA molecules were predominantly observed in the P 0 generation of N2 strain compared to the P 0 AB1 strain. Moreover, fewer ncRNA molecules were expressed when the animals are exposed to an environment change (Fig. 1c, d, and Additional file 1: Figure S3). Our results indicate that ncRNA expression is mainly silenced when the worms are exposed to the tested environmental changes. The fact that a laboratory condition presents an environment change for the wild-isolate may explain why P 0 animals of AB1 showed less ncRNA expression compared to the N2 strain. One particular exception to this pattern is the ncRNA molecules that are commonly expressed in all three generations. We identified that the majority of these ncRNAs (59-66%, n = 322) are conserved among the three experiments (Additional file 2: Dataset S1). This conservation implies that this group of ncRNAs plays important regulatory functions in the cell and their expression is required regardless of the environmental changes. Overrepresented phenotype for these ncRNAs is dendrite development variant (FDR adj. p-value < 0.01, six-fold enrichment) suggesting a potential role for them in generation of neuronal extensions.
We reasoned that if a particular group of genes is only expressed under specific conditions, we might detect transcript evidence to computationally predicted but experimentally unconfirmed genes, especially considering that CeHR is not a common medium used in C. elegans research. Notably, we were able to identify transcript evidence to 21% of the list of unconfirmed genes from WormBase [20] by simply growing the animals in different environmental conditions (Fig. 1e). Unexpectedly, a substantial amount of the transcript evidence was found with animals raised in commonly used OP50 NGM. We compared our list of expressed unconfirmed genes in OP50 NGM-grown N2 to those from a previous study [21] and found that 54% of the unconfirmed genes from our list (n = 292) have also shown expression in the study by Dallaire et al. but were not reported as transcript evidence in WormBase. The remaining 46% of the genes must have been induced due to the subtle differences introduced in laboratory environment (e.g., light and environment). Similar to the confirmed genes, gene expression transmission to the successive generation after environmental change demonstrated a difference between and within the strains (Chi-squared test with Yates correction; P-value < 0.001, P-value = 0.02, respectively). Given that it can robustly provide transcript evidence to unconfirmed genes, our approach and data can be a rich resource for discovery of novel genes.

Environment, genotype, and environment-genotype specific expression of genes
We grouped the genes based on their genotype, environment, or genotype-environment interaction specific expressions ( Fig. 2 and Additional file 3: Dataset S2). For this analysis, a more stringent criterion was used to eliminate any potential noise from the data. That is, we only considered the genes with expression levels higher than pmp-3 of FPKM = 21, and excluded genes with FPKM < 21 from the analysis. This housekeeping gene was selected as a reference because of its stable expression levels found both in our analysis and previous studies [22,23]. Specifically for the ncRNA molecules, we reasoned that expression levels higher than that of a housekeeping gene should provide strong evidence for expression.
We defined the genes that are expressed in only one environmental condition for both the strains as "environment-specific", and found that CeHR-specific gene expression displayed Gene Ontology (GO) enrichment for the neuropeptide signaling pathway. Given that this pathway functions in numerous behavioral activities such as reproduction, locomotion, mechanosensation and chemosensation, it is likely that the neuropeptide signaling pathway plays a crucial role in rapid adaptation to CeHR and is evolutionarily conserved between the strains [24].
The "genotype-specific" expression represents the genes found in only one of the strains under a single environmental condition. All the N2 specific genes showed GO enrichment for body morphogenesis in all the three generations while GO enrichment for reproduction and locomotion were observed in P 0 and F 1 generations (Fig. 2). Previously, we reported that the N2 worms produce less progeny compared to the AB1 animals in CeHR [6]. Our analysis showed that along with the reproduction genes, genes related to hermaphrodite genitalia development were expressed exclusively in the N2 strain in CeHR. This finding may suggest vulva aberrations in the reproductive system. In addition, the exclusive expression of the genes functioning in positive regulation of multicellular organism growth may help the animals readapt to the OP50 NGM condition in the F 2 generation of the N2 animals.
The genes expressed solely under a particular environment and strain were categorized as "genotype-environment interaction" genes ( Fig. 2). This group of genes did not demonstrate an enrichment of GO terms. The proportion of ncRNA molecules in our group of genes showed correlation with our hypothesis that the majority of the ncRNAs are silenced under the tested environmental changes. For example, the N2 worms showed relatively higher proportions of ncRNAs for interaction and genotype-related genes on OP50 NGM, but these proportions were much lower in more stressful CeHR and CeHR reversion (Fig. 2).

N2 strain presents higher differential expression under environment changes
To investigate the differences in transcriptional responses between the strains and among the environmental conditions, we identified the differentially expressed genes (DEGs). We considered genes as differentially expressed if their FPKM is greater than one and their expression is altered over two-folds among the conditions with FDR-adjusted p-values ≤ 0.05. The estimated fold change responses of 20 randomly selected transcripts (ten for OP50 NGM against CeHR and ten for CeHR against CeHR reversion) were confirmed by quantitative RT-PCR (Additional file 1: Figure S4a, b and Additional file 1: Table   Fig. 2 Shown are the number of gene expressions exclusive to environment, genotype, or genotype-environment interactions. The number of the expressed ncRNAs are decreased under stress conditions. The table depicts the enriched gene ontology (GO) terms assigned to environment or genotype specifically expressed genes. Genotype-environment interaction genes did not demonstrate enrichment for GO terms. "None" represents the groups with no GO enrichment. Benjamini Hochberg corrected P-values for the GO terms: * < 0.01; ** < 0.001 S2 and S3). We found the highest number of DEGs between OP50 NGM and CeHR conditions for the N2 strain (Fig. 3). This indicates that the change in diet and physical environment triggers more transcriptional responses for adaptation of the N2 strain. The number of DEGs were much lower for the F 2 generation of the same experiment, suggesting that the majority of the gene expressions started to resemble the P 0 animals. This finding was expected as most of the gene expression patterns return to the steady-state levels after a short period [25]. The change in only the physical condition introduced by the S-Medium apparently did not cause as drastic transcriptional responses as the CeHR since the number of the DEGs were lower (38 and 45% lower number of upregulated and downregulated genes, respectively). However, S-Medium reversion conditions demonstrated higher DEGs compared to OP50 NGM. We note that we are able to maintain the worms in S-Medium for only one generation, potentially because of the toxins released by the bacteria in the liquid environment [26,27]. Therefore, additional genetic mechanisms may be required for re-adaptation of the animals from S-Medium. Transcriptional variations were more modest across the generations for AB1 animals than the N2 strain.
Substantial discrepancies have been reported for the C. elegans gene expressions among axenic media studies [28]. We wished to compare our results to those from CeMM-grown worms as CeMM has a very similar base medium to CeHR. For the 21 upregulated and 26 downregulated genes reported in a previous study [5], we could only identify one common downregulated gene in N2 and AB1 and four and seven common upregulated genes in N2 and AB1, respectively (Additional file 1: Table S4). It is unclear why the DEGs among the axenic media studies mostly disagree, but some potential reasons can be the difference in the technology used and the modifications made in the media. For instance, microarray technology has been reported to introduce Fig. 3 The differential expression of the genes (FPKM > 1, FDR adj. P-value < 0.05, and log2(fold change of the FPKM values) ≥ 2) between the environment conditions and the strains. Red dashed lines represent the two-fold cutoff for differential gene expression. The corresponding genes for the enriched GO terms were represented with the colors. Benjamini Hochberg corrected and FDR adj. P-values < 0.05 for all the enriched GO terms. The GO terms were determined separately for downregulated and upregulated genes more background noise and cross-hybridization and detect a lesser number of DEGs compared to RNA-seq [29]. A more unified study conducted with RNA-seq may help resolve the discrepancy issues among the axenic media experiments.
In our previous study, we have shown that CeHR-grown N2 and AB1 C. elegans display different phenotypical traits on their body morphology [6]. That is, adult animals of both strains demonstrated physical differences in CeHR compared to the OP50 NGM-grown ones, but the CeHR-grown AB1 strain was even longer and thinner than the CeHR-grown N2 strain in average (N2 length is 825.2 ± 74.3 μm and width is 45.9 ± 10.1 μm, and AB1 length is 876.8 ± 67.0 μm and width is 38.3 ± 4.7 μm). In accordance with this finding, the GO enrichment analysis results demonstrated enrichment for body morphogenesis for upregulated genes in CeHR compared to OP50 NGM in the N2 strain, but not the AB1 strain (Fig. 3). Moreover, genes upregulated under the CeHR condition in the AB1 strain compared to the N2 strain were overrepresented for body morphogenesis. The fact that this GO term is not enriched for AB1 may imply that the phenotype is caused by a separate transcriptional mechanism. To test this, we acquired the genes associated with "long" and "thin" phenotypes from WormBase [20] and compared these genes to the DEGs in our list. These phenotypes are generally found by RNAi knockdown of genes. We found that six and nineteen differentially expressed genes functioning in long or thin phenotypes in N2 and AB1, respectively, but only one gene, his-67, was common between the strains (Additional file 1: Table S5). Thus, it is possible that the different gene expression profiles contribute to the body morphology variations between the strains.
Neuropeptide signaling pathway genes were enriched in response to environmental changes. The pathway genes were upregulated in the liquid environments compared to agar in N2 and AB1. The differential expression of this pathway genes in liquid conditions did not maintain their expression patterns in the next generation. This finding supports our hypothesis that neuropeptide signaling pathway can be a crucial mechanism in rapid adaptive responses and it is conserved in the different strains of C. elegans.
Given the component differences in OP50 NGM, CeHR, and S-Medium, it is expected to observe variations in the expressions of metabolic process genes among these conditions. We performed pathway enrichment analysis on DEGs to determine whether dietary sensor pathways, such as insulin/insulin-like growth factor (IGF) and target of rapamycin (TOR), are affected in the tested liquid cultivations. We have not observed enrichment for the renowned dietary sensor pathways but found other metabolic pathways to be affected. For instance, in CeHR, metabolic process in N2 and lysosome and fatty acid metabolism in AB1 demonstrated overrepresentation (Table 1). CeMM-grown animals have been reported not to show pathologies related to starvation [5]. Similarly, we have not detected differential expression for dietary restriction linked insulin-like signaling, pha-4/FoxA, Table 1 Enriched pathways in F1 generations in response to the liquid cultivations. The green and red arrows show upregulation and downregulation compared to the expression in the corresponding P 0 generations. The P-values are Benjamini-Hochberg corrected.
Particular environmental change-triggered DEGs in F 1 maintained their upregulation or downregulation in the next generation. The transmission of the environment-induced differential gene expressions to the F 2 generation from F 1 was between 20 and 34% (Fig. 4a). To assess whether the transmitted differential gene expressions are orchestrated by the same transcriptional regulators (i.e., transcription factors), we first evaluated the promoter regions of these genes for motif enrichments as potential recognition sites. The S-Medium experiment did not present any motif enrichment while the CeHR did for multiple sequence motifs in both the strains. The enriched motifs from both the upregulated and downregulated genes demonstrated differences between the strains, indicating that distinct transcriptional regulators can potentially be in play. We further compared the motifs from the upregulated genes against known transcription factor recognition sites [32] to find if common transcription factors induce gene expressions in CeHR for the strains. One such motif from both the strains is a recognition site for transcription factor MDL-1 (Additional file 1: Figure S5a). This transcription factor has regulatory functions in the inhibition of germline growth and longevity [33]. Other motifs from upregulated genes did not exhibit enrichment for transcription factor recognition. To test the human biological processes related to the putative MDL-1 recognition motifs from both the strains, we performed GO analysis for the DNA motifs against human promoters [32]. Both the motifs were overrepresented with the sensory perception of smell GO term, presumably corresponding to chemosensation functioning in the detection of environmental changes in worms. However, other enriched terms were different for these two motifs. For instance, the motif from the N2 strain had an overrepresentation of DNA damage checkpoint and transcription initiation from RNA polymerase II promoter while the motif from the AB1 strain was enriched for G-protein coupled receptor protein signaling pathway and positive regulation of immune response (Additional file 1: Figure S5a).
We next sought to identify the transcriptional factor genes expressed at higher levels in the F 1 and F 2 generations. We considered the F 1 -induced transcription factor genes as the putative regulators of the environmental change responses and F 2 -induced ones as the putative regulators of re-adaptive responses. The comparison of our list of DEGs to the putative transcription factor genes [34] revealed that the majority of the differentially expressed transcription factors genes are strain and Fig. 4 Different strains and the environment conditions contribute to differences in the maternal gene expression profiles. a The DEGs in the F 1 , F 2 , and both (F 1 F 2 ) generations were identified in comparison to the gene expressions in the P 0 animals. Between 5 to 34% of the DEGs in the F 1 generations were transmitted to the F 2 generations (see the third category -F 1 F 2 ) and these transmissions are higher than expected by chance (hypergeometric test, P-value < 0.001 for all the upregulated and downregulated genes individually). The gene expressions transmitted to the next generation after exposure to CeHR showed distinct enrichments for sequence motifs in their promoter regions for N2 and AB1 strains. The number of DEGs genes in F 1 , F 2 , and F 1 F 2 were represented with blue, green, and purple, respectively. b Top five tissue enrichment of the DEGs under exposure to environmental changes in F 1 (FDR-adj. P-value < 0.05) condition-specific (Additional file 1: Figure S5b). None of the CeHR reversion-induced transcription factor genes were common between the strains indicating different transcriptional re-adaptive responses to agar plates for the N2 and AB1 strains. We identified only five shared CeHR-enriched transcription factor genes between the strains, and six shared genes between the CeHR and S-Medium conditions in N2 (hypergeometric test; p-value < 0.001). Neuron and muscle cell enriched R06C1.6 was commonly induced in CeHR for both the strains and in S-Medium for N2, suggesting a role for R06C1.6 in swimming behavior.
To explore whether the DEGs among the F 1 generations show similar spatial expressions, we performed tissue enrichment analysis (Fig. 4b). We discarded any enrichments on embryonic cells as we conducted all our experiments on the L4 stage worms. The downregulated genes in CeHR demonstrated enrichment for reproductive system-related tissues such as the oocyte, Z2, and Z3 for both the strains. This possibly reflects the changes in the reproduction of the worms in CeHR. However, upregulated genes in CeHR exhibited overrepresentation on different tissues between the strains: DEGs in N2 were mainly enriched for cylindrical hypodermal syncytium in head; DEGs in AB1 were enriched for neuronal tissues. The genes upregulated in response to S-Medium showed overrepresentation for neuronal tissues, but to different ones than the CeHR-exposed AB1 strain genes.
Notably, downregulated genes in CeHR for both the strains and in S-Medium for N2 are enriched in cephalic sheath cell (Fig. 4b). Altogether, these findings suggest that laboratory and wild-isolate strain worms show significant variations in the transcriptomic responses to the environmental alteration.

Functional profiles of the differentially expressed genes show correlation with the observed phenotype
We wondered if the functional data analysis results indeed point to phenotypical alterations in the animals. To test this, we employed microscopy techniques to detect the potential phenotypes. In the N2 strain, the genes functioning in reproduction and oogenesis have been downregulated in CeHR compared to agar and reversion. In AB1, the DEGs among the conditions did not show the same GO enrichment. This pattern is an indicator of phenotypical changes in the reproductive system that can be observed in CeHR but absent in reversion for N2. To test our prediction, we examined the germlines of N2 and AB1 adults in the different growth conditions via confocal microscopy ( Fig. 5a, b, Additional file 1: Figure S8). CeHR grown N2 but not AB1 adults displayed aberrant germlines. Additionally, fewer oocytes were observed in CeHR-grown N2 strain compared to OP50 NGM-grown N2 strain. These findings suggest that the reduced fecundity in CeHR raised N2 animals may have resulted from the inability of animals to make gametes. In agreement with Fig. 5 N2 but not AB1 C. elegans exhibit germline abnormalities in CeHR. a N2 and AB1 germlines exhibit normal germ cell morphologies and numbers in OP50 NGM (aberrant germlines were observed in 21 and 12.5% of the worms, respectively). Germlines of S-Medium grown animals are similar to OP50 NGM grown animals. However, embryos accumulate in 36% of the S-Medium grown adults. As expected, N2 adults grown in CeHR exhibit aberrant germlines with drastically reduced germ cell numbers, which results in an overall smaller germline (86% of the worms). N2 adults in CeHR have a significantly fewer number of oocytes and a high number of internally hatched larvae. These abnormalities are not observed in CeHR-grown AB1 adults (aberrant germlines were observed in 33% of the worms). b Dissected out gonadal arms of N2 C. elegans grown on OP50 NGM and CeHR. Germline size and cellular morphology are severely deformed in animals grown in CeHR our hypothesis, the worms placed back to agar plates did not exhibit the similar phenotypical abnormalities with CeHR grown animals (Additional file 1: Figure S6).
Then, we asked if these germline abnormalities were caused by an increase in physiological germ cell apoptosis. We utilized a transgenic line of C. elegans that expresses functional CED-1::GFP fusion protein in sheath cells [35]. CED-1 functions by initiating a signaling pathway in phagocytic cells that promotes cell corpse engulfment and is commonly used as a marker to detect cells undergoing the apoptotic pathway. In CeHR grown adults, CED-1::GFP was detected in a much higher proportion of germ cells (Fig. 6a, Additional file 1: Figure S8). Why the number of apoptotic germ cells in CeHR grown N2 adults were increased and whether or not this was the sole cause of aberrant germlines will require further investigations.
For the N2 worms, another GO enrichment term was collagen and cuticulin-based cuticle development for upregulated genes in CeHR compared to agar and reversion generations. Hence, we reasoned that the cuticles might show morphological differences between agar and CeHR-grown N2 animals. In adult animals, cuticular structures were difficult to distinguish between N2 and AB1 via scanning electron microscopy (SEM). More pronounced cuticular structures were detected in the earlier larval stages (L3) with SEM. We observed more protruding alae structures in CeHR-grown N2 larvae compared to AB1 larvae, but the annuli were indistinguishable between the two strains (Fig. 6b). Comparing 42 genes annotated with alae morphology variant from WormBase to our data, we identified 15 of these genes (dpy-2, dpy-3, dpy-5, dpy-7, dpy-8, dpy-9, dpy-10, dpy-18, lon-3, qua-1, rol-6, sqt-1, sqt-2, sqt-3, and mlt-10) upregulated under CeHR conditions in N2 but not in AB1 (Additional file 1: Figures  S7 and S8). Taken together, the correlation between the function of the DEGs and the observed phenotype Fig. 6 CeHR-grown C. elegans exhibit higher number of physiological germ cell apoptosis, more pronounced alae structures, and smaller fat stores. a OP50 NGM, S-Medium and CeHR grown MD701 adults were stain with DAPI. One or two germ cells enclosed by CED-1::GFP are detected around the loop region of the gonads in MD701 adults grown in OP50 NGM and S-Medium (84.6 and 100% of the worms, respectively). CED-1::GFP is detected around a large proportion of germ cells in CeHR grown MD701 adults (90.9% of the worms). b Scanning electron microscope (SEM) micrographs of OP50 NGM and CeHR grown N2 and AB1 L3 animals. N2 and AB1 L3 animals display normal alae and annuli with grown on OP50 NGM (81.3 and 85.7% of the worms, respectively). N2 but not AB1 L3 exhibit more pronounced alae in CeHR (72.2 and 12.5% for N2 and AB1, respectively). c Nile red staining of L3 and adult N2 C. elegans grown on OP50 NGM and CeHR. The number of fat stores is indistinguishable between OP50 NGM and CeHR grown N2 animals, but the fat stores significantly larger in OP50 NGM (81.8% large; 18.2% small) than CeHR raised N2 adults (27.3% large; 72.2% small) around the abdominal and tail regions suggests that these genes can potentially be a part of underlying genetic networks causing these phenotypes.
Previously, we demonstrated that the CeHR-grown animals were thinner and longer than the OP50 NGM raised animals [6]. We asked if this change in body morphology was due to a reduction in fat stores that may have resulted from the absence of a bacterial diet. N2 animals grown on OP50 NGM and CeHR were stained with Nile Red [36] to visualize fat stores in all tissues of L3 and adult animals. We observed no significant difference in the number of fat stores between OP50 NGM and CeHR-grown N2 animals; however, larger fat stores were apparent in the abdominal and tail regions of OP50 NGM grown N2 adults (Fig. 6c, Additional file 1: Figure S8). Thus, the lack of a bacterial diet apparently contributes to the morphological differences through the fat storage.

Discussion
In this paper, we demonstrated that environment changes cause significant alterations in the transcriptional responses in a C. elegans model. Different strains of the worms present unique transcriptional characteristics in the adaptive responses, and these characteristics show correlation with the observed phenotypes. For instance, we detected differential expression in the body morphology-related genes along with the phenotypical changes on the body morphology in CeHR-grown N2 strain, but not the AB1 strain. Furthermore, we found that different environments provoke the condition-specific expression of genes, and some of these genes lacked previous experimental EST or cDNA evidence.
Our results revealed that specific groups of genes are expressed only in particular environments. Some of these genes lacked transcript evidence previously (Fig. 1e). Considerable effort has been devoted to annotating the genes in C. elegans, but the standard WormBase [20] annotation still contains thousands of "predicted genes" (about 8%) which lack direct experimental cDNA or EST evidence [37]. The predicted genes in the C. elegans genome have been mainly identified by computational methods. Experimental detection of these genes was missing; mostly because they demonstrate poor expression, weak statistical signals or less conservation across species [38]. Others found that unused genetic information can remain in the genome for many generations [39]. Reactivation of dormant genes has been used for the treatment of particular disorders indicating that the silenced genes can have essential functions in the genome [40]. These studies provide evidence that dormant genes can be available on the genome throughout many generations and they can contain essential genetic information. Thus, we entertain the possibility that these predicted genes simply need a trigger to activate their expression, and this trigger may be elicited when animals are exposed to certain environmental stimuli.
A better understanding of the biological systems can be achieved through mining publicly available data sources [41]. We expanded our knowledge in the affected genetic systems by acquiring the ncRNA profiles from WormBase [20]. Condition-specific expression of ncRNAs implies important regulatory functions for these molecules in adaptive responses. Earlier studies have identified various classes of ncRNAs and their roles in biological systems including RNA splicing, DNA replication, epigenetic regulation of gene expression, and X-chromosome inactivation [42][43][44][45][46]. However, the majority of predicted ncRNAs have properties or functions that have not been identified yet. The list of unclassified ncRNAs, or "ncRNAs," employed in this study was derived from the '7 k-set' that was generated by the modENCODE Consortium [37]. The '7 k-set' was assembled via predictions based on conservation and RNA secondary structure, and therefore, functional genomic studies of these "ncRNAs" are also still lacking. The co-enrichment of the ncRNA expressions with the coding ones suggest an interplay between these two.
In our RNA-seq experiments, we used poly(A) selection method but still observed numerous ncRNA molecules potentially because many eukaryotic ncRNAs are polyadenylated. For example, a poly(A) tail is part of the mature RNA for many long ncRNAs, i.e., Xist that mediates X-chromosome inactivation [47]. We cannot rule out that the ncRNAs sampled in our study resulted from technical artifacts in the RNA-seq. Nevertheless, many ncRNAs had expression levels higher than the housekeeping gene pmp-3. These findings would seem to suggest that the sampling of these RNAs was not an artifact of RNA-seq. The inclusion of non-polyadenylated ncRNA molecules can bring more insights for understanding the ncRNA functions in adaptive responses.
Along with the phenotypical differences, there were considerable differences in RNA expressions between the N2 and AB1 strains which were cultured in the same growth conditions. N2 had been reported to exhibit innate immune response against pathogenic bacteria [48,49]. However, previous studies have determined that the N2 and wild-isolate CB4856 strains display a high variation in gene expressions and many of these strain-specific variations are related to innate immunity genes [50]. Similarly, our GO analysis results showed enrichment in innate immune response for the genes expressed exclusively in the AB1 strain on OP50 NGM (Fig. 2). These genes do not correspond to the previously reported divergent regions between the strains which might have caused the differences in gene expressions [51]. The AB1 animals are challenged more with pathogens in their natural habitat. The AB1-specific expression of innate immune response genes may indicate that the expression of these genes are more readily available in the AB1 strain than N2 to prepare the animals for a potential pathogen encounter. The results seem to suggest that C. elegans in the wild have the natural aptitude to survive not only on different food sources but also many other environmental variables including pathogens or the shifts between the solid and liquid settings. However, years of laboratory domestication of the N2 strain may have given rise to laboratory selections or genetic bottlenecks that may have hindered the ability of the N2 animals to acclimate to changing environments [14]. Alternatively, a potential presence of cis-and trans-eQTL between N2 and AB1 can be the underlying reason for the differences in the enrichment of innate response genes. Follow-up experiments are needed to determine the underlying reason for the expression differences in the innate immunity genes between N2 and the wild-isolates AB1 and CB4856.
For the majority of the C. elegans studies, the worms are cultured on bacteria-seeded agar plates instead of an axenic medium [13,52]. However, a bacterial food source can present confounding effects on certain studies. For example, bacterial by-products can create genetic responses in the animals, and these responses are often overlooked. Axenic media is desirable in space, biochemical, and toxicology studies as it enables automated culturing and experimentation, and it eliminates the potential contamination risks due to a bacterial diet [6,53,54]. Placing the animals in axenic media, however, have profound consequences altering a variety of biological processes. To distinguish whether the genetic responses are from the case study or the axenic media conditions, a separate investigation should be conducted. In light of this, we revealed that C. elegans acquire large variations in gene expressions upon single generation exposure to two types of liquid media, and these variations are more prominent for the N2 strain.

Conclusions
Culturing the worms in axenic medium seems to affect the worms more than only changing their physical environment to liquid in the S-Medium. Our results highlight that caution must be used while studying model organisms as the constant laboratory environments can cause substantial differences in the transcriptomic and phenotypic responses to abrupt environmental changes. We believe our data can provide standard controls for future studies that utilize liquid cultivation of C. elegans for experimentations and is a rich resource for the discovery of genes showing environment-specific expression.

C. elegans strains and growth conditions
Wild-type N2, wild-isolate AB1, and MD701 bcIs39 [lim-7p::ced-1::GFP + lin-15(+)] strains were obtained from the Caenorhabditis Genetics Center (CGC). All animals were grown at 21°C. Embryos were isolated via the bleaching protocol in each step [52] and placed in the following growth media: E. coli OP50 seeded on NGM agar plates (OP50 NGM), 1 mL of S-Medium in 20 mL scintillation vials inoculated with a concentrated E. coli OP50 pellet made in 6 mL of an overnight culture (OP50 S-Medium), and CeHR medium in 20 mL scintillation vials according to Nass and Hamza [3] with minor modifications (we used 50 μM hemin instead of 20 μM and used 30 μM HEPES instead of 10 μM).

RNA isolation, Illumina sequencing
Total RNA from approximately 100,000 synchronized young adult C. elegans (4 h post L4) were isolated with a modified TRIzol protocol and recovered by alcohol precipitation. Total RNA was further purified by PureLink ™ RNA Mini Kit (Life Technologies). Two micrograms of total RNA were used for library preparations using the TruSeq Stranded mRNA LT Sample Prep Kit (Illumina). Libraries were sequenced on an Illumina HiSeq 2500 instrument set to the rapid run mode using single-end, 1 × 51 cycle sequencing reads as per manufacturer's instructions.

qRT-PCR
Quantitative RT-PCR was conducted using 1st strand cDNA synthesized from total RNA and gene-specific primers (Additional file 1: Table S3). Each cDNA sample was amplified using the SYBR® Premix Ex Taq ™ II (Takara Bio) on the ABI 7500 Fast Real-time PCR System (Applied Biosystem) according to the manufacturers. The experiment was performed by three independent experiments with biological triplicates.

Bioinformatics analyses
The Tuxedo pipeline [55] was used to find DEGs with default parameters. The reference genome (WBCel235) was obtained from Ensembl [56]. Data from the study conducted by Dallaire et al. were retrieved from Gene Expression Omnibus (accession number GSE54173) [21] to compare the expression of the unconfirmed genes in OP50 NGM-grown N2, and the same pipeline with the same parameters were used for the analysis of the data. The same pipeline with same parameters were used to analyze RNA-seq data. We considered genes with FPKM > 1, FDR adjusted p-values < 0.05, and log 2 fold change > 2 as differentially expressed. miRNA, piRNA, and rRNA molecules were discarded from the analysis. The ncRNA molecules (WS250) and unconfirmed genes were acquired from WormBase through the ftp site and personal communications, respectively [20]. The motif enrichment analysis was performed with MEME Suite (v4.11.1) [32]. The promoter sequences for the motif enrichment were retrieved from UCSC Genome Browser database for ce11 [57]. Find Individual Motif Occurrences (FIMO) was used to find the enrichment of the detected sequence motifs for the known transcription factor recognition sites [58]. We utilized GOMO (v4.12.0) software to identify the enriched human GO terms for the promoter sequence motifs [59]. GO and pathway enrichment analyses was made with the Database for Annotation, Visualization and Integrated Discovery (DAVID) (v6.8) [60], and Benjamini Hochberg corrected and FDR adj. p-values < 0.05 considered significant. Tissue and phenotype enrichment analyses were made by using WormBase Gene Set Enrichment Analysis tool [61], FDR adj. p-values < 0.05 considered significant. Long, thin, and alae morphology variant phenotype genes were retrieved from WormBase (WS262). R software was used for all statistical analyses and the plots [62].

Microscopy
Nematodes were fixed in cold (− 20°C) methanol for 15 min and stained with SYBR Gold (1:500,000 dilution in PBS) or DAPI (100 ng/mL) for 15 min, or Nile Red (5 ng/ mL) for 30 min, at 21°C with gentle rocking. Worms were destained in PBST, washed three times with M9 and mounted on a 2% agarose pad for microscopy. All micrographs were taken with the Zeiss 710 laser scanning confocal microscope with a 40×/1.2 N.A. water immersion objective. Dissected gonads were stained with SYBR gold as mentioned above. Nematodes were placed onto 0.22 μM filters where the excess liquid was removed to visualize cuticular structures via cryo-SEM. Upon attaching the filter paper onto the sample holder, samples were plunged into liquid nitrogen slush. There a vacuum was pulled allowing sample transfer to the Gatan Alto 2500 cryo chamber at a temperature of − 125°C. Samples were sublimated for 10 min at 90°C followed by cooling to − 125°C. A thin layer of Gold Palladium was sputtered onto the samples. The samples were then transferred into a Hitachi S-4700 Field-Emission Scanning Electron Microscope for imaging. At least 20 animals were screened for each phenotype, and the represented figures were used for the observed phenotypes.