Robustness and variability inCaenorhabditis elegansdauer gene expression

Both plasticity and robustness are pervasive features of developmental programs. The dauer in Caenorhabditis elegans is an arrested, hypometabolic alternative to the third larval stage of the nematode that undergoes dramatic tissue remodeling and gene expression changes compared to conspecifics that continue development. Dauer arrest can be triggered by several adverse environments or genetic mutations that act as independent and parallel inputs into the larval developmental program and is an example of phenotypic plasticity. However, whether gene expression in dauer larvae induced by different genetic or environmental triggers is invariant or varies depending on their route into dauer has not been examined. Here we use RNA-sequencing to characterize gene expression in dauer larvae induced to arrest development in response to different stimuli. By assessing the variance in the expression of all genes and computing the Spearman's rank-order correlation of gene expression within several Gene Ontologies (GO) and gene networks, we find that the expression patterns of most genes are strongly correlated between the different dauer larvae, suggestive of transcriptional robustness. However, we also find that gene expression in specific defense and metabolic pathways varies widely between dauers. We speculate that the transcriptional robustness of core dauer pathways allows for the buffering of variation in the expression of genes involved in the response to the environment, allowing the different dauers to survive in and exploit different niches. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint

Understanding how genotypes map onto phenotypes remains an exciting and unsolved problem in biology1. Phenotypic plasticity, the process by which a genotype can yield several phenotypes, allows organisms to adapt their development, physiology, metabolism, and even morphology to a varying and unpredictable environment2,3 (Fig 1A). However, robustness, the invariant expression of a trait despite genetic or environmental perturbations, is also essential for survival and fitness and is a pervasive feature of developmental programs 4-6(Fig 1B). Dauer formation in the nematode Caenorhabditis elegans is considered an example of extreme phenotypic plasticity7-14. The C. elegans dauer stage is an alternative developmental stage to the third larval stage of the nematode triggered during late larval stage1 (L1)/early larval stage 2 (L2) by environmental stressors such as starvation, crowding, or extreme temperatures. Mutations that downregulate growth-promoting signals that license development also promote dauer arrest and are thought to mimic the environmental triggers 7-14. Thus, the insulin signaling pathway (ILS) whereby insulins released in the presence of food act through the sole insulin-like receptor, DAF- 2 to antagonize the activation of the Forkhead transcription factor DAF-16/FOXO is required for continuous growth; downregulating DAF-2 signaling, as occurs during food scarcity leads to dauer arrest. Likewise, the TGF-β pathway, where the DAF-7 ligand acts through the TGF-β receptors DAF-1/DAF-4 to antagonize the DAF-3/SMAD-DAF-5/Ski transcription factor complex, is also required to permit continuous growth and daf-7 mutants constitutively arrest as dauers under conditions where wild type larvae continue development into reproductive adults 15-17. Recently a cytokine interleukin IL-17 pathway that inhibits the C. elegans p53 ortholog p53/CEP-1 has also been shown to be necessary for continuous growth18. Dauer entry is accompanied by a developmental arrest, gene expression changes 19, morphological changes such as radial shrinkage, pharyngeal constriction, development of a specialized cuticle and buccal plug, physiological changes such as a precipitous decrease in feeding, behavioral changes such as the favoring of nictitation, and metabolic changes such as reduced activities of glycolytic, gluconeogenic, Tricarboxylic Acid Cycle (TCA) cycle, and oxidative phosphorylation pathways 7-10,20. The different environmental and genetic triggers that induce these profound changes are thought to act largely through independent or parallel mechanisms to initiate dauer entry7,8. This has been well studied in the case of the ILS/DAF-2 and TGF-β/DAF- 7 pathways, which regulate parallel and independent inputs into the dauer decision7-10,20, although they share some pathway components: insulin gene daf-28 is regulated by both ILS/DAF-2 and TGF-β/DAF-7 pathways 21. Similarly, although both loss of DAF-2 and ILC-17.1 act genetically upstream of DAF-16/FOXO for dauer arrest, they appear to operate through different mechanisms since ilc-17.1 mutants also requires active CEP-1/p53 to trigger arrest, but daf-2 mutants do not18. All dauer pathways impinge on the steroid hormone signaling pathway, DAF-12, but DAF-12 activity has complex effects on dauer arrest as well as developmental timing, and certain daf-12 mutations cause constitutive dauer arrest while others prevent dauer formation 7-9,11. Thus, it is unclear whether gene expression in the different dauer larvae induced to enter dauer by different environmental or genetic perturbations varies based on the route into dauer, or is invariant and robust, independent of the genetic or environmental trigger that induce dauer. To address this question, we used next-generation RNA sequencing (RNA-seq) to obtain the transcriptional profiles of different dauer larvae. We used dauers initiated by exposure to high temperatures (N2; High Temperature Induced Dauers or HID)22, loss of the insulin signaling (daf- 2)23, downregulated TGF-β signaling ( daf-7)24, loss of the cytokine ILC-17.1 pathway ( ilc-17.1), .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint and activation of the CEP-1/p53 pathway ( cep-1 o/e)18. We then estimated the variance in gene expression across the different dauer larvae and used Spearman's rank-order correlation to compare gene expression within different Gene Ontologies (GO) and functional pathways to also compute an indicator of similarity between dauers. Our analysis shows that while the expression levels of most genes vary widely between the different dauer larvae, expression patterns within most functional pathways are highly similar. Gene expression in a set of stress and defense response pathways does vary and is not correlated between dauer larvae. Our data imply the presence of robust constraints that stabilize gene expression variation. We speculate that gene expression variability in pathways that modulate adaptation to the environment is buffered by the robustness of core dauer pathways, allowing the different dauer larvae to be better suited to survive in and exploit different environmental niches. RESULTS. Dauers induced by different environmental and physiological stimuli utilize similar processes for dauer arrest. To assess the gene expression profiles of dauer larvae generated under different conditions, we conducted RNA sequencing (RNA-seq) of high temperature-induced wild-type dauers (N2; High Temperature Induced Dauers or wild-type HID ; WT HID), daf-2 (e1370) III, daf-7(e1372) III and ilc-17.1 (syb5296) X dauers, and dauers that result from overexpressing the C. elegans p53-like gene cep-1 (cep-1oe) (Fig. 1C, D). Wild-type HID were generated by allowing wild-type larvae to grow at 27°C for 48 hours following hatching. Dauers from the daf-2 (e1370) III, daf-7(e1372) III and ilc-17.1 (syb5296) X and cep-1 oe backgrounds were generated by allowing embryos to grow at 25°C for 48 hours following hatching. In agreement with what has been shown by many previous studies, all dauer larvae displayed the characteristic morphological features of dauer (Fig 1D). In addition, to confirm that the duration and conditions used generated ‘true’ dauers as defined by their resistance to 1% SDS, and to avoid any unknown effects of detergent treatment on our RNA-seq analysis, we conducted a separate pilot experiment and verified that the larvae generated in response to these conditions were SDS-resistant 8. The <1% larvae that escaped dauer arrest were visible as stage 4 larvae (L4) at the time of harvesting and were manually removed prior to mRNA extraction. Thus, the mRNA was highly enriched for dauer-specific genes expressed upon arrest, 48 hours post-hatching. Three independent biological samples were sequenced. In addition, to identify genes that were differentially regulated during dauer arrest, we sequenced late L2/early L3 wild-type larvae that were grown for 32 hours post-hatching at 25°C. These larvae are considered to be at a comparable developmental stage to dauer but do not arrest as dauers and instead are fated to develop into reproductive adults. Their growth at 25°C allowed us to account for the effects of temperature on development. A sample distance matrix produced by comparing expression levels across all genes between samples demonstrated excellent agreement among biological replicates (Fig. 1E; Supplementary Table 1). The mean expression levels (log10 TPM) of all genes were comparable but were modestly lower in the daf-7 and ilc-17.1 larvae perhaps suggestive of global differences in transcription (Supplementary Fig. 1A). Principal Component Analysis (PCA) showed that all dauers separated well from continuously growing larvae, and clustered according to biological replicates (Fig. 1F). Importantly, the different dauer larvae also separated by the trigger that induced dauer entry (genotype or environment), suggesting that they differed from each other, with the separation of wild-type HID from other dauers along PC1 and PC2 being larger than the more modest separation between the daf-2, daf-7, ilc-17.1 or cep-1oe dauers (Fig. 1F). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint As a first step towards evaluating similarities and differences between the dauer larvae, we identified the differentially expressed genes (DEGs) in each of the dauers by comparing their gene expression to that in continuously developing L2/L3 larvae and conducted Gene Ontology (GO) enrichment analysis 25-29. The different dauers displayed similar numbers of DEGs (N2 HID =16302, daf-2 dauers=16885, daf-7 dauers=16557, ilc-17.1 dauers=16747 and cep-1oe dauers= 16565, padj< 0.05), of which 13120 genes were common (Fig 2A; Supplementary Table 2). The DEGs in the different dauers largely represented similar processes, arguably those necessary for their developmental arrest (padj< 0.05; Fig. 2B, C; Supplementary Figs 2 and 3). Thus, all dauer larvae downregulated genes that contributed to the regulation of the cell cycle, DNA metabolic processes, translation, and other anabolic processes. Genes contributing to neuropeptide signaling and ion homeostasis were upregulated in all dauers. In addition to these similarities there were some differences: some GO pathways were enriched in some dauers but not others (Fig. 2B-C; Supplementary Figs 2 and 3; Supplementary Tables 4 and 5), and a small group of between 27 to 109 genes were uniquely altered in each of the dauer larvae (Fig 2A; Fig 2D; Supplementary Table 3). For instance, daf-2 dauers differed from other dauers in the expression of genes related to striated muscle-dense body and contractile fiber, and uniquely altered the expression of genes enriched in “molecular transducer activity GO:0060089” (srb-6, str-220, srg- 2 etc.; Fig 2D; Supplementary Table 6). Similarly, DEGs in larvae that were induced to arrest as dauers due to the loss of the interleukin cytokine gene ilc-17.1 were not enriched for several immune response categories that were enriched in all other dauer larvae, but had an over 20-fold enrichment in GO categories ‘ciliary plasm GO:0097014’ (downregulation of klp-11 and upregulation of H13N06.7) and ‘sodium channel activity GO:0005272’ ; Fig. 2D). Likewise, DEGs in daf-7 dauers also differed in enrichment for a smaller subset of “innate immune response” genes. The GO categories ‘postsynaptic membrane GO:0045211’, ‘regulation of postsynaptic membrane potential GO:0060078l’ and “stabilization of membrane potential GO:0030322” were unique to DEGs in the daf-7 dauers, consistent with what is known regarding DAF-7 expression in neurons, and its effects on the C. elegans nervous system and behavior (Fig. 2D). These data suggested that the different dauer larvae likely utilized similar processes for dauer arrest and other physiological functions, but also displayed differences in their gene expression profiles. Different dauers display high variability in gene expression levels but strong correlation in gene expression patterns. To further assess the extent to which gene expression in the five dauer larvae was similar or differed, we (i) estimated the variance in expression of all expressed genes (mean expression>10), using the coefficient of variation (CV) as an estimate of variance 30-32, and (ii) computed the Spearman's rank-order correlation, rho, between pairs of dauer larvae as an indicator of robustness in their patterns of gene expression33. The coefficient of variation (CV) was computed for each gene by dividing the standard deviation (SD) of its expression across all dauers by its mean expression. The CV across all dauers ranged between 387% and 11%. Somewhat arbitrarily, but based on the genome-wide CV distribution, we defined CV < 30% as low, <30% < CV 50% as high, the latter cut off based on considering that a CV >50% would represent genes whose SD values were over half their mean values. Strikingly, over half the genes (8538/16940, 50.4%) exhibited CVs >50% (Fig 3A; Supplementary Table 7). The relatively high CVs were not simply a consequence of low gene expression, as seen by plotting the SD as a function of the mean and inferred because the CVs were computed after filtering for low mean expression (Supplementary Fig 1B). This suggested that notwithstanding the similarities in the processes upregulated or downregulated by all dauer larvae, the expression .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint levels of genes contributing to these processes varied, and in some cases varied widely, across the different dauers. Surprisingly, however, a pair-wise Spearman's rank-order correlation between individual dauers showed the opposite trend, and gene expression in the different dauers was strongly correlated (Fig 3B). The correlation coefficient, rho, calculated between the log10 gene expression values of expressed transcripts in the five dauer larvae ranged between >0.7, to >0.9, and was significant [we defined rho=0.7-1.0 as strong correlation, rho=0.5-0.7 as moderate, and rho<0.5-0.6 as weak correlation, padj<0.05, according to 1,12]. Thus, the correlation coefficient between the N2 HID dauers and daf-2, daf-7, ilc-17.1 or cep-1oe dauer larvae was rho(16085)=0.81, rho(16085)=0.799, rho(16085)=0.78, and rho(16085)=0.758 respectively (p=2.2e-16 in all cases). Surprisingly, gene expression in daf-2 and daf-7 dauers was also similar [rho (16085) =0.917], even though daf-2 and daf-7 mutations activated separate, independent mechanisms to trigger dauer entry. The more recently discovered dauer pathway triggered by the loss of ilc-17.1, which decreases glucose uptake, resembled daf-2 dauers, which downregulate insulin signaling [rho (16085) =0.948] but differed more from daf-7 dauers. Dauers induced by cep-1 oe most resembled daf-7 dauers and most differed from ilc-17.1 dauers (rho (16085) =0.951 and rho (16085) =0.759 respectively). Thus, together, these data showed that although gene expression varied across the different dauers with over 50% of the genes displaying CVs of >50% there was a strong correlation between all dauers. Notwithstanding the differing levels of resolution provided by CVs and Spearman’s correlation coefficient, these analyses are suggestive of the presence of constraints that stabilize gene expression variation. Gene expression in the core dauer pathway is robust. Previous studies have shown that variation in gene expression is a function of several fundamental biological processes, including the gene regulatory networks and other contexts within which genes operate 1,2,5,6,30,31,33,34. Therefore, to understand the genome-wide CVs, we ranked individual genes according to the CVs computed across all dauers and examined whether genes that displayed relatively low (CV < 30), moderate (30 < CV 50), CVs, were enriched in distinct GO categories that may yield insights into their expression constraints. Indeed, the 10449 most variable (CV>50) genes were enriched in immune and defense pathways. Genes with relatively ‘moderate’ CVs were enriched in pathways related to DNA damage, and the 1490 ‘less’ variable genes were enriched in 15 terms all associated with RNA-related processes (Fig 4A; Supplementary Table 8). These results are consistent with what is known about the enrichment of housekeeping functions amongst genes with low expression variability 6,30,32-34 and suggested that the CVs of genes in the different dauer larvae may indeed be related to the biological pathways within which they functioned. Based on these results, and to better understand the constraints on gene expression variance in the different dauer larvae, we examined the behavior of genes that functioned in the ‘core’ dauer pathways. As ‘core dauer pathway genes’, we chose genes that acted in the ILS (insulins and DAF-16 targets), TGF-β and DAF-12/steroid hormone pathways that are pivotal to trigger the dauer decision, genes that implement the growth arrest, namely those involved in cell cycle regulation (cyclins, the Anaphase Promoting Complex, cyclin-dependent kinases, cyclin inhibitors, etc.), and energy metabolism, and six cuticulin genes that are required for the formation of dauer specific alae8,9. Because daf-2 and daf-7 dauers harbor mutations in the sole insulin receptor and the DAF-7/TGF-β ligand respectively, whereas the other dauers, i.e. wild type HID, ilc-17.1 and cep-1 oe dauers, do not, one prediction was that the CVs of ILS and TGF-β pathway-genes would .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint be relatively ‘high’ when computed across the different dauers, and the pathways would show a weak correlation between daf-2/daf-7 dauers and the others. However, because FOXO/DAF-16 and DAF-3/SMAD-DAF-5/Ski transcription factors responsible for the dauer arrest of daf-2 and daf-7 are also activated in ilcs-17.1 and cep-1 oe larvae during dauer formation and required for their dauer arrest, an alternate possibility was that the activation of these determinative transcription factors, which occurred upon dauer entry through any route, normalized or stabilized gene expression in these pathways. This could result in a strong correlation and, perhaps, low CVs for these genes. The CVs of the ILS, DAF-16, and DAF-12 pathway genes largely supported the former hypothesis, and gene expression in these pathways was more variable than in the cell cycle or energy metabolism pathways. Thus, the majority (45/67, 67%) of cell cycle, glycolysis (8/13, 61.5%), and ETC (49/88, 55.68%) genes had CVs that ranged between moderate to low (20%- 50%), and in the lower half of the genome-wide CV range while the majority of insulins (28/31, 90.32%), DAF-16 targets (56/78, 71.79%) and DAF-12 targets (16/29, 55.17%) were more variable and had relatively high CVs of >50% (Fig 4B; Supplementary Table 9). Surprisingly, the CVs of the TGF-β genes (20/26, 76.9%) were low to moderate and resembled the cell cycle genes (Fig 4B), even though daf-7 dauers had a mutation in the DAF-7/ TGF-β ligand and downregulated TGF-β signaling, suggesting that more needed to be understood regarding the regulation of gene expression in these pathways. A comparison of the mean mRNA expression values (rlog counts), and the log2 fold changes of genes in the core dauer pathway reinforced this interpretation (Supplementary Fig 4A-D). In addition, in the ILS, TGF-β and DAF-12 pathways mRNA expression levels in the different dauers varied in complex ways, across specific genes: for instance, daf-7 mRNA levels were upregulated in all dauer larvae, including in the daf-7 (e1372) dauers, which harbor a point mutation in the daf-7 coding sequence (Supplementary Fig 4C). Expression levels of the six cuticulin genes was more variable (Fig 4B). Despite the high CVs and variability evident in the expression levels of the genes functioning in the ILS, DAF-16, and DAF-12 pathways, the pair-wise correlation coefficients computed between all dauer larvae for the core dauer pathways were uniformly significant and the vast majority of cases, strong [rho >0.7 in all cases; padj<0.05 (Fig 4C; Supplementary Table 10)]. This was true when the expression of individual genes themselves were less variable and had low to moderate CVs, as in the case of the cell cycle and energy metabolism genes, but also when the individual gene expression was variable and displayed high CVs, as in the case of the ILS, DAF-16, and DAF-12 pathways. The few exceptions were in the TGF-β pathway gene expression compared between ilc-17.1 and cep-1 oe dauers, which was significant but moderately correlated (rho=0.59; padj <0.05), and glycolysis genes that differed and were not significantly between cep-1oe dauers and daf-2 or ilc-17.1 (Fig 4C). These data are suggestive of different mechanisms of control over genes that function in the core dauer pathway. Thus, the cell cycle and energy metabolism genes that can be expected to implement the dauer arrest appear to be inherently less variable and their expression levels are more similar across the different dauer larvae. The expression levels of gene in the ILS/TGF- β/DAF-12 pathways that trigger the dauer decision are more variable; however, despite this variation their expression patterns appear scaled, as evidenced by their strong correlation, and they are able, presumably, to reliably enforce a less variable expression of genes that implement dauer arrest. Thus, we posit that together these data show that additional layers of control over gene expression variation renders the core dauer pathway transcriptionally robust. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint Gene expression patterns enriched for a limited set of traits, stress response, immune response and metabolic pathways vary between the different dauer larvae. Our analysis showed that variability in gene expression in the different dauers may be related to specific biological functions. Therefore, to identify pathways that might distinguish the daf-2, daf- 7, lc-17.1 and cep-1 oe dauers, we took advantage of GO enrichment analysis of the DEGs in the different dauer larvae and examined the correlation of gene expression within these GO categories. As shown previously (Fig 2A) a core set of 13120 genes were differentially expressed in all dauers when compared to continuously growing L2/L3 larvae. These shared 13120 DEGs were enriched in 309 GO categories ( Fig 5A; Supplementary Table 11 and 12), and their fold changes in expression compared to the L2/L3 larvae showed modest but distinct differences between dauer larvae (log2 fold change ; Supplementary Fig.1C). Nevertheless, consistent with the genome-wide high correlations between the dauers, gene expression within most of these 309 GO categories was also strongly correlation between all the dauers, with the exception of 21 (Spearman Correlation heat map Fig. 5B ; purple indicated GO categories with low or not significant correlation between any two dauer larvae; rho 0.05). These weakly correlated GO categories were enriched in genes that acted in several different processes such as “response to odorants”, “microtubule depolymerization”, “mRNA destabilization”, etc., and varied between specific dauer larvae, and to different extents (Fig 5C). For instance, daf-2 and ilc-17.1 dauers differed from wild-type HID and cep-1 oe dauers in the expression of genes involved in the “response to odorants”. This GO category (Supplementary Tables 12 and 13) includes hsf-1, tph-1 and nine paralogs of ODR-2 called hot genes (for “homologs of odr two”) that have not been widely studied. ilc-17.1 differed from all other dauers except daf-2, and cep-1 oe dauers differed from all other dauers besides daf-7 in the expression patterns of 15 genes that regulate “microtubule depolymerization”, and cep-1 oe , dauers differed from all other dauers except daf-7 in “nuclear-transcribed mRNA catabolic process” and “mRNA destabilization” (Fig. 5C; Supplementary Tables 12 and 13). Intriguingly cep-1 oe dauers also differed from ilc-17.1 and daf-2 dauers in their expression of genes regulating dendrite development in GO categories “dendrite morphogenesis” and “regulation of dendrite morphogenesis”, and ilc-17.1 dauers differed from wild type HID and cep-1 oe dauers in several aspects of male morphogenesis such as “nematode male tail mating organ morphogenesis”, “male genitalia development”, and “male genitalia morphogenesis” (Fig. 5C; Supplementary Tables 12 and 13). These results suggested that despite the robustness in overall gene expression, the different dauer larvae that arrested development in response to different triggers or stimuli, may also differ in specific traits. We next examined whether the different dauer larvae differed in the expression of genes that are co-expressed to constitute the thirty-four gene expression ‘mountains’, as defined by Kim et al (2001)35. These gene expression maps or ‘mountains’ were defined by analyzing the co- expression of 5361 C. elegans genes across a large number of individual experiments that used wild type and mutant animals at different developmental stages, and fall into categories such as “protein expression,” “histone,” “mitochondria,” “germline” “G protein receptors,” “heat shock” etc., likely representing genes and that act together in a functional and/or spatially coordinated manner. As with the previous analysis, gene expression of all the Kim ‘mountains’ with only six exceptions were also highly correlated between all the dauers ( rho values were > 0.7 for all ‘mountains’ except those purple areas in the heatmap; Fig. 6A, B; Supplementary Tables 14 and 15). The six exceptions which were “heat shock”, “cytochrome p450 (CYP)”, “retinoblastoma complex”, “mechanosensation, “intestine” and “amine oxidases” varied either between all dauers, or between specific dauer larvae (Fig. 6 A, B ; Supplementary Fig. 6 and 7 ; Supplementary Tables .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint 14 and 15). Thus, the expression of the amine oxidases, which in C. elegans consist of a family of seven genes, did not correlate between any two dauer larvae (padj=0.1-1 for all comparisons; Fig. 6A, B). Amine oxidases are heterogenous enzymes that catalyze the oxidative de-amination of polyamines and play essential roles in several biological processes in all organisms 37-39. In C. elegans the putative amine oxidases, amx-2 is thought to be the essential monoamine oxidase in serotonin and dopamine degradation, and is important for lifespan regulation, and modulation of the RAS/MAPK pathway activity during C. elegans vulval development 40. The putative monoamine oxidase gdh-1 is thought to be a mitochondrial gene that plays a role in the innate immune response and susceptibility of C. elegans to Pseudomonas aeruginosa PAO1 pathogenesis41, and the putative monoamine oxidase hpo-15 (hypersensitive to POre-forming toxin 15) 42, has also been implicated in a role in defense response to bacterial pathogens. Consistent with the role of overarching role of the p53 family of genes in the control of the cell cycle in all animals, cep-1 oe dauers differed from all other dauers in the expression pattern of their retinoblastoma complex genes (Fig. 6 A, B ; Supplementary Fig. 6B). Intriguingly, amongst the retinoblastoma complex, the log2 gene expression changes for mdt-22 varied to a larger extent in daf-7 and cep-1 oe larvae. mdt-22 is one of several mediator complex subunits which collaborates with CKI-1, a member of the p27 family of cyclin-dependent kinase inhibitors (CKIs) which themselves are a target of CEP-1/p53, to maintain cell cycle quiescence of vulval precursor cells during larval development43. CYP gene expression varied between wild-type HID dauers and ilc-17.1 dauers, that in the ‘heat shock’ mountain, which consisted mainly of the molecular chaperones, varied between cep-1 of dauers and daf-2 dauers (Fig 6 A, B ; Supplementary Fig. 7), cep-1 of dauers and daf-2 dauers also differed from each other in their expression of specific intestine-related genes and in the genes that affect mechanosensation. Finally, we examined correlation in the expression of 1,799 genes within 85 metabolic pathways described in the iCEL1314 metabolic network model of Nanda et al. (2023) 43, where the genes also clustered based on their co-expression in metabolic pathways. In contrast to the overarching high correlation seen across the GO pathways and the Kim mountains, where only a limited number of pathways differed between the different dauer larvae, gene expression in over a third of the iCEL1314 metabolic pathways (thirty-four of the 85 pathways) showed variability between two or more dauer larvae (we discarded three pathways since they only contained two genes each; Fig 6C, Fig 7A ; Supplementary Tables 16 and 17). Thus, gene expression in the glycine cleavage pathway, histidine degradation, mevalonate metabolism, and propionate degradation pathways differed between all dauer larvae, and did not correlate significantly between practically any pair of dauer larvae [padj.>>0.05; (Fig 7A; Supplementary Figs 8-11). Wild type HID differed from all other dauers in the expression of genes involved in ROS metabolism, folate cycle, glyoxylate and dicarboxylate metabolism, and taurine and hypoxanthine metabolism (Fig 7A ; Supplementary Figs 8-11). Iron metabolism genes, too, differed between wild type HID and the other dauer larvae, but gene expression was significantly and strongly correlated amongst these other dauer larvae. In other pathways too, such molybdenum cofactor biosynthesis pathway and PUFA biosynthesis, the HID differed from all dauer except daf-2 dauers, and in chitin degradation and galactose metabolism, they differed from all other dauer except cep-1 oe dauers, with whom they showed a significant and high correlation of gene expression pattern (Fig 7A; Supplementary Figs 8-11). Similarly, daf-2 and ilc-17.1 dauers differed from all other dauers in the expression of genes in the methionine salvage pathway, while the remaining dauers, wild type HID, cep-1 oe and daf-7 showed a strong and significant correlation in the expression of these genes [rho=0.7-0.8, .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint padj.<0.5; (Fig 7A; Supplementary Tables 16 and 17). The expression of genes involved in fatty acid degradation, Complex II activity and glycine cleavage system were also different in daf-2 and ilc-17.1 dauers compared to all other dauers. There were also differences in gene expression patterns between only specific pairs of dauer larvae. For instance, daf-2 dauers differed from wild-type HID dauers in their expression of genes involved in leucine degradation [padj=0.6, and not significant] and also differed from cep-1 oe dauer larvae in gene expression related to leucine degradation and molybdenum cofactor biosynthesis [padj=0.3 and 0.2, respectively, and not significant]. Wild-type HID dauers and ilc-17.1 dauers differed in the expression of genes involved in the chitin degradation pathway [rho=0.4, padj<0.05], leucine degradation pathway [rho=0.6, padj=0.1 and not significant], molybdenum cofactor biosynthesis [padj=0.1 and not significant], and starch and sucrose metabolism pathways [padj=0.1 and not significant]. The latter difference supported previous observations that ilc-17.1 larvae might be deficient in their ability to take up glucose from their dietary source. Dauers induced by the loss of daf-7 function showed the least differences; they too differed from wild-type HID dauers in the expression of genes in the molybdenum cofactor biosynthesis pathway [padj=0.1 and not significant], but otherwise largely shared gene expression patterns with more than one other type of dauer larva. Gene expression in the Met/SAM cycle differed between cep-1 oe on the one hand and ilc-17.1 and daf-2 on the other (Fig 7A; Supplementary Tables 16 and 17). We summarize all the differences between the dauers (Figure 8). This pervasive variability in the metabolic pathways was also confirmed upon comparing the correlation coefficients for gene expression in all metabolic pathways versus that in all the Kim mountains, where except between wild-type HID dauers and cep-1 oe dauers, the metabolic pathway genes were significantly lesser correlated between all dauers (Fisher r-to-z transformation to calculate the significance of difference between rhos; z=-10.78 to 3.98 in favor of the Kim mountains, padj. <0.05; Supplementary Table 18).

In summary, we show that the expression of individual genes in different C. elegans dauer larvae varies, as indicated by the fact that most genes have >50% CVs. Nevertheless, gene expression patterns between the different dauer larvae are highly correlated, suggestive of transcriptional robustness. We speculate that during dauer entry, as seen during other developmental programs, there are gene regulatory mechanisms that stabilize gene expression variation to generate functional outcomes, which in this case is the developmental arrest of larvae in a hypometabolic dauer state. Intriguingly, our analysis suggests that in general, there appear to be three broad categories of dauer larvae: wild-type HID which differed more from other dauers, dauers generated by the lack of ILS or cytokine ILC-17.1 signaling, both of which regulate the normal response to the bacterial nutrients by controlling insulin signaling, glucose assimilation and, as seen in these studies, innate immunity, and the dauers generated by deficiency in the TGF-β pathway or overactive p53/CEP-1, both of which control cell growth. The presence of these three main types is consistent with the known triggers of dauer arrest 7,8,11,14. We find that despite the presence of what may be three ‘types’ of dauer larvae, gene expression in the core dauer pathway involved in the dauer decision and arrest, i.e. the ILS (insulins and DAF-16 targets), TGF-β, DAF- 12/steroid hormone pathways, cell cycle regulation, glycolysis and ETCs pathways is robust and strongly correlated between dauer larvae irrespective of the stimulus that triggered dauer entry. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint Expression of cuticulins varies both at the level of expression as well as in its correlation; however, whether the alae structure itself varies between dauers was not examined. In addition to robustness in core dauer pathways, a limited number of gene expression pathways that regulate specific morphological features, and specific stress, immune, detoxification and metabolic programs vary between the different dauer larvae. While the rationale for the observed variation, and whether and how it related to the route of dauer entry, is unclear, the existence of these variable pathways suggests that the different dauers possess different suites of traits that may allow for better adaptation to specific niches. For instance, the predicted influence of route of dauer entry on dauer morphology or anatomy is intriguing given that dauer development induces dramatic remodeling of neuron morphology and dendrite arborization 44-47 particularly in the IL2 neurons which are responsible for the dauer-specific behaviors, and the known sex-bias in the ability of larvae to survive dauer 48,49. Similarly, the different dauers varied widely in the expression of genes in metabolic pathways. The reason why metabolic enzymes are more variable in some dauer larvae is not clear. Some of these enzymatic pathways may be inherently more variable, and thus may also vary between dauers. For instance, propionate metabolism has undergone changes in the nematode clade, varying between C. elegans and wild strains, and being lost in several parasitic helminth species 50-53. In addition, since dauers do not feed, they must metabolize stored macronutrients, mainly lipids, but also glycogen and trehalose, and the variability in specific metabolic pathways might be related to the food intake of these larvae prior to dauer arrest. Thus, methionine salvage pathway, leucine degradation and histidine degradation are related to the generation of the essential amino acids histidine, methionine, and leucine, which in C. elegans, are typically obtained from the microbial diet 50-53. Similarly, moco enzymes that are central to the Molybdenum cofactor (Moco) biosynthesis pathway can be obtained from the organism’s microbial diet, or can be synthesized by the animal, and thus food intake, and molybdate biosynthesis by the bacteria prior to dauer arrest might influence the expression patterns of Molybdenum cofactor (Moco) biosynthesis pathway genes 54,55. However, it is also possible that these metabolic pathways act downstream of the sensors that must exist in all dauer larvae to continually monitor the environment and control their entry or exist form dauer arrest. Indeed, of the pathways that varied, such a role has been reported for peroxisomal fatty acid degradation, where deficiency in peroxisomal fatty acid β-oxidation in ASK neurons leads to the premature interruption of the dauer arrest and promotes and untimely exit from dauer to promote continued development even under dauer-inducing conditions such as increased levels of ascarosides 56,57. Such a role can be envisioned for metabolic pathways that are modulated by microbial metabolites, whose presence could signify the advent of favorable conditions. Nematodes fill all trophic levels in the food web 58 and are the most abundant metazoan species on the planet 59-61 . Yet, the mechanisms that facilitate their adaptation to diverse, and often hostile niches remains poorly understood. Almost all nematode species have evolved forms of hypobiosis to adapt their life cycles to variable, unpredictable and harsh environmental conditions 58-65. , of which a specialized form is the developmental diapause ‘dauer’. Indeed, it has been hypothesized that among parasitic nematode species, the dauer stage may have been a prerequisite for the evolution of the diversity of their parasitic lifestyles58-65. Thus, it is tempting to speculate that the transcriptional robustness of core dauer pathways allows for the buffering of variation in the expression of genes involved in their response to the environment. Such differences could be pivotal in allowing the different dauers to be better suited to survive in and colonize diverse niches. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint

1 Wagner, G. P. & Zhang, J. The pleiotropic structure of the genotype-phenotype map: the evolvability of complex organisms. Nat Rev Genet 12, 204-213 (2011). https://doi.org:10.1038/nrg2949 2 Sommer, R. J. Phenotypic Plasticity: From Theory and Genetics to Current and Future Challenges. Genetics 215, 1-13 (2020). https://doi.org:10.1534/genetics.120.303163 3 Felix, M. A. & Barkoulas, M. Pervasive robustness in biological systems. Nat Rev Genet 16, 483-496 (2015). https://doi.org:10.1038/nrg3949 4 Braendle, C. & Felix, M. A. The other side of phenotypic plasticity: a developmental system that generates an invariant phenotype despite environmental variation. J Biosci 34, 543-551 (2009). https://doi.org:10.1007/s12038-009-0073-8 5 Wagner, A. Robustness and evolvability: a paradox resolved. Proc Biol Sci 275, 91-100 (2008). https://doi.org:10.1098/rspb.2007.1137 6 Hill, M. S., Vande Zande, P. & Wittkopp, P. J. Molecular and evolutionary processes generating variation in gene expression. Nat Rev Genet 22, 203-215 (2021). https://doi.org:10.1038/s41576-020-00304-w 7 Fielenbach, N. & Antebi, A. C. elegans dauer formation and the molecular basis of plasticity. Genes Dev 22, 2149-2165 (2008). https://doi.org:10.1101/gad.1701508 8 Hu, P. J. Dauer. WormBook, 1-19 (2007). https://doi.org:10.1895/wormbook.1.144.1 9 Karp, X. Working with dauer larvae. WormBook 2018, 1-19 (2018). https://doi.org:10.1895/wormbook.1.180.1 10 Kiontke, K. & Sudhaus, W. Ecology of Caenorhabditis species. WormBook, 1-14 (2006). https://doi.org:10.1895/wormbook.1.37.1 11 Gems, D. et al. Two pleiotropic classes of daf-2 mutation affect larval arrest, adult behavior, reproduction and longevity in Caenorhabditis elegans. Genetics 150, 129-155 (1998). https://doi.org:10.1093/genetics/150.1.129 12 Golden, J. W. & Riddle, D. L. A pheromone influences larval development in the nematode Caenorhabditis elegans. Science 218, 578-580 (1982). https://doi.org:10.1126/science.6896933 13 Golden, J. W. & Riddle, D. L. The Caenorhabditis elegans dauer larva: developmental effects of pheromone, food, and temperature. Dev Biol 102, 368-378 (1984). https://doi.org:10.1016/0012-1606(84)90201-x 14 Golden, J. W. & Riddle, D. L. A pheromone-induced developmental switch in Caenorhabditis elegans: Temperature-sensitive mutants reveal a wild-type temperature- dependent process. Proc Natl Acad Sci U S A 81, 819-823 (1984). https://doi.org:10.1073/pnas.81.3.819 15 Gunther, C. V., Georgi, L. L. & Riddle, D. L. A Caenorhabditis elegans type I TGF beta receptor can function in the absence of type II kinase to promote larval development. Development 127, 3337-3347 (2000). https://doi.org:10.1242/dev.127.15.3337 16 Riddle, D. L. & Albert, P. S. in C. elegans II (eds D. L. Riddle, T. Blumenthal, B. J. Meyer, & J. R. Priess) (1997). 17 Wolkow, C. A., Munoz, M. J., Riddle, D. L. & Ruvkun, G. Insulin receptor substrate and p55 orthologous adaptor proteins function in the Caenorhabditis elegans daf-2/insulin- like signaling pathway. J Biol Chem 277, 49591-49597 (2002). https://doi.org:10.1074/jbc.M207866200 18 Godthi, A. et al. Neuronal IL-17 controls Caenorhabditis elegans developmental diapause through CEP-1/p53. Proc Natl Acad Sci U S A 121, e2315248121 (2024). https://doi.org:10.1073/pnas.2315248121 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint 19 Wang, J. & Kim, S. K. Global analysis of dauer gene expression in Caenorhabditis elegans. Development 130, 1621-1634 (2003). https://doi.org:10.1242/dev.00363 20 Baugh, L. R. & Hu, P. J. Starvation Responses Throughout the Caenorhabditiselegans Life Cycle. Genetics 216, 837-878 (2020). https://doi.org:10.1534/genetics.120.303565 21 Shaw, W. M., Luo, S., Landis, J., Ashraf, J. & Murphy, C. T. The C. elegans TGF-beta Dauer pathway regulates longevity via insulin signaling. Curr Biol 17, 1635-1645 (2007). https://doi.org:10.1016/j.cub.2007.08.058 22 Ailion, M. & Thomas, J. H. Dauer formation induced by high temperatures in Caenorhabditis elegans. Genetics 156, 1047-1067 (2000). https://doi.org:10.1093/genetics/156.3.1047 23 Vowels, J. J. & Thomas, J. H. Genetic analysis of chemosensory control of dauer formation in Caenorhabditis elegans. Genetics 130, 105-123 (1992). https://doi.org:10.1093/genetics/130.1.105 24 Murakami, M., Koga, M. & Ohshima, Y. DAF-7/TGF-beta expression required for the normal larval development in C. elegans is controlled by a presumed guanylyl cyclase DAF-11. Mech Dev 109, 27-35 (2001). https://doi.org:10.1016/s0925-4773(01)00507-x 25 Ashburner, M. et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet 25, 25-29 (2000). https://doi.org:10.1038/75556 26 Gene Ontology, C. et al. The Gene Ontology knowledgebase in 2023. Genetics 224 (2023). https://doi.org:10.1093/genetics/iyad031 27 Harris, T. W. et al. WormBase: a cross-species database for comparative genomics. Nucleic Acids Res 31, 133-137 (2003). https://doi.org:10.1093/nar/gkg053 28 Harris, T. W. & Stein, L. D. WormBase: methods for data mining and comparative genomics. Methods Mol Biol 351, 31-50 (2006). https://doi.org:10.1385/1-59745-151- 7:31 29 Lee, R. Y . N. et al. WormBase 2017: molting into a new stage. Nucleic Acids Res 46, D869-D874 (2018). https://doi.org:10.1093/nar/gkx998 30 Wolf, S. et al. Characterizing the landscape of gene expression variance in humans. PLoS Genet 19, e1010833 (2023). https://doi.org:10.1371/journal.pgen.1010833 31 de Jong, T. V., Moshkin, Y. M. & Guryev, V. Gene expression variability: the other dimension in transcriptome analysis. Physiol Genomics 51, 145-158 (2019). https://doi.org:10.1152/physiolgenomics.00128.2018 32 Foreman, R. & Wollman, R. Mammalian gene expression variability is explained by underlying cell state. Mol Syst Biol 16, e9146 (2020). https://doi.org:10.15252/msb.20199146 33 Naqvi, S. et al. Precise modulation of transcription factor levels identifies features underlying dosage sensitivity. Nat Genet 55, 841-851 (2023). https://doi.org:10.1038/s41588-023-01366-2 34 Espinosa-Soto, C., Martin, O. C. & Wagner, A. Phenotypic robustness can increase phenotypic variability after nongenetic perturbations in gene regulatory circuits. J Evol Biol 24, 1284-1297 (2011). https://doi.org:10.1111/j.1420-9101.2011.02261.x 35 Kim, S. K. et al. A gene expression map for Caenorhabditis elegans. Science 293, 2087- 2092 (2001). https://doi.org:10.1126/science.1061603 36 Nanda, S. et al. Systems-level transcriptional regulation of Caenorhabditis elegans metabolism. Mol Syst Biol 19, e11443 (2023). https://doi.org:10.15252/msb.202211443 37 Katane, M., Seida, Y., Sekine, M., Furuchi, T. & Homma, H. Caenorhabditis elegans has two genes encoding functional d-aspartate oxidases. FEBS J 274, 137-149 (2007). https://doi.org:10.1111/j.1742-4658.2006.05571.x 38 Wang, C. et al. Identification of amino acids important for the catalytic activity of the collagen glucosyltransferase associated with the multifunctional lysyl hydroxylase 3 (LH3). J Biol Chem 277, 18568-18573 (2002). https://doi.org:10.1074/jbc.M201389200 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint 39 Leung, M. C. et al. Caenorhabditis elegans: an emerging model in biomedical and environmental toxicology. Toxicol Sci 106, 5-28 (2008). https://doi.org:10.1093/toxsci/kfn121 40 Schmid, T. et al. Systemic Regulation of RAS/MAPK Signaling by the Serotonin Metabolite 5-HIAA. PLoS Genet 11, e1005236 (2015). https://doi.org:10.1371/journal.pgen.1005236 41 Balasubramanian, V., Sellegounder, D., Suman, K. & Krishnaswamy, B. Proteome analysis reveals translational inhibition of Caenorhabditis elegans enhances susceptibility to Pseudomonas aeruginosa PAO1 pathogenesis. J Proteomics 145, 141- 152 (2016). https://doi.org:10.1016/j.jprot.2016.04.025 42 Gallotta, I. et al. Extracellular proteostasis prevents aggregation during pathogenic attack. Nature 584, 410-414 (2020). https://doi.org:10.1038/s41586-020-2461-z 43 Buck, S. H., Chiu, D. & Saito, R. M. The cyclin-dependent kinase inhibitors, cki-1 and cki-2, act in overlapping but distinct pathways to control cell cycle quiescence during C. elegans development. Cell Cycle 8, 2613-2620 (2009). https://doi.org:10.4161/cc.8.16.9354 44 Albert, P. S. & Riddle, D. L. Developmental alterations in sensory neuroanatomy of the Caenorhabditis elegans dauer larva. J Comp Neurol 219, 461-481 (1983). https://doi.org:10.1002/cne.902190407 45 Hall, S. E., Beverly, M., Russ, C., Nusbaum, C. & Sengupta, P. A cellular memory of developmental history generates phenotypic diversity in C. elegans. Curr Biol 20, 149- 155 (2010). https://doi.org:10.1016/j.cub.2009.11.035 46 Kyani-Rogers, T. et al. Developmental history modulates adult olfactory behavioral preferences via regulation of chemoreceptor expression in Caenorhabditiselegans. Genetics 222 (2022). https://doi.org:10.1093/genetics/iyac143 47 Yim, H. et al. Comparative connectomics of dauer reveals developmental plasticity. Nat Commun 15, 1546 (2024). https://doi.org:10.1038/s41467-024-45943-3 48 Morran, L. T., Cappy, B. J., Anderson, J. L. & Phillips, P. C. Sexual partners for the stressed: facultative outcrossing in the self-fertilizing nematode Caenorhabditis elegans. Evolution 63, 1473-1482 (2009). https://doi.org:10.1111/j.1558-5646.2009.00652.x 49 Park, J. H. et al. Daumone fed late in life improves survival and reduces hepatic inflammation and fibrosis in mice. Aging Cell 13, 709-718 (2014). https://doi.org:10.1111/acel.12224 50 Chandler, R. J. et al. Propionyl-CoA and adenosylcobalamin metabolism in Caenorhabditis elegans: evidence for a role of methylmalonyl-CoA epimerase in intermediary metabolism. Mol Genet Metab 89, 64-73 (2006). https://doi.org:10.1016/j.ymgme.2006.06.001 51 Fox, B. W. et al. C. elegans as a model for inter-individual variation in metabolism. Nature 607, 571-577 (2022). https://doi.org:10.1038/s41586-022-04951-3 52 Watson, E. et al. Interspecies systems biology uncovers metabolites affecting C. elegans gene expression and life history traits. Cell 156, 759-770 (2014). https://doi.org:10.1016/j.cell.2014.01.047 53 Watson, E. et al. Metabolic network rewiring of propionate flux compensates vitamin B12 deficiency in C. elegans. Elife 5 (2016). https://doi.org:10.7554/eLife.17670 54 Snoozy, J., Breen, P. C., Ruvkun, G. & Warnhoff, K. moc-6/MOCS2A is necessary for molybdenum cofactor synthesis in C. elegans. MicroPubl Biol 2022 (2022). https://doi.org:10.17912/micropub.biology.000531 55 Warnhoff, K. & Ruvkun, G. Molybdenum cofactor transfer from bacteria to nematode mediates sulfite detoxification. Nat Chem Biol 15, 480-488 (2019). https://doi.org:10.1038/s41589-019-0249-y .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint 56 Ludewig, A. H. & Schroeder, F. C. Ascaroside signaling in C. elegans. WormBook, 1-22 (2013). https://doi.org:10.1895/wormbook.1.155.1 57 Park, S. & Paik, Y. K. Genetic deficiency in neuronal peroxisomal fatty acid beta- oxidation causes the interruption of dauer development in Caenorhabditis elegans. Sci Rep 7, 9358 (2017). https://doi.org:10.1038/s41598-017-10020-x 58 Gillan, V. & Devaney, E. Nematode Hsp90: highly conserved but functionally diverse. Parasitology 141, 1203-1215 (2014). https://doi.org:10.1017/S0031182014000304 59 van den Hoogen, J. et al. Soil nematode abundance and functional group composition at a global scale. Nature 572, 194-198 (2019). https://doi.org:10.1038/s41586-019-1418-6 60 van den Hoogen, J. et al. A global database of soil nematode abundance and functional group composition. Sci Data 7, 103 (2020). https://doi.org:10.1038/s41597-020-0437-3 61 Bardgett, R. D. & van der Putten, W. H. Belowground biodiversity and ecosystem functioning. Nature 515, 505-511 (2014). https://doi.org:10.1038/nature13855 62 Crook, M. The dauer hypothesis and the evolution of parasitism: 20 years on and still going strong. Int J Parasitol 44, 1-8 (2014). https://doi.org:10.1016/j.ijpara.2013.08.004 63 Maushe, D. et al. Stress tolerance in entomopathogenic nematodes: Engineering superior nematodes for precision agriculture. J Invertebr Pathol 199, 107953 (2023). https://doi.org:10.1016/j.jip.2023.107953 64 Viney, M. & Diaz, A. Phenotypic plasticity in nematodes: Evolutionary and ecological significance. Worm 1, 98-106 (2012). https://doi.org:10.4161/worm.21086 65 Vlaar, L. E. et al. On the role of dauer in the adaptation of nematodes to a parasitic lifestyle. Parasit Vectors 14, 554 (2021). https://doi.org:10.1186/s13071-021-04953-6 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint

AND MATERIALS Growth conditions for C. elegans strains All strains were grown and maintained at 20°C unless otherwise mentioned. Animals were grown in 20°C incubators (humidity controlled) on 60mm nematode growth media (NGM) plates by passaging 8-15 L4s (depending on the strain) onto a fresh plates. Extra care was taken to ensure equal worm densities across all strains. Animals were fed Escherichia coli OP50 obtained from Caenorhabditis Genetics Center (CGC) that were seeded (OD 600=1.5 and this was strictly maintained throughout the experiments) onto culture plates 2 days before use. The NGM plate thickness was controlled by pouring 8.9ml of autoclaved liquid NGM per 60mm plate. Laboratory temperature was maintained at 20°C and monitored throughout. For all experiments, age- matched day-one hermaphrodites, or larvae timed to reach specific developmental stages as mentioned in the figure legend, were used. C. elegans strains C. elegans strains used in this study are listed in Table 1. Strains were procured from Caenorhabditis Genetics Center (CGC, Twin Cities, MN), generated in the laboratory or generated by Suny Biotech (Suzhou, Jiangsu, China 215028). Table 1 Strain Name Gene name Source Additional information N2, C. elegans var Bristol N2, Wild-type Caenorhabditis Genetics Center CB1370 daf-2 (e1370) III Caenorhabditis Genetics Center VEP032 ilc-17.1 (syb5296) X Prahlad Lab/ SunyBiotech ilc-17.1 deletion, 2173bp deletion, and the 15bp and 127bp sequences were left in the 5' and 3' deletion end, respectively of the 2135 bp ilc-17.1 gene VEP036 unc-119 (ed4); gtIs1 [CEP-1::GFP + unc-119 (+)] Prahlad Lab CEP-1 overexpression (ref: 65) CB1372 daf-7(e1372) III. Caenorhabditis Genetics Center Obtaining larvae and dauers following ‘bleach-hatching’ Populations of 250-300 gravid adults were generated by passaging L4s on NGM plates. These plates were used for obtaining synchronized embryos by bleach-induced solubilization of the adults to then obtain larvae for harvesting mRNA. Specifically, animals were washed off the plates with 1X PBS and pelleted by centrifuging at 2665Xg for 30s. The PBS was removed carefully, and worms were gently vortexed in the presence of bleaching solution [250µl 1N NaOH, 200µl standard (regular) bleach and 550µl sterile water] until all the worm bodies had dissolved (approximately 5-6 minutes), and only eggs were viable. The eggs were pelleted by centrifugation (2665Xg for 45s), bleaching solution was carefully removed and then embryos .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint were washed with sterile water 4-5 times and counted under the microscope. The desired number of embryos were seeded on fresh OP50 plates and allowed to grow at 25°C or 27°C for specific time periods depending on the experimental need. If >5% eggs remained unhatched, these plates were discarded.  32 hrs L2/L3 larvae: Day-one adult worms grown at 20C (I-36LLVL Incubator) were bleach-hatched and ~3200 eggs/genotype (~800 eggs/plate and 4 plates/genotype) were seeded on fresh OP50 plates and allowed to grow for either 15 hours or 32 hours at 25°C in Echotherm incubator IN30-2. Worms were washed with sterile water and total RNA was extracted from biological triplicates using the Direct-zol RNA Miniprep (catalog no. R2050, Zymo Research).  48 hrs dauer larvae: Day-one adult worms grown at 20C (I-36LLVL Incubator) were bleach-hatched, and ~3200 eggs/genotype (~800 eggs/plate and 4 plates/genotype) were seeded on fresh OP50 plates and allowed to grow for 48 hours at 25°C Echotherm incubator IN30-2. All strains other than N2 (wildtype control) showed >99.5% dauers at 48hrs. Any non-dauer worms were picked off the plate before collection to avoid variable staged worms in the RNA prep. N2 (wildtype control) dauers were obtained by incubating the eggs for 48 hours at 27°C in New Brunswick Galaxy 170S Incubator. These plates consisted of >96% dauers and all non-dauer worms were picked off to avoid variable staged worms in the RNA prep. Worms were washed with sterile water, and total RNA was extracted from biological triplicates using the Trizol extraction method. RNA extraction methods Trizol extraction: 300 µl of Trizol (catalog no. 400753, Life Technologies) was added to the samples after collection and snap-frozen immediately in liquid nitrogen. Samples were thawed on ice and then lysed using a Precellys 24 homogenizer (Bertin Corp.). RNA was then purified as detailed with appropriate volumes of reagents modified to 300 µl of Trizol. The RNA pellet was dissolved in 17 µl of RNase-free water. The purified RNA was then treated with deoxyribonuclease using the TURBO DNA-free kit (catalog no. AM1907, Life Technologies) as per the manufacturer’s protocol. cDNA was generated by using the iScript cDNA Synthesis Kit (catalog no. 170–8891, Bio-Rad). qRT-PCR was performed using PowerUp SYBR Green Master Mix (catalog no. A25742, Thermo Fisher Scientific) in QuantStudio 3 Real-Time PCR System (Thermo Fisher Scientific) at a 10 µl sample volume, in a 96-well plate (catalog no. 4346907, Thermo Fisher Scientific). The relative amounts of mRNA were determined using the ΔΔCt method for quantitation. We selected pmp-3 as an appropriate internal control for gene expression analysis in C. elegans. All relative changes of mRNA were normalized to either that of the wild-type control or the control for each genotype (specified in figure legends). Each experiment was repeated a minimum of three times. For qPCR reactions, the amplification of a single product with no primer dimers was confirmed by melt-curve analysis performed at the end of the reaction. Direct-zol RNA Miniprep (catalog no. R2050, Zymo Research): Instructions on the kit were followed. An equal volume of ethanol (95-100%) was added to the samples and mixed thoroughly. The mixture was transferred into a Zymo-Spin™ IICR Column in a Collection Tube and centrifuged at 10,000-16,000 x g for 30 seconds. The column was transferred into a new collection tube, and the flow-through was discarded. The sample was treated with DNase I and incubated at room temperature (20-30°C) for 15 minutes. 400 µl Direct-zol™ RNA PreWash was added to the column and centrifuged at 10,000-16,000 x g for 30 seconds. The flow-through was discarded, and the previous step was repeated. 700 µl RNA Wash Buffer was added to the column and centrifuge for 1 minute to ensure complete removal of the wash buffer. The column was transferred carefully into an RNase-free tube. RNA was eluted by adding 50 µl of .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint DNase/RNase-Free Water directly to the column matrix and centrifuged at 10,000-16,000 x g for 30 seconds. RNA-seq analysis RNA seq data for C. elegans dauer was analyzed and processed using nf-core/rnaseq v3.12.01,and the pipeline was executed with Nextflow v22.10.62. In short, the quality of the reads was assessed with Fastqc v0.11.9, and Trimgalore v0.673 was used to filter low quality reads and remove adapters. The processed reads were to aligned to C. elegans genome (WBcel235 Ensembl release 111 4), using STAR v2.7.9a 5 with default settings. Alignments were quantified with Salmon using the WBCel235 annotation. Quality control of the alignment was performed with Qualimap v2.2.26 and RseQC v3.0.17.  Differential Expression and Normalization The transcript-level estimates were summarized to gene-level counts and to gene-level Transcript per million (TPM) with the tximport package 8. The RNA-seq data from the L2/L3 larvae (WT at 32hours) was analyzed as previously described9. Gene-level counts from this study were merged with the Dauer gene-level counts to perform the normalization and differential expression steps. DESeq210 was used to do differential expression between the samples (Supplementary Table 19). Genes with low read counts (n<10) were removed from the differential expression analysis. Genes with an adjusted p-value of <0.05 (after correction with Benjamini & Hochberg) were considered significant. Changes in expression were presented as log2 fold-change. To obtain the mean gene expression, gene-level TPM, for each sample replicate (dauers only) were transformed log10 (TPM+1) and averaged between the replicates.  Gene expression variability Principal Component Analysis (PCA) and pairwise distance analysis were performed after transforming the raw counts with variance stabilization transformation (VST). PCA was done using the genes with highest variance (top 500), meanwhile the pairwise distances were calculated by determining the Euclidean distance between all the genes and then using hierarchical clustering. Regularized Log Transformation (rlog) transformation for Gene-level counts was used to determine gene variability and expression (Heatmap and meanSdPlot). Coefficient of Variation (CV) of the raw gene-level counts was used to determine the level of gene variability across all samples.  Correlation Analysis Spearman’s correlation and confidence intervals were calculated using the R package psych 11. Global gene correlation was calculated by filtering genes with a TPM value <1, and then plotting them as log10TPM values in scatter plots with the R package ggpubr12. Correlations for Pathways or GO categories between the different type of dauers (WT HID, ilc- 17.1, daf-2, daf-7, cep-1 OE) were calculated by using the log10TPM values of the specified subset of genes, then plotted as dumbbells or heatmap, p-values were adjusted form multiple tests using Bonferroni correction. Rho values <0.6 were consider low, and adjusted p-values <0.05 were considered significant. The significance of the difference between the correlation coefficients was done by performing a Fisher Z transformation with R package TOSTER 13. P-values values <0.05 were considered significant.  Functional Analysis .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint Gene Ontology Analysis was performed using ClusterProfiler14 and Wormbase Enrichment Suite 15. Ontology terms were obtained from the R package org.Ce.eg.db 16 and from Wormbase 17. DEGs and CV ranked genes were used to perform an Over Representation Analysis (ORA). Ontology terms with an adjusted p-value <0.05 (Benjamini & Hochberg) or q-value <0.05 were considered significant. To summarize and visualize the Gene Ontology terms the package rrvgo18 was used. Briefly, the GO terms were grouped by semantic similarity and simplified to their parent terms, using 90% threshold of similarity. Enriched parent terms were plotted using an UMAP projection.  Visualization Venn diagram was generated using the package venn 19. Heatmaps were generated using the package ComplexHeatmap20 and pheatmap21. All the statistical analysis were done in R 22. Data Availability The dauer expression data have been deposited in NCBI's Gene Expression Omnibus 23 with the GEO accession number GSE274872. The WT L2/L3 (WT 32 hours) data was previously published in NCBI’s Gene Expression Omnibus with accession numbers GSE218596 and GSE229132 1 Ewels, P . A. et al. The nf-core framework for community-curated bioinformatics pipelines. Nat Biotechnol 38, 276-278, doi:10.1038/s41587-020-0439-x (2020). 2 Di Tommaso, P . et al. Nextflow enables reproducible computational workflows. Nat Biotechnol 35, 316-319, doi:10.1038/nbt.3820 (2017). 3 Krueger, F. (2021). 4 Yates, A. D. et al. Ensembl Genomes 2022: an expanding genome resource for non- vertebrates. Nucleic Acids Res 50, D996-D1003, doi:10.1093/nar/gkab1007 (2022). 5 Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15-21, doi:10.1093/bioinformatics/bts635 (2013). 6 Garcia-Alcalde, F. et al. Qualimap: evaluating next-generation sequencing alignment data. Bioinformatics 28, 2678-2679, doi:10.1093/bioinformatics/bts503 (2012). 7 Wang, L., Wang, S. & Li, W. RSeQC: quality control of RNA-seq experiments. Bioinformatics 28, 2184-2185, doi:10.1093/bioinformatics/bts356 (2012). 8 Soneson, C., Love, M. I. & Robinson, M. D. Differential analyses for RNA-seq: transcript- level estimates improve gene-level inferences. F1000Res 4, 1521, doi:10.12688/f1000research.7563.2 (2015). 9 Godthi, A. et al. Neuronal IL-17 controls Caenorhabditis elegans developmental diapause through CEP-1/p53. Proceedings of the National Academy of Sciences 121, e2315248121, doi:doi:10.1073/pnas.2315248121 (2024). 10 Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550, doi:10.1186/s13059-014-0550-8 (2014). 11 psych: Procedures for Psychological, Psychometric, and Personality Researc (Northwestern University, 2024). 12 ggpubr: 'ggplot2' Based Publication Ready Plots v. 0.6.0 (2023). 13 Caldwell, A. R. Exploring Equivalence Testing with the Updated TOSTER R Package. PsyArXiv, doi:https://doi.org/10.31234/osf.io/ty8de (2022). 14 Wu, T. et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (N Y) 2, 100141, doi:10.1016/j.xinn.2021.100141 (2021). 15 Angeles-Albores, D., Lee, R., Chan, J. & Sternberg, P. Two new functions in the WormBase Enrichment Suite. MicroPubl Biol 2018, doi:10.17912/W25Q2N (2018). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint 16 Carlson, M. (2019). 17 Davis, P . et al. WormBase in 2022-data, processes, and tools for analyzing Caenorhabditis elegans. Genetics 220, doi:10.1093/genetics/iyac003 (2022). 18 Sayols, S. rrvgo: a Bioconductor package for interpreting lists of Gene Ontology terms. MicroPubl Biol 2023, doi:10.17912/micropub.biology.000811 (2023). 19 Dusa, A. (2024). 20 Gu, Z. Complex heatmap visualization. iMeta 1, e43, doi: https://doi.org/10.1002/imt2.43 (2022). 21 Kolde, R. (2019). 22 R: A Language and Environment for Statistical Computing (2023). 23 Edgar, R., Domrachev, M. & Lash, A. E. Gene Expression Omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res 30, 207-210, doi:10.1093/nar/30.1.207 (2002). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint FIGURE LEGENDS Figure 1. Experimental design. A, B. Models for phenotypic plasticity (A) and phenotypic robustness (B). C. Schematic of experimental design: RNA sequencing (RNA-seq) of high temperature- induced wild-type dauers (N2; High Temperature Induced Dauers or wild type HID or WT HID), daf-2 (e1370) III, daf-7(e1372) III and ilc-17.1 (syb5296) X dauers, and dauers that

from overexpressing cep-1 (cep-1oe). Samples were collected as described in text. Wild-type larvae that were grown for 32 hours post-hatching at 25°C to reach late L2/early L3 stage were used as comparisons. D. Micrographs of dauer larvae as used for RNA-seq. Scale bar=1 mm. E. Pair-wise distance matrix of RNA-seq samples shows the expected clustering of total RNA of the biological triplicates of each strain [Strains used: wild type (N2) L2/L3, wild type (N2) HID, daf-2 (e1370) III, daf-7(e1372) III and ilc-17.1 (syb5296) X and cep-1 (cep-1oe). F. Principal Component Analysis (PCA) of the three repeats of RNA-seq samples. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint Figure 2. Dauers induced by different environmental and physiological stimuli utilize similar processes for dauer arrest. A. Venn Diagram showing overlap between differentially expressed genes (DEGs) in wild- type HID, daf-2, daf-7, ilc-17.1 and cep-1oe dauers compared to wild type L2/L3 larvae. B. Dotplot showing comparison of enrichment between downregulated DEGs in wild type HID, daf-2, daf-7, ilc-17.1 and cep-1oe dauers compared to wild type L2/L3 larvae. Y axis: GO categories (Biological Processes). Color bar: adjusted p-values (Benjamini- Hochberg corrected, p<0.05), lower p-value in red, higher p-value blue. Circle size: Fold Enrichment (Gene Ratio/Background Ratio). C. Dotplot showing comparison of enrichment between upregulated DEGs in wild type HID, daf-2, daf-7, ilc-17.1 or cep-1oe dauers. Y axis: GO categories (Biological Processes). Color bar: adjusted p-values (Benjamini-Hochberg corrected, p<0.05), lower p-value in red, higher p-value blue. Circle size: Fold Enrichment. D. Dotplot showing enrichment of DEGs unique to wild type HID, daf-2, daf-7, ilc-17.1 or cep-1oe dauers. Y axis: GO categories (Wormbase). Color bar: adjusted p-values (Q.value<0.05), lower p-value in red, higher p-value blue. Circle size: Fold Enrichment. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint Figure 3. Dauers display high variability in gene expression levels but strong correlation in gene expression patterns. A. Density plot of the Coefficient of Variation (CV). X-axis: CVs of all expressed genes (mean expression >10 counts) computed by dividing the standard deviation (SD) across all dauers by its mean expression across all dauers. Dotted lines demarcating CVs: low, (CV < 30), moderate (30 < CV 50). Y-axis (left) density, (right): gene count. B. Scatter plots showing pairwise Spearman correlation between the wild type (N2) HID, daf-2 (e1370) III, daf-7(e1372) III and ilc-17.1 (syb5296) X and cep-1 (cep-1oe). Line represents linear regression. TOP: Dauers that are compared, and Spearman rho values shown. p-value is corrected for multiple tests; Benjamini-Hochberg. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint Figure 4. Gene expression in the core dauer pathway is robust. A. Dotplot showing enrichment of less variable (low; CV< 30), moderately variable (moderate; 30<CV 50) genes. Y axis: the GO categories. Color bar show adjusted p-values (Benjamini-Hochberg corrected, p<0.05), lower p-value in red, higher p-value blue. Circle size: Fold Enrichment. B. Density plot of the Coefficient of Variation (CV) of genes in the core dauer pathway. Labels of pathway on top. X-axis: CVs of all genes in pathway computed by dividing the standard deviation (SD) of its expression across all dauers by its mean expression across all dauers. Dotted lines demarcating CVs: low, (CV < 30), moderate (30 < CV 50). Y-axis (left) density. C. Dumbbell Plot showing Spearman’s correlation coefficient (rho) between pairs of dauers in ‘core’ dauer pathways. TOP labels: the pair of dauers compared. Y-axis: pathway. X- axis: rho value and confidence interval (CI): blue dot represents lower CI, yellow dot represents rho value, and red dot represents high CI. p-values, Bonferroni corrected. *p- value <0.05; ns, not-significant. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint Figure 5. Correlation within GO categories enriched by genes differentially expressed by all dauers shows variation in a small subset of categories. A. Scatter plot showing summarized GO categories (Biological Process) enriched from 13210 differentially expressed genes (DEGs) common to all dauers [dauers compared to wild type L2/L3 larvae]. Axes are represented by two UMAP components. Distance between points represent similarity between terms. Size of the points: score in the dissimilarity matrix. B. Heatmap depicting Spearman’s correlation coefficient (rho) for pairwise comparisons between all dauers in each GO category enriched, as in A. Colorbar: Black-white: rho >0.6, purple: rho0.05 (indicates low correlation). C. Dumbbell Plot showing Spearman’s correlation coefficient (rho) between ‘common dauer genes’ compared between pairs of dauers and analyzed within GO categories depicted in A. TOP labels: the pair of dauers compared. Y-axis: GO categories that were collapsed in A. X-axis: rho values and confidence interval (CI): blue dot represents lower CI, yellow dot represents rho value, and red dot depicts high CI. p-values, Bonferroni corrected. *p-value <0.05; ns, not-significant. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint Figure 6. Correlation within categories of co-expressed genes in the C. elegans gene expression map [Kim et al (2001)], and the iCEL1314 metabolic network [Nanda et al. (2023)] shows that dauers vary in stress and metabolic pathways. A. Heatmap depicting Spearman’s correlation coefficient (rho) for pairwise comparisons between all dauers in Kim ‘mountains’ from the C. elegans gene expression map [Kim et al (2001)]. Colorbar: Black-white: rho >0.6, purple: rho0.05 (indicates low correlation). B. Dumbbell Plot showing Spearman’s correlation coefficient (rho) between pairs of dauers compared within ‘Kim mountains’. TOP labels: the pair of dauers compared. Y-axis: pathway. X-axis: rho value and confidence interval (CI): blue dot represents lower CI, yellow dot represents rho value, and red dot depicts high CI. p-values, Bonferroni corrected. *p-value <0.05; ns, not-significant. C. Heatmap depicting Spearman’s correlation coefficient (rho) for pairwise comparisons between all dauers within metabolic pathways in the iCEL1314 metabolic network [Nanda et al. (2023)]. Colorbar: Black-white: rho >0.6, purple: rho0.05 (indicates low correlation). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint Figure 7. Correlation between dauer larvae in metabolic pathway-gene expression reveals broad variation. A. Dumbbell Plot showing Spearman’s correlation coefficient (rho) between pairs of dauers compared within metabolic pathways identified in the iCEL1314 metabolic network [Nanda et al. (2023)]. TOP labels: the pair of dauers compared. Y-axis: pathway. X-axis: rho value and confidence interval (CI): blue dot represents lower CI, yellow dot represents rho value, and red dot depicts high CI. p-values, Bonferroni corrected. *p- value <0.05; ns, not-significant. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint Figure 8. Gene expression differences between dauer larvae induced to arrest development by different stimuli. A limited number of gene expression pathways that regulate specific morphological features, and specific stress, immune, detoxification and metabolic programs vary between the different dauer larvae. These are summarized here. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint

We thank the V.P. laboratory and Dr. Gidalevitz (Drexel University) for comments. Nematode strains were provided by the Caenorhabditis Genetics Center (CGC) (funded by the NIH Infrastructure Programs P40 OD010440). This work was supported by NIH R01 AG060616 (V.P.) and by National Cancer Institute (NCI) grant P30CA016056 involving the use of Roswell Park Comprehensive Cancer Center’s Pathology Network, Genomic, and Biomedical Research Informatics Shared Resources. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint Author Contributions All authors designed the study and performed experiments. J.C.C. and V.P. designed the analyses, J.C.C. conducted the analysis data, and J.C.C. and V.P . drafted the manuscript. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint Figure 1 Phenotype Plasticity Phenotypic RobustnessA C RNA-seq prep B Wild type HID ilc-17.1 daf-7daf-2 cep-1 o/e D E F .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint Figure 2 Upregulated DEGs GO Biological Processes Downregulated DEGs GO Biological Processes A C Wormbase GO (Unique Genes) 109 69 115 87 76 4971 43 70 3577 138272 84494 27 32 72 201 21 43 94395 18 55118 558 178 613 1011 13120 daf-7 daf-2 ilc-17.1 cep-1 OE WT HID DEGs in the dauer larvae B D .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint Figure 3 A B Density CV Scatterplot Gene Expression .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint low high moderate Figure 4 Core dauer pathways A C B CV .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint Figure 5 Spearman correlation of subset of common dauer genes (padj. >0.05 or ρ< .6) A C UMAP GO Biological Process B GO categories .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint Figure 6 Kim Mountains Metabolic pathways A B C .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint Figure 7 A .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint Pathways different from other DauersDauer Amine Oxidase Wild type (N2) HID Folate Cycle Glycine Cleavage System Histidine Degradation Iron Metabolism Mevalonate Metabolism Propionate Degradation Canonical Ros Metabolism Glyoxylate and dicarboxylate metabolism Taurine and Hypoxanthine metabolism PUFA Biosynthesis (except daf-2) Molybdenum Cofactor Biosynthesis (except daf-2) Chitin degradation (except cep-1 oe) PALS (except daf-2) Galactose metabolism (except cep-1 oe) UGT (except cep-1 OE) Protein Complex Oligomerization Amine Oxidase daf-2 (e1370)III Histidine Degradation Methionine Salvage Mevalonate Metabolism Propionate Degradation Canonical Complex II of ETC (except ilc-17.1) Fatty Acid Degradation, other (except ilc-17.1) Glycine Cleavage System (except ilc-17.1) Molecular Transducer Activity Amine Oxidase daf-7 e(1372)III Glycine Cleavage System Histidine Degradation Propionate Degradation Canonical TGFβ (except N2 HID) Mevalonate Metabolism (except cep-1 OE) Methylglyoxyl detoxification (except cep-1 oe) Narrow Pore Channel Activity Chemical Synaptic Transmission Postsynaptic Stabilization Of Membrane Potential Postsynaptic Membrane Regulation Of Postsynaptic Membrane Potential Oxidoreductase Activity Acting On Paired Donors With Incorporation Or Reduction Of Molecular Oxygen Reduced Flavin Or Flavoprotein As One Donor And Incorporation Of One Atom Of Oxygen Steroid Hydroxylase Activity Extracellular Ligand-Gated Monoatomic Ion Channel Activity Passive Transmembrane Transporter Activity Outward Rectifier Potassium Channel Activity Gated Channel Activity Response To Xenobiotic Stimulus Amine Oxidase ilc-17.1 (syb5296)X Histidine Degradation Methionine Salvage Mevalonate Metabolism Propionate Degradation Canonical PUFA Biosynthesis Ubiquinone metabolism (except daf-2) Ros Metabolism (except daf-2) Genitalia Morphogenesis (except daf-2) Microtubule Depolymerization (except daf-2) Glycine Cleavage System (except daf-2) Fatty Acid Degradation, other (except daf-2) Ciliary Plasm Protein Heterodimerization Activity Sodium Channel Activity Complex II of ETC (except daf-2) Structural Constituent Of Chromatin Gated Channel Activity Amine Oxidase cep-1 OE TGFβ Glycine Cleavage System Glyoxylate and dicarboxylate metabolism Pantothenate and COA Biosynthesis Histidine Degradation Propionate Degradation Canonical Galactose metabolism (except N2 HID) Ros Metabolism (except daf-7) PUFA Biosynthesis (except daf-7) Folate Cycle (except daf-7) Mevalonate Metabolism (except daf-7) Methylglyoxyl detoxification (except daf-7) Molybdenum Cofactor Biosynthesis (except daf-7) Lysine degradation (except daf-7) Peroxisomal Fatty Acid Degradation (except daf-7) Nuclear-Transcribed mRNA Catabolic Process, Deadenylation-Dependent Decay (except daf-7) mRNA Destabilization (except daf-7) Microtubule Depolymerization (except daf-7) Figure 8 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 19, 2024. ; https://doi.org/10.1101/2024.08.15.608164doi: bioRxiv preprint